%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Treatment model with Lalonde data - ignoring correlation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tic;  % start stop watch

rand('state',37); % set arbitrary seed for uniform draws
randn('state',37); % set arbitrary seed for normal draws


% load c:\klaus\AAEC6564\mlab\worksp\jtrain2.txt;
% data=jtrain2;
% clear jtrain2;
% save c:\klaus\AAEC6564\mlab\worksp\lalonde data;
% 'ok'
% break
load c:\klaus\AAEC6564\mlab\worksp\lalonde;

% This script uses the job training data from \citeasnoun{lalonde1986}, 
% as described in \citeasnoun{wooldridge2010}, p. 928. 
% The data are labeled ``jtrain2'' and comprise 445 observations on male workers, 
% 118 of whom underwent some job training in the late 1970s. 
% The outcome of interest are real earnings in 1978 (training occurred in 1975-1977). 
% The data set is sorted by treatment (treated first, then untreated), and includes the following variables:


%1  train     =1 if assigned to job training
%2  age       =age in 1977
%3  educ      =years of education
%4  black     =1 if black
%5  hisp      =1 if Hispanic
%6  married   =1 if married
%7  nodegree  =1 if no high school degree
%8  mosinex    =No. months prior to Jan. 78 in experiment
%9  re74       =real earnings, 1974, $1000s
%10 re75       =real earnings, 1975, $1000s
%11 re78       =real earnings, 1978, $1000s
%12 unem74     =1 if unemployed all of 1974
%13 unem75     =1 if unemployed all of 1975
%14 unem78     =1 if unemployed all of 1978
%15 lre74      =log(re74); zero if re74 == 0
%16 lre75      =log(re75); zero if re75 == 0
%17 lre78      =log(re78); zero if re78 == 0
%18 agesq      =age squared
%19 mostrn     =months in training

n=size(data,1);

X=[ones(n,1) data(:,2:6) data(:,9:10) data(:,1)];

%1 intercept
%2  age       =age in 1977
%3  educ      =years of education
%4  black     =1 if black
%5  hisp      =1 if Hispanic
%6  married   =1 if married
%7  re74       =real earnings, 1974, $1000s
%8  re75       =real earnings, 1975, $1000s
%9 train       =1 if assigned to job training (Treatment indicator)

k=size(X,2);

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

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

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

% elements for beta
mu0=zeros(k,1); %diffuse prior for mean of betas
V0=eye(k)*100; % diffuse prior for varcov of beta
betadraw=bols; %use OLS values as starting draws for Gibbs Sampler

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

[betamat,sig2mat]=...
    gs_normal_independent(X,y,n,k,r1,r2,mu0,V0,betadraw,v0,tau0,sig2draw);

'GS done'

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

% Open log file 
[fid]=fopen('c:\Klaus\AAEC6564\mlab\logs\mod9_treatJobTrainNaive.txt','w');
if fid==-1;
   warning('File could not be opened');
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,'\n');

fprintf(fid,'OLS betas \t%6.3f\n',bols); 
fprintf(fid,'OLS s2 \t%6.3f\n',s2); 
fprintf(fid,'\n');
fprintf(fid,'\n');

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


save c:\klaus\AAEC6564\mlab\worksp\mod9_treatJobTrainNaive betamat sig2mat;


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;
