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:

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

HL1996TFIter.gmod

  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.