Applying the Explicit Aggregation Algorithm toHeterogeneous Agent Models in Continuous Time*
Masakazu Emoto

Takeki Sunakawa

First draft: September 2020
This version: April 2021

Abstract

This paper applies the explicit aggregation (XPA) algorithm to the standard heterogeneous agent model with aggregate uncertainty in continuous time. We find that
the XPA algorithm is faster in solving the model than the Krusell–Smith algorithm,
because the XPA algorithm does not rely on simulations within the algorithm to solve
the model. The XPA algorithm is more accurate than the perturbation method whenaggregate uncertainty is large.

Keywords: Continuous Time, Heterogeneous Agent Models, Explicit Aggregation Algorithm.
JEL codes: C63; D52*
We thank Eric Young (the editor), anonymous associate editor and referee, Takashi Kamihigashi, Soyoung Lee,
Tamotsu Nakamura, Galo Nuño, and Haruki Shibuya for comments and suggestions. The programming code used in
the paper is available at https://github.com/Masakazu-Emoto/XPA-in-Continuous-Time.

Graduate School of Economics, Kobe University, 2-1 Rokko-dai, Nada, Kobe, 657-8501; Email:masakazu.emoto@gmail.com
Graduate School of Economics, Hitotsubashi University; Email: takeki.sunakawa@gmail.com.

1

1

Introduction

There is more interest in heterogeneous agent macro models than ever before. Recent studies
such as Achdou et al. (2017) and Ahn et al. (2018) apply newly developed numerical methods
to solve heterogeneous agent models in continuous time. In particular, Ahn et al. (2018) study
heterogeneous agent models with aggregate uncertainty in continuous time and are able to solve
the model quickly. However, there are some challenges to their approach. First, their method solves
the model using a linear approximation and thus does not capture nonlinear effects of aggregate
uncertainty. Second, because of the linearization, the accuracy of solving the model is significantly
compromised when aggregate uncertainty is large and/or the model itself is highly nonlinear.
We present an alternative numerical method to address the issues mentioned above. We introduce the explicit aggregation (XPA) algorithm of den Haan and Rendahl (2010) into the standard
heterogeneous agent model with aggregate shocks of Krusell and Smith (1998) in continuous time.
The XPA algorithm obtains the forecasting rule for aggregate capital from the individual saving
function without relying on simulations within the algorithm.1 Then we compare the XPA algorithm in terms of accuracy and efficiency with the Krusell–Smith (KS) algorithm using simulations
and the Reiter–Ahn (REITER) algorithm using perturbations around the stationary distribution
without aggregate shocks.2
We find that, compared with the KS algorithm, the XPA algorithm is faster than the KS algorithm in solving the standard Krusell–Smith model (5.7 seconds vs. 93.0 seconds in our example)
at the same level of accuracy. Compared with the REITER algorithm, the XPA algorithm is moreaccurate than the REITER algorithm. Although the REITER algorithm can solve the model in a
fraction of second, its accuracy is deteriorated especially when aggregate uncertainty is large.
Our study is closely related to at least two areas of research. One is the literature on the
XPA algorithm. den Haan and Rendahl (2010) is the first paper to apply this method to thestandard heterogeneous agent model in Krusell and Smith (1998). Sunakawa (2020) shows that it
is possible to apply their approach to some other heterogeneous agent models such as Khan and
1

Simulations are a part of the KS algorithm, whereas they are not in the XPA algorithm. See also Section 3.4 fora discussion.
2
Although we focus on these algorithms, some other useful algorithms are also found in den Haan et al. (2010).

2

Thomas (2003, 2008) and Krueger et al. (2016) in discrete time. Therefore, it is straightforward
to apply the XPA algorithm to more empirically plausible models even in continuous time. To thebest of our knowledge, the present paper is the first to apply the XPA algorithm to the standard
heterogeneous agent model with aggregate uncertainty in continuous time. Specifically, we use aforecasting rule of K̇t in continuous time, whereas den Haan and Rendahl (2010); Sunakawa (2020)
use that of Kt+1 in discrete time. One of our contributions is that this alternative forecasting rulein continuous time also works with the XPA algorithm.
The other is the research on methods to solve heterogeneous agent models with aggregate shocks
in continuous time. The pioneering research in this area is Ahn et al. (2018), which is also thefirst to solve the Krusell and Smith (1998) model in continuous time. They adapt the perturbation
method originally developed by Reiter (2009) to heterogeneous agent models with aggregate shocks
in continuous time.3 Fernández-Villaverde et al. (2019a) propose a neural-network algorithm tosolve heterogeneous agent models with aggregate shocks in continuous time. Our algorithm, unlike
Fernández-Villaverde et al. (2019a) and the standard KS algorithm, solves the model without
using simulations. Furthermore, whereas the REITER algorithm solves the model using a linear
approximation, our algorithm solves the model nonlinearly so as to capture nonlinear effects ofaggregate uncertainty.4 This is especially important when we look at, for example, the effect ofaggregate uncertainty on the stochastic steady state (Fernández-Villaverde et al., 2019a).5 As wewill discuss later, the accuracy of solving the model is much higher than what is reported in Ahn
et al. (2018). Our results also hold for different degrees of persistence of the aggregate shock.
The paper consists of the following sections. In Section 2, we apply the XPA algorithm, as wellas the KS and REITER algorithms, to the Krusell and Smith (1998) model in continuous time.
In Section 3, we compare the results of the three algorithms, XPA, KS, and REITER, in terms of3
Reiter (2010a); Winberry (2018) further develop a method to reduce the dimension of the state space by projecting
the distribution onto principal components. Bayer and Luetticke (2020) and Childers (2018) also suggest novel
approaches using linearization.
4
Reiter’s (2010b) backward induction method can also be applied to solve heterogeneous agent models nonlinearly.
The method is applied to the stochastic overlapping generations model with aggregate uncertainty of Khan (2017);
Kim (2018). Okahata (2018) demonstrates that the method can also be merged with the continuous-time methods.
5
Perfect foresight methods (used in e.g., McKay et al. (2016); Kaplan et al. (2018)) are also unable to capture thiseffect. Recently, Schaab (2021) has done such an analysis of continuous time models with aggregate uncertainty byusing a sparse grid method developed by Schaab and Zhang (2021).

3

accuracy and efficiency. Finally, Section 4 concludes.

2

Algorithms

We apply the XPA algorithm first developed by den Haan and Rendahl (2010) to the Krusell andSmith (1998) model in continuous time as studied by Ahn et al. (2018). We choose the Krusell–Smith model as it is known as one of the most popular heterogeneous agent models with aggregate
uncertainty. As the model is well known, we defer the details of the model to Appendix A.

KS and XPA algorithms

Both the XPA and KS algorithms require two types of calculations,the inner loop and the outer loop. The inner loop calculation is common between the XPA and KS
algorithms. In continuous time models, a finite difference method is used to solve the Hamilton–
Jacobi–Bellman (HJB) equation as in Achdou et al. (2017) for the policy function of householdsavings. See Appendix A for details.6
In the outer loop, the policy function s(a, z, Kt , Zt ) obtained in the inner loop is used to obtainthe forecasting rule (perceived law of motion) for the next period’s aggregate capitalK̇t = Γ(Kt , Zt ).

(1)

The XPA and KS algorithms are different in how they represent Equation (1). Let gt (a, z) be thejoint distribution of wealth and productivity. In the XPA algorithm, the forecasting rule is obtainedby aggregating the policy function with the distribution explicitly:

ΓXPA (Kt , Zt ) =

XZ
z

=

s(a, z; Kt , Zt )gt (a, z)da,aX{s(Kt (z), z, Kt , Zt ) + ξ(z)}φ(z),(2)

zwhere Kt (z) is capital conditioned on labor productivity, φ(z) =6Rgt (a, z)da is the proportion of

Fernández-Villaverde et al. (2019b) solve the standard Krusell–Smith model in continuous time using the KS
algorithm. We extend their original programming code especially for the XPA algorithm. See Appendix C and D forthe details of computation.

4

households with z, and ξ(z) is for correcting the biases due to Jensen’s inequality. These objects
can be calculated immediately when we compute the steady state. We just need to evaluate thepolicy function at a = Kt (z). In Appendix B, we explain the details of the XPA algorithm following
den Haan and Rendahl (2010) and Sunakawa (2020).
In the KS algorithm, taking the policy function as given, we simulate the model to obtainthe total factor productivity (TFP) sequence Zt and the mean of the wealth distribution in the
P R
next period K̇t = z a s(a, z; Kt , Zt )gt (a, z)da at each point in time. Then the forecasting rule isobtained by estimating the following forecasting rule using the simulated sequence of {Kt , Zt }.

ΓKS (Kt , Zt ) = β0 + β1 ln Kt + β2 Zt ,(3)

where (β0 , β1 , β2 ) are coefficients estimated by ordinary least squares.
REITER algorithm

The REITER algorithm linearizes the model around the stationary distri-bution obtained from a model without aggregate shocks and uses the system of linearized equations
to solve for the dynamics of the economy in the event of aggregate shocks. As a result, the whole
distribution is used for the perceived law of motion (forecasting rule). See Ahn et al. (2018) for
more details.
The law of motion in the REITER algorithm is obtained as follows. Let gt be a vector of
the wealth–productivity density gt (a, z) discretized by grid points of (a, z) and ĝt = gt − g be its
deviation from the stationary distribution. Then we havedĝt
= B̃g ĝt + B̃Z Zt
dt(4)

where B̃g and B̃Z are the matrices obtained by solving the system of linearized equations.7 GivenRthe sequence of {gt }, the sequence of aggregate capital is obtained by Kt = agt (a, z)dadz.
7

This equation corresponds to Equation (17) in Ahn et al. (2018).

5

3

Numerical results

We compare the numerical results from three algorithms; XPA algorithm, KS algorithm, and
REITER algorithm. The parameter values and detailed settings are found in Appendix F.8

3.1

Forecasting rules and simulation paths

In Figure 1, we show the forecasting rule K̇t (Kt , Zt ) in KS and XPA.9 We consider two cases. One
is with the standard deviation of the disturbance to TFP being equal to σ = 0.007, which is thestandard value. The other is with σ = 0.05, which is very large compared to the benchmark caseand is intended to show how large differences can be in the extreme case. As is clear from thefigure, when σ = 0.007, there is no significant difference in the forecasting rules obtained by eachalgorithm. Each forecasting rule is characterized by a decreasing function with respect to capital
Kt and an increasing function with respect to TFP Zt . Even when σ = 0.05, the forecasting rules
in XPA and KS are very similar to each other, although they are slightly different when the valuesof (Kt , Zt ) are far away from the ones in the stationary distribution.
8
We use the code developed by Ahn et al. (2018) for REITER. All the codes are modified to use the same set ofparameters simply for the purpose of comparison.
9
As in Fernández-Villaverde et al. (2019b), we use K̇(K, Z) = 0 for all the grid points of (K, Z) as an initial guess
for the forecasting rule in both KS and XPA.

6

Figure 1: Forecasting rulesa. σ = 0.007
0.2

0

0.4

0.4

0.2

0.2

0

-0.2

0

-0.2
-0.4

-0.2
30

35

40

30

35

40

30

35

40

30

35

40

b. σ = 0.05
-0.2

-0.4

-0.6

0.4

0.8

0.2

0.6

0

0.4

-0.2

-0.8
30

35

40

0.2
30

35

40

In Figure 2, we show the results of the simulation path derived from the nonlinear model for
10,000 periods in XPA (compared with that in KS), and the simulation path derived from thelinearized model (i.e., Equation (4)) in REITER (compared with that in KS).10 It is clear from theleft panel of Figure 2 that the capital path obtained by each algorithm is very similar to each otherwhen σ = 0.007. In contrast, in the right panel of Figure 2, the capital path obtained by REITER
is very different from those by other algorithms when σ = 0.05.
The capital path in XPA tend to have an upward bias compared with that in KS. The more σincreases, the larger the bias is (about 1% when σ = 0.05). However, each algorithm is accurate interms of the Den Haan errors as shown later in Figure 3.
10

In all the algorithms, the time interval dt = 0.25 used in the simulation is the same as the time interval of
Ahn et al. (2018). In KS and XPA, we use a bilinear interpolation for aggregate capital and TFP to calculate thewealth–productivity density at each point of time. See Appendix D.

7

Figure 2: Simulation pathsa. σ = 0.007
1
0.8
0.6
0.4
0.2
0
-0.2
-0.4
-0.6
XPA-KS (%)
REITER-KS (%)

-0.8
-1
0

2000

4000

6000

8000

10000

b. σ = 0.05
9
8
7
6
5
4
3
2
1
XPA-KS (%)
REITER-KS (%)

0
-1
0

2000

4000

6000

8

8000

10000

3.2

Efficiency

Comparing the results of XPA and KS in Table 1, we can see that XPA solves the model much
faster than KS because it does not use simulations. Comparing the computation times of XPA and
REITER, we can see that REITER is also faster than XPA. However, the difference between thetwo is not that large (XPA: 5.7 seconds vs. REITER: 0.3 seconds).11
Table 1: Computation time
Algorithm
XPA (Our method)
Krusell-Smith
REITER (Ahn et al., 2018)

Computation time
5.714 sec
92.964 sec
0.296 sec

Notes: Computations are done on a workstation with Intel Xeon E5–2699 v4 (2.20Ghz) and 32GB
RAM using MATLAB R2020b. Computation time is the average of 10 runs and excludes the timespent for simulations.

3.3

Accuracy

We check the accuracy of solving the model using the Den haan errors proposed by den Haan (2010).
In the present paper, we simulate two sequences of 10,000 time periods with each algorithm. One
is {K̃t }t∈[0,T ] obtained only from the forecasting rule for each algorithm (Equations (2)–(4)), andthe other is {Kt∗ }t∈[0,T ] obtained from the fully nonlinear model including the the household HJBequation and the Kolmogorov Forward equation.12 Then we use the sequences of {Kt∗ }t∈[0,T ] and{K̃t }t∈[0,T ] to measure the Den haan errors
AXεM≡ 100 · max | ln K̃t − ln Kt∗ |.
DH
t∈[0,T ]Figure 3 summarize the Den haan errors for each algorithm. It is clear that when σ is small,there is not much difference in the Den haan errors. However, when σ is large, the Den haan errors
11

Note that we use MATLAB and no parallelization, so the gap might be smaller when we use a faster language
and/or parallelization.
12
As in Ahn et al. (2018); Bayer and Luetticke (2020), to obtain the sequence of {Kt∗ } in REITER, we solve forthe optimal savings plan s∗t (aj , zi ) at every pair of the wealth and productivity values (aj , zi ) in the histogram underthe equilibrium factor prices (wt , rt ) in every period t. See Appendix E.

9for REITER are considerably larger than those for XPA and KS. Furthermore, the Den Haan errorsof XPA are smaller than those of KS.1314
Figure 3: Comparison of maximum of Den Haan errors
10
9
8
7
6
5
4
3
2
1
0
1

3.4

1.5

2

2.5

3

3.5

4

4.5

5

Discussion

Time for simulation

Simulations are a part of the KS algorithm, whereas they are not in the

XPA algorithm. Specifically, simulations are done in KS to update the forecasting rules in the outerloop. The outer loop in XPA is much faster than in KS, as simulations are unnecessary in XPA.
The computation time in Table 1 is solely for solving the model and includes time for simulations
for the KS algorithm only. Having said that, to get the simulated paths in Figure 2, simulationsin XPA are as slow as in KS. Indeed, taking the converged forecasting rule by each algorithm asgiven, simulation time in XPA is almost identical to that in KS (about 10 seconds to compute asequence of 10,000 time periods in Figure 2).
13

However, there is a possibility of using non-parametric forecasting rules or forecasting rules with higher-orderterms for KS to achieve higher accuracy.
14
We also confirm that our results hold when the persistence of TFP, 1 − µ, is lowered from µ = 0.25 to µ = 0.5
or µ = 0.75. See Appendix F.

10Market clearing It is nontrivial to handle market clearing conditions in XPA or KS. Sunakawa(2020) shows a detailed explanation how to apply XPA in discrete time with market clearing inlumpy investment models developed in Khan and Thomas (2003, 2008).15 XPA can solve the model
faster with a similar degree of accuracy to KS. XPA (or KS) in continuous time should work ina similar way, although it takes more time to solve the model than REITER does. In REITER,
however, market clearing conditions hold only at the stationary distribution.

4

Conclusion

In this paper, we apply the explicit aggregation (XPA) algorithm proposed by den Haan andRendahl (2010) to the standard heterogeneous agent model with aggregate uncertainty (Krusell
and Smith, 1998) in continuous time. Using XPA, we can get an accuracy close to that of KS witha similar speed to that of REITER. Future research will show that XPA can also be useful in models
such as ones with uncertainty shocks (Bloom et al., 2018) where linearization-based methods fail.

References
Achdou, Y., J. Han, J.-M. L. Lasry, P.-L. L. Lions, and B. Moll (2017): “Income and
Wealth Distribution in Macroeconomics: A Continuous-Time Approach,” NBER Working Papers
23732, National Bureau of Economic Research, Inc.
Ahn, S., G. Kaplan, B. Moll, , T. Winberry, and C. Wolf (2018): “When Inequality
Matters for Macro and Macro Matters for Inequality,” NBER Macroeconomics Annual, 1–75.
Bayer, C. and R. Luetticke (2020): “Solving heterogeneous agent models in discrete time withmany idiosyncratic states by perturbation methods,” Quantitative Economics, 11, 1253–1288.
Bloom, N., M. Floetotto, N. Jaimovich, I. Saporta-Eksten, and S. J. Terry (2018):“Really Uncertain Business Cycles,” Econometrica, 86, 1031–1065.
15

In these models, taking the shadow price of consumption as given, we obtain the aggregate consumption andimplied new shadow price. We solve for the fixed point of the shadow price at each grid point in the outer loop in
XPA.

11

Childers, D. (2018): “Solution of Rational Expectations Models with Function Valued States,”
Working Paper.
den Haan, W. J. (2010): “Comparison of Solutions to the Incomplete Markets Model with
Aggregate Uncertainty,” Journal of Economic Dynamics and Control, 34, 4–27.
den Haan, W. J., K. L. Judd, and M. Juillard (2010): “Computational suite of models withheterogeneous agents: Incomplete markets and aggregate uncertainty,” Journal of Economic
Dynamics and Control, 34, 1–3.
den Haan, W. J. and P. Rendahl (2010): “Solving the Incomplete Markets Model with Aggregate Uncertainty Using Explicit Aggregation,” Journal of Economic Dynamics and Control,
34, 69–78.
Fernández-Villaverde, J., S. Hurtado, and G. Nuño (2019a): “Financial Frictions and theWealth Distribution,” NBER Working Papers 26302, National Bureau of Economic Research,Inc.
——— (2019b): “Solving the Krusell-Smith (1998) model,” Manuscript.
Kaplan, G., B. Moll, and G. L. Violante (2018): “Monetary Policy According to HANK,”
American Economic Review, 108, 697–743.
Khan, A. (2017): “Large Recessions in an Overlapping Generations with Unemployment,” 2017
Meeting Papers 1559, Society for Economic Dynamics.
Khan, A. and J. K. Thomas (2003): “Nonconvex Factor Adjustments in Equilibrium Business
Cycle Models: Do Nonlinearities Matter?” Journal of Monetary Economics, 50, 331–360.
——— (2008): “Idiosyncratic Shocks and the Role of Nonconvexities in Plant and Aggregate
Investment Dynamics,” Econometrica, 76, 395–436.
Kim, H. (2018): “Inequality, Portfolio Choice, and the Business Cycle,” Working Paper.

12

Krueger, D., K. Mitman, and F. Perri (2016): “Macroeconomics and Household Heterogeneity,” in Handbook of Macroeconomics, ed. by J. B. Taylor and H. Uhlig, Elsevier, vol. 2,
843–921.
Krusell, P. and A. A. J. Smith (1998): “Income and Wealth Heterogeneity in the Macroeconomy,” Journal of Political Economy, 106, 867–896.
McKay, A., E. Nakamura, and J. Steinsson (2016): “The Power of Forward Guidance Revisited,” American Economic Review, 106, 3133–3158.
Okahata, N. (2018): “An Alternative Solution Method for Continuous-Time Heterogeneous Agent
Models with Aggregate Shocks,” Mimeo.
Reiter, M. (2009): “Solving heterogeneous-agent models by projection and perturbation,” Journalof Economic Dynamics and Control, 33, 649–665.
——— (2010a):“Approximate and Almost-Exact Aggregation in Dynamic Stochastic

Heterogeneous-Agent Models,” Economics Series 258, Institute for Advanced Studies.
——— (2010b): “Solving the incomplete markets model with aggregate uncertainty by backward
induction,” Journal of Economic Dynamics and Control, 34, 28–35.
Schaab, A. (2021): “Micro and Macro Uncertainty,” Mimeo.
Schaab, A. and A. Zhang (2021): “Adaptive Sparse Grid Methods for Differential Equations in
Economics,” Mimeo.
Sunakawa, T. (2020): “Applying the Explicit Aggregation Algorithm to Heterogeneous Macro
Models,” Computational Economics, 55, 845–874.
Winberry, T. (2018): “A Method for Solving and Estimating Heterogeneous Agent Macro Models,” Quantitative Economics, 9, 1123–1151.
Young, E. R. (2005): “Approximate Aggregation,” Mimeo, University of Virginia.

13

Appendix (not for publication)A

The Krusell–Smith model in continuous time

A.1

Environments

Households
There is a continuum of households with a normalized fixed mass indexed by j ∈ [0, 1]. Households
face idiosyncratic uncertainty regarding labor productivity and the borrowing constraint ajt ≥ 0,
where ajt is the value of asset holdings for household j in period t. There are two states oflabor productivity zjt for each household, ze and zu , which follow the Poisson process with arrival
rates λe and λu . ze shows that the household is employed and zu indicates that the household isunemployed. If the household is employed, she/he receives labor income after taxation (1 − τ )wt .
When the household is unemployed, she/he receives unemployment insurance bwt financed by thelabor income tax.
Each household j ∈ [0, 1] chooses their consumption and savings (cjt , ajt ) in each period t ≥ 0
to maximize their expected life-time utility by taking the wage rate wt and the interest rate rt asgiven.
"Z
vj0 = max E0
0

s.t.

∞

#
1−γ
c
jt
e−ρt
dt ,
1−γ

dajt = (rt ajt + (1 − τ )zjt wt + (1 − zjt )bwt − cjt )dt, ajt ≥ 0,zjt ∈ {ze , zu }, ze = 1, zu = 0.

The instantaneous utility function is of constant-relative-risk-aversion form, ρ is the rate of time
preference, and γ is the degree of relative risk aversion. The population of households with a pairof wealth and productivity levels is time-variant and given by gt (a, z). For notational convenience,
we drop time subscripts t from variables at the individual level hereafter.

14

Firm
The representative firm produces the final good Yt using capital Kt and labor Lt . The production
function is Cobb–Douglas form

Yt = eZt Ktα L1−α,twhere α is the capital share. Zt is the logarithm of total factor productivity (TFP) following theOrnstein–Uhlenbeck process

dZt = µ(Z̄ − Zt )dt + σdWt ,Z̄ = 0,where dWt follows a Wiener process. 1 − µ is the persistence of TFP and σ is the volatility of the
TFP. This process is similar to an AR(1) process in discrete time.
The wage rate and the interest rate are obtained from the first-order conditions for the profit
maximization problem as follows:

Zt α−1 1−α
Lt − δ,wt = (1 − α)eZt Ktα L−αt , rt = αe Kt

where δ is the depreciated rate of capital.

Government
The government imposes a tax on labor income to finance unemployment compensation. The
government’s budget is balanced as below

τ wt φ(ze ) = bwt φ(zu ).

That is, the government’s tax revenue from labor income is equal to the government’s expenditureR
to finance unemployment insurance. φ(ze ) = gt (a, ze )da is the share of employment and φ(zu ) =15R

gt (a, zu )da is the share of unemployment in the economy.

Equilibrium
An equilibrium in this economy is consist of a set of prices {wt , rt }, quantities {Kt , Lt } and adensity {gt (a, z)} such that
1. Given wt , rt , the policy functions for consumption and savings, c∗t (a, z) and s∗t (a, z), are thesolution of the the Hamilton–Jacobi–Bellman (HJB) equation

ρvt (a, z) = max u(c) + ∂a vt (a, z)(rt a + (1 − τ )zwt + (1 − z)bwt − c)c+ λz (vt (a, z 0 ) − vt (a, z)) +1
Et [dvt (a, z)],dt

(5)where u(·) is the utility function and vt (·) is the value function in period t, which depends on
a particular realization of the aggregate state (gt (a, z), Zt ).
2. The sequence of {gt (a, z)} is the solution of the Kolmogorov Forward (KF) equation
dgt (a, z)= − ∂a [st (a, z)gt (a, z)] − λz gt (a, z) + λz 0 gt (a, z 0 )dt(6)where st (a, z) = s∗t (a, z) is the optimal saving policy function corresponding to the household
optimization problem.
3. The wage rate and the interest rate are given by
Zt α−1 1−α
wt = (1 − α)eZt Ktα L−α
Lt − δ.
t , rt = αe Kt

4. The government’s budget constraint is satisfied asτ wt φ(ze ) = bwt φ(zu ),where φ(ze ) =Ra gt (a, ze )da,φ(zu ) =Ra gt (a, zu )da are time-invariant distributions of pro-16

ductivity.
5. The capital and the labor markets clear

Kt =

XZ
z

agt (a, z)da, Lt =

a

XZ
z

zgt (a, z)da.

a

HJB equation in the XPA and KS algorithms
In the XPA and KS algorithms, we assume approximate aggregation (Young, 2005) so that thewealth–productivity density {gt (a, z)} as a state variable is approximated by the mean Kt so thatthe aggregate state is (Kt , Zt ). In this case, the value function is denoted by v(a, z, Kt , Zt ) and the
HJB equation can be written asρv(a, z, Kt , Zt ) = maxcc1−γ − 1+ va (a, z, Kt , Zt )ȧ + λz (v(a, z 0 , Kt , Zt ) − v(a, z, Kt , Zt ))1−γ+ vK (a, z, Kt , Zt )K̇t
+ vZ (a, z, Kt , Zt )(−µZt ) +σ2vZZ (a, z, Kt , Zt )2

(7)subject to the budget constraint ȧ = rt a + (1 − τ )zwt + (1 − z)bwt − c and the forecasting ruleK̇t = Γ(Kt , Zt ) for the policy function of household savings, s(a, z, Kt , Zt ). va , vK , and vZ are thefirst order derivatives of the value function in terms of a, K, and Z and vZZ is the second orderderivative of the value function in terms of Z.

B

Details of the XPA algorithm

We provide details of the XPA algorithm following den Haan and Rendahl (2010) and Sunakawa(2020). In contrast with the KS algorithm, the XPA algorithm calculates the forecasting rules
without simulations within the algorithm. First, we rewrite the wealth–productivity density gt (a, z)

17using the conditional probability as

gt (a|z) =R

gt (a, z)
gt (a, z)=φ(z)gt (a, z)da⇔ gt (a, z) = gt (a|z)φ(z)where φ(z) =R

gt (a, z)da is equal to the proportion of households with labor productivity z in theeconomy.16 Then we can rewrite the forecasting rule using the conditional distribution of wealthasXZs(a, z; Kt , Zt )gt (a, z)da

X Z
≈
s
agt (a|z)da, z; Kt , Zt φ(z)X=s(Kt (z), z; Kt , Zt )φ(z)K̇(Kt , Zt ) =where Kt (z) =Ragt (a|z)da is capital conditioned on labor productivity z. We assume that theRhousehold’s policy function s(a, z; Kt , Zt ) is linear at a = Kt (z) so that s(a, z; Kt , Zt )g(a|z)da ≈s (Kt (z), z; Kt , Zt ) holds. We compute Kt (z) by the following equations:
Kss (z)Kss (z)
Kt (z) = ψ(z)Kss , ψ(z) ≡
,
=P R
Kss
z agss (a, z)daRZ
agss (a, z)daKss (z) =agss (a|z)da =,φ(z)where Kss and gss are capital and wealth distribution at the steady state without aggregate uncertainty. We assume that ψ(z) is constant over aggregate fluctuations. Note that the ratio of thecapital conditioned on z to aggregate capital, ψ(z) = Kss (z)/Kss , can be easily obtained in thesteady-state calculation.
Moreover, following den Haan and Rendahl (2010), we conduct bias correction. As we assume
that the policy function is linear at a = Kt (z), there may be bias in the forecasting rule from16
We assume the share of employment φ(ze ) (= 1 − φ(zu )) is time-invariant, although it is straightforward to makethe employment measure be time-variant and depend on aggregate uncertainty as in Krusell and Smith (1998). See
also Sunakawa (2020).

18

Jensen’s inequality. We compute the steady-state counterparts to correct the bias:

ξ(z) = K̇ss (z) − sss (Kss (z), z),RZ
sss (a, z)gss (a, z)daK̇ss (z) =sss (a, z)gss (a|z)da =.
φ(z)Again, ξ(z) can be computed at negligible cost in the steady state. Finally, we can write theforecasting rule as follows

K̇(Kt , Zt ) =X{s(Kt (z), z, Kt , Zt ) + ξ(z)}φ(z).

(8)

zThat is, to obtain the forecasting rule, we just need to evaluate the policy function at a = Kt (z).
The value of Kt (z) may not be on the grid of a, so we use linear interpolation.
We have two important assumptions to derive Equation (8). One is that the policy function islinear at a = Kt (z). The other is that the ratio of capital conditioned on z to aggregate capital,
ψ(z) = Kt (z)/Kt , is constant over aggregate fluctuations. The first assumption is acceptableespecially in the household consumption-saving problem in Krusell and Smith (1998), in which thesaving function is almost linear except for near the origin where poor households face borrowing
constraints. Sunakawa (2020) shows that the algorithm also works with nonlinear (S, s)–type policy
functions such as in Khan and Thomas (2003, 2008).
Also, Sunakawa (2020) shows that we can incorporate more than two states for the idiosyncratic
shock by using the epsilon-indexed aggregation. Note that the original analysis by den Haan andRendahl (2010) has dealt with only two states for the idiosyncratic shock in the Krusell and Smith(1998) model (i.e., employment and unemployment) and they do so by increasing the numberof aggregate states indexed by employment status, K(z) for z ∈ {e, u}.17 Instead, we do notneed to increase the number of aggregate states with the epsilon-indexed approach, as shown inEquation (8). The present paper shows that we can also use the epsilon-indexed aggregation in thecontinuous-time framework as in Sunakawa (2020).
17

The policy function becomes the form of s(k, z, Kt (e), Kt (u), Zt ).

19

Summary of XPA Algorithm
In summary, we perform computations to solve the Krusell–Smith model with the XPA algorithm
as follows. As mentioned above, the XPA algorithm is fast because it does not use any simulations.
1. Compute the stationary distribution without aggregate uncertainty to obtain the conditional
capital ratio ψ(z) and the bias correction term for correcting the forecasting rule ξ(z).
2. (Inner loop) Solve the HJB equation for the policy function taking the forecasting rule asgiven. See Appendix C for details.
3. (Outer loop) Compute the forecasting rule without simulations taking the policy function as
given. The bias correction is also done.
4. Repeat steps 2–3 until the forecasting rule converges.

C

Solving the HJB equation with the upwind scheme

Fernández-Villaverde et al. (2019b) solve the standard Krusell–Smith model in continuous time using the KS algorithm. We extend their original programming code in the following two dimensions:(i) We use the upwind scheme not only individual wealth, a, but also K and Z, which is necessary
for the XPA algorithm.18 (ii) We exclude the direct effect of aggregate variables K and Z on thetransition matrix for simulation.19
We discretize the individual wealth and productivity by grid points (aj , zi ) for i = 1, 2 andj = 1, ..., J so that gt,i,j = gt (aj , zi ) holds. We also discretize the aggregate capital and productivity
by grid points (Kl , Zm ) for l = 1, ..., L and m = 1, ..., M so that vi,j,l,m = v(aj , zi , Kl , Zm ) holds.
n
We are going to solve the HJB equation (7) with an iterative method. Let vi,j,l,m
be the valuefunction at n–th iteration. Then we have (note that we use the implicit method so that t + 1
18

As shown in Appendix D, the upwind scheme for K and Z also improves the accuracy of the KS algorithm.
In the current and next sections, we acknowledge referring to Fernández-Villaverde et al. (2019b) (hereafter
FVHN) for some derivations and notations.
19

20

variables are in the right hand side)n+1nvi,j,l,m− vi,j,l,m∆

n+1+ ρvi,j,l,m=(cni,j,l,m )1−γ − 1
1−γ

n+1
n+1
+ ∂a(F ) vi,j,l,msni,j,l,m,F 1sni,j,l,m,F >0 + ∂a(B) vi,j,l,msni,j,l,m,B 1sni,j,l,m,B <0
(F )

(B)

n+1
n+1
n+1
n+1
+ λi (v−i,j,l,m− vi,j,l,m) + ∂K vi,j,l,mhl,m 1hl,m >0 + ∂K vi,j,l,mhl,m 1hl,m <0(F )σ 2 2 n+1
∂ v
2 ZZ i,j,l,m
(9)

(B)

n+1
n+1
+ ∂Z vi,j,l,mθ(Z̄ − Zm )1Z̄−Zm >0 + ∂Z vi,j,l,mθ(Z̄ − Zm )1Z̄−Zm <0 +

where

n+1
∂a(F ) vi,j,l,m
=

n+1
n+1
vi,j+1,l,m− vi,j,l,m

(F )

∆a
n+1
n+1
vi,j,l+1,m − vi,j,l,m

(F )

∆K
n+1
n+1
− vi,j,l,mvi,j,l,m+1

n+1∂K vi,j,l,m=
n+1∂Z vi,j,l,m
=

,

n+1
∂a(B) vi,j,l,m
=

,

n+1
∂K vi,j,l,m
=

n+1
n+1
vi,j,l,m− vi,j−1,l,m

(B)

∆a
n+1
n+1
vi,j,l,m − vi,j,l−1,m

(B)

∆K
n+1
n+1
− vi,j,l,m−1
vi,j,l,mn+1, ∂Z vi,j,l,m
=
∆Z
n+1
n+1
n+1
− 2vi,j,l,m+ vi,j,l,m−1
vi,j,l,m+1
n+1=∂ZZ vi,j,l,m
,
(∆Z)2

∆Z

,
,
,

and
1/γ


1


sni,j,l,m,F = wl,m zi + rl,m aj −  (F )n∂a vi,j,l,m

1/γ


,

1


sni,j,l,m,B = wl,m zi + rl,m aj −  (B)n∂a vi,j,l,m.

Also, the optimal consumption is set ton)−γcni,j,l,m = (∂a vi,j,l,m

where

n
n
n
n
∂a vi,j,l,m= ∂a(F ) vi,j,l,m
1sni,j,l,m,F >0 + ∂a(B) vi,j,l,m
1sni,j,l,m,B <0 + ∂a v̄i,j,l,m
1sni,j,l,m,F ≤0 1sni,j,l,m,B ≥0 .

n
In the above expression, ∂a v̄i,j,l,m= (c̄ni,j,l,m )−γ where c̄ni,j,l,m is the consumption level such that21sn (ai ) = 0. Equation (9) can be rearranged asn+1
nvi,j,l,m− vi,j,l,m∆

n+1+ ρvi,j,l,m=(cni,j,l,m )1−γ − 1
1−γ
n+1
n+1
n+1
n+1
n
n
n
+ αi,j,l,mvi,j−1,l,m+ βi,j,l,mvi,j,l,m+ ξi,j,l,mvi,j+1,l,m+ λi v−i,j,l,m
n+1
n+1
n
n
+ αK,l,mvi,j,l−1,m+ ξK,l,mvi,j,l+1,m
n+1
n+1
n
n
+ αZ,mvi,j,l,m+1+ ξZ,mvi,j,l,m−1

where

n
αi,j,l,m=n=βi,j,l,m−sni,j,l,m,B 1sni,j,l,m,B <0∆a,sni,j,l,m,B 1sni,j,l,m,B <0 − sni,j,l,m,F 1sni,j,l,m,F >0
∆a

− λi

hl,m (−1hl,m <0 + 1hl,m >0 ) θ(Z̄ − Zm )(−1Z̄−Zm <0 + 1Z̄−Zm >0 )
σ2
−
−
,
∆K
∆Z
(∆Z)2
sni,j,l,m,F 1sni,j,l,m,F >0nξi,j,l,m =,∆ahl,m 1hl,m <0n=−αK,l,m,∆K
hl,m 1hl,m >0n=ξK,l,m,∆K
θ(Z̄ − Zm )1Z̄−Zm <0
σ2
n
αZ,m
=−
+
,
∆Z
2(∆Z)2
θ(Z̄ − Zm )1Z̄−Zm >0
σ2
n
ξZ,m
=
+
.
∆Z
2(∆Z)2−This equation is a system of 2 × J × L × M linear equations and can be stacked into matrix forms
as
1 n+1
(v
− v n ) + ρv n+1 = un + An v n+1
∆

22

(10)

whereAn1 ξn I Z,2 2J×L 02J×Ln
A = −

..

.





02J×L


An1,m



 ξn
 K,2,m I 2J


02J

n
Am = − 

..

.





02J



βn
 1,1,l,m



n I
αZ,1
2J×L

02J×L

An2

n I
αZ,2
2J×L

···

02J×L

n I
ξZ,3
2J×L
..
.

An3
..
.

···
..
.

02J×L
..
.

n
ξZ,M
−1 I 2J×LAnM −1

···

02J×L

n
ξZ,M
I 2J×LnαK,1,m
I 2J

02J

···

02J

An2,m

n
αK,2,m
I 2J

···

02J

n
ξK,3,m
I 2J
..
.

An3,m
..
.

···
..
.

02J
..
.

n
I 2J
ξK,L−1,m

AnL−1

02J

n
ξK,L,m
I 2J

02J×L

02J

···

···

02J×L

02J×L



02J×L



02J×L

,

..

.



n
αZ,M
−1 I 2J×L 

AnM

02J



02J



02J

,

..

.


n
αK,L−1,m I 2J 


n
AL



n
ξ1,1,l,m

0
λ1



n
n

 αn
0
λ1

 1,2,l,m β1,2,l,m ξ1,2,l,m




n
n
n
0
λ1
0
α1,3,l,m β1,3,l,m ξ1,3,l,m






..
..
..
..
.
.
.
.






n
n


0
0
α
β
0
λ
1


1,J,l,m
1,J,l,m
Anl,m = − 
,


n
n
ξ2,1,l,m00β2,1,l,m

 λ2




n
n
n


λ
α
β
ξ
0
2
2,2,l,m2,2,l,m2,2,l,m




n
n
n


λ
0
α
β
ξ
0
2


2,3,l,m
2,3,l,m
2,3,l,m


..
..
..
..


.
.
.
.




n
n
λ2
0
0
α2,J,l,m β2,J,l,m
and

23



un1





v n1







un1,m



v n1,m



















 un 
 vn 
 un 
 vn 
 2  n  2  n
 2,m  n
 2,m nu =  .  , v =  .  , um =  .  , v m =  .  ,
 . 
 . 
 . 
 . 
 . 
 . 
 . 
 . 








unM
v nM
unL,m
v nL,m

1−γ(cn
1,1,l,m )
1−γ
1−γ
(cn
1,2,l,m )
1−γ





n
v1,1,l,m












 vn



 1,2,l,m 





..
..




.
.







 n
 n

 (c1,J,l,m )1−γ  nn=ul,m = ,v
 l,m  v1,J,l,m  .
1−γ






 (cn2,1,l,m )1−γ 



 vn

 2,1,l,m 

1−γ




..
..




.
.







 n
(c2,J,l,m )1−γnv
2,J,l,m1−γGiven v n , the system can be solved for v n+1 . An iterative procedure can be applied to (10) to∗∗∗

obtain a converged v = v n where n∗ is such that v n +1 − v n

is below a very small number.

We also define the transition matrix without the derivatives with regard to K and Znβ̃1,1,l,mn
ξ1,1,l,m

0

λ1







n
n

 αn
0
λ1

 1,2,l,m β̃1,2,l,m ξ1,2,l,m




n
n
n
0
λ1
ξ1,3,l,m0
α1,3,l,m
β̃1,3,l,m






..
..
..
..
.
.
.
.






n
n

0
0
α1,J,l,m β̃1,J,l,m
0
λ1 
n


Ãl,m = − 
,


n
n
λ
0
β̃
ξ
0


2
2,1,l,m2,1,l,m




n
n
n


λ2
α2,2,l,m β̃2,2,l,m ξ2,2,l,m
0




n
n
n


λ2
0
α2,3,l,m β̃2,3,l,m ξ2,3,l,m
0




..
..
..
..

.
.
.
. 




n
n
λ2
0
0
α2,J,l,m
β̃2,J,l,m

24

wherenβ̃i,j,l,m=sni,j,l,m,B 1sni,j,l,m,B <0 − sni,j,l,m,F 1sni,j,l,m,F >0∆a− λi ,

which is used when solving the Kolmogorov Forward equation for simulation.

D

Solving the Kolmogorov Forward equation for simulation

In the outer loop in the KS algorithm and nonlinear simulation in all the algorithms to obtainthe sequence of {Kt∗ }, we solve the Kolmogorov Forward (KF) equation with an infinite differencen∗scheme. In the XPA and KS algorithms, we use a converged Ãl,m = Ãl,m at each grid point of(Kl , Zm ) in the inner loop above for the transition matrix of the wealth–productivity density. Notethat we exclude the direct effect of aggregate variables K and Z (i.e., the derivatives with regard
to K and Z) on the transition matrix for simulation.
The law of motion of the wealth–productivity density is expressed as a form of the KF equation
gt+1,i,j − gt,i,j= − ∂a st (aj , zi )gt+1,i,j − λi gt+1,i,j + λi gt+1,−i,j
∆t

(11)

This equation can be expressed as a matrix form

g t+1 − g t = ∆tΓt g t+1(12)

where g t = [gt,1,1 , gt,1,2 , ..., gt,1,J , gt,2,1 , ..., gt,2,J ]T is a vector of the wealth–productivity density at
grid points of (aj , zi ) and Γt is a transition matrix.
In the XPA and KS algorithms, we have already obtained a converged matrix in the inner loopfor Γt .20 As (Kt , Zt ) in simulation may not be on grid points of (Kl , Zm ) used in the inner loop,T

we use bilinear interpolation of matrices at the nearest four grid points. That is, Γt = Ãt is the20

See FVHN for more details.

25transpose ofÃt =(1 − ωZ,t )((1 − wK,t )Ãl∗ ,m∗ + wK,t Ãl∗ +1,m∗ )+ wZ,t ((1 − wK,t )Ãl∗ ,m∗ +1 + wK,t Ãl∗ +1,m∗ +1 )where wK,t = (Kt − Kl∗ )/(Kl∗ +1 − Kl∗ ) and wZ,t = (Zt − Zm∗ )/(Zm∗ +1 − Zm∗ ) and each of l∗ , m∗satisfies Kt ∈ [Kl∗ , Kl∗ +1 ] or Zt ∈ [Zm∗ , Zm∗ +1 ] at each period t.
Difference from FVHN FVHN solve the HJB equation with the upwind scheme only for a in
the KS algorithm. We find that using the upwind scheme not only for a but also for K and Z isnecessary for the XPA algorithm. We find that it is also important to improve the KS algorithm
in terms of the R–squared and the Den Haan errors as shown in Table 2.
When FVHN solve the KF equation in a form of (12) in the KS algorithm, they use a converged∗

Al,m = Anl,m (without tilde) including the terms of derivatives with K and Z as the transition
matrix of the wealth–productivity density. We exclude these terms as it is more consistent with thenonlinear simulation for the REITER algorithm as explained in Appendix E. However, comparedwith the upwind scheme mentioned above, we find that this has only a marginal effect.
Table 2: Den Haan errors: alternative methods for the KS algorithma. Our approach
Agg Shock σ (%)

0.01

0.1

0.7

1.0

3.0

5.0

AX
εM
DH KS (%)
EANεM
DH KS (%)
R2

0.002

0.015

0.097

0.134

0.493

0.885

0.001

0.011

0.075

0.104

0.286

0.438

0.999

0.999

0.999

0.999

0.999

0.998

b. Alternative approach (as in FVHN)
Agg Shock σ (%)

0.01

0.1

0.7

1.0

3.0

5.0

AX
εM
DH KS (%)
EANεM
DH KS (%)
R2

0.005

0.047

0.322

0.456

1.349

2.463

0.001

0.012

0.084

0.120

0.364

0.627

0.994

0.994

0.994

0.993

0.991

0.988

26

E

Den Haan errors for the REITER algorithm

As in Ahn et al. (2018); Bayer and Luetticke (2020), to obtain the sequence of {Kt∗ } in the REITER algorithm, we solve for the optimal savings plan s∗t (aj , zi ) at every pair of the wealth andproductivity values (aj , zi ) in the histogram under the equilibrium factor prices (wt , rt ) in every
period t.21 To solve for the policy function of savings, we utilize the value function obtained in thesolution with linearization.22 Then we calculate the transition matrix for the law of motion of thewealth–productivity density (i.e., Γt in (12)) from the policy function.
First, we calculate the optimal savings plan s∗t,i,j = s∗t (aj , zi ) at each grid point in period t for
the KF equation (11). Having the current aggregate state {gt,i,j } and Zt , we have

Kt =

XX
i

aj gt,i,j ,jwt = (1 − α)eZt (Kt /Lt )α ,rt = αeZt (Kt /Lt )α−1 − δ.

From the household’s budget constraint, we have savings as st,i,j = st (aj , zi ) = wt zi + rt aj − ct,i,j .
Having the value function from the linearized solution at each grid point at date t, denoted by{ṽt,i,j }, we also have the optimal value of consumption at each (aj , zi ):c∗t,i,j = c∗t (aj , zi ) = (∂a ṽt,i,j )−γ .

(13)

We use the upward scheme to approximate the derivative of the value function as well. Therefore,"s∗t,i,j,F =wt zi + rt aj −

21

#1/γ
,

(F )

∂a ṽt,i,j"s∗t,i,j,B =wt zi + rt aj −

1

1
(B)

#1/γ
,

∂a ṽt,i,j

In general, we need to solve for the equilibrium prices that clear markets. The market clearing condition for
capital is met in the first place in our case with the standard Krusell-Smith model.
22
Bayer and Luetticke (2020) solve the optimization problem by utilizing the expected value function that is
obtained in the solution with linearization.

27

whereṽt,i,j+1 − ṽt,i,j,∆aṽt,i,j − ṽt,i,j−1∂a(B) ṽt,i,j =,∆a∂a(F ) ṽt,i,j =Then, Equation (11) can be rearranged as
gt+1,i,j − gt,i,j=αt,i,j gt+1,i,j+1 + βt,i,j gt+1,i,j + ξt,i,j gt+1,i,j−1 + λ−i gt+1,−i,j
∆t
where

αt,i,j =βt,i,j =ξt,i,j =−s∗t,i,j,B 1s∗t,i,j,B <0∆a,s∗t,i,j,B 1s∗t,i,j,B <0 − s∗t,i,j,F 1s∗t,i,j,F >0∆as∗t,i,j,F 1s∗t,i,j,F >0∆a− λi ,.

This equation can be stacked into a matrix form as Equation (12) such thatβ
ξt,1,10
 t,1,1 αt,1,2 βt,1,2 ξt,1,2
0


 0
αt,1,3 βt,1,3 ξt,1,3



..
..
.
.



 0
0
αt,1,JΓt = 

 λ2



λ2



λ2


..

.



T

λ1
λ1
0
..
.

λ1
..

βt,1,J00βt,2,128

λ1
ξt,2,10αt,2,2 βt,2,2

ξt,2,20αt,2,3 βt,2,3..
.

ξt,2,3..
.

..

αt,2,J

βt,2,J

0

λ2

.

0

0

0
.















 .














F

Further numerical results

Benchmark parameters
We use the parameters in Table 3 in the benchmark case, following those used in Ahn et al. (2018).
Later, we change the volatility and persistence parameters of the TFP to examine its effect on theaccuracy of solving the model between different algorithms. We set the range of grid points for K
as [0.8K̄, 1.2K̄] where K̄ is the value of aggregate capital in the stationary distribution. The rangeof grid points for Z is [−m̄σ, m̄σ] where m̄ is a real number. We set m̄ = 4. We divide the state
space of (Kl , Zm ) by 3 grid points for each.
Table 3: Parameter values
Parameters
γ : Relative Risk Aversion
ρ : Rate of Time Preference
α : Capital Shareδ : Rate of Capital Depreciation
τ : Tax Rate of Labor Incomeb : Rate of Compensation
1 − µ : Persistence of TFP
σ : Volatility of TFP
z : Idiosyncratic of Labor Productivity
λ : Probability of Labor Productivity

Benchmark Value
1.0
0.01
0.36
0.025
0.011
0.15
0.75
0.07
zu = 0, ze = 1
λe = 0.50, λu = 0.03

In Figure 4, we show the results of the simulation path in XPA. The red line {K̃t }t∈[0,T ] inthe figure shows the simulation obtained only from the forecasting rule for 10,000 periods, while
the blue line {Kt∗ }t∈[0,T ] shows the simulation obtained from the nonlinear model including thehousehold HJB equation (with the forecasting rule) and the Kolmogorov Forward equation for
10,000 periods. This is known as the Den haan’s (2010) fundamental plot showing the accuracy ofthe solution. It is clear from the figure that the capital paths resulting from these simulations arevery close.23
23
If the red and blue lines are close, households approximately act according to the correct forecasting rule. If not,households are using the wrong forecasting rule. Therefore, the divergence between the two lines indicates that themodel is not solved correctly based on rational expectations.

29

Figure 4: Simulation path in XPA (Den haan’s fundamental plot)
36.5

36

35.5

35
1000

2000

3000

4000

5000

6000

7000

8000

9000 10000

We also check the robustness of our results with respect to the persistence of TFP. That is,
the persistence of TFP, 1 − µ, is lowered from µ = 0.25 to µ = 0.5 or µ = 0.75. In Table 4, wecompare the mean of the simulation path derived from the nonlinear model for 10,000 periods in
XPA (compared with that in KS) and the mean of the simulation path derived from the linearized
model in REITER (compared with that in KS) for each µ. We see that the gap between XPA and
KS is much smaller than the gap between REITER and KS for alternative values of µ as well.
In Table 5, we show that the Den haan errors of KS, XPA, and REITER for each µ. We use
the sequences of {Kt∗ }t∈[0,T ] and {K̃t }t∈[0,T ] to measure the Den haan errors
AXεM≡ 100 · max | ln K̃t − ln Kt∗ |,DH
t∈[0,T ]∗t∈[0,T ] | ln K̃t − ln Kt |

P
EAN
εM
≡ 100 ·DHT.

It is clear that, regardless of the value of µ, the Den haan errors of XPA are smallest compared tothose of REITER and KS when σ is large.

30

Table 4: Simulation of capital paths
a. Case of µ = 0.25 (benchmark)Agg Shock σ (%)

0.01

0.1

0.7

1.0

3.0

5.0

XPA–KS (%)

0.001

0.013

0.097

0.160

0.529

1.079

REITER–KS (%)

0.002

0.022

0.062

0.084

2.832

6.688

Agg Shock σ (%)

0.01

0.1

0.7

1.0

3.0

5.0

XPA–KS (%)

0.001

0.007

0.055

0.098

0.341

0.711

REITER–KS (%)

0.002

0.015

0.056

0.051

1.524

4.127

Agg Shock σ (%)

0.01

0.1

0.7

1.0

3.0

5.0

XPA–KS (%)

0.000

0.005

0.039

0.074

0.270

0.571

REITER–KS (%)

0.001

0.011

0.047

0.053

0.922

2.908

a. Case of µ = 0.50

b. Case of µ = 0.75

31

Table 5: Den Haan errors: robustness
a. Case of µ = 0.25 (benchmark)Agg Shock σ (%)

0.01

0.1

0.7

1.0

3.0

5.0

AX
εM
DH XP A (%)
EANεMDH XP A (%)
AXεMDH KS (%)
EANεM
DH KS (%)
AXεM
DH REIT ER (%)
EANεMDH REIT ER (%)

0.001

0.011

0.084

0.125

0.336

0.580

0.000

0.004

0.028

0.038

0.091

0.140

0.002

0.015

0.097

0.134

0.493

0.885

0.001

0.011

0.075

0.104

0.286

0.438

0.000

0.003

0.144

0.387

4.379

9.759

0.000

0.002

0.106

0.273

3.381

7.546

Agg Shock σ (%)

0.01

0.1

0.7

1.0

3.0

5.0

AX
εM
DH XP A (%)
EANεMDH XP A (%)
AXεMDH KS (%)
EANεM
DH KS (%)
AXεM
DH REIT ER (%)
EANεMDH REIT ER (%)

0.001

0.008

0.048

0.062

0.157

0.308

0.000

0.003

0.019

0.024

0.053

0.083

0.001

0.009

0.059

0.083

0.278

0.483

0.001

0.007

0.044

0.060

0.155

0.235

0.000

0.001

0.069

0.155

2.431

5.961

0.000

0.001

0.057

0.126

1.989

4.864

Agg Shock σ (%)

0.01

0.1

0.7

1.0

3.0

5.0

AX
εM
DH XP A (%)

0.001

0.006

0.040

0.052

0.110

0.185

EAN
εM
DH XP A (%)
AXεMDH KS (%)
EANεM
DH KS (%)
AXεM
DH REIT ER (%)
EANεMDH REIT ER (%)

0.000

0.003

0.017

0.022

0.041

0.055

0.001

0.006

0.042

0.059

0.200

0.350

0.001

0.005

0.031

0.043

0.107

0.163

0.000

0.001

0.046

0.097

1.593

4.263

0.000

0.001

0.039

0.081

1.312

3.537

a. Case of µ = 0.50

b. Case of µ = 0.75

32