%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Multivariate normal regression model via Bayesian estimation
%  with independent priors
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Application: female earnings regression

% 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\mod2_application.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\labordata.txt;
% data=labordata;
% clear labordata;
% save c:\klaus\AAEC6564\mlab\worksp\labordata data;
% 
% 'ok'
% break


load c:\klaus\AAEC6564\mlab\worksp\labordata;
%We will use  Mroz's (1987) labor data described in Greene (5th edition),
%p. 51

% Contents of Data (columns)

%1       LFP = A dummy variable = 1 if woman worked in 1975, else 0 
%2       WHRS = Wife's hours of work in 1975 
%3       KL6 = Number of children less than 6 years old in household 
%4       K618 = Number of children between ages 6 and 18 in household 
%5       WA = Wife's age 
%6       WE = Wife's educational attainment, in years 
%7       WW = Wife's average hourly earnings, in 1975 dollars 
%8       RPWG = Wife's wage reported at the time of the 1976 interview (not = 1975 estimated wage) 
%9       HHRS = Husband's hours worked in 1975 
%10      HA = Husband's age 
%11      HE = Husband's educational attainment, in years 
%12      HW = Husband's wage, in 1975 dollars 
%13      FAMINC = Family income, in 1975 dollars 
%14      WMED = Wife's mother's educational attainment, in years 
%15      WFED = Wife's father's educational attainment, in years 
%16      UN = Unemployment rate in county of residence, in percentage points. 
%17      CIT = Dummy variable = 1 if live in large city (SMSA), else 0 
%18      AX = Actual years of wife's previous labor market experience 

f=find(data(:,1)==0); % we only include women that worked in 1975
data(f,:)=[];

y=log(data(:,2).*data(:,7)); %dependent variable = log(hrs worked* hourly wage)
% = log(annual earnings)
KL6=data(:,3); %number of children under 6
K618=data(:,4); % Number of children between ages 6 and 18 in household 
age=data(:,5);
age2=age.^2;
edu=data(:,6); %years of education

n=length(y);


X=[ones(n,1) age age2 edu KL6 K618];
k=size(X,2);

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% starting values, priors, and tuners
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% general elements
r1=5000; % 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);

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


%focus on marginal effects of education, K6 and K18
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


edupost=betamat(end-2,:)';
KL6post=betamat(end-1,:)';
K618post=betamat(end,:)';
sig2post=sig2mat';


[f1,x1]=ksdensity(edupost,'kernel','epanechnikov','npoints',300);
[f2,x2]=ksdensity(KL6post,'kernel','epanechnikov','npoints',300);
[f3,x3]=ksdensity(K618post,'kernel','epanechnikov','npoints',300);
[f4,x4]=ksdensity(sig2post,'kernel','epanechnikov','npoints',300);

%figure(1);

subplot(2,2,1);
plot(x1,f1,'-k');
title('posterior for education');
set(gca,'Xlim',[-0.25 0.25]);
set(gca,'Ylim',[0 20]);

subplot(2,2,2);
plot(x2,f2,'-k');
title('posterior for KL6');
set(gca,'Xlim',[-1.5 0]);
set(gca,'Ylim',[0 3]);

subplot(2,2,3);
plot(x3,f3,'-k');
title('posterior for K618');
set(gca,'Xlim',[-0.5 0.2]);
set(gca,'Ylim',[0 10]);

subplot(2,2,4);
plot(x4,f4,'-k');
title('posterior for sig2');
set(gca,'Xlim',[1 2]);
set(gca,'Ylim',[0 5]);

save c:\klaus\AAEC6564\mlab\worksp\mod2application betamat sig2mat X y;

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;










