%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Basic MC3 model
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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.txt','w');
if fid==-1;
warning('File could not be opened');
return
else;
disp('File opened successfully');
end;
load c:\klaus\AAEC6564\mlab\worksp\mod8_SSVS_data;
Z=X;
X=Z(:,2:end); %drop constant
n=size(X,1);
%de-mean:
xbar=mean(X,1);% mean over rows, 1 by k
X=X-repmat(xbar,n,1);
k=size(X,2); %size of beta (without constant)
%note: Here, I will call the dimension of beta k by 1 and the dimension of
%theta = [alpha beta']' k+1 by 1 - it makes the code a
%bit easier to follow
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Preliminaries, Priors, Starting draws
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
r1=50000; % burn-in
r2=10000; % keepers
R = r1+r2;
%Prior for Beta
g=n; %recommended by Fernandez et al (2001)
del0=0.5; %prior probability that a regressor is in model
% Note: in that case all model priors are
% equal and drop out of the MH formula. So del0 is never actually used in
% the Gibbs Sampler code
%Starting vector for psi
%%%%%%%%%%%%%%%%%%%%%%%%%
gamstart=ones(k,1);%all variables in model
% 'ok'
% break
[amat,betamat,sig2mat,gammat,acceptcount]=...
gs_MC3(X,y,k,r1,r2,gamstart,g);
'GS done'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Capture key results
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Run diagnostics on beta and sigma
allmat=[amat;betamat;sig2mat];
kdiag=klausdiagnostics_greater0(allmat);
fprintf(fid,'total number of iterations =\t%6.0f \n',R);
fprintf(fid,'burn-in iterations =\t%6.0f \n',r1);
fprintf(fid,'\n');
fprintf(fid,'true betas \t%6.3f\n',btrue);
fprintf(fid,'true sig2 \t%6.3f\n',sig2true);
fprintf(fid,'\n');
fprintf(fid,'\n');
% beta stuff
out=kdiag(1:k+1,:)';
fprintf(fid,'Output table for betas \n\n');
fprintf(fid,'mean\t\tstd\t\tp(>0)\t\tnse\t\tIEF\t\tm*\t\tcd\n\n');
fprintf(fid,'%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\n',out);
fprintf(fid,'\n');
%sig2 stuff
out=kdiag(k+2,:)';
fprintf(fid,'Output table for sig2 \n\n');
fprintf(fid,'mean\t\tstd\t\tp(>0)\t\tnse\t\tIEF\t\tm*\t\tcd\n\n');
fprintf(fid,'%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\n',out);
fprintf(fid,'\n');
%gamma stuff
out=[mean(gammat,2)]; %inclusion probabilities for each beta
fprintf(fid,'Inclusion probabilities \n\n');
fprintf(fid,'%6.3f\n',out);
fprintf(fid,'\n');
finish = toc/60;
fprintf(fid,'Time elapsed in minutes \n\n');
fprintf(fid,'%6.3f\n',finish);
save c:\klaus\AAEC6564\mlab\worksp\m8_MC3 betamat sig2mat gammat X y;
st=fclose(fid);
if st==0;
disp('File closed successfully');
else;
warning('Problem with closing file');
end;