LifeCycleModel16.m

%% Life-Cycle Model 16 Consumption and Borrowing Constraints 2
% Households would like to consumption smooth
% Exogenous shocks, if enough bad shocks in a row occour, can be another
% reason for households to run up against borrowing constraints (we already
% saw one reason in Life-Cycle Model 15).

% We use Life-Cycle Model 9, where exogenous shocks were an idiosyncratic
% AR(1) process on labor productivity units. We will use just 5 points to
% discretize the AR(1) process, which makes a poor approximation but makes
% it easier for us to look at how a series of bad shocks can lead to the
% borrowing constraint binding.

% We cannot use life-cycle profiles to look at role of the borrowing 
% constraint as these pick up what happens to households in
% general, while we are interested in households that have a series of bad
% shocks and as a result run up against the borrowing constraint. The
% 'mean' of households will not be running into borrowing constraints due
% to shocks.

% We will instead simulate a large panel data set and then finding some
% households in that with a series of bad shocks. You could just
% simulate a single household given a specific series of shocks and while
% this is not difficult to code nor is it a standard command.

% From the panel we then look for households that get good shocks for the
% first fifteen periods and mostly bad shocks for the next ten periods (see
% code below for exact definition of good/bad). In figure 3 we plot the
% time series for ten such households.

% In the bottom panel of Figure 3 you can see how most of our candidate households 
% manage to save some small amount of assets in the first fifteen periods (note 
% that the mean life-cycle profile is also low), and then during the bad shocks between
% periods 16 to 25 their assets go to zero and the borrowing constraints
% become binding (actually plot next period assets as it is these to which 
% the borrowing constraint applies).

% The top and middle panels of Figure 3 show that when the bad shocks hit during periods 
% 16-25 notice how the consumption  of our candidate households falls below the mean 
% life-cycle profile, while the marginal utility of consumption increases above the mean
% life-cycle profile. This is the impact of the borrowing constraints binding.

% Note: LifeCycleModel13 explained simulating panel data

% Note: I use just 5 states for the exogenous shock z, this is not accurate
% for an AR(1), but means it is easier define how to find simulations with runs of good vs bad
% shocks.

%% How does VFI Toolkit think about this?
%
% One decision variable: h, labour hours worked
% One endogenous state variable: a, assets (total household savings)
% One stochastic exogenous state variable: z, an AR(1) process (in logs), idiosyncratic shock to labor productivity units
% Age: j

%% Begin setting up to use VFI Toolkit to solve
% Lets model agents from age 20 to age 100, so 81 periods

Params.agejshifter=19; % Age 20 minus one. Makes keeping track of actual age easy in terms of model age
Params.J=100-Params.agejshifter; % =81, Number of period in life-cycle

% Grid sizes to use
n_d=51; % Endogenous labour choice (fraction of time worked)
n_a=201; % Endogenous asset holdings
n_z=5; % Exogenous labor productivity units shock
N_j=Params.J; % Number of periods in finite horizon

%% Parameters

% Discount rate
Params.beta = 0.96;
% Preferences
Params.sigma = 2; % Coeff of relative risk aversion (curvature of consumption)
Params.eta = 1.5; % Curvature of leisure (This will end up being 1/Frisch elasty)
Params.psi = 10; % Weight on leisure

% Prices
Params.w=1; % Wage
Params.r=0.05; % Interest rate (0.05 is 5%)

% Demographics
Params.agej=1:1:Params.J; % Is a vector of all the agej: 1,2,3,...,J
Params.Jr=46;

% Pensions
Params.pension=0.8;

% Age-dependent labor productivity units
Params.kappa_j=[linspace(1.5,2,Params.Jr-15),linspace(2,1,14),zeros(1,Params.J-Params.Jr+1)];
% Exogenous shock process: AR1 on labor productivity units
Params.rho_z=0.9;
Params.sigma_epsilon_z=0.03;

% Conditional survival probabilities: sj is the probability of surviving to be age j+1, given alive at age j
% Most countries have calculations of these (as they are used by the government departments that oversee pensions)
% In fact I will here get data on the conditional death probabilities, and then survival is just 1-death.
% Here I just use them for the US, taken from "National Vital Statistics Report, volume 58, number 10, March 2010."
% I took them from first column (qx) of Table 1 (Total Population)
% Conditional death probabilities
Params.dj=[0.006879, 0.000463, 0.000307, 0.000220, 0.000184, 0.000172, 0.000160, 0.000149, 0.000133, 0.000114, 0.000100, 0.000105, 0.000143, 0.000221, 0.000329, 0.000449, 0.000563, 0.000667, 0.000753, 0.000823,...
    0.000894, 0.000962, 0.001005, 0.001016, 0.001003, 0.000983, 0.000967, 0.000960, 0.000970, 0.000994, 0.001027, 0.001065, 0.001115, 0.001154, 0.001209, 0.001271, 0.001351, 0.001460, 0.001603, 0.001769, 0.001943, 0.002120, 0.002311, 0.002520, 0.002747, 0.002989, 0.003242, 0.003512, 0.003803, 0.004118, 0.004464, 0.004837, 0.005217, 0.005591, 0.005963, 0.006346, 0.006768, 0.007261, 0.007866, 0.008596, 0.009473, 0.010450, 0.011456, 0.012407, 0.013320, 0.014299, 0.015323,...
    0.016558, 0.018029, 0.019723, 0.021607, 0.023723, 0.026143, 0.028892, 0.031988, 0.035476, 0.039238, 0.043382, 0.047941, 0.052953, 0.058457, 0.064494,...
    0.071107, 0.078342, 0.086244, 0.094861, 0.104242, 0.114432, 0.125479, 0.137427, 0.150317, 0.164187, 0.179066, 0.194979, 0.211941, 0.229957, 0.249020, 0.269112, 0.290198, 0.312231, 1.000000]; 
% dj covers Ages 0 to 100
Params.sj=1-Params.dj(21:101); % Conditional survival probabilities
Params.sj(end)=0; % In the present model the last period (j=J) value of sj is actually irrelevant

% Warm glow of bequest
Params.warmglow1=0.3; % (relative) importance of bequests
Params.warmglow2=3; % bliss point of bequests (essentially, the target amount)
Params.warmglow3=Params.sigma; % By using the same curvature as the utility of consumption it makes it much easier to guess appropraite parameter values for the warm glow

%% Grids
% The ^3 means that there are more points near 0 and near 10. We know from
% theory that the value function will be more 'curved' near zero assets,
% and putting more points near curvature (where the derivative changes the most) increases accuracy of results.
a_grid=10*(linspace(0,1,n_a).^3)'; % The ^3 means most points are near zero, which is where the derivative of the value fn changes most.

% First, the AR(1) process z
[z_grid,pi_z]=discretizeAR1_FarmerToda(0,Params.rho_z,Params.sigma_epsilon_z,n_z);
z_grid=exp(z_grid); % Take exponential of the grid
[mean_z,~,~,statdist_z]=MarkovChainMoments(z_grid,pi_z); % Calculate the mean of the grid so as can normalise it
z_grid=z_grid./mean_z; % Normalise the grid on z (so that the mean of z is exactly 1)

% Grid for labour choice
h_grid=linspace(0,1,n_d)'; % Notice that it is imposing the 0<=h<=1 condition implicitly
% Switch into toolkit notation
d_grid=h_grid;

%% Now, create the return function 
DiscountFactorParamNames={'beta','sj'};

% Notice we still use 'LifeCycleModel8_ReturnFn'
ReturnFn=@(h,aprime,a,z,w,sigma,psi,eta,agej,Jr,pension,r,kappa_j,warmglow1,warmglow2,warmglow3,beta,sj) LifeCycleModel8_ReturnFn(h,aprime,a,z,w,sigma,psi,eta,agej,Jr,pension,r,kappa_j,warmglow1,warmglow2,warmglow3,beta,sj)

%% Now solve the value function iteration problem, just to check that things are working before we go to General Equilbrium
disp('Test ValueFnIter')
vfoptions=struct(); % Just using the defaults.
tic;
[V, Policy]=ValueFnIter_Case1_FHorz(n_d,n_a,n_z,N_j, d_grid, a_grid, z_grid, pi_z, ReturnFn, Params, DiscountFactorParamNames, [], vfoptions);
toc

% Because we just want to simulate panel data we have no need for the
% stationary distribution. We will need to create a distribution from which
% the initial characteristics of households are drawn for our panel data.
% We will make them all households of agej=1 (age 20 years old), and with
% zero assets. Note that this is just jequaloneDist, the same as previously.

%% Initial distribution of agents at birth (j=1)
% Before we plot the life-cycle profiles we have to define how agents are
% at age j=1. We will give them all zero assets. We use the stationary
% distribution of z.
jequaloneDist=zeros([n_a,n_z],'gpuArray'); % Put no households anywhere on grid
jequaloneDist(1,:)=statdist_z; % All agents start with zero assets, and the median shock

%% We now compute the 'stationary distribution' of households
% Start with a mass of one at initial age, use the conditional survival
% probabilities sj to calculate the mass of those who survive to next
% period, repeat. Once done for all ages, normalize to one
Params.mewj=ones(1,Params.J); % Marginal distribution of households over age
for jj=2:length(Params.mewj)
    Params.mewj(jj)=Params.sj(jj-1)*Params.mewj(jj-1);
end
Params.mewj=Params.mewj./sum(Params.mewj); % Normalize to one
AgeWeightsParamNames={'mewj'}; % So VFI Toolkit knows which parameter is the mass of agents of each age
simoptions=struct(); % Use the default options
StationaryDist=StationaryDist_FHorz_Case1(jequaloneDist,AgeWeightsParamNames,Policy,n_d,n_a,n_z,N_j,pi_z,Params,simoptions);
% Again, we will explain in a later model what the stationary distribution
% is, it is not important for our current goal of graphing the life-cycle profile

%% FnsToEvaluate are how we say what we want to graph the life-cycles of
% Like with return function, we have to include (h,aprime,a,z) as first inputs, then just any relevant parameters.
FnsToEvaluate.fractiontimeworked=@(h,aprime,a,z) h; % h is fraction of time worked
FnsToEvaluate.earnings=@(h,aprime,a,z,w,kappa_j) w*kappa_j*z*h; % w*kappa_j*z*h is the labor earnings (note: h will be zero when z is zero, so could just use w*kappa_j*h)
FnsToEvaluate.assets=@(h,aprime,a,z) a; % a is the current asset holdings
FnsToEvaluate.consumption=@(h,aprime,a,z,agej,Jr,w,kappa_j,r,pension) (agej<Jr)*(w*kappa_j*z*h+(1+r)*a-aprime)+(agej>=Jr)*(pension+(1+r)*a-aprime);
FnsToEvaluate.marginalutilityofcons=@(h,aprime,a,z,agej,Jr,w,kappa_j,r,pension,sigma) ((agej<Jr)*(w*kappa_j*z*h+(1+r)*a-aprime)+(agej>=Jr)*(pension+(1+r)*a-aprime))^(-sigma); % u(c)=(c^(1-sigma))/(1-sigma), therefore u'(c)=c^(-sigma); note that we are using a seperable utility fn
FnsToEvaluate.z=@(h,aprime,a,z) z; % Need this to identify the households with a long run of 'bad shocks'
FnsToEvaluate.totalsavings=@(h,aprime,a,z) aprime; % aprime>=0 is the borrowing constraint

%% Calculate the life-cycle profiles
AgeConditionalStats=LifeCycleProfiles_FHorz_Case1(StationaryDist,Policy,FnsToEvaluate,Params,[],n_d,n_a,n_z,N_j,d_grid,a_grid,z_grid,simoptions);
% Plot the life cycle profiles
figure(1)
subplot(3,2,1); plot(1:1:Params.J,AgeConditionalStats.fractiontimeworked.Mean)
title('Life Cycle Profile: Fraction Time Worked (h)')
subplot(3,2,2); plot(1:1:Params.J,AgeConditionalStats.earnings.Mean)
title('Life Cycle Profile: Labor Earnings (w kappa_j h)')
subplot(3,2,3); plot(1:1:Params.J,AgeConditionalStats.assets.Mean)
title('Life Cycle Profile: Assets (a)')
subplot(3,2,4); plot(1:1:Params.J,AgeConditionalStats.consumption.Mean)
title('Life Cycle Profile: Consumption (c)')
subplot(3,2,5); plot(1:1:Params.J,AgeConditionalStats.marginalutilityofcons.Mean)
title('Life Cycle Profile: Marginal Utility of Consumption (u''(c))')
subplot(3,2,6); plot(1:1:Params.J,AgeConditionalStats.z.Mean)
title('Life Cycle Profile: AR(1) (in logs) on labor productivity units (z)')


%% Simulate panel data

% To simulate panel data we will set the number of time periods
simoptions.simperiods=N_j; % N_j is the default value
simoptions.numbersims=10^4; % 10^3 is the default value
% To simulate panel data you have to define 'where' an individual household
% simulation starts from, we will use the jequaloneDist (from which
% starting points will be drawn randomly)
InitialDist=jequaloneDist;

SimPanelValues=SimPanelValues_FHorz_Case1(InitialDist,Policy,FnsToEvaluate,Params,[],n_d,n_a,n_z,N_j,d_grid,a_grid,z_grid,pi_z, simoptions);
% Simulates a panel based on PolicyIndexes of 'numbersims' agents of length

% For example
% SimPanelValues.earnings
% is the simulate panel of FnsToEvaluate.earnings
% size(SimPanelValues.earnings) is [simperiods,numbersims] 
% (what econometric theory on panel data would typically call T-by-N)

%% Find households with and good 'start' and then lots of 'bad shocks'

% I want a household that does not have a bad shock in first 15 periods
nobadshocksatfirst=(prod(SimPanelValues.z(1:15,:)>z_grid(2),1)==1);
[sum(nobadshocksatfirst==0),  sum(nobadshocksatfirst==1)] % Second of these is number of households satisfying this criterion
% Note: look for households that have median or better shock in all of first 15 periods

% I want a household that has seven or more bad shocks in next 10 periods
sevenormorebadinnext10=(sum(SimPanelValues.z(16:25,:)==z_grid(1),1)>=7);
[sum(sevenormorebadinnext10==0),  sum(sevenormorebadinnext10==1)] % Second of these is number of households satisfying this criterion
% Note: look for households that have the worst shock in at 7 or more of the next 10 periods

% Both at once
likely=(nobadshocksatfirst==0).*(sevenormorebadinnext10==1);
sum(likely) % Couple of hundred canditate households

% The indexes for those which satisfy both criteria
temp=1:1:simoptions.numbersims;
likelyindex=temp(logical(likely));
% Note: length(likelyindex) is same as sum(likely), so appear to work as intended
% Below I arbitrarily use the first ten of these candidates

%% Plot that compares the consumption and marginal utility of consumption for our households likely to hit
% borrowing constraints against the mean life-cycle profiles. 

% Plot just the first 30 periods, as this it is this part we have chosen our candidates from.
% Remember we choose candidates who have good shocks in first 15 periods,
% then mostly bad shocks for periods 16-25, after that we make no selection.

% Arbitrarily plot the first ten of the candidate households
figure(3)
subplot(3,1,1);  plot(1:1:30,AgeConditionalStats.consumption.Mean(1:30))
hold on 
plot(1:1:30,SimPanelValues.consumption(1:30,likelyindex(1:10)))
hold off
legend('Mean')
title('Consumption (c)')
subplot(3,1,2);  plot(1:1:30,AgeConditionalStats.marginalutilityofcons.Mean(1:30))
hold on 
plot(1:1:30,SimPanelValues.marginalutilityofcons(1:30,likelyindex(1:10)))
hold off
legend('Mean')
title('Marginal Utility of Consumption (u''(c))')
subplot(3,1,3);  plot(1:1:30,AgeConditionalStats.totalsavings.Mean(1:30))
hold on 
plot(1:1:30,SimPanelValues.totalsavings(1:30,likelyindex(1:10)))
hold off
legend('Mean')
title('Total Savings (next period assets, aprime)')
% You can see how most of our candidate households manage to save some
% small amount of assets in the first fifteen periods (note that the mean
% life-cycle profile is also low), and then during the bad shocks between
% periods 16 to 25 their assets go to zero and the borrowing constraints
% become binding (actually plot next period assets as it is these to which 
% the borrowing constraint applies).

% When the bad shocks hit during periods 16-25 notice how the consumption
% of our candidate households falls below the mean life-cycle profile,
% while the marginal utility of consumption increases above the mean
% life-cycle profile. This is the impact of the borrowing constraints
% binding.


%% Note: The literature on this has more commonly focused on regressions to look at consumption volatility
% among people near a budget constraint. In part this is because they want to compare findings to the 
% data and the marginal utility of consumption is of course unobservable in the data.