# 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:

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

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 54 55 56 | ```
#for i=1:N
var_interp w1_future_#i;
initial w1_future_#i 0.5; % A place holder
#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

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

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

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.