Heaton and Lucas (1996): Incomplete Markets with Portfolio Choices

The benchmark model in Heaton and Lucas (1996) is a good starting point to demonstrate the capability of the current framework in dealing with endogenous state variables with implicit law of motions. The model encompasses many ingredients that appear in recent macroeconomic studies, such as incomplete markets, portfolio choices, occasionally binding constraint, non-stationary shock process, and asset pricing with non-trivial market-clearing conditions. We show how the model can be solved with wealth share or consumption share as the endogenous state, which are the two prominent approaches in the literature and naturally fit in our toolbox framework.

The Model

This is an incomplete-markets model with two representative agents $i\in\mathcal{I}=\{1,2\}$ who trade in equity shares and bonds. The aggregate state $z\in\boldsymbol{Z}$, which consists of capital income share, agents’ income share, and aggregate endowment growth, follows a first-order Markov process. $p_{t}^{s}(z^t)$ and $p_{t}^{b}(z^t)$ denote share price and bond price at time $t$ and in shock history $z^t=\{z_0,z_1,\dots,z_t\}$. To simplify the notations, we omit the explicit dependence on shock history.

Agent $i$ takes the share and bond prices as given and maximizes her inter-temporal expected utility

\[\mathcal{U}_{t}^{i}=\mathbb{E}_{t}\left[\sum_{\tau=0}^{\infty}\beta^{\tau}\frac{\left(c_{t+\tau}^{i}\right)^{1-\gamma}}{1-\gamma}\right]\]

subject to

\[c_{t}^{i}+p_{t}^{s}s_{t+1}^{i}+p_{t}^{b}b_{t+1}^{i}\leq(p_{t}^{s}+d_{t})s_{t}^{i}+b_{t}^{i}+Y_{t}^{i}\]

and

\[\begin{split}s_{t+1}^{i} & \geq0, \\ b_{t+1}^{i} & \geq K^b_t,\end{split}\]

where $Y^a_t$ denotes the aggregate income. $d_t = \delta_t Y^a_t$ is total dividend (capital income) and $Y^i_t = \eta^i_t Y^a_t$ is labor income of agent $i$. Aggregate income grows at a stochastic rate $\gamma^a_t = \frac{Y^a_t}{Y^a_{t-1}}$. $z_t = \{\gamma^a_t,\delta_t,\eta^1_t\}$ follows a first-order Markov process estimated using U.S. data. The borrowing limit is set to be a constant fraction of per capita income, i.e., $K^b_t = \bar{K}^b Y^a_t$.

In equilibrium, prices are determined such that markets clear in each shock history:

\[\begin{split}& s_{t}^{1}+s_{t}^{2}=1,\\ & b_{t}^{1}+b_{t}^{2}=0.\end{split}\]

We use the financial wealth share

\[\omega_{t}^{i}=\frac{(p_{t}^{s}+d_{t})s_{t}^{i}+b_{t}^{i}}{p_{t}^{s}+d_{t}}\]

as an endogenous state variable. In equilibrium, the market clearing conditions imply that $\omega^1_t + \omega^2_t = 1$.

For any variable $x_t$, let $\hat{x}_t$ denote the normalized variable: $\hat{x}_t=\frac{x_t}{Y^a_t}$ (except $b^i_t$ for which $\hat{b}^i_t = \frac{b^i_t}{Y^a_{t-1}}$). Using this normalization, agent i’s budget constraint can be rewritten as

\[\hat{c}_{t}^{i}+\hat{p}_{t}^{s}s_{t+1}^{i}+p_{t}^{b}\hat{b}_{t+1}^{i}\leq\left(\hat{p}_{t}^{s}+\hat{d}_{t}\right)\omega_{t}^{i}+\hat{Y}_{t}^{i}.\]

The financial wealth share is rewritten as

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

The optimality of agent i’s consumption and asset choices is captured by the first-order conditions in $s^i_{t+1}$ and $b^i_{t+1}$:

\[\begin{split}1& =\beta\mathbb{E}_{t}\left[\left(\frac{\hat{c}_{t+1}^{i}}{\hat{c}^i_t}\right)^{-\gamma}\left(\gamma_{t+1}^{a}\right)^{1-\gamma}\frac{\hat{p}_{t+1}^{s}+\hat{d}_{t+1}}{\hat{p}_{t}^{s}}\right]+\hat{\mu}^{i,s}_t\\ 1& =\beta\mathbb{E}_{t}\left[\left(\frac{\hat{c}_{t+1}^{i}}{\hat{c}^i_t}\right)^{-\gamma}\left(\gamma_{t+1}^{a}\right)^{-\gamma}\frac{1}{p_{t}^{b}}\right]+\hat{\mu}^{i,b}_t,\end{split}\]

where $\hat{\mu}^{i,s}_t$ and $\hat{\mu}^{i,b}_t$ are the Lagrangian multipliers on agent i’s no short sale constraint and borrowing constraint, respectively. The multipliers and portfolio choices satisfy the complementary-slackness conditions:

\[\begin{split}0 & = \hat{\mu}^{i,s}_t s^i_{t+1}, \\ 0 & = \hat{\mu}^{i,b}_t (\hat{b}^i_{t+1} - \bar{K}^b).\end{split}\]

Evaluate the Accuracy of Solutions

The converged policy iterations deliver both the policy functions and the state transition functions, which can be used conveniently to evaluate the accuracy of the solutions by e.g., inspecting the Euler equation errors. Define the unit-free Euler equation errors for shares and bonds as

\[\begin{split}& \mathcal{E}^{s,i}_t = -1 + \beta\mathbb{E}_{t}\left[\left(\frac{\hat{c}_{t+1}^{i}}{\hat{c}^i_t}\right)^{-\gamma}\left(\gamma_{t+1}^{a}\right)^{1-\gamma}\frac{\hat{p}_{t+1}^{s}+\hat{d}_{t+1}}{\hat{p}_{t}^{s}}\right]+\hat{\mu}^{i,s}_t\\ & \mathcal{E}^{b,i}_t = -1 + \beta\mathbb{E}_{t}\left[\left(\frac{\hat{c}_{t+1}^{i}}{c^i_t}\right)^{-\gamma}\left(\gamma_{t+1}^{a}\right)^{-\gamma}\frac{1}{p_{t}^{b}}\right]+\hat{\mu}^{i,b}_t.\end{split}\]

And we are to evaluate these errors starting from a simulated ergodic set of states. Due to symmetry, we focus on Agent 1 only. This can be done by simulating, starting from the ergodic set, for one period forward, and evaluating these errors according to the equations defined above using the simulated paths. In particular, to calculate the expectation of the objects in the equations along the simulated path, we should use the exact transition path for integration to eliminate sampling errors.

The MATLAB code that implements the above procedure is listed below (or download eval_euler_errors.m), which should be self-explanatory with the documentation contained.

% Extract the ergodic set
NUM_PERIODS = 1000;
w1 = reshape(SimuRslt.w1(:,end-NUM_PERIODS:end),1,[]);
shock = reshape(SimuRslt.shock(:,end-NUM_PERIODS:end),1,[]);
% Replicate the sample to accommodate future shock realizations
w1 = repmat(w1,IterRslt.shock_num,1);
shock1 = repmat(shock,IterRslt.shock_num,1);
shock2 = repmat([1:IterRslt.shock_num]',1,size(w1,2));

% Simulate forward for one period
simuOptions = struct;
simuOptions.init.w1 = w1(:);
simuOptions.init.shock = [shock1(:),shock2(:)];
% The following line states that the first two-period shock indexes are
% supplied and not regenerated
simuOptions.GEN_SHOCK_START_PERIOD = 2;
simuOptions.num_samples = numel(w1);
simuOptions.num_periods = 2;    % One-period forecasting error
% The following line simulates one period forward by starting from w1 and
% shock contained in simuOptions
simuForward = simulate_HL1996(IterRslt,simuOptions);

% Calculate Errors
beta = IterRslt.params.beta;
gamma = IterRslt.params.gamma;
c1 = simuForward.c1(:,1); c1n = simuForward.c1(:,2);
ps = simuForward.ps(:,1); psn = simuForward.ps(:,2);
pb = simuForward.pb(:,1);
ms1 = simuForward.ms1(:,1); mb1 = simuForward.mb1(:,1);
gn = IterRslt.var_shock.g(shock2(:))';
dn = IterRslt.var_shock.d(shock2(:))';
es1_error = -1 + beta*gn.^(1-gamma).*(c1n./c1).^(-gamma).*(psn+dn)./ps + ms1;
eb1_error = -1 + beta*gn.^(-gamma).*(c1n./c1).^(-gamma)./pb + mb1;
% Calculate expectation errors, integrating using the transition matrix
shock_trans = IterRslt.shock_trans(shock,:)';
shock_num = IterRslt.shock_num;
es1_expect_error = sum(shock_trans.*reshape(es1_error,shock_num,[]),1);
max_abs_es1_error = max(abs(es1_expect_error))
mean_abs_es1_error = mean(abs(es1_expect_error))
eb1_expect_error = sum(shock_trans.*reshape(eb1_error,shock_num,[]),1);
max_abs_eb1_error = max(abs(eb1_expect_error))
mean_abs_eb1_error = mean(abs(eb1_expect_error))

Running the code produces:

max_abs_es1_error =

    0.0057

mean_abs_es1_error =

2.5290e-05

max_abs_eb1_error =

    0.0036

mean_abs_eb1_error =

2.1279e-05

which says that the max and mean absolute errors (across states in the ergodic distribution) for stock Euler equations are 0.0057 and 2.5290e-05 respectively, and the max and mean absolute errors for bond Euler equations are 0.0036 and 2.1279e-05 respectively.

These are relatively errors in marginal utility, to convert them into relative errors in consumption, we just need to multiply them by $1/\gamma$. Therefore, the max and mean absolute errors in stock holding decisions are $38 and $0.169 per $10,000 in consumption. The max and mean absolute errors in bond holding decisions are $24 and $0.142 per $10,000 in consumption.

Increasing the number of grid points reduces the errors. For example, increasing the number of grid points to 1000 reduces the max absolute (Euler equation) errors below 1e-3 and mean absolute errors below 1e-5. However, a more effective approach is to use the adaptive-grid interpolation method included in the toolbox. Applying to the current model, the method reduces the max absolute errors below 1e-5 and mean absolute errors below 1e-6, with the number of grid points smaller than 1000. The user only needs to specify a one-line option to enable the adaptive-grid method, but does need to initialize var_interp with a more flexible model_init block. See example Bianchi (2011) for how to define a model_init block and apply the adaptive-grid method.

Consumption Share as the Endogenous State

The model can be solved using consumption share as the endogenous state. In this case, the budget constraint

\[\hat{c}_{t+1}^i=s_{t+1}^i (\hat{p}_{t+1}^s + \hat{d}_{t+1})+ \frac{\hat{b}_{t+1}^i }{g_{t+1}}+ \underbrace{\eta_{t+1}^i - \hat{p}_{t+1}^s s_{t+2}^i-p_{t+1}^b \hat{b}_{t+2}^i}_{\text{Financial Wealth}_{t+1}}\]

is a natural consistency equation for the transition of consumption share $\hat{c}^1$. Specifically, with consumption share, the recursive equilibrium can be defined as ${s^i}'(z,\hat{c}^1),{\hat{b}^i}'(z,\hat{c}^1), \hat{p}^s(z,\hat{c}^1),p^b(z,\hat{c}^1), {\hat{c}^1}'(z';z,\hat{c}^1)$ such that

\[\begin{split}-1+\beta \mathbb{E}_t \Big[\gamma^{1-\gamma}_{t+1}\frac{[\hat{c}_{t+1}^i]^{-\gamma}}{[\hat{c}_t^i]^{-\gamma} } \frac{\hat{p}_{t+1}^s + \hat{d}_{t+1}}{\hat{p}_t^s}] + \hat{\mu}^{i,s}_t=0, \forall i=1,2 \\ -1+\beta \mathbb{E}_t \Big[\gamma_{t+1}^{-\gamma}\frac{[\hat{c}_{t+1}^i]^{-\gamma}}{[\hat{c}_t^i]^{-\gamma} } \frac{1}{{p}_t^b}\Big] + \hat{\mu}^{i,b}_t=0, \forall i=1,2 \\ \hat{b}_{t+1}^1+\hat{b}_{t+1}^2=0 \\ s_{t+1}^1+s_{t+1}^2=1 \\ \hat{c}^1_{t+1}=s_{t+1}^i (\hat{p}_{t+1}^s + \hat{d}_{t+1})+ \frac{\hat{b}_{t+1}^i }{\gamma_{t+1}}+ \underbrace{\eta_{t+1}^i - \hat{p}_{t+1}^s s_{t+2}^i-p_{t+1}^b \hat{b}_{t+2}^i}_{\text{Financial Wealth}_{t+1}}, \forall z_{t+1}\end{split}\]

where $\hat{c}^2$ (and $(\hat{c}^2)'$) can be trivially inferred from the goods market clearing condition $\hat{c}^1+\hat{c}^2=1 +\hat{d}$, and, hence does not need to be defined as extra unknowns when evaluating the equation system. The “Financial Wealth” is a function of future endogenous states, and can be part of the policy functions, returned as one of var_aux at each iteration.

The gmod file that implements the recursive system is (HL1996_consumption_share.gmod)

% 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 = 8;
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 c1;  % consumption
c1 = linspace(0.2,0.8,101);
% Policy variables and bounds that enter the equations
var_policy s1p nb1p nb2p ms1 ms2 mb1 mb2 ps pb c1n[8];
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 2 adaptive(1.5);
inbound pb 0 2 adaptive(1.5);
inbound c1n 0.0 1.0;
% Other policy variables
var_aux equity_premium c2 w1 flow;
% Interpolation variables for policy and state transitions
var_interp ps_future flow_future;
initial ps_future 0.0;
initial flow_future eta1;
ps_future = ps;
flow_future = flow;

model;
  % Interpolation
  [psn',flow_future'] = GDSGE_INTERP_VEC'(c1n');
  % Goods market clear
  c2n' = 1+d'-c1n';
  c2 = 1+d-c1;
  % Expectations in Euler Equations
  es1 = GDSGE_EXPECT{g'^(1-gamma)*(c1n'/c1)^(-gamma)*(psn'+d')/ps};
  es2 = GDSGE_EXPECT{g'^(1-gamma)*(c2n'/c2)^(-gamma)*(psn'+d')/ps};
  eb1 = GDSGE_EXPECT{g'^(-gamma)*(c1n'/c1)^(-gamma)/pb};
  eb2 = GDSGE_EXPECT{g'^(-gamma)*(c2n'/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
  c1_consis' = s1p*(psn'+d') + b1p/g' + flow_future' - c1n';
  % Other policy variables
  w1 = (c1 + ps*s1p + pb*b1p - eta1) / (ps + d);
  flow = eta1 - ps*s1p - pb*b1p;
  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;
    c1_consis';
  end;
end;

simulate;
  num_periods = 10000;
  num_samples = 24;
  initial c1 0.5;
  initial shock 1;
  var_simu w1 c2 ps pb equity_premium;
  c1' = c1n';
end;

As shown, compared to the one with wealth share as the endogenous state, the new implementation is made possible by declaring ${\hat{c}_1}'(z')$ to be c1n in

27var_policy s1p nb1p nb2p ms1 ms2 mb1 mb2 ps pb c1n[8];

by defining the “Financial Wealth” by flow as var_interp

41var_interp ps_future flow_future;

by defining the consistency equations for $\hat{c}_1'$ in

67  c1_consis' = s1p*(psn'+d') + b1p/g' + flow_future' - c1n';

and including them as part of the equation system.

Finally we compare the solutions solved with wealth share as the endogenous state and consumption share as the endogenous state. This can be done by projecting the solutions to the same endogenous state. For example

69  w1 = (c1 + ps*s1p + pb*b1p - eta1) / (ps + d);

constructs the wealth share from the budget constraint of Agent 1.

../../_images/policy_premium_overlapped.png

As shown, the two solutions (solid lines for wealth share as endogenous state and markers “X” for consumption share as endogenous state) are not visually distinguishable.

Comparison with Alternative Algorithms

The Simultaneous Transition and Policy Function Iteration (STPFI) algorithm and toolbox especially demonstrates its strength at solving models with endogenous state variables that feature implicit laws of motions, like the current one (here, wealth share is such a state).

We also compare the performance of the algorithm with alternative methods that can solve the current model. The results are summarized in the table below.

As shown, Algorithms (1)-(3) use the STPFI algorithm, and show that higher accuracies can be achieved by using a denser grid at costs that increase linearly in the number of grids, or, using an adaptive grid at lower costs.

Algorithm (4) is the transition function iteration (TFI) method developed by Elenev, et al. (2021). This algorithm is also based on time iterations, but unlike ours which solves the transition and policy functions jointly in each iteration, theirs fixes the state transition function implied by the policy functions solved from the previous iteration. We implement this alternative algorithm with the same numerical routines (including the equation solver and function approximation procedure) in C++, starting from the same initial transition function, and fine-tuned for maximum numerical efficiency— e.g., we pre-compute the expectation terms when they do not depend on current policy variables under the fixed state transition functions —so the performance of the two algorithms can be compared fairly. Using dampened updating with a damping parameter of 0.25 (i.e., with weight 0.75 assigned to updated transition functions), Algorithm (4) can converge successfully, and the Euler equation errors are comparable to STPFI. However, the speed is significantly slower than STPFI, mainly because their algorithm fixes the state transition function and does not allow the future endogenous state (in the current example, future wealth share) to respond to current policy variables (share and bond choices etc.), which renders finding a solution to the equilibrium system more difficult. A more detailed discussion and the implementation of TFI using the toolbox can be found at here.

The second alternative algorithm, (5) in the table, is the nested fixed point algorithm. The algorithm proceeds in a nested loop, with the inner loop solving for converged policy functions fixing the state transition function, and the outer loop solving for the state transition function that is consistent with the one implied by the inner loop. We find that such an algorithm is very unstable, and can only start to converge with a sufficient large dampening parameter for updating the transition function, especially for early iterations. Even with a large dampening parameter (0.75) the algorithm can still only converge to a level of 5e-5 (measured by the metric between policy functions of adjacent iterations), although the Euler equation errors are close to those of the STPFI benchmark. Due to its nested nature, the speed is also significantly slower than STPFI. Guvenen (2009) uses a similar nested fixed point algorithm for solving an asset pricing model with a richer setup. We demonstrate how to solve the model in the toolbox and discuss the performance comparison in more detail here.

The third alternative algorithm, (6) in the table, is the one developed by Kubler and Schmedders (2003) and is labeled as NFP-KS. Like STPFI, their algorithm solves policy variables and future state variables at each time step that respect all short-run equilibrium conditions but does this in a nested manner: in the inner nest the future state variables are solved to respect the consistency equations with current policy variables fixed, and in the outer nest the policy variables are searched to satisfy the remaining equilibrium conditions. As shown in the table, the algorithm is slower than STPFI, but once converged under the same criterion, delivers solutions with similar accuracies. NFP-KS is slower because it incurs substantially higher number of evaluations of the consistency equations, which outweigh the benefits of solving a smaller system of equations each time.

To understand what drives the difference in the performance across algorithms, we have produced a detailed breakdown of computation costs of all the listed algorithms. To summarize, the difference is mainly due to that other algorithms economize on evaluations of a subset of equilibrium conditions but require substantially more evaluations of other conditions, whereas the STPFI algorithm achieves a balance in evaluating all equilibrium conditions.

In more detail, all algorithms are based on time iterations, and at each time step solve an equilibrium system of equations which involves evaluating equilibrium conditions repeatedly. The equilibrium conditions can be broken into three blocks: intra-temporal evaluations (i.e., market clearing conditions, budget constraints etc. which involve only current equilibrium variables), inter-temporal evaluations (i.e., conditions which involve future variables—in particular—Euler equations), and consistency equations; different algorithms have comparative advantages in reducing the number of evaluations for a certain block. As shown in Table R1, the STPFI algorithm solves the whole short-term equilibrium system simultaneously, so the number of evaluations of the three categories are equal. This offers a useful benchmark to understand the costs of each computation block for a typical economic model like the current one: the computation costs of intra-temporal and inter-temporal evaluations, in terms of number of evaluations, are of the same order of magnitude (the computation time of evaluating inter-temporal conditions is still higher since they involve evaluating function approximations and integrations when constructing the expected values of future policies), whereas the costs of evaluating consistency equations are much smaller. This observation makes it clear why the STPFI algorithm can be potentially efficient: by including extra consistency equations in the equation system, it enables the whole short-term equilibrium system to be solved in one step, and minimizes the number of evaluations of intra-temporal and inter-temporal conditions.

To demonstrate this feature of the STPFI algorithm more clearly, let us compare the baseline algorithm (Algorithm (1)) and the Transition Function Iteration algorithm (TFI, Algorithm (4)). As described in detail in the paper, the TFI algorithm fixes the state transition function implied by the previous time step instead of solving a consistent one, and only updates the state transition function at the end of a time step. Such a design minimizes the number of evaluations for inter-temporal evaluations, since the values of future state variables do not depend on policy variables within a time step, and thus all the expectation terms that enter the inter-temporal conditions are only evaluated once and then fixed when solving the system of intra-temporal equilibrium equations; also, the costs of evaluating consistency equations are minimal as they are only required at the end of each time step when updating the state transition function. However, such a design is at the cost of increasing the number of intra-temporal evaluations, since by fixing the value of future states instead of allowing them to respond to current policy variables, the model is distant from its true dynamics, which makes solving the remaining equilibrium system rather difficult. This extra burden imposed on intra-temporal evaluations is huge: as shown in the table, the TFI algorithm requires intra-temporal evaluations one order of magnitude more than the baseline, rendering it a less effective algorithm for the current problem.

Similarly, the nested fixed point algorithm of Kubler and Schmedders (2003) (Algorithm (6)), which solves consistency equations state by state as inner nested problems and the remaining system of equations as the outer problem, is designed to reduce the system of equations to be solved simultaneously at a time. However, due to the nested nature of the algorithm, the number of evaluations of consistency equations is one order of magnitude more than that required by the STPFI algorithm; besides, evaluating consistency equations for the NFP-KS algorithm is also more costly, as it involves extra evaluations of off-equilibrium future share prices when searching for consistent future wealth shares in the inner loop—such evaluations of future share prices when using the STPFI algorithm do not incur extra costs, because future share prices are also used in evaluating inter-temporal conditions, of which the computation costs are counted there. This explains why NFP-KS is a slower algorithm for the current problem.

The nested fixed point algorithm of Guvenen (2009) (Algorithm (5)) is the slowest as it involves solving an inner converged policy functions, holding the state transition functions at the outer loop that are potentially distant from the equilibrium ones. As shown, the required number of policy function iterations is much larger, despite that later policy function iterations can converge in only a few steps using the “warm-up” solutions. Besides, the nested iterative algorithm is not amenable to an automated design, and the toolbox needs to be combined with MATLAB code that is responsible for conducting the outer loops for transition function iterations. This is why there is substantially more amount of time spent in the “other” category, as loops in MATLAB do not achieve the same performance as the toolbox, of which the main computations are delegated to compiled C++ binaries.

In summary, the STPFI algorithm is more efficient for the current model since it balances the number of evaluations of the three computation blocks, and for the current model, the costs of evaluating inter-temporal and intra-temporal evaluations are of the same order of magnitude. In other models where inter-temporal evaluations are more costly, such as models involving many Euler equations for multiple asset choices or with less complex intra-temporal conditions, or in models where the number of policy variables is large so reducing the system of equations is of first-order importance, the other algorithms may show their strength and likely to be superior.

Toolbox enabled by STPFI

Lastly, compared with nested algorithms such as NFP-KS (Algorithm (6)) and NFP-Guvenen (Algorithm (5)), the non-nested feature of STPFI also makes it more amenable to an automated design, for which the current toolbox can be developed. To demonstrate such an advantage, below we list the code that implements the NFP-KS algorithm as a comparison. The implementation can actually be done using the toolbox, by manually inserting C++ codes that are used to solve the inner problem (as a reminder, at each time step, the NFP-KS algorithm solves for consistent future endogenous states in the inner nest fixing values of policy variables, and solves for policy variables in the outer nest.)

% Parameters
cinclude('rtnewt.h');
cinclude('stdlib.h');
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 = 8;
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.05,1.05,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.05 1.0;
inbound c2 0.05 1.0;
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 3 adaptive(1.5);
inbound pb 0 3 adaptive(1.5);
% Other policy variables
var_aux equity_premium w1n[8];
% Interpolation variables for policy and state transitions
var_interp ps_future c1_future c2_future;
initial ps_future 0.0;
initial c1_future w1.*d+eta1;
initial c2_future (1-w1).*d+1-eta1;
ps_future = ps;
c1_future = c1;
c2_future = c2;

model;
  % Transform bond back
  b1p = nb1p + Kb;
  b2p = nb2p + Kb;

  vector w1nn[8];
  cxx;
  for (int i_shock=1; i_shock<=8; i_shock++)
  {
      auto eval_consis_resid = [&] (adouble w1_next)
      {
          auto ps_next = ps_future_adouble(i_shock, w1_next);
          auto total_share_value = ps_next + d_GRID(i_shock);
          auto consis_resid = (s1p*total_share_value + b1p/g_GRID(i_shock)) / total_share_value - w1_next;
          return consis_resid;
      };
      
      
      const double delta = 1e-6;
      const double tol_eq = 1e-6;
      auto funcd_eval_consis_resid = [&] (adouble w1_next, adouble* df)
      {
          auto consis_resid = eval_consis_resid(w1_next);
          auto w1_next_plus = w1_next+delta;
          auto consis_resid_plus = eval_consis_resid(w1_next_plus);
          *df = (consis_resid_plus - consis_resid) / (delta);
          return consis_resid;
      };
      
      const double w1n_max = 1.5;
      const double w1n_min = -0.5;
      adouble w1n0 = 0.5;
      w1nn(i_shock) = rtnewt(funcd_eval_consis_resid, w1n_min, w1n_max, w1n0, tol_eq);
  }
  endcxx;
  w1n' = w1nn'; % Return the solution

  % Evaluate interpolation
  [psn',c1n',c2n'] = GDSGE_INTERP_VEC'(w1n');
  % Calculate expectations that enter the Euler Equations
  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)};
  es1 = es1_future*c1^gamma/ps;
  es2 = es2_future*c2^gamma/ps;
  eb1 = eb1_future*c1^gamma/pb;
  eb2 = eb2_future*c2^gamma/pb;

  % 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
  w1_consis' = (s1p*(psn'+d') + b1p/g')/(psn'+d') - w1n';
  % 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;        % Normalized by total budget
    budget_2/(1-w1);
  end;
end;

simulate;
  num_periods = 10000;
  num_samples = 6;
  initial w1 0.5;
  initial shock 1;
  var_simu c1 c2 ps pb equity_premium ms1 mb1;
  w1' = w1n';
end;

As shown, Line 57-87 of the modified code is used to solve the inner problem, and codes for the equation solver need to be “included” (by Line 1 and 2, download the included file here). The cxx; block is essentially a hack implemented by the toolbox to insert a block of raw C++ codes. The block here defines the inner problem (the inner lambda function eval_consis_resid defined starting at Line 61), which is a single-unknown equation for the future state variable given current policies, and calls the solver to solve the problem, for each future state. The Newton solver, used here, achieves high performance but requires a delicate choice of the initial point for the inner problem. The choice of the initial point of 0.5, guarantees that the solutions can be found across the state space and time iterations, with values of policy variables not too distant from the equilibrium ones.

Although doable, inserting add-hoc C++ codes like this tends to be inconvenient and err-prone, and thus not recommended to regular researchers. This example makes it clear why designing an interface that accommodates add-hoc nested algorithms is difficult, whereas it can be made with the STPFI algorithm.

Finally, the existence of the nested problem breaks the efficient automatic differentiation for calculating the Jacobian Matrix of the outer problem, since the nested problem needs to use an iterative solver that is not amenable to efficient automatic differentiation. As shown in Table R1 and discussed above, the inner problem requires substantially more number of evaluations of consistency equation (eval_consis_resid), which require repeatedly calling the interpolation procedure to evaluate future share prices at off-equilibrium wealth shares (Line 63). This makes the algorithm not as performanent as the STPFI.

What’s Next?

Through this simple example, you understand the power of the toolbox and all the essential ingredients to solve a modern macro model.

For the time-iteration algorithm to work robustly, a crucial step is to define the starting point of the iteration properly. A candidate that delivers good theoretical property and proves to be numerically stable is to start from a last-period problem, so the algorithm can be viewed as taking the limit of the solution from finite-horizon iterations (Duffie et al (1994); Magill and Quinzi (1994); Cao (2020)).

The last-period problem has been so far trivial in the RBC model and Heaton and Lucas (1996) examples, but could turn out to be more complex and requires to define a different system of equations than the main model; block. Also, at the boundary of the state space, the equilibrium may be characterized by a different system of equations, and such boundary conditions turn out to be necessary to solve many models robustly (for example, consumption might be exactly zero at the boundary, violating the Inada condition). To see how these issues are addressed very conveniently in the toolbox, see example Cao and Nie (2017).

Or you can directly proceed to the toolbox’s API.