Exogenous shock both age and type dependent

Consider a life-cycle model with two exogenous shocks (\eta,h). The shock \eta \in \left\lbrace \eta_1, \eta_2,..,\eta_{n_{\eta}}\right\rbrace is a standard Markov shock with transition matrix \pi_{\eta}(\eta,\eta). The transition probabilities of the shock h \in \{0,1 \} instead depends both on age j=1,..,J and on a permanent type denoted as \theta \in \{0,1 \}. We can therefore write the transition matrix for h as \pi_h(h,h'|j,\theta).

Possible applications. Imagine \eta is a labor productivity shock discretized with Tauchen or Farmer-Toda’s method. h is a health shock, where e.g. h=0 denotes unhealthy, whose evolution over the life-cycle is stochastic, depends on age (older people are more likely to be unhealthy), and it also depends on a fixed education type.

How do we set up this in the toolkit?
First, we discretize the standard Markov shock \eta getting a grid eta_grid (size [n_eta,1]) and a transition matrix pi_eta (size [n_eta,n_eta]).
Second, the grid for h does not change with j nor with \theta, so it is just h_grid=[0,1]'. The transition matrix is actually a structure with two fields, pi_h.theta0 (size [n_h,n_h,J]), pi_h.theta1 (size [n_h,n_h,J]).

The final step is to combine the two shocks into a single shock, as z=(\eta,h), with \eta ordered first, and n_z=[n_eta,n_h]. The big transition matrix \pi_z is actually pi_z_J.theta0 (with size [n_eta*n_h,n_eta*n_h,J]) and pi_z_J.theta1 (with size [n_eta*n_h,n_eta*n_h,J]) and are defined as

\pi_z.\theta_0(:,:,j) = \pi_h.\theta_0(:,:,j) \otimes \pi_{\eta}

for all j=1,..,J, and

\pi_z.\theta_1(:,:,j) = \pi_h.\theta_1(:,:,j) \otimes \pi_{\eta}

for all j=1,..,J. Note the kron in reverse order: \eta is ordered first, but comes second in the kron product (if you don’t see why, check here)

I coded the above as follows

	n_z    = [n_eta,n_h]; %row vector
	z_grid = [eta_grid;health_grid]; % column vector
	
    % Names_i  = {'th0','th1'}
	pi_z_J.th0 = zeros(prod(n_z),prod(n_z),N_j);
	pi_z_J.th1 = zeros(prod(n_z),prod(n_z),N_j);
	
	for j_c=1:N_j
		% Note reverse order
		pi_z_J.th0(:,:,j_c) = kron(pi_health.th0(:,:,j_c),pi_eta);
		pi_z_J.th1(:,:,j_c) = kron(pi_health.th1(:,:,j_c),pi_eta);
	end

It would be great to have @robertdkirkby’s blessing that this setup is correct :slight_smile:
I think it could be useful for researchers working on models with health risk. From a more selfish angle, I’m using this setup myself, but I’ve noticed that some of the age-conditional stats don’t look quite right. I’ll outline the possible bug in a follow-up question, but for now I just wanted to double-check that the overall approach is sound.

Note: I’ve tried out the equivalent approach with ExogShockFn and it gives exactly the same (possible wrong :face_with_peeking_eye:) results.

This looks correct to me.

The only thing I am not certain about is what is going to happen when z_grid is not dependent on age and ptype, while pi_z is dependent on age and ptype. I think this will be fine, but I have never actually run one with this kind of combo (my test examples all had the same dependences in z_grid and pi_z).

Easy test that occurs to me to find out if they worked correctly, looking at StationaryDist but using sum() to remove all the dimensions except for the (eta,h,j). You can then easily manually double check the first two or three periods. This would test the transition probabilities. You could then calculate AgeConditionalStats for ‘h’ and ‘eta’, and manually compute them as just ‘sum(h_grid.*StationaryDist)’ [skipping over the details in terms of which dimensions of StationaryDist here as that is likely obvious], again can just check the first two or three periods. This would test the grids being handled correctly.

1 Like

Thanks for your suggestion, I will compute the evolution of the exogenous states (\eta,h,j) “manually” and check with the corresponding elements of the agents distribution computed by the toolkit.

Another quick check is just to replicate z_grid along the age and type dimensions and verify the results are the same, e.g.

z_grid.th0 = repmat(z_grid,[1,N_j]);
z_grid.th1 = repmat(z_grid,[1,N_j]);

so that each z_grid.PTYPE is an array with size [sum(n_z),N_j].

Update

I have checked with

z_grid_vec = [eta_grid;health_grid];

% The next two lines may be redundant
z_grid.th0 = repmat(z_grid_vec,[1,N_j]);
z_grid.th1 = repmat(z_grid_vec,[1,N_j]);

and results are identical.

Another quick check is just to replicate z_grid along the age and type dimensions and verify the results are the same, e.g.

z_grid.th0 = repmat(z_grid,[1,N_j]);
z_grid.th1 = repmat(z_grid,[1,N_j]);

so that each z_grid.PTYPE is an array with size [sum(n_z),N_j].

This is fine, but it is checking your own code, not the toolkit code (still a good idea, just want to be clear).

1 Like