%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Sample selection model
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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\select_data;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% prepare data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

y=y2; %latent dep variable for outcome
y(f)=0;  % zero out observations associated with a negative entry in the selection equation
clear y1 y2;

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=15000; % burn-ins
r2=5000; % keepers
R=r1+r2;

% generic OLS
bols1=(X1'*X1)\X1'*y;
bols2=(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(y==0); % zero amount
g=find(y~=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

y1draw = y; 
y2draw = y;
ydraw=[y1draw y2draw];


[betamat,Emat]=gs_select(X,k1,k,m,n,r1,r2,mu0,V0,betadraw,v0,S0,Edraw,...
    f,g,nf,ng,l0,r0,lpos,rpos,ydraw);

'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_select.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,'true beta and gamma \t%6.3f\n',[beta1true;beta2true]); 
fprintf(fid,'\n');
fprintf(fid,'\n');


fprintf(fid,'true E \t%6.3f\n',var_vec(Etrue)); 
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');

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


save c:\klaus\AAEC6564\mlab\worksp\mod9_select 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;
