%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Savage-Dickey Density Ratio (using SUR model as an example)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Set random number generator, start stop watch, open output file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
rand('state',37); % set arbitrary seed for uniform draws
randn('state',37); % set arbitrary seed for normal draws

tic;  % start stop watch

[fid]=fopen('c:\klaus\AAEC6564\mlab\logs\mod4_sddr.txt','w');
if fid==-1;
   warning('File could not be opened');
   return
else;
   disp('File opened successfully');
end;

% load SUR results
load c:\klaus\AAEC6564\mlab\worksp\mod4_sur betamat Emat;

%load & prepare data& define priors
load c:\klaus\AAEC6564\mlab\worksp\sur_data;
n=length(y1); %since all equations have same # obs's

k1=size(X1,2);
k2=size(X2,2);
k3=size(X3,2);

k=k1+k2+k3;
M=3; % number of equations

% create individual-specific panels - useful for Gibbs Sampler
Xcell=cell(n,1);
ycell=cell(n,1);
for i=1:n
    Xcell{i}=[X1(i,:) zeros(1,k2+k3); zeros(1,k1) X2(i,:) zeros(1,k3);...
        zeros(1,k1+k2) X3(i,:)];
    ycell{i}=[y1(i);y2(i);y3(i)];
end
X=cell2mat(Xcell); %convert cell array to full matrix
y=cell2mat(ycell);

% elements for beta
mu0=zeros(k,1); %diffuse prior for mean of betas
V0=eye(k)*100; % diffuse prior for varcov of beta

% elements for E
% E will be 2 by 2
v0=M+1; % number of equations
S0=eye(M); % So these are the diffusest priors possible

% test restriction that first beta = 1.2 (the true beta1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%evaluate prior density under restriction
denom=normpdf(1.2,mu0(1),sqrt(V0(1,1))); %since the marginal prior of beta1 is normal(0,100)
% -see sur script; normpdf requires the standard deviation, i.e. 10 in this
% case

%evaluate conditional posterior at all possible values of E

R=size(betamat,2);
nummat=zeros(R,1);

for i=1:R
   Edraw=vec_var(Emat(:,i));
   %re-construct VCOV E from vector of elemnets captured in Emat
   int1=0; 
   int2=0;
   j=1;
   for j=1:n
       Xi=cell2mat(Xcell(j)); % note: round brackets needed for this
       yi=cell2mat(ycell(j));
       int1=int1+Xi'*inv(Edraw)*Xi;
       int2=int2+Xi'*inv(Edraw)*yi;
   end
   V1=inv(inv(V0)+int1);
   mu1=V1*(inv(V0)*mu0+int2);
   nummat(i)=normpdf(1.2,mu1(1),sqrt(V1(1,1))); 
   %evaluate marginal posterior at b1=1.2
  end

num=mean(nummat);

logSDDR=log(num)-log(denom);
fprintf(fid,'logSDDR for b1=1.2 =\t%6.4f \n',logSDDR);


% test restriction that first three betas are 1.2, 0.4, and 0.8,
% respectively (true betas)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%evaluate prior density under restriction
denom=mvnpdf([1.2 0.4 0.8]',mu0(1:3),V0(1:3,1:3)); 

%evaluate conditional posterior at all possible values of E
R=size(betamat,2);
nummat=zeros(R,1);

for i=1:R
   Edraw=vec_var(Emat(:,i));
   int1=0; 
   int2=0;
   j=1;
   for j=1:n
       Xi=cell2mat(Xcell(j)); 
       yi=cell2mat(ycell(j));
       int1=int1+Xi'*inv(Edraw)*Xi;
       int2=int2+Xi'*inv(Edraw)*yi;
   end
   V1=inv(inv(V0)+int1);
   mu1=V1*(inv(V0)*mu0+int2);
   nummat(i)=mvnpdf([1.2 0.4 0.8]',mu1(1:3),V1(1:3,1:3)); 
end

num=mean(nummat);

logSDDR=log(num)-log(denom);

fprintf(fid,'logSDDR for b1=1.2, b2=0.4, b3=0.8 =\t%6.4f \n',logSDDR);


st=fclose(fid);
if st==0;
   disp('File closed successfully');
else;
   warning('Problem with closing file');
end;










