%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Multivariate normal regression model - false model
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

rand('state',37); % set arbitrary seed for uniform draws
randn('state',37); % set arbitrary seed for normal draws

tic;  % start stop watch - always fun to know your run time.

% Open log file - for your text output
[fid]=fopen('c:\klaus\AAEC6564\mlab\logs\mod1s2b.txt','w');
if fid==-1;
   warning('File could not be opened');
   return
else;
   disp('File opened successfully');
end;

load c:\klaus\AAEC6564\mlab\worksp\mod1s2a y X btrue sig2true betaPriordraws ...
    betaPostdraws logpy;

%%%%%%%%%%%%%%%%%%%%%%%%%%
%Use "wrong" X matrix
%%%%%%%%%%%%%%%%%%%%%%%%%%
X(:,2)=[]; %drop second variable
k=size(X,2);

btrue(2)=[]; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Conjugate estimation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% define priors - same as before
%%%%%%%%%%%%%%%%
mu0=zeros(k,1); %Prior mean for betas
c=10;
V0=c*eye(k); %Prior variance component for betas 

v0=2.1; % prior shape for ig density
tau0=20; %prior scale for ig density

% get OLS results for convenient expressions of posterior
bols=inv(X'*X)*(X'*y); 
%faster and more accurate, but less transparent:
bols=(X'*X)\(X'*y); 
e=y-X*bols;% Get residuals. 
SSE=e'*e; %get the regression error (estimated variance of "eps").

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% derive marginal posterior moments
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%sig2
%%%%%%%%%%%%%%
V1wig=inv(inv(V0)+X'*X);
v1=(2*v0+n)/2;
tau1=tau0+(1/2)*(SSE +(bols-mu0)'*V1wig*(bols-mu0));

%beta
%%%%%%%%%%%%%%
mu2=V1wig*(inv(V0)*mu0 + X'*y);
V2=(tau1/v1)*V1wig;
v2=2*v1;

%Send prior and posterior moments to output file
fprintf(fid,'$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$\n');
fprintf(fid,'Prior and posterior moments (ananlytical) \n'); 
fprintf(fid,'$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$\n'); 
fprintf(fid,'\n');
fprintf(fid,'\n');

out=[btrue mu0 mu2]; 
names=cell(k,1); %define k empty cells; note the round brackets
names{1}='constant'; 
names{2}='x2';
names{3}='x3';

fprintf(fid,'Mean for beta \n\n'); 
fprintf(fid,'%15s','variable','true','prior','post');
fprintf(fid,'\n');
for i=1:k
    fprintf(fid,'%15s',names{i}); %print ith variable name
    fprintf(fid,'%15.4f',out(i,:));%print results for variable i
    fprintf(fid,'\n');
end
fprintf(fid,'\n');
fprintf(fid,'\n');

out=[sqrt(diag(V0)) sqrt(diag(V2))];
fprintf(fid,'Std for beta \n\n'); 
fprintf(fid,'%15s','variable','prior','post');
fprintf(fid,'\n');
for i=1:k
    fprintf(fid,'%15s',names{i}); %print ith variable name
    fprintf(fid,'%15.4f',out(i,:));%print results for variable i
    fprintf(fid,'\n');
end
fprintf(fid,'\n');
fprintf(fid,'\n');


out=[sig2true (tau0/(v0-1)) (tau1/(v1-1))]; 
names{1}='sig2'; 
fprintf(fid,'Output table for sigma \n\n'); 
fprintf(fid,'%15s','variable','true','prior','post');
fprintf(fid,'\n');
for i=1:1
    fprintf(fid,'%15s',names{i}); 
    fprintf(fid,'%15.4f',out(i,:));
    fprintf(fid,'\n');
end
fprintf(fid,'\n');
fprintf(fid,'\n');

%plot prior vs. posterior
%%%%%%%%%%%%%%%%%%%%%%%%%
%Draw from the prior - straightforward:
r1=1000;
r2=10;
R=r1*r2;

%draw r1 sigma2's from its prior density
sig2vec=1./gamrnd(v0,1/tau0,r1,1);  

%for each of those, draw r2 beta vectors from their prior density

betaPriorcell=cell(r1,1);
for i=1:r1
    sig2=sig2vec(i);
    betaPriorcell{i}= mvnrnd(mu0,sig2*V0,r2);%r2 by k
end

betaPriordraws=cell2mat(betaPriorcell); %stack cells into a (r1*r2)xk matrix

%Drawing from the multivariate t is a bit tedious
%Matlab's canned mvtrnd expects a correlation matrix, while we 
%prefer to work with a generic scale matrix.
%The following works:

%We want to draw vectors of beta from t(mu2,V2,v2).

[Rho,stdV] = corrcov(V2);%stdV=std's
betaPostdraws = repmat(mu2',R,1) + mvtrnd(Rho,v2,R)*diag(stdV); %R by k

%Lets plot this for beta3:
b3Prior=betaPriordraws(:,2); %now beta3 is the second column
b3Post=betaPostdraws(:,2);

[f1,x1]=ksdensity(b3Prior,'kernel','epanechnikov','npoints',200);
[f2,x2]=ksdensity(b3Post,'kernel','epanechnikov','npoints',200);

figure(2)
plot(x1,f1,'-b',x2,f2,'-k','LineWidth',1);
title('prior and posterior distributions for $\beta_3$',...
    'interpreter','latex','fontsize',14);
xlabel('$\beta_3$','interpreter','latex','fontsize',12);
ylabel('density');
h=legend('$p\left(\beta_3\right)$','$p\left(\beta_3|y,X\right)$');
set(h,'interpreter','latex','fontsize',12);
set(gca,'Xlim',[-5 5]);
set(gca,'XTick',-5:0.2:5);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Marginal likelihood
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

logpyF=v0*log(tau0)+gammaln(v1)-...
((n/2)*log(2*pi)+gammaln(v0)) +...
(-1/2)*log(det(V0))+(1/2)*log(det(V1wig))+...
(-v1)*log(tau1);

fprintf(fid,'Log marg. LH for false model\n');
fprintf(fid,'log(mLH) =\t%6.3f \n',logpyF);
fprintf(fid,'\n');
fprintf(fid,'\n');

%%%%%%%%%
% log BF in favor of correct model
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

logBF=logpy-logpyF;
fprintf(fid,'Log BF in favor of correct model\n');
fprintf(fid,'log(mLH) =\t%6.3f \n',logBF);
fprintf(fid,'\n');
fprintf(fid,'\n');

%stop timer
finish = toc/60;
fprintf(fid,'Time elapsed in minutes \n\n');
fprintf(fid,'%6.3f\n',finish);

%save your numerical output
save c:\klaus\AAEC6564\mlab\worksp\mod1s2b y X btrue sig2true betaPriordraws ...
    betaPostdraws logpy;

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



