# Heaton and Lucas (1996): Transition Function Iteration

In the Heaton and Lucas (1996) example,
one of the challenges in solving the model
is the fact that the endogenous state variable, agents 1’s wealth share (or consumption share), has an implicit law of motion.
We tackle this issue by including future values of the state variables among the unknown policy variables
to solve for in each iteration. We also include the *consistency equations* in the system of equilibrium conditions.
These equations require that future realizations of wealth share are consistent with current portfolio choices
and future realizations of stock price, dividend, and aggregate growth rate:

An alternative approach to the implicit law of motion is, in each policy function iteration, to start from a conjectured transition function
for the endogenous state variable and use it to forecast future policy variables (consumption, stock price, etc.).
These forecasts can then be used to solve for current policy variables using the same system of equilibrium conditions (excluding the consistency equations).
At the end of each iteration, the conjectured transition function is updated using the solution of the current policy variables.
In a way, the approach iterates over both policy and transition functions.
Therefore, to differentiate from the standard policy function iteration, we call this method *transition function iteration*.

There are two advantages of this alternative approach compared to our main policy function iteration approach with consistency equations. First, there are a fewer equations to be solved in each iteration. Second, the forecasts for future policy variables can be computed at the beginning of each iteration, independent of equation solving during the iteration. A potential drawback is that the solution might not be as accurate as the original method because the consistency equations might not be satisfied with high precision. However, one might start with this alternative to obtain solution with reasonable accuracy quickly then switch back to the our main approach to achieve higher accuracy.

The alternative approach is developed by Vadim Elenev (Johns Hopkins University), Tim Landvoigt (University Pennsylvania), and Stijn Van Nieuwerburgh (Columbia University), to solve their model with rich empirical content, featuring both financial and non-financial sectors and a government. (The endogenous state variables with implicit laws of motion include aggregate producers’ leverage and aggregate banks’ net-worth.) The model captures salient features of macro-finance dynamics during the Great Recession and has important implications on bank capital regulation.

The authors showed us that this approach can also be used to solve the model in Heaton and Lucas (1996)
with reasonable accuracy and have kindly contributed this example to this website. Our GDSGE toolbox can implement this approach and its gmod file for
Heaton and Lucas (1996)’s model is given below.
We design two features of the toolbox, macros for loop and *pre-model* block, to
faciliate the implementation of the approach and to make the most of its computation economizing property.
We discuss these features in more detail below.

Notice that this approach requires some dampening in the update of policy and transition functions:

```
57damp = 0.25;
58ps_future = ps_future*damp + (1-damp)*ps;
59c1_future = c1_future*damp + (1-damp)*c1;
60c2_future = c2_future*damp + (1-damp)*c2;
61index_current = @(vec,current_idx) squeeze(vec(:,current_idx,:));
62#for i=1:N
63 w1_future_#i = w1_future_#i*damp + (1-damp)*index_current(reshape(w1n,[shock_num GDSGE_SIZE]),#i);
64#end
```

As also noted above, over the iterations, the consistency equations are not required to be satisfied so they are excluded from the system of equilibrium conditions in the gmod file:

```
99 equations;
100 -1+beta*es1+ms1;
101 -1+beta*es2+ms2;
102 -1+beta*eb1+mb1;
103 -1+beta*eb2+mb2;
104 ms1*s1p;
105 ms2*s2p;
106 mb1*nb1p;
107 mb2*nb2p;
108 b1p+b2p;
109 budget_1 / w1;
110 budget_2 / (1-w1);
111 end;
```

Therefore, once the algorithm converges for policy functions, one can switch back to the original gmod file in Heaton and Lucas (1996) with the consistency equations in order to make sure that these equations are satisfied with higher precision. However, the new converged solution would not differ significantly from the old one.

## The gmod File

```
1#define N 8;
2SaveFreq = inf;
3% Parameters
4parameters beta gamma Kb;
5beta = 0.95; % discount factor
6gamma = 1.5; % CRRA coefficient
7Kb = -0.05; % borrowing limit in ratio of aggregate output
8% Exogenous state variables
9var_shock g d eta1;
10% Enumerate exogenous states and transition matrix
11shock_num = N;
12g = [.9904 1.0470 .9904 1.0470 .9904 1.0470 .9904 1.0470];
13d = [.1402 .1437 .1561 .1599 .1402 .1437 .1561 .1599];
14eta1 = [.3772 .3772 .3772 .3772 .6228 .6228 .6228 .6228];
15shock_trans = [
16 0.3932 0.2245 0.0793 0.0453 0.1365 0.0779 0.0275 0.0158
17 0.3044 0.3470 0.0425 0.0484 0.1057 0.1205 0.0147 0.0168
18 0.0484 0.0425 0.3470 0.3044 0.0168 0.0147 0.1205 0.1057
19 0.0453 0.0793 0.2245 0.3932 0.0157 0.0275 0.0779 0.1366
20 0.1366 0.0779 0.0275 0.0157 0.3932 0.2245 0.0793 0.0453
21 0.1057 0.1205 0.0147 0.0168 0.3044 0.3470 0.0425 0.0484
22 0.0168 0.0147 0.1205 0.1057 0.0484 0.0425 0.3470 0.3044
23 0.0158 0.0275 0.0779 0.1365 0.0453 0.0793 0.2245 0.3932
24 ];
25% Endogenous state variables
26var_state w1; % wealth share
27w1 = linspace(-0.04,1.04,201);
28% Policy variables and bounds that enter the equations
29var_policy c1 c2 s1p nb1p nb2p ms1 ms2 mb1 mb2 ps pb;
30inbound c1 0.1 1;
31inbound c2 0.1 1;
32inbound s1p 0.0 1.0;
33inbound nb1p 0.0 1.0; % nb1p=b1p-Kb
34inbound nb2p 0.0 1.0;
35inbound ms1 0 1; % Multipliers for constraints
36inbound ms2 0 1;
37inbound mb1 0 1;
38inbound mb2 0 1;
39inbound ps 0 10 adaptive(1.5);
40inbound pb 0 10 adaptive(1.5);
41% Other policy variables
42var_aux equity_premium w1n[N];
43% Interpolation variables for policy and state transitions
44var_interp ps_future c1_future c2_future;
45ergodic = shock_trans^50;
46ergodic = ergodic(1,:);
47deltabar = ergodic * d';
48etabar = ergodic * eta1';
49gammabar = ergodic * g';
50initial ps_future beta*deltabar*gammabar^(1-gamma)/(1-beta*gammabar^(1-gamma));;
51initial c1_future (1-deltabar)*etabar+0.5*deltabar;
52initial c2_future (1-deltabar)*(1-etabar)+0.5*deltabar;
53#for i=1:N
54 var_interp w1_future_#i;
55 initial w1_future_#i 0.5; % A place holder
56#end
57damp = 0.25;
58ps_future = ps_future*damp + (1-damp)*ps;
59c1_future = c1_future*damp + (1-damp)*c1;
60c2_future = c2_future*damp + (1-damp)*c2;
61index_current = @(vec,current_idx) squeeze(vec(:,current_idx,:));
62#for i=1:N
63 w1_future_#i = w1_future_#i*damp + (1-damp)*index_current(reshape(w1n,[shock_num GDSGE_SIZE]),#i);
64#end
65
66pre_model;
67 % Evaluate interpolation
68 w1n_fixed' = 0;
69 #for i=1:N
70 if shock==#i
71 w1n_fixed' = w1_future_#i'(w1);
72 end
73 #end
74 [psn',c1n',c2n'] = GDSGE_INTERP_VEC'[1:3](w1n_fixed');
75 es1_future = GDSGE_EXPECT{g'^(1-gamma)*(c1n')^(-gamma)*(psn'+d')};
76 es2_future = GDSGE_EXPECT{g'^(1-gamma)*(c2n')^(-gamma)*(psn'+d')};
77 eb1_future = GDSGE_EXPECT{g'^(-gamma)*(c1n')^(-gamma)};
78 eb2_future = GDSGE_EXPECT{g'^(-gamma)*(c2n')^(-gamma)};
79end;
80
81model;
82 % Calculate expectations that enter the Euler Equations
83 es1 = es1_future*c1^gamma/ps;
84 es2 = es2_future*c2^gamma/ps;
85 eb1 = eb1_future*c1^gamma/pb;
86 eb2 = eb2_future*c2^gamma/pb;
87 % Transform bond back
88 b1p = nb1p + Kb;
89 b2p = nb2p + Kb;
90 % Market clearing of shares
91 s2p = 1-s1p;
92 % Budget constraints
93 budget_1 = w1*(ps+d)+eta1 - c1 - ps*s1p - pb*b1p;
94 budget_2 = (1-w1)*(ps+d)+(1-eta1) - c2 - ps*s2p - pb*b2p;
95 % Consistency equations
96 w1n' = (s1p*(psn'+d') + b1p/g')/(psn'+d');
97 % Other policy variables
98 equity_premium = GDSGE_EXPECT{(psn'+d')/ps*g'} - 1/pb;
99 equations;
100 -1+beta*es1+ms1;
101 -1+beta*es2+ms2;
102 -1+beta*eb1+mb1;
103 -1+beta*eb2+mb2;
104 ms1*s1p;
105 ms2*s2p;
106 mb1*nb1p;
107 mb2*nb2p;
108 b1p+b2p;
109 budget_1 / w1;
110 budget_2 / (1-w1);
111 end;
112end;
113
114simulate;
115 num_periods = 10000;
116 num_samples = 24;
117 initial w1 0.5;
118 initial shock 1;
119 var_simu c1 c2 ps pb equity_premium ms1 mb1;
120 w1' = w1n';
121end;
```

Notice that, to implement this alternative algorithm, the gmod code uses two features of the toolbox on top of the original example with consistency equations.
First, since this alternative algorithm needs to evaluate the state transition function for the current state
to get the future endogenous state variables that corresponds to each of the future realizations, one needs to include the state transition function
*for each of the current state* as *var_interp*. For the current example, this means 8 state transition functions corresponding to each of the current state.
Accordingly, these state transition functions need to be initialized, updated, and evaluated in the model block. This can be done compactly by using a macro for loop
in the gmod file, replicated below

```
53#for i=1:N
54 var_interp w1_future_#i;
55 initial w1_future_#i 0.5; % A place holder
56#end
```

This block expands to

```
var_interp w1_future_1;
initial w1_future_1 0.5; % A place holder
var_interp w1_future_2;
initial w1_future_2 0.5; % A place holder
...
var_interp w1_future_8;
initial w1_future_8 0.5; % A place holder
```

In particular, the #for … #end block will repeat the code in the block, and replace the iterator with each value in the value list in the header line starting with #for. And the above code block thus declares
w1_future_1 to w1_future_8 as *var_interp* and initialize their values. Correspondingly, the following block specifies the updates for each of these state transition functions
after a time step of the policy iteration

```
61index_current = @(vec,current_idx) squeeze(vec(:,current_idx,:));
62#for i=1:N
63 w1_future_#i = w1_future_#i*damp + (1-damp)*index_current(reshape(w1n,[shock_num GDSGE_SIZE]),#i);
64#end
```

in which *index_current* is a helper function to look up the corresponding current shock index from the returned *var_aux* w1n after each iteration. Therefore, *w1_future_1* indeed corresponds to
the state transition function for the current exogenous state indexed by 1, and so on.

The other feature that maximizes the efficiency of the algorithm is by noticing that the future endogenous states are directly evaluated with the state transition function from the last iteration.
Then the future endogenous states are used to evaluate necessary future variables using the policy function from the last iteration.
Therefore, given current
endogenous states and exogenous states, the future endogenous states and future variables do not depend on current policy variables, and thus can be evaluated before solving the equation.
This is enabled by defining the *pre_model* block in the gmod code, replicated below

```
66pre_model;
67 % Evaluate interpolation
68 w1n_fixed' = 0;
69 #for i=1:N
70 if shock==#i
71 w1n_fixed' = w1_future_#i'(w1);
72 end
73 #end
74 [psn',c1n',c2n'] = GDSGE_INTERP_VEC'[1:3](w1n_fixed');
75 es1_future = GDSGE_EXPECT{g'^(1-gamma)*(c1n')^(-gamma)*(psn'+d')};
76 es2_future = GDSGE_EXPECT{g'^(1-gamma)*(c2n')^(-gamma)*(psn'+d')};
77 eb1_future = GDSGE_EXPECT{g'^(-gamma)*(c1n')^(-gamma)};
78 eb2_future = GDSGE_EXPECT{g'^(-gamma)*(c2n')^(-gamma)};
79end;
```

As shown, the *pre_model* block essentially evaluates the future endogenous states, using the state transition function that corresponds to the index of the current exogenous state,
and then evaluates the future variables that are necessary in evaluating the current equation, by feeding the future endogenous states into the policy functions from the last iteration.
Actually, the *pre_model* block does a bit more by also calculating the expectations that enter the inter-temporal equations, as they do not depend on current policy variables either.
These expectation terms, combined with current policy variables, arrive at the inter-temporal equations (the Euler equations for bonds and shares) that enter the system, defined in the model block. It should be clear that any variables defined in the *pre_model* block can be directly used in the *model* block, and for obvious reasons, the *pre_model* block can not use any
policy variables declared as *var_policy*. Evaluating the state transition functions and expectations beforehand has greatly improved the efficiency of this alternative algorithm (for the current example, it reduces the computation time by half), by avoiding redundant evaluations of objects that do not vary with *var_policy*.

Generally, this alternative algorithm, using the most efficient implementation enabled by procedures above, delivers better performance in computation speed for a time step compared to the algorithm with consistency equations
in Heaton and Lucas (1996).
The tradeoff in computation time is that the algorithm with consistency equations solves a larger system of equations
(by including the consistency equations and endogenous states as unknowns)
and does not allow interpolating for future endogenous states and computing expectation beforehand,
but this alternative algorithm constructs more *var_interp* and needs evaluate interpolation twice
– first to get future endogenous states, then future policy variables.
On the theoretical ground, the algorithm with consistency equations finds solutions as the limit of a finite-horizon economy
and thus can also be used to solve stochastic transition paths.
Besides, it converges more robustly without resorting to ad-hoc dampening.