%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Get model probabilities and model-averaged results
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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\mod8_MC3_growth_convTest.txt','w');
if fid==-1;
   warning('File could not be opened');
   return
else;
   disp('File opened successfully');
end;

load c:\klaus\AAEC6564\mlab\worksp\m8_MC3_growth betamat sig2mat gammat X y;
k=size(X,2);
kf = k; %total number of "flexible coefficients"
n=size(X,1);
g=k^2;
ybar=mean(y);

load c:\klaus\AAEC6564\mlab\worksp\m8_MC3_growth_modProb;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Capture empirical model probabilities
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
postMprobE=relcount; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Derive analytical marg. LH for the same visited models
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
U=size(relcount,1);
logmlvec=zeros(U,1);
m=1;

for jj=1:U
   gamnew=Mmat(jj,:)'; %pick one of the visited models 
   f=find(gamnew==0);
   Xn=X;
   Xn(:,f)=[];%eliminate corresponding columns in X
   kn=size(Xn,2);
   vn=(n-1)/2;
   taun=(1/2)*((g/(1+g))*(y'*y - (y'*Xn)*((Xn'*Xn)\(Xn'*y))) + ...
       (1/(1+g))*(y-ybar)'*(y-ybar)-(g/(1+g))*n*ybar^2);
   %now we need to full expression of the marginal LH, not just the kernel
   logmlnew= -((n-1)/2)*log(2*pi)-(1/2)*log(n)+gammaln(vn)+...
       (kn/2)*log(1/(1+g))-((n-1)/2)*log(taun);
    
   logmlvec(jj)=logmlnew;
    
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Compute analytical posterior model probabilities using the approach
%outlined in Koop, ch. 11 (p. 276-277)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%subtract logML of top model from all others to avoid "zero-sum" problem

int=logmlvec-max(logmlvec);
postMprobA=exp(int)/sum(exp(int));

%compare empirical and analytical model frequencies for UP TO top 10 models,
%a la Koop ch. 11, p. 276-277
l=length(postMprobE);
out=[postMprobE(1:min(l,10)) postMprobA(1:min(l,10))]
fprintf(fid,'posterior probabilities for top ten models\n\n');
fprintf(fid,'empirical\t\tanalytical\n\n');
fprintf(fid,'%6.8f\t%6.8f\n',out');
fprintf(fid,'\n');

%Compute correlation between empirical and analytical posterior model
%probabilities a la Fernandez et al, 2001, p. 398

rho=corr([postMprobE postMprobA])

fprintf(fid,'Correlation emp. Mprob, anal. Mprob., all visited models \n\n');
fprintf(fid,'%6.5f\n',rho(1,2));
fprintf(fid,'\n');


save c:\klaus\AAEC6564\mlab\worksp\mod8_MC3_growth_convTest logmlvec ...
    postMprobA postMprobE;



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




