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:

\[\omega_{t+1}^{1}=\frac{(\hat{p}_{t+1}^{s}+\hat{d}_{t+1})s_{t+1}^{1}+\frac{\hat{b}_{t+1}^{1}}{\gamma^a_{t+1}}}{\hat{p}_{t+1}^{s}+\hat{d}_{t+1}}.\]

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

HL1996TFIter.gmod

  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.