%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% AR1 model with ic MH
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
rand('state',37); % set arbitrary seed for uniform draws
randn('state',37); % set arbitrary seed for normal draws
tic;  % start stop watch

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

T=size(y,1);
k=size(X,2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% starting values, priors, and tuners
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
r1=10000; % burn-in
r2=10000; % keepers
R = r1+r2;

% generic OLS (ignoring nonlinear part, just for s2)
bols=(X'*X)\X'*y;
res=y-X*bols;
s2=(res'*res)/(T-k);

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

% elements for sig2
sig2draw=s2; % use OLS variance as starting draw for GS
nu0=1/2;
tau0=1/2; % diffuse prior shape and scale

% elements for rho
murho=0;
Vrho=2;
tc=1; %tuner for MH
tv=10;
rhodraw=tnormrnd(murho,sqrt(Vrho),-1,1);
rhodraw=0.5;

% 'ok'
% return

[betamat,sig2mat,rhomat,acceptmat]=...
gs_AR1_ic(X,y,k,T,r1,r2,...
mu0,V0,betadraw,...
nu0,tau0,sig2draw,...
murho,Vrho,tc,tv,rhodraw);


'GS done'

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

fprintf(fid,'total number of iterations\t%6.0f \n',R);
fprintf(fid,'burn-in iterations\t%6.0f\n',r1);
fprintf(fid,'acceptance rates in MH part for rho \t%6.3f\n',acceptmat);

fprintf(fid,'True betas \n\n');
fprintf(fid,'%6.3f\n',btrue);
fprintf(fid,'\n');

fprintf(fid,'True sig2 \n\n');
fprintf(fid,'%6.3f\n',sig2true);
fprintf(fid,'\n');

fprintf(fid,'True rho \n\n');
fprintf(fid,'%6.3f\n',rhotrue);
fprintf(fid,'\n');


% put all draws together & get diagnostics
allmat=[betamat;sig2mat;rhomat];%
kdiag=klausdiagnostics_greater0(allmat);

% beta stuff
out=kdiag(1:k,:)';
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+1,:)';
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');

%rho stuff
out=kdiag(k+2,:)';
fprintf(fid,'Output table for rho \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');



save c:\klaus\AAEC6564\mlab\worksp\AR1ic betamat sig2mat rhomat;


fprintf(fid,'Elapsed time in minutes\n\n');
fprintf(fid,'%6.0f\t%6.8f\n',toc/60);

st=fclose(fid);

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


