Counterfactual3

%% Problem Set III

%% COUNTERFACTUAL I


%The parameters obtained are:

xo=x0;

%Loading the data
a=load('PS3_data.csv');
AGE=a(:,2);
L=a(:,3);
W=a(:,4);
Y=a(:,5);
S=a(:,7);
A=a(:,6);
B=a(:,8);

%% Loading parameters


gamma0=xo(1);
gamma1=xo(2);
gamma2=xo(3);
gamma3=xo(4);
gamma4=xo(5);
alpha0=xo(6);
alpha1=xo(7);
stheta=xo(8);
sxi=xo(9);

%% We will find the three cutoff values and the 2 emax functions
%  for every period and for every individual:

%Also, note that I will have to generate the shock and the wages

%Setting the seed
rng(2581633)

%Generate the shocks: (This should be done 50 times, think about how to
%better accomplish this.

SHOCKS=normrnd(0,sxi,1000,15,50);

%In this vector I will store the history of labor force participation
%For each individual, time and for each simulation. 
LABOR3=zeros(1000,15,50);



%These two functions will help us for the truncated log-normal
%expectations when double truncation or just one:


%Double truncation
TREXP = @(xiR,xiT,sxi)(normcdf(sxi-xiR/sxi)-normcdf(sxi-xiT/sxi)/(normcdf(xiT/sxi)-normcdf(xiR/sxi)));


%Single truncation
TREXP2=@(ximax,sxi)(normcdf(sxi-ximax/sxi)/(1-normcdf(ximax/sxi)));

%Predicted wages (Non-random component)
WAGESNR=gamma0+gamma1*S+gamma2*A-gamma3*A.^2-gamma4*B;

%% For the first part I just need the same as the simulations and I will 
% get simply one of the simulations out of the 50 for each individual

%% I have already stored the cutoff levels for each possible experience 
%  time and individual in the previous files. Then I don't need to re-run
%  simulations here, just to compute for the adequate level of experience
LABOR2=zeros(1000,15);

%% Begin with the loop:   

%Getting accumulated experience until period 8:

ACEXP=zeros(1000,1);
for i=1:1:1000
    dd=(i-1)*15+1;
    ACEXP(i)=sum(L(dd:dd+6));
end


%Once I get the accumulated experience, what I do is to 
%compare the shocks witht the corresponding threshold levels

LABOR3=zeros(1000,15);

for i=1:1:1000
    tt=7;
    ii=15*(i-1)+tt;
    b=ii-(tt-1);
    %Difference between initial and actual experience
    uu=ACEXP(i)+1;
    for tt=8:1:15
        if SHOCKS(i,tt,rr)>XIMAT(i,1,uu,tt) %This should come from Counterfactuals1.m
            LABOR3(i,tt)=1;
            uu=uu+1;
        end
    end
end

AVEXP=mean(LABOR3,1);

%I fill the rest of LABOR3 vector with the participation predicted in the
%simulated data:

for i=1:1:1000
    for t=1:1:7
        ii=15*(i-1)+t;
        LABOR3(i,t)=LABOR(i,t);
    end
end

AVEXP=mean(LABOR3,1);



%Predicted wages
%Total accumulated experience:

EXP=zeros(15000,1);
PW= zeros(15000,1);
for i=1:1:1000
    for t=1:1:15
         ii=15*(i-1)+t;
         b=ii-(t-1);
         EXP(ii)=A(b)+LABOR(i,t);
    end
end


%Number of periods worked in the 50 simulations
WO=sum(LABOR,3);

%Predicted wage in each observation
PREDW=zeros(15000,1);
for i=1:1:1000
    for t=1:1:15
        %For each period-individual a vector of predicted wage
        PW=zeros(50,1);
        
        %Indicator of location for time-individual
        ii=15*(i-1)+t;
        
        %Defining non random component of wages
        WNRT = (gamma0+gamma1*S(ii)+gamma2*EXP(ii)-gamma3*EXP(ii)^2-gamma4*B(ii));
        
        for rr=1:1:50
            %Predicted wage including shocks-Taking into account only
            %workers
            PW(rr)=exp(WNRT+SHOCKS(i,t,rr))*(LABOR(i,t,rr)==1);
        end
        PREDW(ii)=sum(PW)/WO(i,t);
    end
end



PREDWB=zeros(15,1);

for i=8:1:15
   PREDWB(i,1)=nanmean(PREDW(AGE==i));
end