# 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

subject to

and

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:

We use the financial wealth share

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

The financial wealth share is rewritten as

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}\):

where \(\hat{\mu}^{i,s}_t\) and \(\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:

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

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.

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 | ```
% 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.

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