% In this model, labor supply is introduced. 


%In this model, a stochasitic death is considered.
%No private annuity is allowed. 
% For a model that allows annuity, please look at tran70_stochas_annuity
% This program uses find_ss_ke1 and super_ss_ke6
% In this program, probability distribution of death can be different
% before and after. 
clear all
close all
clear all

global time_pref alpha beta theta t_1 t_2 t_3  coef tech precision gamma x_age TT_1 TT_2 ...
   TT_3 deprec  accuracy1 accuracy2
    


t_1=3
t_2=6
t_3=t_1+t_2
x_age=t_1
TT_1=200

SS_T=160

TT_2=200;
% TT_1 is the period to calculate the transitional path
%TT_2 is the perios to calculate the steady state situation

TT_3=200
% TT_3 is the period from which  the capital stock per capita is being the steady 
% state in the transition. Thus, ke(TT_3)=af_ke . 


deprec=0.047
time_pref=0.019
%time_pref=-0.0109
beta=1/(1+time_pref);
alpha=0.99
gamma=2%zeta=0.0195
%zeta=0.0194
%zeta=0.0194
zeta=0


tech=0.949



%alpha is the parameter of the cobb-dougls utility function when
%leisure(labo supply) is endogenous
%Note that we assume that at the age i (1<=i<t_1)
% U_i=(1/(1-gamma)) *[(c^alpha)*(1-lbsply)^(1-alpha)]^(1-gamma)
% beta is the discount factor. I use the same notation as Storeletten for
% notational simplicity
theta=0.4
% theta is the capital share in the cobb-douglas production function
eta=0.01
eta_bf=eta

dump=1

% eta_int(t) is the populationo growth rate of age 1 at period t. That is 
% eta_int(t)=tmp_pop(1,t)/tmp_pop(1,t-1) for 2<= t<= SS_T. 
% eta_int must be calculated from the data. For the moment, I put some
% temporary number. 
for t=2:SS_T
    eta_int(t)=eta
end




eta_af=eta

dynamic_inefficiency_border(1:TT_1)=eta+zeta+deprec;

precision=-11
%precision is criteria used in the main program. 
accuracy1=-11
accuracy2=-11
%accuracy1 is the criteria used in find_ss_ke1. 
rep_ratio_bf=0.35
rep_ratio_af=0.35
tau_k_bf=0
tau_k_af=0

%1/(1+time_pref) is the discount rate,i.e. u(c_1)+(1/(1+time_pref))u(c_2)+...
%theta is the income share of capital
%eta is the population growth rate
% gamma is the parameter of the CRRA, i.e, (1/(1-gamma))*C^(1-gamma)
%tech is technological parameter. 
% rep_rate_bf and rep_rate_af are the replacement of the social security. 
% tau_k_af and tau_k_bf are the capital income tax rate

%zeta is consumption growth rate and also income per capita growth rate. 

% x_age is the age that an agent receives the unintended bequest.
% It should be greater than 1 and less than or equal to t_1


%TT_1 is the very large number within which the economy will reach the steady-state.
% t_1 is the lengh of working and t_2 is the lengh of retired period
%SS_T is the year that intrinsic economic data such as death probability rearches the
%steady state
% i is tte index of age
% j is the index of cohort
% t is the index of the period. 
% deprec is the depreciation rate
% l_bound_r is the lower bound of interest rate to avoid negative interest
% rate during the computation. (not  at the equilibrium)

initial_debt_ratio_bf(1)=0.0;
for t=2:TT_1
debt_ratio_bf(t)=0;
end
for t=2:TT_1
    debt_ratio_af(t)=0;
end
% initial_debt_rate is the total debt level divided by total saving at
% period 0. 
%initial_debt_rate_af is determined inside of the model
% debt_ratio_bf(t) and debt_ratio_af(t) are the share of the value of governmet
% debt in the total saving(t-1). That is
% debt_ratio(t)=(total_debt(t)/ total_saving(t-1)

initial_debt_ratio_af=debt_ratio_bf(TT_1);

for t=1:TT_1
gov_exp_bf(t)=0.0;
end
for t=1:TT_1
gov_exp_af(t)=0.0;
end



% gov_exp is the government expenditure per labor income. That is
% gov_exp(t)=total_government_expenditure(t)/(w(t)*LF(t))


shock_age1=t_1-1;
% This is the parameter to draw the graph. we draw a graph who experience a
% shock at age t_1-1. 






% prob_bf can be calculated from the data using period 0 and period 1 and
% assuming that the economy is on the steady state at period 0 and period 1. 

%prob_af can be calculated from the data using period TT_1-1 and perod TT_1
%assuming that the economy reaches new steady state at TT_1-1 and TT_1. 

% prob=csvread('prob.csv');
% 
% for i=1:t_3
%     prob_bf(i)=prob(i);
% end
% 
%  
% for i=1:t_3
%     prob_af(i)=prob(i);
% end

for i=1:t_1
    prob_bf(i)=1
end
for i=t_1+1:t_3-3
    prob_bf(i)=1
end

for i=t_3-2:t_3
    prob_bf(i)=0.2
end


prob_af=prob_bf;
% % 

% Let prob_int(j,i) be the
%probability that an agent of a cohort j be alive at the age i given that he is alive at
%age i-1 . By definition, prob_int(j,1) is always equal to 1. 


%for cohort j<=t_3-1and i<t_3-j, it must be that prob_int(j,i) should be equal to prob_bf(i). 

for j=1:t_3-1
    for i=1:t_3-j+1
    prob_int(j,i)=prob_bf(i);
    end
end
    % the following must be calculatd from the data using the data from period
    % 1, 2, 3 ....
    % For the moment I set this is equal to prob_af(i)
    for j=2:t_3-1
    for i=t_3-j+2:t_3
    prob_int(j,i)=prob_af(i);
    end
end


% note that SS_T is the perod that the economic data on survival probabiity
% reaches the steady state. This implies that SS_T must be striclty smaller
% than TT_1. 
% Then, oldest cohort that reaches SS_T period  is the cohort SS_T since
% cohort SS_T lives until the period
% SS_T. 


%The following must be calculated from the data. For the moment I set this
%equal to prob_bf. But it must be calculated from the population data. 
for j=t_3:SS_T
for i=1:t_3
    prob_int(j,i)=prob_af(i);
end
end

% The following part must be calculated from the data. 

% The cohort j=SS_T+1:SS_T+t_3-1 is the cohort where some of the lives
% overlaps with the stedy state situation since the economic probabiity data reaches the
% steady state at SS_T. For example, the cohort SS_T+t_3-2 will be born at
% SS_T-1 at age one. From the period SS_T-1 to SS_T+1( from age to age 2),  he will follow the real data .  From age 2 to age 3 and so on, he will have
% prob_af(i). 


for j=SS_T+1:SS_T+t_3-2
for i=1:t_3-(j-SS_T)
prob_int(j,i)=prob_af(i);
end



for i=t_3-(j-SS_T)+1:t_3
prob_int(j,i)=prob_af(i);
end
end


%The following must be consistent with prob_af. Thus it is set equal to
%prob_af
for j=SS_T+t_3-1:TT_1+t_3-1
for i=1:t_3
    prob_int(j,i)=prob_af(i);
end
end










%For the cohort j<t_3, those probability is not observed from the data.
%Thus, it is reasonable to assume that before 1, it is steady state. 

%Let initial_prob(i) be the probability that can be calculated from the
%data in period 1 and period 2. We can assume that this probability
%distribution govern the cohort j<t_3. 

%The following must be calcuated from the data at perid 1 and period 2.



% Also define, for i, acc_prob(i)=p(1)*p(2)*p(3)...*p(i). That is, acc_prob(i) is the
% probability that an agent will be alive at age i. 

for j=1:TT_1+t_3-1
for i=1:t_3
    accm_prob(j,i)=prod(prob_int(j,1:i));
end
end

% From here, I calculate the coefficient of the uility function
coef(1:TT_1+t_3-1,1:t_3)=zeros(TT_1+t_3-1,t_3);
for j=1:t_3-1
    for i=t_3-j+1:t_3
   coef(j,i)=beta^(i-1)*accm_prob(j,i);
    end
end


for j=t_3:TT_1+t_3-1
    for i=1:t_3
    coef(j,i)=beta^(i-1)*accm_prob(j,i);
    end
end

% From here, probability path of alive is calculated. This path is used for
% the consumption function and income function. Prob_path2(j,i) is the
% probablity that the cohort of j is alive at the age i.  
prob_path2=zeros(TT_1+t_3-1,t_3);
for j=t_3:TT_1+t_3-1
prob_path2(j,1:t_3)=accm_prob(j,1:t_3);
end

for j=1:t_3-1
prob_path2(j,t_3-j+1:t_3)=accm_prob(j,t_3-j+1:t_3)./accm_prob(j,t_3-j+1) ;
end

        
        













% Do not change the above value. This must be small. If it is too small, interest rate becomes too high.
% and there is zero division which stops the program. 
stepsize=0.05;
scale=0.5;
% stepsize is the distance of each grid search. 
%scale is the parameter. It should be greater than 0 and striclty less than
%1. 


 bf_saving0=ones(1,t_3)*2;
 binding_saving=ones(1,TT_1+t_3-1)*(t_3);
 bf_ke0=3;
 
 for j=1:t_3-1
    start_age(j)=t_3-j+1;
end

for j=t_3:TT_1+t_3-1
    start_age(j)=1;
end

for j=t_2+1:TT_1+t_3-1
    re_age(j)=t_1+1;
end


 



[bf_ke1 bf_saving1 bf_conpath initial_debt_ratio1 all_lbsply1 re_age1 binding_saving1]= find_ss_ke1( bf_saving0,bf_ke0, ...
    rep_ratio_bf,eta,zeta,initial_debt_ratio_bf,debt_ratio_bf,...
    gov_exp_bf, tau_k_bf,prob_bf, re_age, binding_saving, stepsize,scale);


 %binding2(1:TT_2+t_3-1)=binding1(1:TT_2+t_3-1)
 %binding2(TT_2+t_3:TT_1+t_3-1)=ones(1,TT_1-TT_2)*binding1(TT_2+t_3-1);
 %bf_ke_path=bf_ke1*ones(1,TT_1);
  %[bf_ke bf_saving bf_conpath initial_debt_ratio1,binding]=super_ss_ke6( bf_saving1,bf_ke_path,  rep_ratio_bf,eta,zeta,initial_debt_ratio_bf,debt_ratio_bf, gov_exp_bf,tau_k_bf,prob_bf,binding1 )
bf_ke=bf_ke1;
af_ke=bf_ke1;
af_saving=bf_saving1;
bf_saving=bf_saving1;
binding_saving1(1:TT_1+t_3-1)=binding_saving1(TT_1+t_3-1)*ones(1,TT_1+t_3-1);

for j=1:TT_1+t_3-1
    if binding_saving1(j)<=start_age(j)
        binding_saving1(j)=start_age(j)
    end
end

binding_saving=binding_saving1
  %[af_ke af_saving bf_conpath initial_debt_ratio1 binding_af]= find_ss_ke1( bf_saving0,bf_ke0,  rep_ratio_af,eta,zeta,initial_debt_ratio_af,debt_ratio_af, gov_exp_af, tau_k_af,prob_af, binding, stepsize,scale)


  
  
  % bf_ke_path=bf_ke1*ones(1,TT_1);
%   [bf_ke1 bf_saving bf_conpath initial_debt_ratio1,binding]=super_ss_ke6( bf_saving,bf_ke_path,  rep_ratio_bf,adj_eta,zero_zeta,initial_debt_ratio_bf,debt_ratio_bf, gov_exp_bf,tau_k_bf,prob_bf,binding1 )

%  bf_ke2=bf_ke1*ones(1,TT_1);
 %[bf_ke3 bf_saving2 bf_conpath initial_debt_ratio1,binding3]=super_ss_ke6( bf_saving,bf_ke_path,  rep_ratio_bf,eta_bf,initial_debt_ratio_bf,debt_ratio_bf, gov_exp_bf,tau_k_bf,prob_bf,binding)
% 
%  
 
 
 %[af_ke af_saving af_conpath initial_debt_ratio_af binding4]= find_ss_ke1( bf_saving, bf_ke,  rep_ratio_af,eta_af, initial_debt_ratio1, debt_ratio_af, gov_exp_af,tau_k_af,prob_af,binding1,stepsize,scale)




initial_debt_ratio=0;

% period is a function that calculat the period from the cohort index and
% age index


for j=1:TT_1+t_3-1
  for  i=1:t_3
    period(j,i)=j+i-t_3;
end
end



for t=1:TT_1
   for i=1:t_3
      chrt(i,t)=(t_3+1-i)+(t-1);
    end
end


    

% base_pop(j) is the population of cohort j at age 1. 
for j=1:t_3
    base_pop(j)=(1+eta_bf)^(j-1);
end
for j=t_3+1:SS_T+t_3-1
    base_pop(j)=base_pop(j-1)*(1+eta_int(j-t_3+1));
end

for j=SS_T+t_3:TT_1+t_3-1
    base_pop(j)=base_pop(j-1)*(1+ eta_af);
end


for t=1:TT_1
    for i=1:t_3
        pop(i,t)=base_pop(chrt(i,t)) *accm_prob(chrt(i,t),i) ;
    end 
end



 

for t=1:TT_1
sum_pop(t)=sum(pop(1:t_3,t));
end
% period is a function to calculate the period when age i and cohort j are
% given. 

num_accd_death=zeros(t_3,TT_1);
for t=1:TT_1
  for i=2:t_3
      num_accd_death(i,t)= base_pop(chrt(i,t) )*accm_prob(chrt(i,t),i-1)*(1-prob_int(chrt(i,t),i));
  end
      num_accd_death(1,t)=0;
end
%num_accd_death(i,t) is the number of death that people lives until i-1 th age and dies at age i in period t.   


for i=1:t_3-1
pop_0(i)=pop(i+1,1)/prob_bf(i+1);
end
pop_0(t_3)=(pop(1,1)/(1+eta_bf))*accm_prob(1,t_3) ;

%pop_0 is the hypothetical population at the period 0 that is consisttent with
%the population at the perido 1 and the prob(j,i). 




for t=1:TT_1+t_3-1
    agg_effcy(t)=(1+zeta)^(t-1);
end


  

  %Define the labor supply vector here. 
  % all_lbsply(t,i)  is the labor supply of individual at period t and at the age i. 
  %  lbsply_chrt(j,i) is the labor supply of cohort j at age i.  
  %We assume that the time endowment is one
  %The following number is temporary.  
  
  
  
  
%   
%   lbsply_chrt(t_3:TT_1+t_3-1,1:t_1)=0.5*ones(TT_1,t_1);
%   lbsply_chrt(t_3:TT_1+t_3-1,t_1+1:t_3)=zeros(TT_1,t_2);
%   
%   lsr_chrt(t_3:TT_1+t_3-1,1:t_3)=ones(TT_1,t_3)-lbsply_chrt(t_3:TT_1+t_3-1,1:t_3);
%   
%   
%   for j=t_2+1:t_3-1
%     for i=start_age(j):t_1
%     lbsply_chrt(j,i)=0.5;
%     end
%     for i=t_1+1:t_3
%         lbsply_chrt(j,i)=0;
%     end
%     for i=start_age(j):t_3
%         lsr_chrt(j,i)=1-lbsply_chrt(j,i);
%     end
% 
%   end
%    
% 



for t=1:TT_1
 all_lbsply(t,1:t_1)=all_lbsply1(TT_1,1:t_1);
end

for t=1:TT_1
 all_eff_lbsply(t,1:t_1)=all_lbsply(TT_1,1:t_1).(*hcpf(1:t_1)');
end



for j=1:t_2
for i=start_age(j):t_3
    lsr_chrt(j,i)=1;
end
end

for j=t_2+1:t_3-1
    for i=start_age(j):t_1
lbsply_chrt(j,i)=all_lbsply(period(j,i),i);
lsr_chrt(j,i)=1-lbsply_chrt(j,i);
    end
    for i=t_1+1:t_3
        lbsply_chrt(j,i)=0;
        lsr_chrt(j,i)=1;
    end
end
  
for j=t_3:TT_1
   for i=1:t_1
    lbsply_chrt(j,i)=all_lbsply(period(j,i),i);
    lsr_chrt(j,i)=1-lbsply_chrt(j,i);
   end
   for i=t_1+1:t_3
    lbsply_chrt(j,i)=0;
    lsr_chrt(j,i)=1;
end
end

for j=TT_1+1:TT_1+t_3-1
    lbsply_chrt(j,1:t_3)=lbsply_chrt(TT_1,1:t_3);
    lsr_chrt(j,1:t_3)=lsr_chrt(TT_1,1:t_3);
end



  for t=1:TT_1
      ELF(t)=sum( agg_effcy(t)* (all_eff_lbsply(t,1:t_1))*pop(1:t_1,t) );
  end
  for t=1:TT_1
      LF(t)=sum( pop(1:t_1,t) );
  end
  
for t=1:TT_1

    all_saving(t, 1:t_3)=((1+zeta)^(t))*af_saving(1:t_3);
end
% The above all_saving is initial value for finding the steady state.
% Thefore it is needed to be to use af_saving. 


K1(1,1)=0;
for i=2:t_3
   K1(i,1)=bf_saving(i-1);
end
% The above must be bf_saving because it is pre-deermined. 


KST1(1)=(1-initial_debt_ratio)*( (pop_0(1:t_3-1) )*K1(2:t_3,1));



for t=2:TT_1
    KST1(t)=(1-debt_ratio_af(t))*(all_saving(t-1,1:t_3))*pop(1:t_3,t-1);   
end


ke(1:TT_1)=ones(1,TT_1)*bf_ke;

k0(1:TT_1)=ones(1,TT_1);


 %This is the initial value for calculation. 

 initial_debt_ratio1=0;
counter=0
 
 
 
  while max(abs(k0-ke))+max(abs(binding_saving-binding_saving1))>10^(precision)
 %bf_saving=bf_saving1;
 binding_saving=binding_saving1;
  [max_dif_ke_k0 indx]=max(abs(k0-ke))
 dif_k0_ke_at_t_70=abs(k0(70)-ke(70)) ;
   dif_k0_ke_at_t_50=abs(k0(50)-ke(50)); 
   dif_k0_ke_at_t_2=abs(k0(2)-ke(2))  ;
 
 %k0(2:TT_3-1)=k0(2:TT_3-1)+dump*( ke(2:TT_3-1)-k0(2:TT_3-1));
k0(1:TT_3)=k0(1:TT_3)+dump*(ke(1:TT_3)-k0(1:TT_3)); 
ke(TT_3+1:TT_1)=ke(TT_3);
k0(TT_3+1:TT_1)=ke(TT_3);
%pop(i,t) is the population of  age i at time t.
%LF(t) is the labor force at period t.


%K1 is the mana of the economy




   for t=1:TT_3   
  
   r(t)=(1-tau_k_af)*( theta*tech*((k0(t))^(theta-1))-deprec );
   efw(t)=(1-theta)*tech*((k0(t))^theta);
   w(t)=efw(t)*agg_effcy(t);
   %efw(t) is the wage rate per efficiency unit person
   % w(t) is the wage rate that the person receives for those whose human
   % capital is equal to one. 
end
for t=TT_3:TT_1+t_3-1
    r(t)=(1-tau_k_af)*( theta*tech*(k0(TT_3))^(theta-1)-deprec );
    w(t)=agg_effcy(TT_3)*((1-theta)*tech*(k0(TT_3)^theta))*(1+zeta)^(t-TT_1);
 end

 
compndr(1)=1+r(1);
for t=2:TT_1+t_3-1
   compndr(t)=compndr(t-1)*(1+r(t));
end

%from here, I define the pricepath2(j,1:t_3), which is the consumer price
%of cohort j , to calculate the consumption path

pricepath2=zeros(TT_1+t_3-1,t_3);
for j=t_3:TT_1+t_3-1
start_y=j-t_3+1;
pricepath2(j,1:t_3)=compndr(start_y)./compndr(start_y:start_y+t_3 -1);
end

for j=1:t_3-1
   pricepath2(j,t_3-j+1:t_3)=compndr(1)./compndr(1:j);
end
    
   
% from here, I define the wagepath2, which is the wage path of cohort j. 



wagepath2=zeros(TT_1+t_3-1,t_3);
for j=t_3:TT_1+t_3-1
wagepath2(j,1:t_1)=w(1,j-t_3+1:j-t_3+t_1).*hcpf(1:t_1);
wagepath2(j,t_1+1:t_3)=zeros(1,t_2);




end

if (t_2+1)<=( t_3-1)
    for j=t_2+1:t_3-1
      wagepath2(j,t_1-j+t_2+1:t_1)=w(1,1:j-t_2).*hcpf(start_age(j):t_1);
      % for a cohort j which is <t_3, he starts to work at the age
      % t_1-j+t_2+1. As a result, the length of his work is j-t_2. 
      
    end
end
% In the case that t_2+1>t_3-1, then  1=t_1.thus worker works only one period.
% In such a case, a cohort before t_3 do not work at all. Thus, the wage
% profile can be defined only from t_3 work. 








%from here, I calculate the newly issued bond
debt_ratio_af(1)=initial_debt_ratio1;

new_bond_ratio(1)=(debt_ratio_af(2)*( (all_saving(1,1:t_3))*pop(1:t_3,1)      )-debt_ratio_af(1)* ( bf_saving(1:t_3-1)*pop_0(1:t_3-1 )'  ))/LF(1);

for t=2:TT_1-1
    new_bond_ratio(t)=(debt_ratio_af(t+1)*( (all_saving(t,1:t_3))*pop(1:t_3,t)        )-debt_ratio_af(t)*(  (all_saving(t-1,1:t_3))*pop(1:t_3,t-1)   ) )/LF(t);
end

for t=TT_1:TT_1+t_3-1
    new_bond_ratio(t)=new_bond_ratio(TT_1-1)*((1+zeta)^(TT_1-t));
end

%new bond ratio is newly issued bond per efficiency unit person






















for t = 1:TT_1 
benefit2( t ) = rep_ratio_af * w( t ) *(all_lbsply(t,1:t_1)*pop(1:t_1,t) )/LF(t);
end
%replacement ratio is defined as labor income per person of the young at
%period t

 for t=TT_1+1:TT_1+t_3-1
     benefit2(t)=benefit2(TT_1)*(1+zeta)^(t-TT_1);
 end
benefits_path=zeros(TT_1+t_3-1,t_3);
for  j=1:t_2
    benefits_path(j,t_3-j+1:t_3)=benefit2(1:j);
    
end
for j= t_2+1:t_3 -1
    
    benefits_path(j,t_1+1:t_3)=benefit2(j-t_2 +1:j) ;    
    
end
for j=t_3:TT_1+t_3-1
    benefits_path(j,t_1+1:t_3)=benefit2(j-t_2+1:j) ;
end    
%for  j=TT_1+1:TT_1 +t_3-1
%    benefits_path(j,t_1+1:t_3)=ones(1,t_2)*benefit2(TT_1) ;
%end

 


    sum_benefit(1)=benefit2(1)*sum(pop( t_1+1:t_3,1 ) );
    tax(1)=(sum_benefit(1)+debt_ratio_af(1)*( bf_saving(1:t_3-1)*pop_0(1:t_3-1)' )*r(1)-r(1)*KST1(1)*tau_k_af/(1-tau_k_af)-new_bond_ratio(1)*LF(1)+gov_exp_af(1)*LF(1)*w(1) )/((all_lbsply(1,1:t_1)*pop(1:t_1,1))*w(1)    );
    net_tax(1)=1-tax(1);




for t=2:TT_1
    sum_benefit(t)=benefit2(t)*sum(pop( t_1+1:t_3,t ) );
    tax(t)=(sum_benefit(t)+debt_ratio_af(t)*(  (all_saving(t-1,1:t_3))*pop(1:t_3,t-1)    )*r(t)-r(t)*KST1(t)*tau_k_af/(1-tau_k_af)-new_bond_ratio(t)*LF(t)+gov_exp_af(t)*LF(t)*w(t) )/((all_lbsply(t,1:t_1)*pop(1:t_1,t))*w(t) );
    net_tax(t)=1-tax(t);
end 

for t=1:TT_1
    if isnan(net_tax(t))==1
        net_tax(t)=1;
    end
end

for t=1:TT_1
    if net_tax(t)<=0
        net_tax(t)=0.2;
    end
end
%This part is  necessary to prevent the after tax
%wage rate to be negative durig calculationg. (not at the equilibrium).




%Note that debt_ratio_af(t)*LF(t)*w(t) comes from the total debt level is
%debt_ratio_af*w*LF/r . Thus, interest payment due to the debt is
% debt_ratio_af*w*LF

 for t=TT_1+1:TT_1+t_3-1
    net_tax(t)=net_tax(TT_1);
end


 %from here,I define the net_tax path
 
 net_taxpath=zeros(TT_1+t_3-1,t_3);
for j=t_3:TT_1+t_3-1
net_taxpath(j,1:t_1)=net_tax(j-t_3+1:j-t_3+t_1);
end

for j=t_2+1:t_3-1
      net_taxpath(j,t_1-j+t_2+1:t_1)=net_tax(1:j-t_2);
      
end
 

for j=1:TT_1+t_3-1
 net_wagepath(j,1:t_3)=net_taxpath(j,1:t_3).*wagepath2(j,1:t_3);
end
 
 
 
 
 
 
 

 

 
 
% From here, gift as unintended bequest is calculated.



gift(1)=(  (bf_saving(1:t_3-1))*num_accd_death(2:t_3,1)*(1+r(1)) )/ ( pop(x_age,1));



for t=2:TT_1
    gift(t)= ( (all_saving(t-1,1:t_3-1))* num_accd_death(2:t_3,t)*(1+r(t)) )/(pop(x_age,t));
end

for t=TT_1+1:TT_1+t_3-1
    gift(t)=gift(TT_1)*((1+zeta)^(t-TT_1));
end

gift_path=zeros(TT_1+t_3-1,t_3);
for j=1:t_3-1
    if x_age<t_3-j+1
gift_path(j,1:t_3)=zeros(1,t_3);
    elseif t_3-j+1<=x_age
        gift_path(j,x_age)=gift(period(j,x_age));
    end
end        
% note that t_3-j+1 is the starting age for cohort j=1 to t_3-1
for j=t_3:TT_1+t_3-1
    gift_path(j,x_age)=gift(period(j,x_age));
    
end









%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%Calculate the life time wealth of jth cohort up to the point where saving
%constraint becomes binding
%The following 
for j=1:t_2
bd=binding_saving(j);
 lftwlt(j)=K1(t_3+1-j,1)*(1+r(1)) +( benefits_path(j,1:bd) + gift_path(j,1:bd) )*(pricepath2(j,1:bd)');
end
% In the above, K1(i,1) is the asset of age i at the period 1.


  for j=t_2+1:t_3-1
      bd=binding_saving(j);
   lftwlt(j)=  K1(t_3+1-j,1)*(1+r(1))+(net_wagepath(j,1:bd))*(pricepath2(j,1:bd)') +...
          ( benefits_path(j,1:bd) + gift_path(j,1:bd)) *(pricepath2(j,1:bd)');
   end
% Remember that wagepth(j,w) is the wage path of jth cohort

for j=t_3:TT_1+t_3-1
    bd=binding_saving(j);
   lftwlt(j)=( net_wagepath(j,1:bd) )*(pricepath2(j,1:bd)')+ ...
       ( benefits_path(j,1:bd)+ gift_path(j,1:bd)) *(pricepath2(j,1:bd)');

end





%From here, the program calculates the consumption function under CRRA. 
%For derivation, see my handout. 

% we consider the problem in two steps optization. The first step is to
% allocate the expenditure(sum of consumption plus leisure) for each
% period. Then allocate pure consumption and leisure. Note that once
% leisure is determined, then labor supply is automatically determined. 

% indlut(j,i) is the indirect utility of the cobb-douglas
% utility(consumption and leisure) of the cohort j at the age i.

% re_age(j) is the retirement age of cohort j. 
% start_age(j) is the starting age of cohort j. For j>=t_3 start_age(j)=1. 
% for the cohort j<=t_3-1, start_age(j)=t_3-j+1

ind_coef=( (alpha^alpha)*((1-alpha)^(1-alpha))  )^(1-gamma);
alpha_gamma=1-alpha*(1-gamma);
coefvector_alive=zeros(TT_1+t_3-1,t_3);

for j=t_3:TT_1+t_3-1
    for i=1:t_1
    wagehat(j,i)=net_wagepath(j,i)^((1-alpha)*(gamma-1));
    end
end
for j=t_2+1:t_3-1
    for i=start_age(j):t_1
    wagehat(j,i)=net_wagepath(j,i)^((1-alpha)*(gamma-1));
    end
end

for j=t_3:TT_1+t_3-1
    bd=binding_saving(j);
 if re_age(j)==t_1+1     
    for i=1:t_1
    con_coef(j,i)=( (coef(j,i)*wagehat(j,i)*pricepath2(j,1))/(coef(j,1)*wagehat(j,1)*pricepath2(j,i)) )^(1/gamma) ;
    end
    
    for i=t_1+1:t_3
        con_coef(j,i)=( (coef(j,i)*alpha*pricepath2(j,1))/(coef(j,1)*ind_coef*wagehat(j,1)*pricepath2(j,i)) )^(1/alpha_gamma); 
    end
% from here, calculate the coeffcient of equation
    

    % calculate the polynominal equation coefficient
    
    poly(j,1)=sum(con_coef(j,1:t_1)*   (pricepath2(j,1:t_1)')  );
    poly(j,2)=sum(con_coef(j,t_1+1:bd)*(pricepath2(j,t_1+1:bd)')  );
    poly(j,3)=-lftwlt(j);
   my_fun =@(m_1) m_1*poly(j,1)+poly(j,2)*m_1^(gamma/alpha_gamma) +poly(j,3);

   m_1=fzero(my_fun,[0 lftwlt(j)]);
   
   for i=1:t_1
       exp_chrt_alive(j,i)=m_1*con_coef(j,i);
       c_chrt_alive(j,i)=exp_chrt_alive(j,i)*alpha;
       lsr_chrt(j,i)=exp_chrt_alive(j,i)*(1-alpha)/net_wagepath(j,i);
       lbsply_chrt(j,i)=1-lsr_chrt(j,i);
   end
   if bd==t_3
   for i=t_1+1:t_3
       exp_chrt_alive(j,i)=(m_1^(gamma/alpha_gamma)) *con_coef(j,i);
       c_chrt_alive(j,i)=exp_chrt_alive(j,i);
       lsr_chrt(j,i)=1;
   end
   elseif bd<=t_3-1 & bd>=t_1+1
   for i=t_1+1:bd
       exp_chrt_alive(j,i)=(m_1^(gamma/alpha_gamma)) *con_coef(j,i);
       c_chrt_alive(j,i)=exp_chrt_alive(j,i);
       lsr_chrt(j,i)=1;
   end
   for i=bd+1:t_3
       exp_chrt_alive(j,i)=benefits_path(j,i);
       c_chrt_alive(j,i)=exp_chrt_alive(j,i);
       lsr_chrt(j,i)=1;
   end
   end
       
       
       
 elseif re_age(j)<=t_1 
for i=1:re_age(j)-1
        con_coef(j,i)=( (coef(j,i)*wagehat(j,i)*pricepath2(j,1))/(coef(j,1)*wagehat(j,1)*pricepath2(j,i)) )^(1/gamma) ; ;
    end
    
    for i=re_age(j):t_1
 con_coef(j,i)=( (coef(j,i)*alpha*pricepath2(j,1))/(coef(j,1)*ind_coef*wagehat(j,1)*pricepath2(j,i)) )^(1/alpha_gamma); 
    end
    for i=t_1+1:t_3
   con_coef(j,i)=( (coef(j,i)*alpha*pricepath2(j,1))/(coef(j,1)*ind_coef*wagehat(j,1)*pricepath2(j,i)) )^(1/alpha_gamma);  
    end
    
% from here, calculate the coeffcient of equation


    % calculate the polynominal equation coefficient
    poly(j,1)=sum(con_coef(j,1:re_age(j)-1)*   (pricepath2(j,1:re_age(j)-1)')  );
    poly(j,2)=sum(con_coef(j,re_age(j):bd)*(pricepath2(j,re_age(j):bd)')  );
    poly(j,3)=net_wagepath(j,re_age(j):t_1)*(pricepath2(j,re_age(j):t_1)')-lftwlt(j);
   my_fun =@(m_1) m_1*poly(j,1)+poly(j,2)*m_1^(gamma/alpha_gamma) +poly(j,3);

   m_1=fzero(my_fun,[0 lftwlt(j)]);

    
    
    
    
    

   

    for i=1:re_age(j)-1
       exp_chrt_alive(j,i)=m_1*con_coef(j,i);
       c_chrt_alive(j,i)=alpha*exp_chrt_alive(j,i);
       lsr_chrt(j,i)=(1-alpha)*exp_chrt_alive(j,i)/net_wagepath(j,i);
       lbsply_chrt(j,i)=1-lsr_chrt(j,i);
   end
   for i=re_age(j):t_1
       exp_chrt_alive(j,i)=net_wagepath(j,i)+(m_1^(gamma/alpha_gamma)) *con_coef(j,i);
       c_chrt_alive(j,i)=exp_chrt_alive(j,i)-net_wagepath(j,i); 
       lsr_chrt(j,i)=1;
        lbsply_chrt(j,i)=0;
   end
   if bd==t_3
   for i=t_1+1:t_3
       exp_chrt_alive(j,i)=(m_1^(gamma/alpha_gamma)) *con_coef(j,i);
       c_chrt_alive(j,i)=exp_chrt_alive(j,i);
       lsr_chrt(j,i)=1;
   end
   elseif bd<=t_3-1
       for i=t_1+1:bd
       exp_chrt_alive(j,i)=(m_1^(gamma/alpha_gamma)) *con_coef(j,i);
       c_chrt_alive(j,i)=exp_chrt_alive(j,i);
       lsr_chrt(j,i)=1;
       end
       for i=bd+1:t_3
       exp_chrt_alive(j,i)=benefits_path(j,i);
       c_chrt_alive(j,i)=exp_chrt_alive(j,i);
       lsr_chrt(j,i)=1;
       end
   end
       


 end
end

% from here calculate the expenditure and labor supply for cohort
% 1<=j<=t_3-1

for j=1:t_2
bd=binding_saving(j);

    for i=start_age(j):t_3

con_coef(j,i)=(coef(j,i)/pricepath2(j,i));
    end
    for i=start_age(j):t_3
        eqn_coef(j,i)=pricepath2(j,i)*( (con_coef(j,i)/con_coef(j,start_age(j)))^(1/alpha_gamma));
    end
    
    m_1=lftwlt(j)/sum(eqn_coef(j,start_age(j):bd));
    if bd==t_3
    for i=start_age(j):t_3
        exp_chrt_alive(j,i)=m_1*(con_coef(j,i)/con_coef(j,start_age(j)))^(1/alpha_gamma);
        c_chrt_alive(j,i)=exp_chrt_alive(j,i);
        lsr_chrt(j,i)=1;
    end
    elseif bd<=t_3-1 
        for i=start_age(j):bd
         exp_chrt_alive(j,i)=m_1*(con_coef(j,i)/con_coef(j,start_age(j)))^(1/alpha_gamma);
         c_chrt_alive(j,i)=exp_chrt_alive(j,i);   
         lsr_chrt(j,i)=1;
        end
        for i=bd+1:t_3
            exp_chrt_alive(j,i)=benefits_path(j,i);
           c_chrt_alive(j,i)=exp_chrt_alive(j,i); 
           lsr_chrt(j,i)=1;
        end
    end
end


% from here, I consider the cohort j=t_2+1:t_3-1
for j=t_2+1:t_3-1
    bd=binding_saving(j);
  % This is the cohort that some peeriod they can work. 
  if re_age(j)>=start_age(j)+1 
      % This is the case when the agent works  least at the first age. 
    for i=start_age(j):re_age(j)-1
        con_coef(j,i)=( (coef(j,i)*wagehat(j,i)*pricepath2(j,start_age(j)))/(coef(j,start_age(j))*wagehat(j,start_age(j))*pricepath2(j,i)) )^(1/gamma) ; ;

    end
    
    
        for i=re_age(j):t_1
 con_coef(j,i)=( (coef(j,i)*alpha*pricepath2(j,start_age(j)))/(coef(j,start_age(j))*ind_coef*wagehat(j,start_age(j))*pricepath2(j,i)) )^(1/alpha_gamma); 
    end
    for i=t_1+1:t_3
   con_coef(j,i)=( (coef(j,i)*alpha*pricepath2(j,start_age(j)))/(coef(j,start_age(j))*ind_coef*wagehat(j,start_age(j))*pricepath2(j,i)) )^(1/alpha_gamma);  
    end
    
        
 
 

    % calculate the polynominal equation coefficient
    
    poly(j,1)=sum(con_coef(j,start_age(j):re_age(j)-1)*   (pricepath2(j,start_age(j):re_age(j)-1)')  );
    poly(j,2)=sum(con_coef(j,re_age(j):bd)*(pricepath2(j,re_age(j):bd)')  );
    poly(j,3)=net_wagepath(j,re_age(j):t_1)*(pricepath2(j,re_age(j):t_1)')-lftwlt(j);
   my_fun =@(m_1) m_1*poly(j,1)+poly(j,2)*m_1^(gamma/alpha_gamma) +poly(j,3);

    

   my_fun =@(m_1) m_1*poly(j,1)+poly(j,2)*m_1^(gamma/alpha_gamma)+poly(j,3);
   
      m_1=fzero(my_fun,[0 lftwlt(j)]);
    for i=start_age(j):re_age(j)-1
       exp_chrt_alive(j,i)=m_1*con_coef(j,i);
       c_chrt_alive(j,i)=alpha*exp_chrt_alive(j,i);
       lsr_chrt(j,i)=(1-alpha)*exp_chrt_alive(j,i)/net_wagepath(j,i);
       lbsply_chrt(j,i)=1-lsr_chrt(j,i);
   end
   for i=re_age(j):t_1
       exp_chrt_alive(j,i)=net_wagepath(j,i)+(m_1^(gamma/alpha_gamma)) *con_coef(j,i);
       c_chrt_alive(j,i)=exp_chrt_alive(j,i)-net_wagepath(j,i); 
       lsr_chrt(j,i)=1;
        lbsply_chrt(j,i)=0;
   end
    if bd==t_3
   for i=t_1+1:t_3
       exp_chrt_alive(j,i)=(m_1^(gamma/alpha_gamma)) *con_coef(j,i);
       c_chrt_alive(j,i)=exp_chrt_alive(j,i);
       lsr_chrt(j,i)=1;
       
   end
    elseif bd<=t_3-1
        for i=t_1+1:bd
       exp_chrt_alive(j,i)=(m_1^(gamma/alpha_gamma)) *con_coef(j,i);
       c_chrt_alive(j,i)=exp_chrt_alive(j,i);
       lsr_chrt(j,i)=1;
        end
        for i=bd+1:t_3
        exp_chrt_alive(j,i)=benefits_path(j,i);
       c_chrt_alive(j,i)=exp_chrt_alive(j,i);    
       lsr_chrt(j,i)=1;
        end
    end
    
   
 elseif re_age(j)==start_age(j)
    for i=start_age(j):t_3

        con_coef(j,i)=(coef(j,i)/pricepath2(j,i));
    end
    for i=start_age(j):t_3
        eqn_coef(j,i)=pricepath2(j,i)*( (con_coef(j,i)/con_coef(j,start_age(j)))^(1/alpha_gamma));
    end
    
    m_1_wage= (lftwlt(j)-net_wagepath(j,start_age(j):t_1)*(pricepath2(j,start_age(j):t_1)'))/sum(eqn_coef(j,start_age(j):bd));
% m_1_wage is m_x -w_x where x is the start_age
    for i=start_age(j):t_1
        exp_chrt_alive(j,i)=net_wagepath(j,start_age(j))+m_1_wage*(con_coef(j,i)/con_coef(j,start_age(j)))^(1/alpha_gamma);
        c_chrt_alive(j,i)=exp_chrt_alive(j,i)-net_wagepath(j,i);
        lsr_chrt(j,i)=1
    end
    if bd==t_3
    for i=t_1+1:t_3
        exp_chrt_alive(j,i)=m_1_wage*(con_coef(j,i)/con_coef(j,start_age(j)))^(1/alpha_gamma);
        c_chrt_alive(j,i)=exp_chrt_alive(j,i);
        lsr_chrt(j,i)=1;
    end
    elseif bd<=t_3-1
    for i=t_1+1:bd
        exp_chrt_alive(j,i)=m_1_wage*(con_coef(j,i)/con_coef(j,start_age(j)))^(1/alpha_gamma);
        c_chrt_alive(j,i)=exp_chrt_alive(j,i);
        lsr_chrt(j,i)=1;
    end    
        for i=bd+1:t_3
        exp_chrt_alive(j,i)=benefits_path(j,i);
        c_chrt_alive(j,i)=exp_chrt_alive(j,i);
        lsr_chrt(j,i)=1
        end
    end
end
end




   
   
   
   
% From here, I consider the non-negative constraint of labor supply. 
% change_re_age(1:TT_1+t_3-1)=zeros(1,TT_1+t_3-1);
% for j=t_2+1:TT_1+t_3-1
%    for i=1:re_age(j)-1
%        if lbsply_chrt(j,i)<- 10^(accuracy2) 
%            lbsply_chrt(j,i)=0;
%            if re_age(j)>=start_age(j)+1
%            re_age(j)=re_age(j)-1;
%            change_re_age(j)=1;
%            end
%        end
%     end
% end
           
           


% From here, I consider the non-working condition of leisure.
% When re_age <t_1+1 and net_wage*MU_con -MU_leisure>0,  , then they should decrease the consumtion of lesiure
% and should deter the retirement. 

% for j=t_2+1:TT_1+t_3-1
% if re_age(j)<=t_1 & change_re_age(j)==0;
%     
%   if  net_wagepath(j,re_age(j))*alpha*exp_chrt_alive(j,re_age(j))^(alpha-1)...
%         -(exp_chrt_alive(j,re_age(j))^alpha)*(1-alpha)>10^accuracy2
%  re_age(j)=re_age(j)+1;
%   end
% end
% end
% 
%  
 
 


    


%saving block




saving_chrt(1,1:t_3)=zeros(1,t_3);






% 
% 
% for j=2:TT_1+t_3-1
%     bd=binding(j);
%     if bd==t_3
%     exp_chrt_alive(j,1:t_3)=coefvector_alive(j,1:t_3).*lambda_hat(j);
%     end
%     if bd<=t_3-1 & t_1+1<=bd+1
%        
%     exp_chrt_alive(j,1:bd)=coefvector_alive(j,1:bd).*lambda_hat(j);
%     for i=bd+1:t_3
%     exp_chrt_alive(j,i)=benefit2(period(j,i));
%     end
% 
%     elseif bd+1<=t_1 & bd<t_3-1
%        exp_chrt_alive(j,1:bd)=coefvector_alive(j,1:bd).*lambda_hat(j);  
%         
%         for i=bd+1:t_1
%         exp_chrt_alive(j,i)=net_tax(period(j,i))*w(period(j,i));
%         end
%         for i=t_1+1:t_3
%             exp_chrt_alive(j,i)=benefit2(period(j,i));
%         end
%     end 
% end
% 
% 




for j=1:t_2
    %this is the cohort who do not work at all
    for i=start_age(j):t_3
      if i==start_age(j) 
    dispin_chrt(j,i)=K1(i,1)*(1+r(1))+benefits_path(j,i)+gift_path(j,i); 
    saving_chrt(j,i)=dispin_chrt(j,i)-c_chrt_alive(j,i);   
        
      elseif i>=start_age(j)+1
        
 dispin_chrt(j,i)=saving_chrt(j,i-1)*(1+r(period(j,i)))+benefits_path(j,i)+gift_path(j,i);
 saving_chrt(j,i)=dispin_chrt(j,i)-c_chrt_alive(j,i);   

 
    end        
    end
end


for j=t_2+1:t_3-1
    % this is the cohort where at some period they can work. 
for i=start_age(j):t_1
    if i==start_age(j)
     dispin_chrt(j,i)=K1(i,1)*(1+r(1))+net_wagepath(j,i)*lbsply_chrt(j,i)+gift_path(j,i);
     saving_chrt(j,i)=dispin_chrt(j,i)-c_chrt_alive(j,i);   
    elseif i>=start_age(j)+1
               dispin_chrt(j,i)=saving_chrt(j,i-1)*(1+r(period(j,i)))+net_wagepath(j,i)*lbsply_chrt(j,i)+gift_path(j,i);
               saving_chrt(j,i)=dispin_chrt(j,i)-c_chrt_alive(j,i);   
        
    end
end
for i= t_1+1:t_3
    
   
   dispin_chrt(j,i)=saving_chrt(j,i-1)*(1+r(period(j,i)))+benefits_path(j,i);
   saving_chrt(j,i)=dispin_chrt(j,i)-c_chrt_alive(j,i);
end
end


for j=t_3:TT_1+t_3-1
% this is the cohort who live the whole life
    
    for i=1:1
            dispin_chrt(j,i)=net_wagepath(j,i)*lbsply_chrt(j,i)+gift_path(j,i);        
            saving_chrt(j,i)=dispin_chrt(j,i)-c_chrt_alive(j,i);
        end
            for i=2:t_1
                dispin_chrt(j,i)=saving_chrt(j,i-1)*(1+r( period(j,i) ) )...
                    +net_wagepath(j,i)*lbsply_chrt(j,i)+gift_path(j,i);
                saving_chrt(j,i)=dispin_chrt(j,i)-c_chrt_alive(j,i);
            end
            for i=t_1+1:t_3
                dispin_chrt(j,i)=saving_chrt(j,i-1)*(1+ r(period(j,i)))+benefits_path(j,i);
                saving_chrt(j,i)=dispin_chrt(j,i)-c_chrt_alive(j,i);
            
            end
end

    
   
    
    
    
    
 
    
 change_binding_saving(1:TT_1+t_3-1)=zeros(1,TT_1+t_3-1);    
 for j=2:(t_3-1)
     for i=t_3-j+1:t_3
 if   ( saving_chrt(j,i)<=0  & t_1+1<=i)
% the above condition implies that the binding constraint move left if the
% zero saving age is lower than binding_saving(j) and it is greater or
% equal to t_1+1
     saving_chrt(j,i:t_3)=zeros(1,t_3-i+1);
     if i<binding_saving(j)
     binding_saving(j)=i;
change_binding_saving(j)=1;
     end
 end
    end
end

for j=t_3:TT_1+t_3-1
    for i=1:t_3
        if (saving_chrt(j,i)<= 0 &  t_1+1<=i)
           saving_chrt(j,i:t_3)=zeros(1,t_3-i+1);
          if i<binding_saving(j) 
              binding_saving(j)=i;
            change_binding_saving(j)=1;
          end
        end
    end
end





% From here, I consider the non-negative constraint of the saving.
% MRS> 1+interest rate
if max(change_binding_saving)==0
for j=2:TT_1+t_3-1
       
    
    if binding_saving(j)<=t_3-1 & start_age(j)<=binding_saving(j) 
            bd=binding_saving(j);
           if ( -alpha*(c_chrt_alive(j,bd)^(-alpha_gamma))* lsr_chrt(j,bd)^((1-alpha)*(1-gamma)) ...
                +(1+r(period(j,bd+1)) )*alpha*(c_chrt_alive(j,bd+1)^(-alpha_gamma))* lsr_chrt(j,bd+1)^((1-alpha)*(1-gamma)) >10^accuracy2 & ...
               saving_chrt(j,bd)==0 )
               binding_saving(j)=binding_saving(j)+1;
           else
            binding_saving(j)=binding_saving(j);
            
           end
    end
end    

end



    
    
    
for t=1:TT_1
for    i=1:t_3
all_saving(t,i)=saving_chrt(chrt(i,t), i)  ;
end
end




for t=1:TT_1
    for i=1:t_3
 all_exp(t,i)=exp_chrt_alive(chrt(i,t),i);
 

 % Note  that an agent who is i years old at period t will dies at the
 % t+t_3-i. 
 % This means that this agent is t+t_3-i cohort. 
    end
end


for t=1:TT_1      
 for i=1:t_3
     if isnan(all_saving(t,i))==1
         all_saving(t,i)=((1+zeta)^t)*af_saving(i);
     end
 end
end


    




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
%The following code is very important
 %K1(2:t_3)=all_saving(1:t_3-1,TT_1);
%bf_saving(1:t_3-1)=all_saving(1:t_3-1,TT_1);

% KST1(1)=(1-initial_debt_ratio)*( (pop(2:t_3,1)')./prob(2:t_3) )*K1(2:t_3,1);


for t=2:TT_1
    KST1(t)=(1-debt_ratio_af(t))*(all_saving(t-1,1:t_3)*pop(1:t_3,t-1));
  
    
end

% trasnforming labor supply from cohort-based to period based
for t=1:TT_1
    for i=1:t_1

 all_lbsply(t,i)=lbsply_chrt(chrt(i,t),i);
    end
end
for t=1:TT_1
 all_eff_lbsply(t,1:t_1)=all_lbsply(TT_1,1:t_1).(*hcpf(1:t_1)');
end
  


  
  for t=1:TT_1
      ELF(t)=sum( agg_effcy(t)* (all_eff_lbsply(t,1:t_1))*pop(1:t_1,t) );
  end







% Define ke here







   ke=KST1./ELF;
   
   for t=1:TT_1
       if isnan(ke(t))==1
           ke(t)=0.1;
       end
   end
  
   
   
   if max(ke)>500
   for t=1:TT_1
      if ke(t)>500;
         ke(t)=500;
      end
  end
  
   end
   
   
   if min(ke)<=0
   for t=1:TT_1
      if ke(t)<=0;
         ke(t)=0.01;
      end
  end
  
   end
   
   
   
   
counter=counter+1



binding_saving1=binding_saving;

end %This is the end of while loop.


all_saving(TT_1-t_3, 1:t_3-1)./all_saving(TT_1-t_3-1,1:t_3-1)

 exp_chrt_alive(TT_1-10,1:t_3)./exp_chrt_alive(TT_1-11,1:t_3)

%Pure consumption block. I use the word pure consumption to distinguish the
%expenditure which is the sum of the pure consumption and the consumption
%of leisure

   


    

for t=1:TT_1
    for i=1:t_3
all_c(t,i)=c_chrt_alive(chrt(i,t),i);       
        
    end
end















%Calculate the life time wealth of jth cohort
%The following 
for j=1:t_2

 total_lftwlt(j)=K1(t_3+1-j,1)*(1+r(1)) +( benefits_path(j,1:t_3) + gift_path(j,1:t_3) )*(pricepath2(j,1:t_3)');
end
% In the above, K1(i,1) is the asset of age i at the period 1.


  for j=t_2+1:t_3-1
   total_lftwlt(j)=K1(t_3+1-j,1)*(1+r(1))+(net_taxpath(j,1:t_3).*wagepath2(j,1:t_3))*(pricepath2(j,1:t_3)') +...
          ( benefits_path(j,1:t_3) + gift_path(j,1:t_3)) *(pricepath2(j,1:t_3)');
   end
% Remember that wagepth2(j,i) is the wage path of jth cohort

for j=t_3:TT_1+t_3-1
   total_lftwlt(j)=( net_taxpath(j,1:t_3).*wagepath2(j,1:t_3) )*(pricepath2(j,1:t_3)')+ ...
       ( benefits_path(j,1:t_3)+ gift_path(j,1:t_3)) *(pricepath2(j,1:t_3)');

end





%The following code will check Walras law. 
for t=1:TT_1-1
   Y(t)=tech*(KST1(t)^theta)*ELF(t)^(1-theta);
   Exp(t)=all_c(t,1:t_3)*pop(1:t_3,t)+KST1(t+1)+gov_exp_af(t)*w(t)*LF(t);
   Exp2(t)=all_c(t,1:t_3)*pop(1:t_3,t)+KST1(t+1)+gov_exp_af(t)*w(t)*LF(t);
   Inc(t)=(all_lbsply(t,1:t_1)*pop(1:t_1,t))*w(t)+KST1(t)*(1+r(t)/(1-tau_k_af));
   
end



        

for j=t_3:TT_1+t_3-1
start_y=j-t_3+1;
pricepath_base(j,1:t_3)=compndr(start_y)./compndr(start_y:start_y+t_3 -1);
end

for j=1:t_3-1
   pricepath_base(j,1:j) =compndr(1)./compndr(1:j);
end



for j=t_3:TT_1+t_3-1
exp_price_alive(j,1:t_3)=pricepath_base(j,1:t_3);
end


for j=1:t_3-1
    exp_price_alive(j,1:j)=pricepath_base(j,1:j);
end









for j=t_3:TT_1
   
   expend_alive(j)=exp_chrt_alive(j,1:t_3)*exp_price_alive(j,1:t_3)';
   total_exp(j)=exp_chrt_alive(j,1:t_3)*(pricepath2(j,1:t_3)');
end




for j=1:t_3-1
    
   expend_alive(j)=exp_chrt_alive(j,t_3-j+1:t_3)*exp_price_alive(j,1:j)';

   total_exp(j)=exp_chrt_alive(j,t_3-j+1:t_3)*(pricepath2(j,t_3-j+1:t_3)');
end


   


 
for j=1:TT_1   
Error40(j)=total_lftwlt(j)-total_exp(j);


end




 


                                                                                                                                                                                

time3=1:length(Error40);
%Y(t) is the total output at time t.
%Exp(t) is the sum of consumption and saving at time t.
%Inc(t) is the total income at time t. 

for t=1:TT_1-1
Error1(t)=Y(t)+(1-deprec)*KST1(t)-Exp(t);
Error2(t)=Inc(t)-Exp(t);
Error3(t)=Y(t)+(1-deprec)*KST1(t)-Inc(t);


%Error 3 is based on eulers theorem.

end
Error5=all_saving(1:TT_1,t_3)';
%consumption path of jth cohort at age i
time1=1:TT_1;
time2=1:TT_1-1;

size_bf_chrt=20;
all_ke= [bf_ke1*ones(1,size_bf_chrt) ke];
length_allke=length(all_ke);


for t=1:length_allke
all_wage(t)= (1-theta)*tech*((all_ke(t))^theta);
end

for t=1:length_allke
 all_interest_rate(t)=(1-tau_k_af)*(theta*tech*((all_ke(t))^(theta-1))-deprec);
end

for t=1:length_allke
trans_time(t)=-size_bf_chrt+ t;
end


for j=t_3:TT_1+t_3-1
    all_conpath(j,1:t_3)=c_chrt_alive(j,1:t_3);
end

for j=1:TT_1+t_3-1
utility_int(j)=lftutility(all_conpath(j,1:t_3),j);
end



trans_utility=[ utility_int(1:TT_1+t_3-1)];


%expenditure_int=expd_fun_int( tau_k_bf, bf_ke, utility_int);



%trans_expenditure=[expenditure_int];

% for j=t_3:TT_1+t_3-1
% dif_lftwlt_expd(j)=lftwlt(j)-expenditure_int(j);
% end

for t=1:length(trans_utility)
    cohort(t)=-size_bf_chrt+t;
end

for j=1:t_3-1
cohort_saving(j,1:t_3-j)=bf_saving(1:t_3-j);
for i=t_3-j+1:t_3
    cohort_saving(j,i)=all_saving(period(j,i),i);
end
end

for j=t_3:TT_1
    for i=1:t_3
    cohort_saving(j,i)=all_saving(period(j,i),i);
    end
end

af_saving=((1+zeta)^(-TT_1))*all_saving(TT_1,1:t_3);

for t=1:TT_1-1
saving_growth(t,1:t_3-1)=all_saving(t+1,1:t_3-1)./all_saving(t,1:t_3-1);
end



for t=1:TT_1
MPK(t)= theta*tech*((ke(t))^(theta-1));
end

for t=1:TT_1-1
output_capital_ratio(t)=KST1(t)/Y(t);
end


figure(1)
plot(trans_time,all_ke,'b-');
title('Transitional Path of K/L ') ;
xlabel('Time Period');
ylabel('K/L');

figure(2)
plot(trans_time,all_wage,'b*-');
title('Transitional Path of wage rate ') ;
xlabel('Time Period');
ylabel('wage rate');

figure(3)
plot(trans_time,all_interest_rate,'b*-');
title('Transitional Path of interest rate ') ;
xlabel('Time Period');
ylabel('interest rate');


figure(4)


plot(cohort,trans_utility,'b*-');
title('Transitional Path of life time utility') ;
xlabel('cohort');
ylabel('life-time utility');



% figure(5)
% 
% 
% plot(cohort,trans_expenditure,'b*-');
% title('Transitional Path of expenditure function') ;
% xlabel('cohort');
% ylabel('expenditure given the utility');

% figure(6)
% 
% 
% plot(1:length(dif_lftwlt_expd),dif_lftwlt_expd,'b*-');
% title('Checking expenditure minus lifetime income') ;
% xlabel('cohort minus t_3 +1');
% ylabel('difference of expendituer minus lftwlt');






figure(7)
plot(1:length(Error1),Error1,'r-.');
title('Checking Walras Law1: Output+capital -exp');
xlabel('Time Period');
ylabel('Error1');

figure(8)
plot(time2,Error2,'b-.');
title('Checking Walras Law2');
xlabel('Time Period');
ylabel('Error2');

figure(9)
plot(time2,Error3,'b-');
title('Eulers theorem ');
xlabel('Time Period');
ylabel('Error3');

figure(10)
plot(time3,Error40,'b-');
title('Checking income equal to expenditure condition');
xlabel('cohort');
ylabel('Error40');


figure(11)
plot(time1,Error5,'b-');
title('Checking the balance of saving at the death');
xlabel('cohort');
ylabel('Error5');

figure(12)
plot(1:t_3,bf_saving,'b-*', 1:t_3,af_saving,'r-+');
xlabel('age');
title('saving path at the SS: Before the shock (blue *) and after  the schok(adjusted) (red +)  ') ;
ylabel('saving balance');





figure(13)
plot(1:TT_1-1,saving_growth( 1:TT_1-1,t_1),'b-');
xlabel('time');
title('saving growth ') ;
ylabel('saving growth rate');


figure(14)
plot(1:t_3,all_conpath(t_3-shock_age1,1:t_3),'b*-');
xlabel('age');
title('life-time consumption path who experiece a shock at the age t_1-1  ') ;
ylabel('consumption ');

figure(15)
plot(1:t_3,all_conpath(TT_1+t_3-1,1:t_3),'b*-');
xlabel('age');
title('life-time consumption path at the ex-post steady state  ') ;
ylabel('consumption ');



figure(16)
plot(1:TT_1,MPK,'b-',1:TT_1,dynamic_inefficiency_border(1:TT_1), 'r-');
xlabel('time');
title(' Checking dynamic inefficiency: MPK(blue) and golden rule(red)');
ylabel('Marginal product of capital and golden rule ');

figure(17)
plot(1:TT_1-1,output_capital_ratio,'b-');
xlabel('time');
title('capital output ratio: at the begining it should be around 3.2 for the US ');
ylabel('capital output ratio ');


figure(18)
plot(t_3:TT_1+t_3-1,lbsply_chrt(t_3:TT_1+t_3-1,t_1),'b-');
xlabel('cohort');
title('Labor Supply at the t_1 th Age of Each Cohort ');
ylabel('Labor supply ');

figure(19)
plot(1:t_3,lsr_chrt(TT_1+t_3-1,1:t_3),'b*-');
xlabel('age');
title('Consumpiton of lesisure  ') ;
ylabel('Leisure ');