Modeling Uncertain Retirement Income

Hi Robert,

I am trying to model a small probability that retirement income drops drastically. A household may consider such disaster due:

i.) Death of a spouse: If one member of the couple passes away, the household’s average income drops significantly and does not recover until the death of the surviving spouse.
ii.) Health expenses due to illness or entry into a retirement home: Expenses due to illness may be temporary or permanent, while expenses related to entering a retirement home are likely to be permanent.

I aim to introduce uncertain retirement income into the portfolio model, which already works well and accounts for persistent (z1) and transitory (z2) shocks to labor income, as well as the possibility of a stock return disaster. Now, I want to consider two scenarios separately:

a.) A tail event resulting in a temporary drop in retirement income (e.g., temporary illness expenses).
b.) A tail event resulting in a permanent drop in retirement income (e.g., permanent illness expenses or/and death of a spouse or/and entry into a retirement home).

I believe I conceptually understand how to model the disastrous event in scenario a). Unfortunately, I am struggling to set up the more interesting and important concept in scenario b).

Finally, my issues are:

1.) I receive an error message when I start modeling scenario a), and I am unsure why. Specifically, I am introducing a temporary shock using e and connecting it to retirement income in the return function. After that, I planned to introduce a small probability of disaster with
p and a normal process with 1−p (similar to the stock disaster case modeled on the stochastic error u).

Setup:

n_e=3;

Params.sigma_epsilon_e=0.1

[e_grid,pi_e]=discretizeAR1_FarmerToda(0,0,Params.sigma_epsilon_e,n_e);
e_grid=exp(e_grid);
pi_e=pi_e(1,:)'; 
mean_e=pi_e'*e_grid; 
e_grid=e_grid./mean_e; 

vfoptions.n_e=n_e;
vfoptions.e_grid=e_grid;
vfoptions.pi_e=pi_e;
simoptions.n_e=vfoptions.n_e;
simoptions.e_grid=vfoptions.e_grid;
simoptions.pi_e=vfoptions.pi_e;

% + I add e in return function …

Error:

Unrecognized function or variable 'e_gridvals_J'.

Error in ValueFnIter_Case1_FHorz_RiskyAsset (line 69)
      [VKron, PolicyKron]=ValueFnIter_Case1_FHorz_RiskyAsset_noa1_e_raw(n_d,n_a2,n_z,vfoptions.n_e,n_u, N_j, d_grid, a2_grid, z_gridvals_J, e_gridvals_J, u_grid, pi_z_J, pi_e_J, pi_u, ReturnFn, aprimeFn, Parameters, DiscountFactorParamNames, ReturnFnParamNames, aprimeFnParamNames, vfoptions);
  
Error in ValueFnIter_Case1_FHorz (line 596)
[V,Policy]=ValueFnIter_Case1_FHorz_RiskyAsset(n_d,n_a1,n_a2,n_z,vfoptions.n_u, N_j, d_grid, a1_grid, a2_grid, z_gridvals_J, vfoptions.u_grid, pi_z_J, vfoptions.pi_u, ReturnFn, vfoptions.aprimeFn, Parameters, DiscountFactorParamNames, ReturnFnParamNames, vfoptions);
  
Error in LifeCycleModel31_11A_Stock_disaster_ret (line 196)
[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);

2.) Any example, advice or instructions on how to set up the code to model scenario b) would be very helpful to me.

Thank you so much!

1 Like

I will take a look at issue 1 and let you know in the next few days.

For issue 2, how to set up an uncertain permanent loss of income. The best way is probably as an exogenous markov (z) state. Set z_grid=[1;0.5]; (0.5 here is going to represent losing 50% of income). And then set pi_z=[0.9, 0.1; 0, 1]; (so there is a 10% chance you lose half your income, and this is permanent; notice that the second row shows that once you move to this second value you stay there permanently). You likely want to set up pi_z_J so that the probability of being hit by this permanent shock can vary by age (makes more sense given it is about death of a spouse or health problems). You can then just make retirement income equal to what you previously had it set at times the value of z. (You can of course do something with two different size permanent shocks with z_grid=[1;0.6;0.3] and pi_z=[0.8,0.1,0.1; 0,1,0; 0,0,1].)

2 Likes

Thank you for your help, Robert.

I will carefully go through your suggestions. I will let you know if I get stuck somewhere.

1 Like

Issue 1 (e in risky asset model). Should work now, just pushed fix.

2 Likes

Thank you very much!

I cannot work on modelling until tomorrow evening. I will start running codes again on Wednesday.

1 Like

Issue 1 (e in risky asset model): A tail event resulting in a temporary drop in retirement income (e.g., temporary illness expenses).

Thanks, now it goes through without errors. Just to check if I understand the code correctly:

a1) Setting new transitory shock e:

Params.sigma_epsilon_e=0.1

[e_grid,pi_e]=discretizeAR1_FarmerToda(0,0,Params.sigma_epsilon_e,n_e;
e_grid=exp(e_grid);
pi_e=pi_e(1,:)';  
mean_e=pi_e'*e_grid;
e_grid=e_grid./mean_e;

% + other modification of code as new process e is incorporated

This is just the usual temporary shock, which I connected to retirement income in the return function. At this point, it would be equivalent to use z2 instead of e in the return function as z2 is also a temporary iid shock with equal settings. Are these lines of code set correctly?

a2) Generating disaster within e process:

Params.sigma_epsilon_e=0.1

[e_grid,pi_e]=discretizeAR1_FarmerToda(0,0,Params.sigma_epsilon_e,n_e-1);
e_grid=exp(e_grid);
pi_e=pi_e(1,:)';
e_grid=[0.5; e_grid];
pi_e=[0.1; (1-0.1)*pi_e];
mean_e=pi_e'*e_grid;
e_grid=e_grid./mean_e;

% + other modification of code as new process e is incorporated

I added two lines and changed n_e to n_e-1 as the value function complained about an incorrect shape of e. I am not sure if these lines are actually correct and if they introduce a tail event with a 10% probability that retirement income can drop by 50% at each age. I also noticed that the life-cycle profile of assets and the share of assets invested in risky assets do not change much, which seems somewhat weird to me. Should the code be written differently?

1 Like

Issue 2 (e in risky model): A tail event resulting in a permanent drop in retirement income (e.g., permanent illness expenses or/and death of a spouse or/and entry into a retirement home)

Thanks!

Unfortunately, I am not able to implement your suggestion. In fact, I can implement it in the model where this shock is only shock. For example, if I modify model 8 by replacing the unemployment shock (z_grid=[1;0] and pi_z=[0.7, 0.3; 0.5, 0.5]) with a retirement shock (z_grid=[1;0.5] and pi_z=[0.9, 0.1; 0, 1]) and connect z with retirement income in return function, it works as expected. Conversely, in model where I have already set z1 and z2 as persistent and transitory shocks, I cannot set a retirement shock as markov shock e. When I try to set it (e_grid=[1;0.5] and pi_e=[0.9, 0.1; 0, 1]), I receive an error:

Error using ReturnFnParamNamesFn (line 56)
Cannot find the parameter e in the Parameters structure (it is needed as an input to the ReturnFn)

Error in ValueFnIter_Case1_FHorz (line 335)
    ReturnFnParamNames=ReturnFnParamNamesFn(ReturnFn,n_d,n_a,n_z,N_j,vfoptions,Parameters);

If I set parameter e = 1 then I receive error:

Test ValueFnIter
Error using CreateReturnFnMatrix_Case2_Disc_Par2 (line 132)
gpuArray/arrayfun encountered an issue while compiling 'LifeCycleModel31_11A_Stock_disaster_ReturnFn (line 1)'.

Error in ValueFnIter_Case1_FHorz_RiskyAsset_noa1_raw (line 33)
        ReturnMatrix=CreateReturnFnMatrix_Case2_Disc_Par2(ReturnFn, n_d, n_a, n_z, d_grid, a_grid, z_gridvals_J(:,:,N_j), ReturnFnParamsVec);

Error in ValueFnIter_Case1_FHorz_RiskyAsset (line 83)
            [VKron, PolicyKron]=ValueFnIter_Case1_FHorz_RiskyAsset_noa1_raw(n_d,n_a2,n_z,n_u, N_j, d_grid, a2_grid, z_gridvals_J, u_grid, pi_z_J, pi_u, ReturnFn, aprimeFn, Parameters, DiscountFactorParamNames, ReturnFnParamNames, aprimeFnParamNames, vfoptions);

Error in ValueFnIter_Case1_FHorz (line 596)
    [V,Policy]=ValueFnIter_Case1_FHorz_RiskyAsset(n_d,n_a1,n_a2,n_z,vfoptions.n_u, N_j, d_grid, a1_grid, a2_grid, z_gridvals_J, vfoptions.u_grid, pi_z_J, vfoptions.pi_u, ReturnFn, vfoptions.aprimeFn, Parameters, DiscountFactorParamNames, ReturnFnParamNames, vfoptions);

Caused by:
    LifeCycleModel31_11A_Stock_disaster_ReturnFn takes at most 11 inputs. 12 were supplied.

The way you do:

e_grid=[0.5; e_grid];
pi_e=[0.1; (1-0.1)*pi_e];

is correct to “introduce a tail event with a 10% probability that retirement income can drop by 50% at each age.”
For this to cause a 50% drop in retirement income though will obviously depend on how you use the value of e in the return fn (and elsewhere). The code so far is just setting a 10% probability of a shock with a value of 0.5, that the 0.5 corresponds to a 50% fall in retirement income depends on how you use the e shocks.

1 Like

For what you say about Issue 2 with a permanent shock. It sounds like you are not setting up multiple exogenous shocks correctly.

When you use multiple markov z shocks you need to set up n_z=[n_z1, n_z2]; z_grid=[z1_grid; z2_grid]; and pi_z=kron(pi_z2, pi_z1); (note reverse order in the kron, also the pi_z is here assuming the two shocks are uncorrelated).

When you use multiple iid e shocks you need to set up n_e=[n_e1, n_e2]; e_grid=[e1_grid; e2_grid]; pi_e=kron(pi_e2, pi_e1); (note that pi_e is now a column vector, whereas pi_z was a matrix, hence the difference; again, pi_e is here assuming the two shocks are independent/uncorrelated).

In the ReturnFn and in FnsToEvaluate you would have (…,z,e,…). If you don’t use z or e then you just omit it. If you use multiple z and multiple e you just put (…,z1,z2,e1,e2,…).

My impression of your error is that you set up two z and one e, but you didn’t modify ReturnFn and FnsToEvaluate to have inputs (…,z1,z2,e,…).

2 Likes

A quick note: if you use multiple iid variables “e”, how do you input the big probability vector pi_e?

I would do

pi_e = kron(pi_e2, pi_e1);

1 Like

@aledinola you are quite right, corrected the expression for pi_e in my above post. Thanks!

2 Likes

Thanks, I think Issue 1 works fine now for me.

At this point, I model a 10% possibility that retirement income drops by half in each age of retirement. Otherwise (with the remaining 90% possibility), it follows an iid process with sigma_epsilon_e = 0.2. More specifically, I write the critical parts of the code below. Please let me know if you see anything weird.

n_e=5;
Params.rho_e=0;
Params.sigma_epsilon_e=0.2;

[e_grid,pi_e]=discretizeAR1_FarmerToda(0,0,Params.sigma_epsilon_e,n_e-1);
e_grid=exp(e_grid);
pi_e=pi_e(1,:)';
e_grid=[0.5; e_grid];
pi_e=[0.1; (1-0.1)*pi_e];
mean_e=pi_e'*e_grid;
e_grid=e_grid./mean_e;

vfoptions.n_e=n_e;
vfoptions.e_grid=e_grid;
vfoptions.pi_e=pi_e;
simoptions.n_e=vfoptions.n_e;
simoptions.e_grid=vfoptions.e_grid;
simoptions.pi_e=vfoptions.pi_e;

% + I add e in returm function after z1, z2

% ***

% Return function

if agej<Jr % If working age
c=w*kappa_j*z1*z2+a-savings; % Add z here
else % Retirement
c=pension*kappa_j*e+a-savings;

% ***

% Plotting

% 1)

eind=floor(n_e+1)/2;

figure(1)

subplot(2,1,1); surf(a_grid*ones(1,Params.J),ones(n_a,1)*(1:1:Params.J),reshape(V(:,zind_1,zind_2,eind,:),[n_a,Params.J]))
title('Value function: median value of z')
xlabel('Age j')
ylabel('Assets (a)')
subplot(2,1,2); surf(a_grid*ones(1,Params.J),ones(n_a,1)*(Params.agejshifter+(1:1:Params.J)),reshape(V(:,zind_1,zind_2,eind,:),[n_a,Params.J]))
title('Value function: median value of z')
xlabel('Age in Years')
ylabel('Assets (a)')

% + I add eind accordingly on other places to generate other graphs prepared in model 31.

% 2)

jequaloneDist=zeros([n_a,n_z,n_e],'gpuArray');
jequaloneDist(1,floor((n_z(1)+1)/2),floor((n_z(2)+1)/2),floor((n_e+1)/2))=1

FnsToEvaluate.riskyshare=@(savings,riskyshare,a,z1,z2,e) riskyshare;
FnsToEvaluate.stockmarketparticpation=@(savings,riskyshare,a,z1,z2,e) (savings>0)*(riskyshare>0);
FnsToEvaluate.earnings=@(savings,riskyshare,a,z1,z2,e,w,kappa_j) w*kappa_j*z1*z2*e;
FnsToEvaluate.assets=@(savings,riskyshare,a,z1,z2,e) a;
1 Like

Thank you, Robert and Aledinola!

Robert, you are right, I saw in the morning that I did not change the name of the return function, and consequently, the program looked into the wrong return function where shock e was not specified. I corrected this and added Params.e = 1 as the value iteration function required (this seems weird to me as code in model 8 does not require specifying parameter z …). Then the code goes through without error; however, nothing changes in the graphs. I think I should modify something in the plotting section but do not know where and how. :frowning: I tried a similar modification as in Issue 1 (I posted some lines of code a few hours ago) but cannot solve the problem in the same fashion.

I created a GitHub profile and added (Code) in case you immediately see what is going wrong and can advise me on how to solve the problem under Issue 2. (I just modified model 31 and added iid shocks z1 and z2 and exogenous shock e. Without shock e, the code works fine.)

1 Like

Looks fine, presumably you have a reason for wanting the pension to be pension*kappa_j*e rather than just pension*e (just checking as could be right or wrong, depends on what you want).

1 Like

[this posts misuses the term state space in a way that is completely incorrect with it’s technical meaning, but I cannot think of a better one and hopefully you get the idea]

You do NOT want Params.e = 1. VFI Toolkit uses the info in things like n_a and n_z to figure out how many of the first few inputs to ReturnFn are the state space [e.g., if n_d=3, n_a=5, and n_z=4, then the first inputs are (d,aprime,a,z,…) and so it expects the first four will be the state space. Whereas, e.g., if n_d=0, n_a=5, and n_z=[3,4,2] then the first inputs are (aprime,a,z1,z2,z3,…) and so it expects the first five will be the state space). In your case it is using riskyasset, which means that aprime is not part of the state space inputs.

Everytime the toolkit uses ReturnFn or FnsToEvaluate, it counts the number of state space inputs, and then assumes everything after this is a parameter and goes looking for it in the parameters structure.

In your case, e is supposed to be part of the state space (it is an i.i.d. beginning-of-period shock). Your code is not set up properly in terms of n_d, n_a, etc., and so the toolkit is not counting the right number of state space inputs (it is coming up one short). Thus it asks for e which it thinks is one of the parameters, but is in fact part of the state space.

Looking at your code your return function has first inputs (savings,riskyshare,a,z1,z2,e,…). This corresponds correctly to your
n_d=[201,51];
n_a=201;
n_z=[21,5];
n_e=2;
n_u=5;
[as you use riskyasset should have (d1,d2,a,z,e,…) as first return fn inputs]

The error is because you need to put n_e, e_grid and pi_e into both vfoptions and simoptions (currently when you run the value fn command, it is unaware of e, so it thinks the state space is one less than it should be and hence it asks for a parameter e which is actually the last state space entry). [Note that n_d, n_a and n_z are inputs to the value function iteration command, everything else (e,u,riskyasset,expasset, etc.) has to be mentioned using vfoptions, most of these are also relevant to simoptions and so need to be mentioned there too.]

2 Likes

Thank you! (Just want to test something else.)

1 Like

Thank you, Rober!

I set:

vfoptions.n_e = n_e;
vfoptions.e_grid = e_grid;
vfoptions.pi_e = pi_e;
simoptions.n_e = n_e;
simoptions.e_grid = e_grid;
simoptions.pi_e = pi_e;

Unfortunately, I receive:

[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
Test ValueFnIter
Error using ValueFnIter_Case1_FHorz (line 187)
vfoptions.pi_e is not the correct shape (should be of size N_e-by-1 or N_e-by-N_j)

Sorry for so many interactions.

1 Like

Is pi_e a row or a column vector? I think the function you are calling wants pi_e to be a column vector

2 Likes

Thanks, Aledinola. Pi_e seems like a matrix of size 2x2 (see, code).

[I just mimicked model 8 and added n_e, e_grid, and pi_e into both vfoptions and simoptions as Robert suggested (however, model 8 does not add n_z, z_grid, and pi_z into vfoptions and simoptions). Maybe [V, Policy] is not prepared to take exogenous Markov e analogously to exogenous Markov z in model 8. This might be because I already have exogenous iid shocks z1 and z2. Just thinking, but I do not know for sure :(]

1 Like

pi_e should NOT be a matrix, it should be a (column) vector, and this is the error (as it says, pi_e is the wrong size).

The shock is i.i.d, so there should be one probability for each value (each grid point). [Possibly you should have this for each model period j, if you want the probabilities to depend on age.]

PS. Life-Cycle Model 11B shows how to use i.i.d. (e) shocks.

2 Likes