%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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;