%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Treatment model with Lalonde data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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);
% Treatment equation
%%%%%%%%%%%%%%%%%%%%%
X1=[ones(n,1) data(:,2:7) data(:,9) data(:,12)];
%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 nodegree =1 if no high school degree
%8 re74 =real earnings, 1974, $1000s
%9 unem74 =1 if unemployed all of 1974
% nodegree and unem74 are the identification variables that are not in X2
X2=[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)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% prepare data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
m=2; %number of equations
y=data(:,11); %equation 2 outcome - real earnings, 1978, $1000s
T=data(:,1); %0/1 treatment dummy
k1=size(X1,2);
k2=size(X2,2);
k=k1+k2;
%generate Xi and collect as separate matrices
X=cell(n,1);
for i=1:n
X{i}=[X1(i,:) zeros(1,k2); zeros(1,k1) X2(i,:)];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% starting values, priors, and tuners
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% general elements
r1=20000; % burn-ins
r2=10000; % keepers
R=r1+r2;
% generic OLS
bols1=inv(X1'*X1)*X1'*y;
bols2=inv(X2'*X2)*X2'*y;
clear X1 X2;
% elements for beta
% note: no need to distinguish between coefficients across the two
% equations at this point
mu0=zeros(k,1); %diffuse prior for mean of betas
V0=eye(k)*100; % diffuse prior for varcov of beta
betadraw=[bols1;bols2];
% elements for E
% E will be 2 by 2
v0=m; % number of equations
S0=eye(m); % So these are the diffusest priors possible for the IW
Edraw=iw_sig11(v0,inv(S0),1); % use Nobile's (2000) contrained IW draws-routine
% elements for y1star (participation) and y2star (outcome)
f=find(T==0); % zero amount
g=find(T~=0); % positive amount
nf=length(f); % number of zero cases
ng=length(g); % number of >0 cases
l0=-inf*ones(nf,1); %lower truncation point for zero cases
r0=zeros(nf,1);% upper truncation point for zero cases
lpos=zeros(ng,1);% lower truncation point for positive amounts
rpos=inf*ones(ng,1);%upper truncation point for positive amount
Tdraw = T;
[betamat,Emat]=gs_treat(X,k1,k,m,n,r1,r2,mu0,V0,betadraw,v0,S0,Edraw,...
f,g,nf,ng,l0,r0,lpos,rpos,Tdraw,y);
'GS done'
% put all draws together & run diagnostics
allmat=[betamat;Emat];
kdiag=klausdiagnostics_greater0(allmat);
% Open log file
[fid]=fopen('c:\Klaus\AAEC6564\mlab\logs\mod9_treatJobTrain.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');
% beta stuff
out=kdiag(1:k1,:)';
fprintf(fid,'Output table for betas, treatment equation \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');
out=kdiag(k1+1:k,:)';
fprintf(fid,'Output table for betas, outcome equation \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');
out=kdiag(k+1:end,:)';
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');
%Correlation coefficient
%Estuff as correlations
%output in correlation form
rho12=Emat(1,:)./sqrt(Emat(2,:));
kdiag=klausdiagnostics_greater0(rho12);
out=kdiag';
fprintf(fid,'Output table for correlation\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_treatJobTrain betamat Emat;
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;