%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Hierarchical normal regression model
% v1 version uses fully augmented setup; draws bf, br
% conditional on bri
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Set random number generator, start stop watch, open output file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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\mod6_HNRM_v1.txt','w');
if fid==-1;
   warning('File could not be opened');
   return
else;
   disp('File opened successfully');
end;

% Load & prepare data
%%%%%%%%%%%%%%%%%%%%%%
load c:\klaus\AAEC6564\mlab\worksp\mod6_data;
kf=2;
kr=2;
n=200;
J=6;
obs=n*J;

X=[Xf Xr];
k=kf+kr;

% break y, Xf, and Xr into individual-specific cells
yimat=mat2cell(y,J*ones(n,1),1);
Ximat=mat2cell(X,J*ones(n,1),k);

% Estimation
%%%%%%%%%%%%%%%%%%%%%%%%%

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

% generic OLS
bols=inv(X'*X)*X'*y;
res=y-X*bols;
s2=(res'*res)/(obs-k);

% elements for betaf (fixed coefficients)
muf0=zeros(kf,1); %diffuse prior for mean of fixed coefficients
Vf0=eye(kf)*100; % diffuse prior for varcov of fixed coefficients
% since we're drawing beta first, won't need starting draw

% elements for betar (hierarchical means of random coefficients
mur0=zeros(kr,1); 
Vr0=eye(kr)*100; 
betardraw=bols(kf+1:end);

% elements for E -hierarchical VCOV
v0=kr+2; 
S0=eye(kr);
Edraw=iwishrnd(S0,v0);

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

[betafmat,betarmat,Emat,sig2mat,betarimat]...
    =gs_HNRM_v1(y,Xf,Xr,Ximat,yimat,k,kf,kr,n,J,obs,r1,r2,muf0,Vf0,...
    mur0,Vr0,betardraw,v0,S0,Edraw,eta0,tau0,sig2draw);


'GS done'
% put all draws together & run diagnostics
allmat=[betafmat;betarmat;Emat;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 betafs \t%6.3f\n',bftrue); 
fprintf(fid,'\n');
fprintf(fid,'true betars \t%6.3f\n',brtrue); 
fprintf(fid,'\n');
fprintf(fid,'true E \t%6.3f\n',var_vec(Etrue)); 
fprintf(fid,'\n');
fprintf(fid,'true sig2 \t%6.3f\n',sig2true); 
fprintf(fid,'\n');
fprintf(fid,'\n');

% betaf stuff
out=kdiag(1:kf,:)';
fprintf(fid,'Output table for betafs \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');

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

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


save c:\klaus\AAEC6564\mlab\worksp\mod6_HNRM_v1 betafmat betarmat Emat sig2mat betarimat;

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


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










