Applying the Explicit Aggregation Algorithm to Heterogeneous
Macro Models∗
Takeki Sunakawa†
July 2019

Abstract
This paper applies the explicit aggregation (Xpa) algorithm developed by den Haan andRendahl [2010. Solving the incomplete markets model with aggregate uncertainty using explicit
aggregation. Journal of Economic Dynamics and Control, 69–78] to the heterogeneous–firm
models of Khan and Thomas (2003, 2008) and the heterogeneous–household models of Krueger
et al. (2016). We find the Xpa algorithm is an order of magnitude faster for solving these models
than the standard Krusell–Smith (KS) algorithm because it does not need to simulate the distribution of individual capital and productivity when updating the aggregate forecasting rules.
However, the simulation results in the Xpa and KS algorithms in terms of both the micro– andmacro–level moments are almost identical, even in levels, because of bias correction. The Xpa
algorithm also exhibits accuracy comparable with the KS algorithm.

∗

This paper is originally a part of my dissertation at the Ohio State University. The author is especially gratefulto Aubhik Khan and Julia Thomas for their valuable guidance and support, and to Pontus Rendahl for sharing
his code for bias correction. The author also appreciates discussions with anonymous referees, Jie April Cai, Grey
Gordon, Youngsoo Jang, Heejeong Kim, Timothy Kam, Toshihiko Mukoyama, Taisuke Nakata, Nobuhide Okahata,
Tatsuro Senga, Shuhei Takahashi, Stephen Terry, Tomoaki Yamada, and Minchul Yum, as well as comments byseminar participants at Kyoto University, the Ohio State University, the Australian National University, and theShanghai University of Finance and Economics, and attendees at the 2013 Conference on Computing in Economics andFinance, the 10th Dynare Conference, the 2019 Spring Midwest Macroeconomics Meetings, and the 2019 Conference
on Computing in Economics and Finance. The usual disclaimer applies. This research is financially supported by aJapan Society for the Promotion of Science Grant–in–Aid for Scientific Research (KAKENHI) for Young Scientists(B) Project No. 15K20918 and the Seimeikai Foundation. The programming code used in the paper is publicallyavailable at https://github.com/tkksnk/Xpa.
†
Kobe University Center for Social Systems Innovation. E–mail: takeki.sunakawa@gmail.com.

1

Keywords: Projection methods, Heterogeneous agents, Explicit aggregation.
JEL codes: C63; D52

2

1

Introduction

Solving heterogeneous–agent macro models accurately and efficiently is more than ever an importantresearch agenda. Some recent studies revisit the importance of understanding the effect of micro–
level heterogeneity on aggregate macro dynamics (Ahn et al., 2018; Krueger et al., 2016). For
this reason, as more micro–level panel data consistent with macro aggregates become available,
structural estimations of heterogeneous–agent models using both micro– and macro–level data areof increasing interest. Accordingly, a faster and accurate computation method for heterogeneous–
agent models is necessary. den Haan and Rendahl (2010) develop the explicit aggregation (Xpa)
algorithm to solve the heterogeneous–household model studied in Krusell and Smith (1998). A
special issue in the Journal of Economic Dynamics and Control (den Haan et al., 2010) considers
different algorithms and shows that the Xpa algorithm can solve the model more efficiently than
the standard Krusell–Smith (KS) algorithm (Maliar et al., 2010).
In this paper, the Xpa algorithm originally proposed by den Haan and Rendahl (2010) is
applied to the heterogeneous–firm models studied in Khan and Thomas (2003, 2008) (hereafter KT)
and to the heterogeneous–household models examined in Krueger et al. (2016) (hereafter KMP).
The latter application is especially interesting as they demonstrate that micro–level heterogeneity
matters for the aggregate macro dynamics when the model is matched with the micro–level wealthdistribution, which is not the case in Krusell and Smith (1998). Indeed, the Xpa algorithm is anorder of magnitude faster than the KS algorithm for solving these models because it does not needto simulate the distribution of individual capital and productivity when updating the aggregate
forecasting rules. The Xpa algorithm is also easy to implement as we merely evaluate the individual
decision rules as if the amount of individual capital is equal to that of the aggregate capital, e.g.,k = Kt , especially in the case without idiosyncratic shocks.
We first show how to apply the Xpa algorithm to the KT (2003) model without idiosyncratic
productivity shocks (see also Thomas, 2002). Given Jensen’s inequality, bias correction using
steady–state information is the key to accuracy in terms of forecasting errors such as the den Haan
(2010) statistics. We find the simulation results using the Xpa and KS algorithms are almost

3

identical in terms of both the micro– and macro–level moments, even in levels, because of the bias
correction. Further, the Xpa algorithm is nearly equally as accurate as the KS algorithm. We alsopresent how the Xpa algorithm can apply to the KT (2008) model with idiosyncratic productivity
shocks, where the Xpa algorithm has den Haan errors at most only twice as large as those using
the Krusell–Smith algorithm. We accomplish this developing a novel approach we call –indexed
aggregation to deal with the idiosyncratic shocks.1
In terms of related work, Terry (2017) applies different algorithms to solve the KT (2008)
model and asserts that the Xpa algorithm is indeed faster, but less accurate, than other algorithms
employing projection methods. However, that implementation of the Xpa algorithm applies a naive
aggregation.2 Instead, in the present paper, we argue that the –indexed aggregation provides much
more accurate solutions by a factor of three or four more than what Terry (2017) reports for theden Haan errors.
We then show that the Xpa algorithm also works in a KMP (2016) environment in which heterogeneity matters. The –indexed aggregation is modified so that we can deal with exogenously
time–varying employment (Krusell and Smith, 1998). This is a non–trivial extension as the dimension of idiosyncratic shocks is made much larger to include not only employment status, but alsoearnings risk and discount factor shocks, in KMP’s benchmark model.3 Once again, the simulation
results using the Xpa and KS algorithms are very close and the accuracy in the forecasting errors
only slightly worsens with the Xpa algorithm. In our computation, we solve the benchmark modelin KMP in seconds with the Xpa algorithm using Fortran and OpenMP.
Krusell and Smith (1998) first develop a simulation–based method for solving an incomplete–market heterogeneous–household model with aggregate uncertainty. A computational difficulty1The –indexed aggregation evaluates the individual decision rules at k = Kt (i ), where Kt (i ) is the total amountof capital held by individuals that have a common idiosyncratic shocksof Rhouseholds
R i . Given that the measurePµt (i , k) can be decomposed into µkt (k|i )φ(i ), we can use Kt (i ) = kµkt (k|i )dk instead of Kt = i kµt (i , k)to evaluate the individual decision rules and take a weighted average of the evaluated functions with φ(i ) being theweight. As we know only the values of Kt as grid points, we approximate Kt (i ) ≈ ψ(i )Kt where ψ(i ) = Kss (i )/Kssby using steady–state information.
2
The Xpa algorithm is first adapted and applied to the heterogeneous–firm models by Sunakawa (2012), which
also uses naive aggregation.
3
den Haan and Rendahl (2010) include each –indexed capital in the indivisual decision rule as a state variable.
This approach is very costly or infeasible as the number of grid points for i increases because we need more dimensionsof the state space.

4

arises from the fact that the state of the economy is the whole distribution of individual capitaland productivity. To address this, they assume that approximate aggregation (Krusell and Smith,
1998; Young, 2005) holds, such that the entire state of the economy is approximated by the meanof the distribution. They then update the forecasting rules for individual households to infer thestate of the economy by simulating the model for a long time, which is time–consuming.4
In contrast, Reiter (2009) proposes a way to handle the whole distribution without assuming theaggregate approximation.5 The distribution is approximated by the non–stochastic method withhistogram (Young, 2010). Then, the equilibrium dynamics with aggregate shocks are approximatedat the first order by linearizing the model around the deterministic steady state. The method canthen handle hundreds of state variables because of such linearization.6 Although Reiter’s (2009)
method is also much faster than the standard KS algorithm, as is the Xpa algorithm, the method
is based on first–order perturbation and cannot consider the higher–order effects of aggregate
uncertainty as certainty equivalence holds.7
The reminder of the paper proceeds as follows. Section 2 presents the heterogeneous–firm models in KT (2003; 2008). Section 3 explains the KS and Xpa algorithms for each of the cases without
and with idiosyncratic shocks. Especially in the latter, we introduce –indexed aggregation todeal with the idiosyncratic shocks. We also show how to adjust the bias in levels from Jensen’s
inequality by following den Haan and Rendahl (2010). Section 4 analyzes the numerical results of
the heterogeneous–firm models with the calibrated parameters in KT (2003; 2008) and Section 5
presents the heterogeneous–household models in KMP (2016) and provides the calibration and numerical results. We find that in both cases of the heterogeneous–firm and heterogeneous–household
models, the Xpa algorithm is much faster and comparably accurate to the KS algorithm. Section
4

Alternatively, Algan et al. (2008, 2010) and Winberry (2018) show how to solve the model by parametrizing the
distribution. Gordon (2011) proposes a way to approximate the whole distribution with hundreds of state variables(i.e., bins in a histogram) with sparse grid.
5
Reiter’s (2010b) backward induction method can also be applied to models in which approximate aggregation does
not hold. The method is applied to stochastic overlapping generations models with aggregate uncertainty by Khan
(2017); Kim (2018). Okahata (2018) demonstrates the method can also merge with the continuous–time methods in
Ahn et al. (2018).
6
Reiter (2010a); Ahn et al. (2018) further develop a method to reduce the dimension of the state space by projecting
the distribution onto principal components. See also Bayer and Luetticke (2018). Childers (2018) suggests a distinct
approach using linearization.
7
Terry (2017) compares Reiter’s method and the standard KS method in solving variants of the KT (2008) model.

56 concludes. Additional details of the computation are in the Appendix.

2

Khan–Thomas models

In this section, we briefly explain the heterogeneous–firm models in KT (2003; 2008) and refer theinterested readers to the original articles for more detailed discussion. In brief, KT (2008) provides
two extensions to KT (2003), namely, idiosyncratic productivity shocks and constrained investment
without fixed costs, so as to match the micro–level moments of the investment–capital ratio in themodel and the data.
In this economy, there is a continuum of firms. Each firm produces output using the production
function: y = i zt F (k, n), where k is the capital stock, n is the labor input, i ∈ {1 , ..., N } is theidiosyncratic firm–specific productivity and zt ∈ {z1 , ..., zNz } is the aggregate productivity common
across firms. Each of i and zt follows a Markov chain. F (k, n) is a common production functionthat exhibits decreasing returns to scale. A random fixed cost ξ is incurred on adjusting individual
capital in each period.
A firm is defined by the individual state variables (i , k, ξ). The aggregate state of the economy
is given by (zt , µt ), where µt is the distribution of individual capital and productivity (i , k) heldby each firm, which is an infinitesimal object.8 Taking as given the shadow price of the household’sutility pt and the wage rate wt , the firm chooses labor input n and the next period’s capital k 0 .
The firm’s choice on labor is static. The firm chooses n so as to maximize the current profit

π(i , k; zt , µt ) = max {i zt F (k, n) − wt n} + (1 − δ)k,n(1)

where wt = Γw (zt , µt ) is the wage rate, which is a function of the aggregate state variables. Notethat the profit includes undepreciated capital (1 − δ)k with the depreciation rate of δ ∈ (0, 1] fromthe previous period.
The firm’s choice on capital is dynamic and discrete. That is, if the firm pays a fixed cost
(proportional to the wage rate) ξwt , the firm can choose any level of capital in the next period, k 0 .
8

We use µt and µt (i , k) interchangeably for the distribution, which is also the measure of firms with (i , k).

6

If not, the firm can choose k 0 within a limited range of Ω(k).9 The firm’s value by choosing thenext period’s capital k 0 isṽ0 (i , k, ξ; zt , µt )= π(i , k; zt , µt ) + max {−ξwt




X
X
p
t+1
0 0
0
0
ṽ(
,
k
;
z
,
µ
)
,
π
(z
|z
)
+ max
−γk
+
β
π
(
|
)
t+1 t+1
z t+1 ti
k0 >0 pt0
 |i
zt+1 |zt




X
X
pt+1 0 0πz (zt+1 |zt )max−γk 0 + βπ (0 |i )ṽ( , k ; zt+1 , µt+1 ),pt
k0 ∈Ω(k) 
0
 |i

(2)

zt+1 |zt

where pt = Γp (zt , µt ) is the (shadow) price and µt+1 = Γµ (zt , µt ) is the forecasted distribution in the next period. π (0 |i ) (πz (zt+1 |zt )) is the conditional probability of the next period’s
idiosyncratic (aggregate) productivity given the current period’s productivity. ṽ(i , k; zt , µt ) =R ξ̄
0 ṽ0 (i , k, ξ; zt , µt )dξ is the ex ante value of the firms that commonly have (i , k) before they drawa random fixed cost ξ. Each firm discounts the future by a common stochastic discount factor
βpt+1 /pt , whereas β ∈ (0, 1) is the deterministic discount factor. γ is the deterministic trend of theeconomy.
We define v(i , k; zt , µt ) ≡ pt ṽ(i , k; zt , µt ). By multiplying pt on both sides of (2) and integrating
the equation over ξ, we have a transformed Bellman equation that is recursive in v(i , k; zt , µt ):v(i , k; zt , µt ) = pt π(i , k; zt , µt )Z ξ̂(i ,k;zt ,µt )− pt
ξG0 (ξ)dξwt
0+ α(i , k; zt , µt )E0 (i ; zt , µt ) + (1 − α(i , k; zt , µt )) E1 (i , k; zt , µt ),9

If Ω(k) = (1 − δ)k is a point set, firms hold the current level of capital after depreciation.

7(3)

where




X
X
0 0
0
0
πz (zt+1 |zt )v( , k ; zt+1 , µt+1 ) , (4)−γpt k + βπ ( |i )E0 (i ; zt , µt ) = max

k0 >0 0 |izt+1 |zt




X
X
πz (zt+1 |zt )v(0 , k 0 ; zt+1 , µt+1 ) . (5)E1 (i , k; zt , µt ) = max−γpt k 0 + βπ (0 |i )k0 ∈Ω(k) 0
 |izt+1 |zt

If the firm chooses to pay the fixed cost, it selects k ∗ (i ; zt , µt ) regardless of current capital k andthe value of adjustment is E0 (i ; zt , µt ). If the plant chooses not to pay the fixed cost, the nextˆ i , k; zt , µt ) =period’s capital is k c (i , k; zt , µt ) and the value of non–adjustment is E1 (i , k; zt , µt ). ξ((E0 (i ; zt , µt ) − E1 (i , k; zt , µt )) / (wt pt ) is the threshold level of the fixed cost at or below which
the firm chooses to pay the fixed cost as E0 (i ; zt , µt ) ≥ E1 (i , k; zt , µt ) holds. G is the cumulative¯ → [0, 1] and α(i , k; zt , µt ) = G ξ(ˆ i , k; zt , µt ) is thedensity function of the distribution G : [0, ξ]
probability of adjusting capital.
Given the decision rules for adjusters and non–adjusters, the ex post firm–level decision rule forcapital isk 0 = K(i , k, ξ; zt , µt ) =k ∗ ( ; z , µ ),
i

t

t

¯ ξ(
ˆ i , k; zt , µt )},if ξ ≤ min{ξ,



¯ ξ(
ˆ i , k; zt , µt )}.
k c (i , k; zt , µt ), if ξ > min{ξ,
By convexifying K(i , k, ξ; zt , µt ) over ξ among the firms that commonly have k, we obtain the exante firm–level decision rule for capital:
gk (i , k; zt , µt ) = α(i , k; zt , µt )k ∗ (i ; zt , µt ) + (1 − α(i , k; zt , µt ))k c (i , k; zt , µt ).

(6)

In addition, the firm–level decision rule for output is given bygy (i , k; z, µ) = zt F (i , k, n∗ (i , k; zt , µt )),(7)

where n∗ (i , k; zt , µt ) is the optimal level of labor input that the firms choose by solving the intratemporal problem (1).

8The individual decision rules are aggregated as

Yt =

XZ

gy (i , k; zt , µt )µt (i , k)dk,(8)

gk (i , k; zt , µt )µt (i , k)dk,(9)i

Kt+1 =

XZ
iwhere µt (i , k) is the measure of firms with (i , k). We normalize the distribution so thatP R
iµt (i , k)dk =1 holds. Having Yt and Kt+1 at hand, the resource constraint implies

Ct = Yt + (1 − δ)Kt − Kt+1

where Ct is the aggregate consumption and Kt =P R
ikµt (i , k)dk is the aggregate capital. The

labor market clears such that

Nt =

XZ

"
n∗ (i , k, zt , µt ) +Z ξ̂(i ,k;z,µ)#ξG0 (ξ)dξ µt (i , k)dk

0iholds, where Nt is the aggregate hours worked.
Finally, the prices and wages are functions of the aggregate state of the economy because thehousehold’s consumption–saving and working decisions have effects on the prices and wages in a
general equilibrium. Let U (Ct , 1 − Nt ) = log Ct + η(1 − Nt ) be the household’s utility function
period by period. Then the prices and wages are given by pt = D1 U (Ct , 1 − Nt ) = 1/Ct and
wt = D2 U (Ct , 1 − Nt )/D1 U (Ct , 1 − Nt ) = η/pt as a result of the household’s optimization.
A recursive competitive equilibrium is defined so as to satisfy (i) the household’s optimality,(ii) each firm’s optimality, (iii) feasibility (i.e., markets clearing), and (iv) consistency between theindividual decision rules and the aggregate forecasting rules.
Regarding the feasibility stated in the definition above, the individual decision rules gy (i , k; zt , µt )and gk (i , k; zt , µt ) are implicit functions of the prices pt and wages wt = η/pt . For a conjectured
price p̃t , (4) and (5) are solved for the individual decision rules. Then (8) and (9) yield the aggregate
variables Yt , Kt+1 , and the updated price pt = Ct−1 . There is a mapping denoted by pt = F (p̃t )9from the conjectured price to the updated price. To obtain the market–clearing price, a fixed point
problem p∗t = F (p∗t ) is solved.10
Regarding consistency, each of the Xpa or KS algorithms achieves this in a different manner,as shown in the following section.

3

Two algorithms

In this section, we explain the two algorithms used in this study, namely, the Xpa algorithm
developed by den Haan and Rendahl (2010) and the standard KS algorithm. The individual
decision rules are consistent with the forecasting rules for the aggregate state of the economy. The
Xpa and KS algorithms achieve consistency in a different manner. The only difference in thesealgorithms is their way of aggregating the individual decision rules to obtain the forecasting rules.
For the most part, the Xpa algorithm is easier to implement as it does not require simulations, but
it is not straight–forward to apply the method to models with idiosyncratic shocks. For this reason,
we present a novel approach we call –indexed aggregation to deal with idiosyncratic shocks.

3.1

Without idiosyncratic shocks

We first consider the case without idiosyncratic shocks as in KT (2003). We remove idiosyncratic
shocks i and constrain investment by setting Ω(k) = (1−δ)k (i.e., no adjustment is allowed without
paying the fixed cost) in the model presented in the previous section. We also assume that theapproximate aggregation of Krusell and Smith (1998) and Young (2005) holds. That is, the meanRof the distribution at the beginning of period t, Kt = kµt (k)dk, is sufficient for individuals toforecast the state of the economy. Therefore, Kt enters into the individual decision rules insteadof µt . Taking the individual decision rules with the aggregate approximation y = gy (k; zt , Kt ) andk 0 = gk (k; zt , Kt ) as given, we obtain the aggregate output and the next period’s aggregate capital10
Ct and Nt are unpredetermined variables. Therefore, the market–clearing prices are solved to ensure the feasibility. This is also the case in Chang and Kim (2007, 2014); Takahashi (2014); Jang et al. (2018). In contrast, theprices and wages are functions of only predetermined variables in Krusell and Smith (1998).

10

by
Z
Yt =

gy (k; zt , Kt )µt (k)dk,(10)

gk (k; zt , Kt )µt (k)dk.

(11)

Z
Kt+1 =We also obtain Ct = Yt + (1 − δ)Kt − Kt+1 and pt = Ct−1 . As discussed earlier, we impose market
clearing.
In the KS algorithm, we simulate the model for a long time to obtain the sequence of {zt , Kt , pt }Tt=1
taking the sequence of exogenous aggregate productivity shock {zt } and the initial value of K1 asgiven. Using the sequence obtained, we construct a subset of the sequence sorted by the realizationof zt for each zi ∈ {z1 , ..., znz } and estimate the forecasting rule by ordinary least squares in theform oflog Kt+1 = bK,0 (zi ) + bK,1 (zi ) log Kt ,log pt = bp,0 (zi ) + bp,1 (zi ) log Kt ,KS
for each i = 1, ..., nz . We denote them by Kt+1 = ΓKS
K (Kt , zt ) and pt = Γp (Kt , zt ).

In the Xpa algorithm, instead of simulating the model, we just move the integral over k into
the individual decision ruleZ
Kt+1 =

gk (k; zt , Kt )µt (k)dk,
Z≈ gk ( kµt (k)dk; zt , Kt ),= gk (Kt ; zt , Kt ).

Note that when we solve for the individual decision rule, Kt is given as a grid point instead of theintegral of the distribution. We immediately obtain the forecasting rule for aggregate capital in
the Xpa algorithm Kt+1 = ΓXpa (Kt , zt ) ≡ g(Kt ; zt , Kt ) by evaluating the individual decision rulegk (k; zt , Kt ) at k = Kt . We also obtain Yt = gy (Kt ; zt , Kt ), which implies the forecasting rule for11the shadow price of the household’s utility by Ct = Yt + (1 − δ)Kt − Kt+1 and pt = Ct−1 . We again
consider market clearing.

Bias correction Because of Jensen’s inequality, the Xpa algorithm yields biases in Yt and Kt+1 .
The individual decision rules gy (k; zt , Kt ) and gk (k; zt , Kt ) must be linear at k = Kt for the Xpa
algorithm to have the same aggregate variables as in (10) and (11) obtained in the KS algorithm.
den Haan and Rendahl (2010) suggested a way to correct these biases in Yt and Kt+1 . Let Kss =Rḡk (k)µ̄(k)dk be the level of aggregate capital, which is obtained by the steady–state individualR
decision rule ḡk (k) and the stationary distribution µ̄(k). Also let K̃ss = ḡk ( k µ̄(k)dk) = ḡk (Kss )be the decision rule ḡk (k) evaluated at k = Kss . Then a time–invariant number for bias correction
is given by ζk ≡ Kss − K̃ss = Kss − ḡk (Kss ) , which measures Jensen’s inequality in the steadystate. We do the same for output and obtain ζy ≡ Yss − Ỹss = Yss − ḡy (Kss ). We then add thebias–correction terms to the forecasting rules in the Xpa algorithm, that is,

Kt+1 = gk (Kt ; zt , Kt ) + ζk ,Yt = gy (Kt ; zt , Kt ) + ζy .

Figure 1 illustrates the shape of the individual decision rule gk (k; zt , Kt ) evaluated at k = Kt .
The figure is a typical example of (S, s)–type capital adjustment. That is, only when the current
level of individual capital is sufficiently far from the target level will the firm choose to adjust capital;
otherwise, the firm will continue holding the current level of capital (after capital depreciation).
Even after convexifying the ex post firm–level decision rule over random fixed costs, the ex antefirm–level decision rule exhibits considerable nonlinearity. This is quite different from the individual
decision rules in the household consumption–saving problem in Krusell and Smith (1998), which
are almost linear save points near the origin where poor households face borrowing constraints.
Even though the individual decision rule gk (k; zt , Kt ) is nonlinear in k, when it is evaluated at
k = Kt , the resulting function gk (Kt ; zt , Kt ) is linear in Kt .11 The red line in the figure depicts
11

This result may suggest that the Xpa algorithm works when the forecasting rule obtained by the KS algorithm
exhibits small nonlinearity, e.g., a high R–squared.

12this relationship. Given Jensen’s inequality, there is a downward bias in gk (Kt ; zt , Kt ) comparedwith ΓK (zt , Kt ). Note that gk (k; zt , Kt ) is convex in k taking (zt , Kt ) as fixed. We can correct
for this bias using the steady–state counterparts, ḡk (k). After bias correction, we can see that theintersection of ΓK (zt , Kt ) and the 45–degree line (in black) corresponds to the deterministic steadystate of Kss (the vertical gray line).12
12

Terry (2017) uses a different method for correcting biases. The intercepts in the forecasting rules, log Kt+1 =

b

(z )

K,0 i
bK,0 (zi ) + bK,1 (zi ) log Kt , are adjusted so that the steady state implied by the forecasting rules exp 1−b
K,1 (zi )coincides with the deterministic steady state Kss for each level of zi . In the proposed method, we undertake such a
correction using the individual decision rules in the steady state.

13

Figure 1: Individual decision rules and forecasting rules.
g k (k;z t ,K t )
2
1.8

k'

1.6
1.4
1.2
1

1.4
1.2
1

Kt

1.5

1

0.5

2.5

2

3

k

Notes: The x–axis is the current period’s individual capital k, the y–axis is the aggregate capital Kt , andthe z–axis is the next period’s individual capital k 0 . zt = 1.0.
1.5

1.4

Kt+1

1.3

1.2

1.1

g (K ;z ,K )

1

k

t

t

t

(z ,K )
K t t

1

1.1

1.2

1.3

1.4

1.5

Kt

Notes: The x–axis is the current period’s aggregate capital Kt and the y–axis is the next period’s aggregate
capital Kt+1 .

14

3.2

With idiosyncratic shocks

When we allow for idiosyncratic shocks, aggregate capital in the next period (we omit aggregate
output, but the same discussion applies) is given by

Kt+1 =XZ

gk (i , k; zt , Kt )µt (i , k)dk,i

The measure of (i , k) can be decomposed intoµt (i , k) = µkt (k|i )φ(i ),where φ(i ) =Rµt (i , k)dk is the marginal distribution of i assumed to be a time–invariant measure(i ,k)of i and µkt (k|i ) = R µµtt(is the conditional distribution of µt (i , k).13 Then, the aggregation,k)dki

in the Xpa algorithm becomes

Kt+1 =XZ

gk (i , k; zt , Kt )µkt (k|i )φ(i )dk,

i

≈

X

Z
gk (i ,kµkt (k|i )dk; zt , Kt )φ(i ),

i

=

X

gk (i , Kt (i ); zt , Kt )φ(i ),iwhere Kt (i ) ≡Rkµkt (k|i )dk is the amount of capital indexed by i at the beginning of period t.
P RPNote that Kt = i kµt (k, i )dk = i Kt (i )φ(i ) holds.
–indexed aggregation

We wish to use Kt (i ) to evaluate the individual decision rule. However,

we know the value of Kt only as a grid point. Therefore, we apply so–called –indexed aggregation as
follows. Kt (i ) is not on grid points and we evaluate the individual decision rule at Kt (i ) = ψ(i )Kt
whereψ(i ) =13

Kss (i )Kss (i )=P,
Kss
i Kss (i )φ(i )Note that its integral over k is normalized to one, i.e.,

R

15

R

µ ( ,k)dkµkt (k|i )dk = R µtt (ii ,k)dk = 1.

Kss (i ) =Rk µ̄k (k|i )dk, and Kss =P Rik µ̄(i , k)dk =P

i Kss (i )φ(i ) (µ̄(i , k) is the stationary

distribution of (k, i )).14
Bias correction

When we allow for idiosyncratic shocks, we perform bias correction for eachRlevel of idiosyncratic shocks. Let Kss (i ) = ḡk (i , k)µ̄k (k|i )dk be the amount of i –indexed
capital in the steady state, which is obtained by the decision rule ḡk (i , k) and the stationaryR
distribution µ̄(i , k). Also let K̃ss (i ) = ḡk (i , k µ̄k (k|i )dk) = ḡk (i , Kss (i )) as the decision rule
evaluated at Kss (i ). Then the bias–correction term depends on each i , ζk (i ) ≡ Kss (i )− K̃ss (i ) =Kss (i ) − ḡk (i , Kss (i )), and we haveKt+1 =

X[gk (i , Kt (i ); zt , Kt ) + ζk (i )] φ(i ).

i

Discussion

Terry (2017) and Sunakawa (2012) undertake a naive explicit aggregation. Thatis, they evaluate the policy function gk (i , k; zt , Kt ) at k = Kt instead of k = Kt (i ). However,this aggregation may lead to large forecasting errors, as shown in Terry (2017). den Haan andRendahl (2010) assume that the individual decision rule is a function of each –indexed capital.
For example, in the standard Krusell and Smith (1998) model, the idiosyncratic shock takes twovalues, i ∈ {u, e}, and the policy function becomes g(i , k; zt , Kt (u), Kt (e)). We evaluate this at
k = Kt (i ) for each i as we know the value of –indexed capital as a grid point. This approachis very costly or infeasible as the number of grid points for i increases because we need moredimensions of the state space.
We use steady–state information for both the bias correction and the –indexed aggregation.
Therefore, an accurate steady–state solution is necessary. We find that non–stochastic methods(Young, 2010) augmented with the use of eigenvalue and eigenvector decomposition often yield aprecise approximation to the stationary distribution. See Appendix A.1 for details. In contrast,
stochastic methods may give a poor approximation to the stationary equilibrium that also worsens
14

The assumption that ψ(i ) = Kt (i )/Kt is time–invariant may not be valid with large aggregate shocks when theshare of each –indexed capital significantly changes over time. Later we relax this assumption with an exogenously
time–varying employment measure in Section 5.

16

the precision of the Xpa algorithm.15

4

Numerical examples

In this section, we present the numerical results from applying the KS and Xpa algorithms to
the heterogeneous–firm models in KT (2003; 2008). We examine both cases without idiosyncratic
shocks (KT 2003, aka the traditional model) and with idiosyncratic shocks (KT 2008, aka theextended model), although we mainly focus on the latter in the main text.

4.1

Parameterization

The parameter values used in this paper closely follow those in KT (2003; 2008). The functional
forms are U (C, 1 − N ) = log C − ηN and F (k, n) = k α nν , where η > 0, α ∈ (0, 1), ν ∈ (0, 1), and¯ → [0, 1].
α + ν < 1. The distribution of the random fixed costs drawn by firms is uniform G : [0, ξ]
The Markov chain for aggregate– and firm–level productivity approximates an AR(1) process by
Tauchen’s (1986) method; z 0 = ρz z + ηz0 and ε0 = ρε ε + ηε0 for each, where ηz0 ∼ N (0, σzη ) and
i
h
1−δ+b
k,
k
, whereηε0 ∼ N (0, σεη ). The range of constrained investment is given by Ω(k) = 1−δ+aγγ−a = b > 0.
The parameter values and targeted moments are in Tables 1 and 2, respectively. The unit oftime is a year. The mean growth rate of technology γ implies a 1.6 percent average annual growthrate. The discount factor β is from an average annual real interest rate βγ −1 − 1 of 4 percent. The
annual depreciation rate is 6.9 percent to match the average investment–capital ratio in the fixed
asset tables. Labor’s share ν is 0.64 as in Prescott (1986).
We internally calibrate the other parameters in the model. The capital share α is set to match acapital–output ratio of 2.353. η is a normalization parameter and set so that the steady state labor¯ b, σξη ) are chosen to match the data of the average firm–level investment
is 1/3. The values of (ξ,rates in Cooper and Haltiwanger (2006) as shown in Table 2. Our steady–state results show that
we match these internal calibration targets quite well. In the traditional model of KT (2003), firm–
15
Rios-Rull (1999) mentions that non–stochastic methods provide a better approximation of the steady state thanstochastic methods.

17

level productivity is constant and there is no constrained investment, i.e., b = σξη = 0. ξ¯ = 0.014
is set to match the average investment rate of positive spikes (i/k > 0.20). In the extended modelof KT (2008), the other moments are reasonably matched between the model and the data. Notethat we set ρε = ρz as there is little agreement about the persistence of the firm–level idiosyncratic
productivity.
Table 1: Parameters.
a. Externally calibrated
β

γ

δ

ν

ρz

σzη

ρε

Traditional

0.977

1.016

0.069

0.640

0.859

0.014

—

Extended

0.977

1.016

0.069

0.640

0.859

0.014

0.859

α

η

ξ¯

b

σεη

Traditional

0.256

2.400

0.014

0.0

—

Extended

0.256

2.400

0.0083

0.011

0.022

b. Internally calibrated

Notes: The values of (ξ, b, σεη ) are chosen by following KT (2008).

Table 2: Targeted moments.
a. Aggregate variables (without aggregate uncertainty)
K/YN

Traditional

2.34870.3337

Extended

2.35150.3338b. Firm–level investment rates (without aggregate uncertainty)
Mean

Std.dev.

Inaction

Spike+

Spike-

Invest+Invest-Traditional

0.10460.02140.21130.12840.00000.78870.0000Extended

0.11580.07160.00970.15950.01300.81270.1779Census data

0.122

—

0.081

0.186

0.018

0.815

0.104

Notes: Census data are from Cooper and Haltiwanger (2006). Inaction: |i/k| < 0.01, Spike+:
i/k ≥ 0.2, Spike-: i/k ≤ −0.2, Invest+: i/k ≥ 0.01, and Invest-: i/k ≤ −0.01.
18

4.2

Model statistics

To start, we detail the unconditional business cycle statistics in Table 3. We simulate the modelfor 2,500 periods and discard the first 500 periods. We take logs of the simulated data and use
the Hodrick–Prescott filter to remove the low frequency components. We can then see that after
being detrended, the outcomes of the two algorithms are very similar. This is also true for thehigher–order moments, including skewness and excess kurtosis, of the aggregate investment–capitalratio, although estimates of these higher–order moments are usually unstable.
As we can see from the plot of the simulated time series of the aggregate variables in Figure
2, not only the second– and higher–order moments, but also the levels are very close to eachother for the two algorithms. This is through the bias–correction procedure explained in Section
3. Compared with previous studies that apply the Xpa algorithm (Terry, 2017; Sunakawa, 2012),
the –indexed aggregation with bias correction yields almost indistinguishable results in both levels
and higher–order moments from those with the KS algorithm.
We also consider the conditional moments by computing the impulse responses to a positive
shock to the aggregate productivity. We simulate the economy with no shocks for a long period
so that the economy reaches the stochastic steady state (Coeurdacier et al., 2011). Then we hit
the economy with a one–time persistent positive productivity shock. As depicted in Figure 3, we
can see that the outcomes of the two algorithms are again quite similar, although the responses ofaggregate capital as computed by the Xpa algorithm are slightly weaker.

19

Table 3: Business cycle statistics, extended model.

Standard deviationKS
Xpa
Output correlation
KS
Xpa

Y

C

I

N

K

Z

(2.5329)(2.5100)

0.43440.44013.61093.57670.66520.65640.51380.50970.59720.60271.00001.00000.87320.8802

0.9731
0.97340.93640.93600.04590.04490.99960.9995

Persist.

Std.dev.

Skewness

Exc.kur.

0.64990.65250.01070.01050.07140.0575-0.0943-0.0915

Agg. investment rate, I/K
KS
Xpa

Notes: We take logs of the simulated data and use the Hodrick–Prescott filter (with a filteringparameter λ = 100) to remove the low–frequency components.

Figure 2: Unconditional business cycle simulation, extended model.

Output

Log

−0.2

Xpa
KS

−0.4

0.24

Log

0.22
0.20
0.18
0.16

−1.050Labor−1.075−2.2−1.100

−2.4

−0.6
−0.8

Investment−1.125−1.150−2.6
1 10 20 30 40 50
1 10 20 30 40 50
1 10 20 30 40 50

Capital−0.84

Consumption

TFP
0.02

−0.86

0.00

−0.88

−0.02

−0.90

−0.04

1 10 20 30 40 50
Year

1 10 20 30 40 50
Year

20

1 10 20 30 40 50
YearFigure 3: Impulse responses, extended model.

OutputXpa
KS

Percent2

Percent

6

0.5
0.0

0
5 10 15 20

0CapitalConsumption

1.5

1.00

1.0

0.75

5 10 15 20

5 10 15 20
Year

0.00

5 10 15 20

1.0
0.5

0.25
0

0

TFP

0.50

0.5
0.0

1.0

2
0

Labor

1.5

4

1
0

Investment

8

0

5 10 15 20
Year

0.0

0

5 10 15 20
Year

Table 4: Level differences between KS and Xpa, extended model.
a. Aggregate variables
K/YN

p = C −1KS

2.35060.33382.3666Xpa

2.35380.33392.3656SS

2.35150.33382.3663Mean

Std.dev.

Inaction

Spike+

Spike-

Invest+Invest-KS0.11580.07170.00980.15950.01300.81260.1779Xpa

0.11580.07170.00970.15910.01300.81260.1779SS

0.11580.07160.00970.15950.01300.81270.1779b. Firm–level investment rates

Notes: Inaction: |i/k| < 0.01, Spike+: i/k ≥ 0.2, Spike-: i/k ≤ −0.2, Invest+: i/k ≥ 0.01, andInvest-: i/k ≤ −0.01.
Table 4 displays the aggregate and disaggregate variables in levels at the stochastic steady state.
Interestingly, the levels of the capital–output ratio and the hours worked are almost the same as21

in the deterministic steady state. This is also true when we look at the distribution of firm–level
investment rates. This means that the higher–order effects of aggregate uncertainty on the levels of
aggregate and disaggregate variables are very limited, as opposed to the heterogeneous–household
models of Krusell and Smith (1998) in which households are risk averse.
In the traditional KT model (2003), we also observe a strong similarity in the unconditionaland conditional moments between the outcomes of the Xpa and KS algorithms. The corresponding
tables and figures are in Appendix A.3.

4.3

Accuracy and efficiency

To evaluate accuracy in terms of forecasting errors, we employ not only the R–squared and the root
mean squared errors (RMSE), but also the den Haan (DH) statistics based on the forecasting rules
ΓK (Kt , zt ) and Γp (Kt , zt ) as in den Haan (2010). We calculate both the static and dynamic forecast
errors as in Terry (2017). To calculate the dynamic forecast errors, taking as given the sequenceof exogenous shock {zt }Tt=1 and the initial value of K1 , we use the forecasting rules recursively togenerate the sequence of the next period’s aggregate capital and the shadow price, {K̃t+1 , p̃t }Tt=1 .
The sequence is compared with the simulated data from the model, {Kt+1 , pt }. Then we haveuK,t = K̃t+1 − Kt+1 ,up,t = |p̃t − pt | ,for t = 1, ..., T . The static forecast errors are one–step ahead forecast errors, i.e., eK,t = |ΓK (Kt , zt ) − Kt+1 |and ep,t = |Γp (Kt , zt ) − pt |. We sort the forecast errors by the realization of zt for each zi .16 The
dynamic forecast errors are more robust measures of accuracy as errors accumulate over time, andwe report the max and mean DH errors based on the dynamic forecast errors.
Figure 4 is the so–called den Haan fundamental plot for aggregate capital and consumption(the inverse of the shadow price). We compare the simulated data with the dynamic and staticqP
Ti
2
Using the static forecast errors sorted for each zi with Ti samples, we also have RMSEs as
t=1 eK,t
qP
P
P
P
P
Ti
Ti
Ti
2
2
2
i
i
and
and 1 − Tt=1
e2p,t / Tt=1
(pt − p̄)2 wheret=1 ep,t and R–squares as 1 −t=1 eK,t /t=1 (Kt+1 − K̄)PTiPTiK̄ = (1/Ti ) t=1 Kt+1 and p̄ = (1/Ti ) t=1 pt .
16

22

forecasts. We can see that the differences between the simulated data and the dynamic or static
forecasts (i.e., uKt and upt or eKt and ept ) are very small throughout. Specifically, by using theextended model, we report in Table 5 that the maximum values of the DH statistics in the Xpa
algorithm are at most twice those in the KS algorithm (0.792 vs. 0.337 for uKt and 0.263 vs. 0.119
for upt ). Indeed, they are much smaller (by a factor of three or four times) than what Terry (2017)
reports for the max DH statistics (3.457 for uKt and 1.210 for upt ). In the traditional model, thedifference is even smaller as the Xpa algorithm is almost equally as accurate as the KS algorithm.
Figure 4: den Haan fundamental plots, extended model.

KS:Capital

0.4

−0.80
Dynamic
StaticActual

0.0
1

10

20

30

40

Log

Log

0.2

−0.2

−0.85
−0.90

50

Xpa:Capital

0.4

−0.95

1

10

20

30

40

50

Xpa:Consumption

−0.75
−0.80
Log

Log

0.2
0.0
−0.2

KS:Consumption

−0.75

−0.85
−0.90

1

10

20

30

40

50

23

−0.95

1

10

20

30

40

50

Table 5: DH statistics, RMSE, and R–squared.
a. Extended model

KS

Xpa

z1
z2
z3
z4
z5
z1
z2
z3
z4
z5

DHmax, %K0p0.3095 0.1174
0.3386 0.1189
0.2605 0.1122
0.3064 0.1062
0.2297 0.0743
0.7921 0.1858
0.6629 0.2535
0.5470 0.2631
0.3941 0.2341
0.2227 0.1892

DHmean, %K0p0.0055 0.0021
0.0166 0.0063
0.0178 0.0063
0.0149 0.0052
0.0040 0.0016
0.0404 0.0055
0.0831 0.0197
0.0558 0.0228
0.0165 0.0145
0.0032 0.0057

RMSE, %
K0p
0.0256 0.0064
0.0274 0.0060
0.0270 0.0056
0.0256 0.0064
0.0184 0.0056
0.1660 0.1574
0.0820 0.0701
0.0379 0.0243
0.0196 0.0312
0.0243 0.0600

R–squaredK0p0.9999 1.0000
0.9999 1.0000
0.9999 1.0000
0.9999 1.0000
0.9999 1.0000
0.9966 0.9873
0.9993 0.9978
0.9999 0.9997
1.0000 0.9995
0.9999 0.9973

DHmax, %K0
p0.3255 0.1490
0.3675 0.1528
0.2945 0.1407
0.3389 0.1271
0.2371 0.0866
0.4755 0.1374
0.4457 0.1624
0.4125 0.1769
0.3425 0.1779
0.2747 0.1706

DHmean, %K0p0.0055 0.0023
0.0181 0.0075
0.0197 0.0080
0.0159 0.0062
0.0042 0.0019
0.0256 0.0046
0.0628 0.0158
0.0559 0.0214
0.0277 0.0172
0.0072 0.0063

RMSE, %
K0p
0.0293 0.0189
0.0314 0.0184
0.0307 0.0170
0.0276 0.0160
0.0195 0.0134
0.0941 0.0708
0.0620 0.0364
0.0366 0.0171
0.0216 0.0268
0.0149 0.0405

R–squaredK0p0.9999 0.9998
0.9999 0.9999
0.9999 0.9999
0.9999 0.9999
0.9999 0.9999
0.9989 0.9975
0.9996 0.9994
0.9999 0.9999
0.9999 0.9997
1.0000 0.9988

b. Traditional model

KS

Xpa

z1
z2
z3
z4
z5
z1
z2
z3
z4
z5

i
i
Notes: Taking the dynamic and static forecast errors (sorted for each zi ) {uy,t }Tt=1
and {ey,t }Tt=1forTiiy ∈ {K, p}qas given, DHmax is the maximum value of {uy,t }t=1 , DHmean is the average of {uy,t }Tt=1,
PTi 2
PTi 2 PTi
P
T
2i
RMSE ist=1 ey,t , and R–squared is 1 −t=1 ey,t /
t=1 (yt − ȳ) where ȳ = (1/T )t=1 yt .

Finally, in Table 6, we compare the two algorithms in terms of computation time. The programming code is written in Fortran with OpenMP directives so that all the cores on a computer
are utilized. We find the Xpa algorithm is about 20–100 times faster than the KS algorithm. This
result is comparable with that obtained by den Haan and Rendahl (2010) for the Krusell and Smith(1998) model. Specifically, the speed gain is mainly in the outer loop to obtain the forecasting rules,24

as we do not need to simulate the model to update the forecasting rules in the Xpa algorithm. Notethat the two algorithms share the same subroutines for calculating the steady state, inner loops,
and unconditional/conditional simulations.
Table 6: Computational time.

ExtendedKS
Xpa
Xpa/KS
Traditional
KS
Xpa
Xpa/KS

Inner

Outer

# of loop

Total

0.750
0.745
0.993

122.607
0.032<0.00110.00011.0001.100

1233.570
8.547
0.007

0.386
0.386
1.000

6.532
0.007
0.001

6.000
6.000
1.000

41.5052.355
0.057

Notes: Computation using a workstation with Xeon E5–2696v4 (2.2 Ghz, 44 cores).

5

Krueger–Mitman–Perri models

Krueger, Mitman and Perri (2016; hereafter KMP) extend the analysis in Krusell and Smith (1998)
and show that micro–level heterogeneity matters for the aggregate dynamics.17 In this section, weshow that the Xpa algorithm can be applied to the models in KMP.

5.1

Model

In the economy, there are the representative firm, heterogeneous households, and the government.
The continuum of households differs in their asset holdings, denoted as k, and the other idiosyncratic
shocks, represented by a vector xi ≡ (ei , yi , βi ), which include the following.
• Employment status ei : ei ∈ {1, ρ}, where ρ is the replacement ratio. That is, households haveunemployment insurance.
17

We use the model in a working paper version of Krueger et al. (2015). There are three major differences from
Krueger et al. (2016): (i) earnings risk is an AR(1) process instead of a joint process combining persistent andtransitory components, (ii) there is no retirement and social insurance, which influences (iii) the calibration of β.
Nonetheless, introducing retirement and social insurance and/or a more complicated stochastic process for earnings
risk is relatively straight–forward.

25

• Labor productivity yi : Households face idiosyncratic earnings risk,
log y 0 = φ log y + η 0 ,where η 0 ∼ N (0, ση2 ). We transform the continuous variables in the AR(1) process into thediscretized grid points of a Markov chain. That is, yi ∈ {y1 , ..., yNy }.
• Discount factor βi : Households are born with βi ∈ [β̄ − , β̄ + ] and live with β their entire
life. They die with probability (1 − θ) ∈ [0, 1]. A fraction θ of households is born in every
period with no assets.
A household is defined by the individual state variables (xi , k). The aggregate state of the economy
is given by (zt , µt ), where zt is aggregate productivity and µt is the distribution of individual assets
and the other characteristics held by each household. Taking as given the wage rate wt and thereal interest rate rt , the household (xi , k) chooses consumption c and the next period’s capital k 0 .
The household’s life–time utility from choosing (c, k 0 ) is

max {u(c)XXX+ θβi
πz (zt+1 |zt )πe (e0 |ei , zt+1 , zt )πy (y 0 |yi )v(xi , k; zt , µt ) =c,k0 ≥0

zt+1 e00 0y0×v(e , y , βi , k 0 ; zt+1 , µt+1 ) ,(12)subject to the budget constraint and forecasting rule for the aggregate state in the next periodc + k 0 = (1 − τ (zt ; ρ))wt yi ei + (1 + rt − δ)k/θ,µt+1 = Γ(zt , µt ),where each of πz (zt+1 |zt ), πe (e0 |ei , zt+1 , zt ), and πy (y 0 |yi ) is the conditional probability of the next
period’s aggregate or idiosyncratic state given the current period’s state. τ (zt ; ρ) is a flat income tax
with which the government finances unemployment insurance. By solving the household’s problem,26

we have the individual decision rule for saving
k 0 = g(xi , k; zt , µt ),(13)

which is aggregated as

Kt+1 =XZ

g(xi , k; zt , µt )θµt (xi , k)dk,

iwhere µt (xi , k) is the measure of households with the amount of individual asset k and otherpeculiarities xi and θ is the survival probability from period t to period t + 1. Note that the
replacing newborns have no assets.
The representative firm has access to a Cobb–Douglas production technology Yt = zt Ktα L1−αtand maximizes profit Yt − (r + δ)Kt − wt Lt in each period. The first–order necessary conditions are
wt = (1 − α)zt Ktα L−αt ,rt = αzt Ktα−1 L1−α− δ,twhere Kt =P R
ikµt (xi , k)dk is the aggregate capital and Lt =P R
i

I(ei = 1)µt (k, xi )dk is theamount of aggregate hours worked.
A recursive competitive equilibrium is defined so as to satisfy (i) each household’s optimality,(ii) the firm’s optimality, (iii) feasibility: Ct = Yt + (1 − δ)Kt − Kt+1 , and (iv) consistency between
the individual decision rules and the aggregate forecasting rules.

–indexed aggregation with time–variant employment measure

We assume that approx-

imate aggregation holds. That is, v(xi , k; zt , µt ) and g(xi , k; zt , µt ) are replaced by v(xi , k; zt , Kt )and g(xi , k; zt , Kt ) and the aggregate forecasting rule is given by Kt+1 = Γ(zt , Kt ). The aggregate
capital in the next period is given by

Kt+1 =XZ

g(xi , k; zt , Kt )θµt (xi , k)dk,i

27

The measure of (xi , k) can be decomposed intoµt (xi , k) = µkt (k|xi )φ(xi ; zt ),where φ(xi ; zt ) =Rµt (xi , k)dk is the marginal distribution of xi that is assumed to be a time–

variant measure of xi depending on zt .18 This implies exogenously time–varying hours worked

Lt =

XZ

I(ei = 1)µt (k, xi )dk

i

=

X

I(ei = 1)φ(xi ; zt )i

= 1 − u(zt )With the time–varying employment measure φ(xi ; zt ), the aggregation in the Xpa algorithm
becomes

Kt+1 =XZ

g(xi , k; zt , Kt )θµkt (k|xi )φ(xi ; zt )dk,

i

≈

X

Z
g(xi ,

kµkt (k|xi )dk; zt , Kt )θφ(xi ; zt ),

i

=

X

g(xi , Kt (xi ; zt ); zt , Kt )θφ(xi ; zt ),iwhere Kt (xi ; zt ) ≡Rkµkt (k|xi )dk is the amount of capital indexed by xi at the beginning of periodt. As before, Kt (xi ; zt ) is not on the grid points and we evaluate the individual decision rule atR
Kt (xi ; zt ) = ψ(xi ; zt )Kt where ψ(xi ; zt ) = Kss (xi ; zt )/Kss and Kss (xi ; zt ) = k µ̄k (k|xi ; zt )dk.
We also calculate the bias–correction terms

ζ(xi ) = Kss (xi ) − ḡ(Kss (xi )),where Kss (xi ) =Rḡk (i , k)µ̄k (k|xi )dk. Note that the bias–correction terms do not depend on zt .19(xi ,k)µkt (k|xi ) = R µµtt(xis the conditional distribution of µt (xi , k), which is also time–variant.
i ,k)dk19
There are two different decompositions of the stationary distribution here. One is µ̄(xi , k) = µ̄k (k|xi )φ(xi ) and18

28

Then we haveKt+1 =

X[g(xi , Kt (xi ; zt ); zt , Kt )θ + ζ(xi )] φ(xi ; zt ).

i

5.2

Parameterization

We calibrate the model to quarterly US data, where the capital share is α = 0.36 and the depreciation rate is δ = 0.025. A two–state Markov process is used for aggregate productivity withtransition matrixπz (z 0 |z) = ρl
1 − ρh1 − ρl .
ρhWe have two calibration targets to determine the persistence parameters (pl , ph ): the average
length of a severe recession ELl = 1/(1 − ρl ) and the fraction of time in a severe recession Πl =(1 − ρh )/(2 − ρl − ρh ). With ELl = 22 quarters and Πl = 16.48% obtained from the sample periods,
we have (ρl , ρh ) = (0.9545, 0.9910).20
To identify the levels of aggregate technology at each state, we also target Yl /Yh = 0.9298,
corresponding to a fall in GDP per capita during the great recession of about 7%. We also havethe average unemployment rates u(zl ) = 8.39% and u(zh ) = 5.33% from the sample periods.
Then we have zl /zh = (Yl /Ll )/(Yh /Lh ) = 0.9608 assuming Ll /Lh = Kl /Kh . A normalization
Πl zl + (1 − Πl )zh = 1 provides (zl , zh ) = (0.9676, 1.0064). For the employment status, Krueger
et al. (2016) use the job finding and separation rates from CPS and obtain transition matrices
the other is µ̄(xi , k) = µ̄k (k|xi ; zt )φ(xi ; zt ).
20
Krueger et al. (2016) define a recession to be one where the unemployment rate rises above 9% for one quarterand continues to be above 7%. During 1948.I–2014.III, two recession periods, 1980.II–1986.II and 2009.I–2013.III are
duly identified.

29

depending on the current and next period’s aggregate technology z, z 0 ∈ {zl , zh }: 0.3378πe (e0 |e, zl , zl ) = 
0.0606 0.3382
πe (e0 |e, zh , zl ) = 
0.0696

0.6622  0.18900 , πe (e |e, zh , zh ) = 
0.93940.0457
0.6618  0.22200 , πe (e |e, zl , zh ) = 
0.93040.0378
0.8110 ,
0.9543

0.7780 .
0.9622

Note that these matrices are consistent with the time–varying unemployment rates u(zl ) andu(zh ).21 For the labor productivity, Krueger et al. (2016) estimate (φ̂, σ̂y2 ) from the annual PSID
data and translate them for the quarterly model, (φ, σy2 ) = (0.9457, 0.0359). Rouwenhorst’s (1995)
method instead of Tauchen (1986) is used for discretization because of the long persistence φ. The
set of discount factors B = [β̄ − , β̄ + ] is divided into grid points where (β̄, ) = (0.98349, 0.01004).
θ = 0.99375 so that the probability of death is 1/160. For the replacement ratio, ρ = 0.5 is usedt )ρholds from the government’s balanced
for the benchmark case. Note that τ (zt ; ρ) = 1−u(zu(zt )+u(zt )ρ

budget. The parameter descriptions and values are summarized in Table 7.
21
The stationary distribution of employment status implied by πe (e0 |e, zl , zl ) (πe (e0 |e, zh , zh )) coincides with theemployment rate 1 − u(zl ) (1 − u(zh )).

30

Table 7: Parameters.
Parameters
Capital share
Depreciation rate
Risk aversion
Agg. productivity
Unemployment rateαδσ(zl , zh )(ul , uh )Transition for agg. prod.

πz (z 0 |z)

Replacement ratio
Transition for emp. statusρπe (e0 |e, zh , zh )πe (e0 |e, zh , zl )πe (e0 |e, zl , zh )πe (e0 |e, zl , zl )Persistence of idio. prod.
Std. dev. of idio. prod.
Discount factor
Discount factor shock
Survival rate

φ
σy2
β̄

θ

0.36
0.025
1.0
(0.9676, 1.0064)
 (0.0839, 0.0533) 
0.9545 0.0455
0.0090 0.9910
0.50


0.1890 0.8110
 0.0457 0.9543 
0.3382 0.6618
 0.0696 0.9304 
0.2220 0.7780
 0.0378 0.9622 
0.3378 0.6622
0.0606 0.9394
0.94570.03590.98349
0.01004
0.99375

Source: Krueger, Mitman, and Perri (2015, working paper).
In Table 8, we replicate the stationary distribution results in KMP. We consider three calibrations. The baseline calibration is denoted by KMP. In the KS calibration, we omit the earnings
risk and discount factor heterogeneity. We also set θ = 1 and ρ = 0.01. That is, we have nostochastic death and unemployment insurance (we need a positive ρ > 0 so that the amount ofconsumption is positive). We add the earnings risk only to the KS calibration, denoted +σ(y). The
KMP calibration shows that there are many households that have no net wealth. In contrast, thetop 20% of households hold most wealth (82.7% in the 2006 PSID vs. 82.9% in the model).22 Themodel with the baseline calibration also has a Gini coefficient of 0.78, which is comparable withthe data.
The KS calibration has low wealth inequality and a small Gini coefficient (0.32). In our com22

Unfortunately, the model cannot generate the concentration in the top 1% as observed in the data.

31

putation, the +σ(y) calibration lies in between the KMP and KS calibrations in terms of wealthinequality. The three calibrations have a similar capital–output ratio at the aggregate level, although the KS calibration has a slightly low value.
Table 8: Net wealth distribution (without aggregate uncertainty).

% of Share:
Q1
Q2
Q3
Q4
Q5
90–95
95–99
T1
Wealth Gini
K/Y

Data
PSID, 06 SCF, 07
-0.9
-0.2
0.8
1.2
4.4
4.6
13.0
11.9
82.7
82.5
13.7
11.1
22.8
25.3
30.9
33.5
0.77
0.78
—
—

KMP
0.11820.8796
3.5852
12.5377
82.879319.2195
29.3285
16.01370.784311.0821

ModelsKS
7.369412.338117.1474
23.7850
39.3601
9.91659.99753.47000.320510.2674+σ(y)1.76196.432213.6520
24.9265
53.2275
13.4301
14.58315.44050.5162
11.5596

Notes: The figures in columns “Data” (2006 PSID and 2007 SCF) are from Krueger et al. (2016).
KS: Remove the earnings risk and discount factor heterogeneity and set θ = 1 and µ = 0.01. +σ(y):Remove the discount factor heterogeneity and set θ = 1 and ρ = 0.01.

5.3

Model statistics

To start, we compute the unconditional business cycle statistics for each calibration in Table 9. We
can see that the two algorithms yield very similar results in all three specifications. The outcomesare also similar in the second–order moments as well as in the levels. Specifically, we have ahigh correlation between output and consumption, volatile investment, and persistent output. The
baseline KMP calibration has the largest procyclicality of consumption, less volatile investment,
and less persistent output. Figure 5 also depicts that the unconditional simulation paths for eachalgorithm are close to each other in the KMP calibration. This is also the case for the other two
calibrations. See Appendix B.3.

32

Table 9: Business cycle statistics.
Model
KMP
KS
Xpa
KS
KS
Xpa
+σ(y)
KS
Xpa

Mean(K/Y )Corr(Y, C)

Std.dev.(I)

Corr(Y, Y−4 )11.0627
11.06240.94050.93910.04900.04940.88510.885610.2887
10.2873

0.92080.9226

0.05170.05120.89470.894011.4663
11.46440.91390.91680.0567
0.05670.89110.8913

Notes: All variables unfiltered and in levels as in Krusell and Smith (1998).

Figure 5: Unconditional business cycle simulation, KMP calibration.

CapitalConsumption

1.07

3.79

1.06

3.78

1.05

Log

Log

3.80

3.77

1.04
1.03

3.76

Xpa
KS

1.02
1

40

80

120 160 200

1

40

1.38

0.10

1.36

0.05

1.34

0.00

1.32

−0.05
1

40

80

120 160 200

Investment

Log

Log

Output80

120 160 200

1

40

80

120 160 200

Table 10 provides the micro–level moments at the stochastic steady state as well as the aggregate
capital–output ratio and the fall in aggregate consumption from a negative aggregate productivity
shock. Once again, the results from the two algorithms are almost indistinguishable. In bothalgorithms, the capital–output ratio is slightly higher than in the deterministic steady state becauseof precautionary saving in the presence of aggregate uncertainty (see also Table 8). This cannot be
captured by Reiter’s (2009) linearization method. The consumption decline in response to a one–33

time negative shock is also similar. We confirm this in Figure 6 by plotting the impulse responses
for the KMP calibration in which heterogeneity matters.
Interestingly, in the KS calibration, the micro–level moments with aggregate uncertainty areslightly different from those in the stationary distribution. This is also because of precautionarysaving. With aggregate uncertainty, the poor tend to save more, which makes the rich dissave as
the interest rate falls (K/Y becomes higher compared with the deterministic steady state). As
a result, the wealth Gini is even smaller in the stochastic steady state than in the deterministic
steady state.
Table 10: Net wealth distribution (with aggregate uncertainty) and consumption decline.

% of Share:
Q1
Q2
Q3
Q4
Q5
90–95
95–99
T1
Wealth Gini
K/Y
∆C

KMP
KS
Xpa
0.11590.11580.88390.88333.63593.634012.6631 12.6591
82.7012 82.7077
19.2587 19.2582
29.2337 29.2371
15.7834 15.7900
0.78270.782811.0678 11.0672-2.2794 -2.3179

ModelsKS
KS
Xpa
8.35038.301713.4778 13.4081
18.0413 17.9760
23.8262 23.8071
36.3043 36.5071
9.14909.19858.92018.99642.97873.01640.28110.283510.3031 10.2993-1.7865 -1.8358

34

+σ(y)
KS
Xpa
1.74691.74676.48246.482113.7234 13.7233
24.9657 24.9658
53.0816 53.0820
13.3950 13.3951
14.5243 14.5244
5.40925.40930.51490.514911.4788 11.4781-1.7581 -1.8120

Figure 6: Impulse responses, KMP calibration.

TFP

−0.5

−1
−2
−3
−4

−1.5

0

1

2

3

4

5

Xpa
KS
0

Output

−2
−4

1

2

3

4

5

4

5

Investment

0
Percent

Percent

5.4

−1.0
−2.0

0

−6

Consumption

0.0
Percent

Percent

0

−5
−10
−15

0

1

2

3

4

5

0

1

2

3

Accuracy and efficiency

Table 11 reports the forecasting errors for the dynamic and static forecasts as in Section 4.3. We
can see that the forecasting rules in the Xpa algorithm are very accurate and comparable with the
KS algorithm. For example, the max DH error is 0.345 in the Xpa algorithm and 0.219 in the KS
algorithm for the KMP calibration. We see that the mean DH error tends to be higher in goodtimes (i.e., zt = zh ) in both the Xpa and KS algorithms.
Finally, we display the computation time in Table 12. The Xpa algorithm is a factor of 100–200
times faster than the KS algorithm. While the number of outer loops tends to be higher in the
Xpa algorithm, it takes less than 1/100th of a second per loop to calculate.
In our computation, we solve the KMP model within seconds with the Xpa algorithm using
Fortran and OpenMP. Note that most of the computation time is in the inner loop.23
23
We can solve the inner loop even faster by using the endogenous grid point method (Carroll, 2006; Barillas andFernández-Villaverde, 2007; Fella, 2014), for example.

35

Table 11: DH statistics, RMSE, and R–squared.

KMP
KS
Xpa
KS
KS
Xpa
+σ(y)
KS
Xpa

DHmax, %DHmean, %RMSE, %

R–squared

zh
zl
zh
zl

0.20690.21930.34040.34450.04610.0062
0.09880.01680.00330.00430.00500.01181.00001.00001.00001.0000

zh
zl
zh
zl

0.04650.05020.16590.16110.01280.00210.10510.01010.00100.00170.00470.00591.00001.00001.00001.0000

zh
zl
zh
zl

0.03850.03710.29940.26250.00670.0016
0.11420.00790.00070.00160.00730.01261.00001.00001.00001.0000ii
Notes: Taking the dynamic and static forecast errors (sorted for each zi ) {uK,t }Tt=1
and {eK,t }Tt=1
Ti
Ti
as q
given, DHmax is the maximum value of {uK,t }t=1 , DHmean is the average of {uK,t }t=1 , RMSE
PTi
PTi 2 PTi
PTi 2
2
is
t=1 Kt+1 .
t=1 (Kt+1 − K̄) where K̄ = (1/T )t=1 eK,t /t=1 eK,t , and R–squared is 1 −Table 12: Computational time.

KMP
KS
Xpa
Xpa/KSKSKS
Xpa
Xpa/KS+σ(y)KS
Xpa
Xpa/KS

Inner

Outer

# of loop

Total

0.455
0.360
0.792

68.2280.003<0.00114.00014.0001.000

961.568
5.078
0.005

0.203
0.301
1.484

30.9390.000<0.00110.00015.0001.500

311.421
4.520
0.015

0.243
0.365
1.503

37.0820.001<0.001

12.00013.0001.083

447.897
4.754
0.011

Notes: Computation using a workstation with Xeon E5–2696v4 (2.2 Ghz, 44 cores).

36

6Concluding remarks

In this paper, the Xpa algorithm originally proposed by den Haan and Rendahl (2010) is applied tothe heterogeneous–firm models in Khan and Thomas (2003, 2008) and the heterogeneous–household
models in Krueger et al. (2016). We demonstrate that the Xpa algorithm is much faster than the
KS algorithm and has a similar accuracy in terms of errors in the forecasting rules. Future research
may include applying the Xpa algorithm to other classes of heterogeneous–agent models, such as
HANK models (e.g., Bayer et al., 2018; Gornemann et al., 2016; Kaplan et al., 2018). As we
can solve the models using a projection method that preserves nonlinearity at the aggregate level,
structural estimations by matching the micro– and macro–level moments in these models by using
the Xpa algorithm may also be of great interest.24

24See, e.g., Mongey and Williams (2017); Winberry (2018); Williams (2017) for recent attempts at structural
estimations using Reiter’s (2009) perturbation methods. Note that none of these studies can investigate the higher–
order effects of aggregate uncertainty as certainty equivalence holds.

37

References
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.
Algan, Y., O. Allais, and W. J. den Haan (2008): “Solving heterogeneous-agent models
with parameterized cross-sectional distributions,” Journal of Economic Dynamics and Control,
32, 875–908.
——— (2010): “Solving the incomplete markets model with aggregate uncertainty using parameterized cross-sectional distributions,” Journal of Economic Dynamics and Control, 34, 59–68.
Barillas, F. and J. Fernández-Villaverde (2007): “A generalization of the endogenous grid
method,” Journal of Economic Dynamics and Control, 31, 2698–2712.
Bayer, C. and R. Luetticke (2018): “Solving heterogeneous agent models in discrete time withmany idiosyncratic states by perturbation methods,” Working Paper.
Bayer, C., R. Luetticke, L. Pham-Dao, and V. Tjaden (2018): “Precautionary Savings,
Illiquid Assets, and the Aggregate Consequences of Shocks to Household Income Risk,” Econometrica, forthcoming.
Carroll, C. (2006): “The method of endogenous gridpoints for solving dynamic stochastic optimization problems,” Economics Letters, 91, 312–320.
Chang, Y. and S.-B. Kim (2007): “Heterogeneity and Aggregation: Implications for LaborMarket Fluctuations,” American Economic Review, 97, 1939–1956.
——— (2014): “Heterogeneity and Aggregation: Implications for Labor-Market Fluctuations: Reply,” American Economic Review, 104, 1461–1466.
Childers, D. (2018): “Solution of Rational Expectations Models with Function Valued States,”
Working paper.

38

Coeurdacier, N., H. Rey, and P. Winant (2011): “The Risky Steady State,” American
Economic Review, 101, 398–401.
Cooper, R. W. and J. C. Haltiwanger (2006): “On the Nature of Capital Adjustment Costs,”
Review of Economic Studies, 73, 611–633.
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.
Fella, G. (2014): “A generalized endogenous grid method for non-smooth and non-concave problems,” Review of Economic Dynamics, 17, 329–344.
Gordon, G. (2011): “Computing Dynamic Heterogeneous-Agent Economies: Tracking the Distribution,” PIER Working Paper 11-018.
Gornemann, N. M., K. Kuester, and M. Nakajima (2016): “Doves for the Rich, Hawks for
the Poor? Distributional Consequences of Monetary Policy,” Board of Governors of the Federal
Reserve System International Finance Discussion Papers 1167.
Jang, Y., T. Sunakawa, and M. Yum (2018): “Is Household Heterogeneity Important for
Business Cycles?” Working Paper.
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.
39

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.
Krueger, D., K. Mitman, and F. Perri (2015): “Macroeconomics and Household Heterogeneity,” Working paper.
——— (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.
Maliar, L., S. Maliar, and F. Valli (2010): “Solving the Incomplete Markets Model with
Aggregate Uncertainty Using the Krusell-Smith Algorithm,” Journal of Economic Dynamics
and Control, 34, 42–49.
Mongey, S. and J. Williams (2017): “Firm dispersion and business cycles: Estimating aggregate
shocks using panel data,” Mimeo.
Okahata, N. (2018): “An Alternative Solution Method for Continuous-Time Heterogeneous Agent
Models with Aggregate Shocks,” Mimeo.
Prescott, E. C. (1986): “Theory Ahead of Business Cycle Measurement,” Federal Reserve Bankof Minneapolis Quarterly Review, 10, 9–22.
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.
40

——— (2010b): “Solving the incomplete markets model with aggregate uncertainty by backward
induction,” Journal of Economic Dynamics and Control, 34, 28–35.
Rios-Rull, V. (1999): “Computation of Equilibria in Heterogenous Agent Models,” in Computational Methods for the Study of Dynamic Economies, Oxford University Press.
Rouwenhorst, K. G. (1995): “Asset Pricing Implications of Equilibrium Business Cycle Models,”
in Frontiers of Business Cycle Research, ed. by T. Cooley, New Jersey: Princeton University
Press, 294–330.
Sunakawa, T. (2012):“Applying the Explicit Aggregation Algorithm to Discrete Choice

Economies: With an Application to Estimating the Aggregate Technology Shock Process,”
Mimeo.
Takahashi, S. (2014): “Heterogeneity and Aggregation: Implications for Labor-Market Fluctuations: Comment,” American Economic Review, 104, 1446–1460.
Tauchen, G. (1986): “Finite state markov-chain approximations to univariate and vector autoregressions,” Economics Letters, 20, 177–181.
Terry, S. J. (2017): “Alternative Methods for Solving Heterogeneous Firm Models,” Journal of
Money, Credit and Banking, 49, 1081–1111.
Thomas, J. K. (2002): “Is Lumpy Investment Relevant for the Business Cycle?” Journal of
Political Economy, 110, 508–534.
Williams, J. (2017): “Bayesian estimation of DSGE models with heterogeneous agents,” Mimeo.
Winberry, T. (2018): “A Method for Solving and Estimating Heterogeneous Agent Macro Models,” Quantitative Economics, forthcoming.
Young, E. R. (2005): “Approximate Aggregation,” Mimeo, University of Virginia.

41

——— (2010): “Solving the incomplete markets model with aggregate uncertainty using the
Krusell-Smith algorithm and non-stochastic simulations,” Journal of Economic Dynamics andControl, 34, 36–41.

42

Appendix (not for publication)
A

Details of heterogeneous–firm models

A.1

Stationary equilibrium

We use steady–state information for both the bias correction and the –indexed aggregation in the
Xpa algorithm. As a precise solution of the stationary equilibrium is necessary, we explain how tocompute the stationary equilibrium in the following steps.
1. Make an initial guess of w and p = η/w because of the household’s utility with indivisible
labor.
2. Taking p and w as given, solve for the value function and the individual decision rules(a) Make an initial guess of the value function v (0) (i , k).25(b) Taking v (l−1) (i , k) as given (l is an index for iteration in the inner loop), approximatePthe expected value function 0 |i π (0 |i )v (l) (0 , k 0 ) and solve




X
0
0
(l) 0 0E0 (i ) = max
−γpk
+
β
π
(
|
)v
(
,
k
)
,

i

k0 >0 
0
 |i




X
−γpk 0 + βπ (0 |i )v (l) (0 , k 0 ) ,E1 (i , k) = max
k0 ∈Ω(k) 
0

(14)

(15)

 |i

for k ∗ (i ) and k c (i , k). The individual decision rules are given byḡk (i , k) = α(i , k)k ∗ (i ) + (1 − α(i , k))k c (i , k),ḡy (i , k) = i F (k, n∗ (i , k)),where n∗ (i , k) = arg maxn {i F (k, n) − wn}.
25

k is a continuous variable and we discretize it on a computer by grid points. We suppose that k takes a valueof grid points kj ∈ {k1 , ..., knk } and use unidimensional piecewise cubic splines to interpolate the values between thegrid points. We omit the index for k for the sake of exposition.

43

(c) Update the value function v (l) (i , k) byv (l) (i , k) = p max {i F (k, n) − wn} + p(1 − δ)knZ ξ̂(i ,k)− η
ξG0 (ξ)dξ0+ α(i , k)E0 (i ) + (1 − α(i , k)) E1 (i , k),

(16)

where
ˆ i , k) = (E0 (i ) − E1 (i , k)) /η,ξ(ˆ i , k) .
α(i , k) = G ξ((d) Iterate Steps (b) and (c) until v (l) − v (l−1) is sufficiently small.
3. Solve for the stationary distribution µ̄(i , k) implied by the individual decision rules (seebelow).
4. Update the prices: Given the stationary distribution, the individual decision rules are aggregated as

Y =

XZ

ḡy (i , k)µ̄(i , k)dk,(17)ḡk (i , k)µ̄(i , k)dk,

(18)

i

K=

XZ
i

C = Y − δK and p = C −1 .
5. Iterate Steps 2–4 until convergence.
In Step 3, we approximate the dynamics of µ̄(i , k) using the non–stochastic method with histogram(Young, 2010). That is, µ̄(i , k) is a continuous object with regard to k and we suppose that ktakes a value of grid points kj ∈ {k1 , ..., kñk }.26 Given the current period’s distribution µ(i , kj )26

We have finer grid points for approximating the distribution than interpolating the value function in Step 2, i.e.,ñk > nk .

44

for i = 1, ..., ne and j = 1, ..., ñk , the next period’s distribution µ(i0 , kj 0 ) for i0 = 1, ..., ne andj 0 = 1, ..., ñk is obtained from the transition probabilities from the current period’s index (i, j) tothe next period’s index (i0 , j 0 ). There are two cases of transition depending on the value of a fixedˆ i , kj ) withcost drawn by the firms indexed by (i, j). (i) When a firm pays the fixed cost ξ ≤ ξ(probability α(i , kj ): for each i, there is an index j00 such that k 0 = k ∗ (i ) ∈ [kj00 , kj00 +1 ].27 Then thetransition probabilities from (i, j) to (i0 , j00 ) and (i0 , j00 + 1) areP(i0 , j00 |i, j) =π (i0 |i )ω ∗ (i )α(i , kj ),P(i0 , j00 + 1|i, j) =π (i0 |i ) (1 − ω ∗ (i )) α(i , kj ),where π (i0 |i ) is the transition probability from i to i0 and ω ∗ (i ) =kj 0 +1 −k∗ (i )0

kj 0 +1 −kj 00. (ii) When a0

ˆ i , kj ) with probability 1 − α(i , k): for each (i, j), there is
firm does not pay the fixed cost ξ > ξ(an index j10 such that k 0 = k c (i , kj ) ∈ [kj10 , kj10 +1 ]. Then the transition probabilities from (i, j) to(i0 , j10 ) and (i0 , j10 + 1) areP(i0 , j10 |i, j) =π (i0 |i )ω c (i , kj ) (1 − α(i , kj )) ,P(i0 , j10 + 1|i, j) =π (i0 |i ) (1 − ω c (i , kj )) (1 − α(i , kj )) ,where ω c (i , kj ) =kj 0 +1 −kc (i ,kj )1

kj 0 +1 −kj 0
1. Given the transition probabilities for all (i, j), we have a transition

1

matrix P such that µ~ = P~µ holds, where µ~ is a stacked (ne ñk × 1) vector of the stationary
distribution. µ~ is obtained by an iterative method with an initial distribution or eigenvalue andeigenvector decomposition of the matrix P.28
As shown in Figure 7, the distribution takes the form of a discrete function, especially inthe traditional model of KT (2003) without idiosyncratic productivity shocks and constrained
investment. Therefore, we need fine grid points to approximate such a function.
27

Note that all the firms indexed by i choose the same level of capital regardless of kj .
P is a sparse matrix. To compute the eigenvalue and eigenvector decomposition, Reiter (2009) suggests using
MATLAB’s command eigs. We implement its functionality in Fortran by using the libraries arpack and Sparse BLAS.
28

45

Figure 7: Stationary distributions.
a. Traditional model
0.14
0.12

Density

0.10
0.08
0.06
0.04
0.02
0.00

1

2

Capital

3

4

5

b. Extended model
ε1

ε2

ε3

Density

0.020
0.002

0.010

0.001
0.000

0.04

0.015

0.02

0.005
2

4

0.000

2

ε4

Density

0.06

ε5

0.006

0.04

0.004

0.02

0.002

0.00

4

2
4
Capital0.000

24Capital

46

0.00

2

4

A.2

Xpa and KS algorithms

In both the Xpa and KS algorithms, we solve for the individual decision rule g and the aggregate
forecasting rule Γ by the following nested loops. In the inner loop, given the forecasting ruleΓ(n−1) (n is an index for iteration), we solve for the individual decision rule g (n) . In the outerloop, given the individual decision rule obtained in the inner loop, we solve for the forecasting ruleΓ(n) . The inner and outer loops are iterated until both g and Γ converge. At the convergence,
the individual decision rule is consistent with the forecasting rule. Each KS and Xpa algorithm
achieves consistency in a different manner. The inner loop is common to both algorithms and only
the outer loop differs.

Inner loop(n−1)

Taking as given the forecasting rules Kt+1 = ΓK

(n−1)

(zt , Kt ) and pt = Γp(zt , Kt ),the common inner loop for the two algorithms is as follows.
1. Make an initial guess of the value function v (0) (i , k; zt , Kt ).29
(n−1)

2. Taking v (l−1) (i , k; zt , Kt ) and Kt+1 = ΓK(zt , Kt ) as given (l is an index for iteration in

the inner loop),(a) Approximate the expected value function in the next periodh(l−1)(k 0 , Kt ) '
i zt

X
0

π (0 |i )

X

(n−1)

πz (zt+1 |zt )v (l−1) (0 , k 0 ; zt+1 , ΓK(zt , Kt ))zt+1

for each i and zt .
(b) Solveno0(l−1) 0E0 (i ; zt , Kt ) = max
−γp
k
+
βh
(k
,
K
)
,
t
t
i zt
k0 >0
n
o
0
E1 (i , k; zt , Kt ) = max −γpt k 0 + βh(l−1)
(k
,
K
)
,
t
i zt
k0 ∈Ω(k)

29

(19)
(20)

(k, Kt ) are continuous variables and we discretize them on a computer by a rectangle of grid points. We supposethat k takes a value of kj ∈ {k1 , ..., knk } and Kt takes a value of Km ∈ {K1 , ..., KnK } and use bidimensional piecewise
cubic splines to interpolate the values between the grid points.

47

for k ∗ (i ; zt , Kt ) and k c (i , k; zt , Kt ). The individual decision rules are given by
gk (i , k; zt , Kt ) = α(i , k; zt , µt )k ∗ (i ; zt , Kt )+ (1 − α(i , k; zt , Kt ))k c (i , k; zt , Kt ),gy (i , k; zt , Kt ) = i zt F (k, n∗ (i , k, zt , Kt )),where n∗ (i , k, zt , Kt ) = arg maxn {i zt F (k, n) − wt n}.
3. Update the value function v (l) (i , k; zt , Kt ) byv (l) (i , k; zt , Kt ) = pt max {i zF (k, n) − wt n} + pt (1 − δ)knZ ξ̂(i ,k;zt ,Kt )− η
ξG0 (ξ)dξ + α(i , k; zt , Kt )E0 (i ; zt , Kt )0

+ (1 − α(i , k; zt , Kt )) E1 (i , k; zt , Kt ),whereˆ i , k; zt , Kt ) = (E0 (i ; zt , Kt ) − E1 (i , k; zt , Kt )) / (wt pt ) ,ξ(ˆ i , k; zt , Kt ) .
α(i , k; zt , Kt ) = G ξ(4. Iterate Steps 2 and 3 until v (l) − v (l−1) is sufficiently small.
Outer loop in the KS algorithm In the outer loop in the KS algorithm, the individual decisionrules aggregate with µt (i , kj ) in each period. The model is simulated for long enough to obtainthe sequence of aggregate variables.
Taking the value function v(i , k; zt , Kt ) obtained in the inner loop as given, the outer loop in
the KS algorithm is as follows.
1. Make a sequence of {zt } for t = 1, ..., T and an initial distribution µ1 (i , kj ) approximated by
P P
a histogram at each grid point for i = 1, ..., ne and j = 1, ..., ñk so that i j µ1 (i , kj ) = 1
holds.
48

2. In each period t = 1, ..., T , zt and the current period’s distribution µt (i , kj ) are given.
(a) The aggregate capital is given by Kt =(n−1)

Kt+1 = ΓKP Pij kj µt (i , kj ).

Taking v(i , kj ; zt , Kt ) and(zt , Kt ) as given, make a conjecture of p̃t (which is not necessarily theforecasted price) and solve for gy (i , kj ; zt , Kt ) and gk (i , kj ; zt , Kt ) by the same procedure
as in Step 2 for the inner loop. Then aggregate the individual decision rules with the
distribution to obtain
XX

Yt =

i

XX

Kt+1 =gy (i , kj ; zt , Kt )µt (i , kj ),ji

gk (i , kj ; zt , Kt )µt (i , kj ),jCt = Yt − γKt+1 + (1 − δ)Kt , and pt = Ct−1 . That is, there is a mapping from theconjectured price to the new price and it is denoted by pt = F (p̃t ).
(b) Iterate Step 2(b) to find the market–clearing price p∗t satisfying p∗t = F (p∗t ).
(c) Update the next period’s distribution µt+1 (i , kj ) by µ
~ t+1 = P~µt where µ~ t is a stacked
vector of the distribution in period t and P is a transition matrix with the transition
probabilities from the current period’s index (i, j) to the next period’s index (i0 , j 0 ) andˆ i , k; zt , Kt ),(i0 , j 0 + 1). There are two cases. When ξ ≤ ξ(P(i0 , j00 |i, j) =π (i0 |i )ω ∗ (i ; zt , Kt )α(i , kj ; zt , Kt ),P(i0 , j00 + 1|i, j) =π (i0 |i ) (1 − ω ∗ (i ; zt , Kt )) α(i , kj ; zt , Kt ),where ω ∗ (i ; zt , Kt ) =kj 0 +1 −k∗ (i ;zt ,Kt )0

kj 0 +1 −kj 00ˆ i , k; zt , Kt ),and when ξ > ξ(0

P(i0 , j10 |i, j) =π (i0 |i )ω c (i , kj ; zt , Kt ) (1 − α(i , kj ; zt , Kt )) ,P(i0 , j10 + 1|i, j) =π (i0 |i ) (1 − ω c (i , kj ; zt , Kt )) (1 − α(i , kj ; zt , Kt )) ,where ω c (i , kj ; zt , Kt ) =kj 0 +1 −kc (i ,kj ;zt ,Kt )1

kj 0 +1 −kj 0
1

1

49

.

(n)

(n)

3. Update the forecasting rules Kt+1 = ΓK (zt , Kt ) and pt = Γp (zt , Kt ) by using the sequenceof {zt , Kt , pt }Tt=t0 +1 obtained in Step 2 for t = 1, ..., T (the first t0 periods are discarded).
Specifically, using the sequence obtained, we construct a subset of the sequence sorted by therealization of zt for each zi ∈ {z1 , ..., znz } and estimate the forecasting rule by ordinary leastsquares in the form oflog Kt+1 = bK,0 (zi ) + bK,1 (zi ) log Kt ,log pt = bp,0 (zi ) + bp,1 (zi ) log Kt ,for each i = 1, ..., nz .
Outer loop in the Xpa algorithm In the Xpa algorithm, aggregate technology zm and aggregate capital Km are given as a grid point for m = 1, ..., ne nK .30 Taking the expected valuefunctions as given, the fixed point problem pm = F (pm ) is solved as in the KS algorithm but at
each grid point. The market clearing price is solved for nz nK grid points in the Xpa algorithm,
whereas it is solved in each period of simulations for T periods in the KS algorithm. The Xpa
algorithm is faster than the KS algorithm as nz nK << T holds.
Taking the value function v(i , k; z, K) obtained in the inner loop as given, the outer loop in
the Xpa algorithm is as follows.
1. Calculate the bias–correction terms

ζk (i ) = Kss (i ) − ḡk (Kss (i )),ζy (i ) = Yss (i ) − ḡy (Kss (i )),where Kss (i ) =P

j kj µ̄k (kj |i ) =P

j kj µ̄(i , kj )/φ(i ).

Also compute ψ(i ) = Kss (i )/Kssfor the –indexed aggregation.
30

We explicitly consider the aggregate state variables (zm , Km ) as grid points and use subscripts m instead oft. The next period’s aggregate capital is denoted by K 0 instead of Kt+1 . The other aggregate variables are also
computed at each grid point, so they are indexed by m rather than time t.

50

2. At each grid point (zm , Km ) indexed by m = 1, ..., nz nK ,(a) The aggregate capital Km is given as a grid point.
(n−1)

0 = Γ
(b) Taking v(i , k; zm , Km ) and KmK(zm , Km ) as given, make a conjecture of p̃m(which is not necessarily the forecasted price) and solve for gy (i , k; zm , Km ) and gk (i , k; zm , Km )by the same procedure as in Step 2 for the inner loop. Then evaluate the individual decision rules at k = K(i ) ≈ ψ(i )Km for each i and calculate the weighted average toobtain

Ym =

X

0
Km

X

[gy (i , ψ(i )Km ; zm , Km ) + ζy (i )] φ(i ),i=[gk (i , ψ(i )Km ; zm , Km ) + ζk (i )] φ(i ),i
0 + (1 − δ)K , and p = C −1 . That is, there is a mapping denoted by
Cm = Ym − γKmmmmpm = F (p̃m ).
(c) Iterate Step 2(b) to find the market–clearing price p∗m satisfying p∗m = F (p∗m ).
(n)

(n)

0 , p ) obtained
3. Update the forecasting rules K 0 = ΓK (z, K) and p = Γp (z, K) by using (Kmmin Step 2 at each grid point (zm , Km ).
Numerical setup

The individual capital k ranges between [0.1,5.0] and is spaced into nk = 101

grid points. The aggregate capital K ranges between [0.75Kss ,1.25Kss ] around the steady statevalue of Kss and is evenly spaced into nK = 5 grid points. The distribution of k is approximatedby ñk = 2001 grid points by the non–stochastic method with histogram (Young, 2010). For theidiosyncratic and aggregate shocks, ne = 5 and nz = 5 grid points are used with Tauchen’s (1986)
method to approximate the AR(1) processes. We set the grid points to cover 99% (95%) of the
distribution implied by the original AR(1) process for the idiosyncratic (aggregate) productivity
shock.31 A two–dimensional cubic spline is used to approximate the expected value function. To
solve for k ∗ and k c in (19) and (20), a variant of Newton’s method or the golden section search
31

That is, ne = −1 = 2.58σ /pp
1 − ρ2 and znz = −z1 = 1.96σz / 1 − ρ2z .

51

method is used. Convergence criteria are 10−4 for the outer loop, 10−5 for the inner loop, 10−6 formarket clearing in the outer loop, and 10−10 for the other optimizations. The initial forecasting
rules are from the social planner’s problem. In the KS algorithm, the model is simulated for 2,500
periods and the first 500 periods are discarded. The same sequence of {zt } throughout all the outer
loops is used to avoid chatter, which prevents the algorithm from converging within a finite sample.

A.3

More numerical results
Table 13: Business cycle statistics, traditional model.

Standard deviationKS
Xpa
Output correlation
KS
Xpa

Y

C

I

N

K

Z

(2.5293)(2.5134)

0.43290.43703.61093.59010.66800.66190.51690.51420.59810.60191.00001.00000.87620.88070.97400.97400.93790.93740.0451
0.04470.99950.9995

Persist.

Std.dev.

Skewness

Exc.kur.

0.65710.65900.01080.01060.0747
0.0619-0.0968-0.0946

Agg. investment rate, I/K
KS
Xpa

52

Figure 8: Unconditional business cycle simulation, traditional model.

Output
Xpa
KS

−0.2
Log

Investment−1.125−1.150−2.6
1 10 20 30 40 50
1 10 20 30 40 50
1 10 20 30 40 50

Capital

0.22
0.20
Log

−1.100

−2.4

−0.6
−0.8

−1.075

−2.2

−0.4

Labor

−1.050

0.18
0.16
0.14

Consumption

TFP

−0.86

0.02

−0.88

0.00

−0.90

−0.02

−0.92

−0.04

1 10 20 30 40 50
Year

1 10 20 30 40 50
Year

1 10 20 30 40 50
Year

Figure 9: Impulse responses, traditional model.

OutputXpa
KS

Percent2

Percent

60.0

05 10 15 20

0CapitalConsumption

1.5

1.00

1.0

0.75

5 10 15 20

5 10 15 20
Year

0.00

5 10 15 20

TFP

0.5

0.25
0

0

1.0

0.50

0.5
0.0

0.5

2
0

Labor

1.5
1.0

4

1
0

Investment

8

0

5 10 15 20
Year

53

0.0

0

5 10 15 20
Year

Table 14: Level differences between KS and Xpa, traditional model.
a. Aggregate variables
K/YN

p = C −1KS

2.34760.33372.4002Xpa

2.34890.33372.3998SS

2.34870.33372.3998Mean

Std.dev.

Inaction

Spike+

Spike-

Invest+Invest-KS0.10460.02140.21150.12790.00000.78850.0000Xpa

0.10460.02140.21150.12800.00000.78850.0000SS

0.10460.02140.21130.12840.00000.78870.0000b. Firm–level investment rates

Figure 10: den Haan fundamental plots, traditional model.

KS:Capital

0.4

−0.80
Dynamic
StaticActual

0.0
1

10

20

30

40

Log

Log

0.2

−0.2

−0.85
−0.90

50

Xpa:Capital

0.4

−0.95

1

10

20

30

40

50

Xpa:Consumption

−0.75
−0.80
Log

Log

0.2
0.0
−0.2

KS:Consumption

−0.75

−0.85
−0.90

1

10

20

30

40

50

54

−0.95

1

10

20

30

40

50

B

Details of heterogeneous–household models

B.1

Stationary equilibrium

We explain how to compute the stationary equilibrium in the following steps.
1. Make an initial guess of K and obtain w = (1 − α)(K/L)−α and r = α(K/L)α−1 − δ where
L = 1 − (Πl u(zl ) + (1 − Πl )u(zh )).
2. Taking K, w and r as given, solve for the value function and the individual decision rules(a) Make an initial guess of the value function v (0) (xi , k).
(b) Taking v (l−1) (i , k) as given (l is an index for iteration in the inner loop), approximate
the expected value function and solve0v (l) (xi , k) = max
u
(1
−
τ
)wy
e
+
(1
+
r)k/θ
−
k
i
i
0
k ≥0

+

θβi

XX
e0

πe (e0 |ei )πy (y 0 |yi )v (l−1) (e0 , y 0 , βi , k 0 )




,



y0

for the individual decision rule k 0 = ḡ(xi , k).
(c) Iterate Step (b) until v (l) − v (l−1) is sufficiently small.
3. Solve for the stationary distribution µ̄(xi , k) implied by the individual decision rules.
4. Given the stationary distribution, the individual decision rule is aggregated as

K=

XZḡ(xi , k)θµ̄(xi , k)dk,

iwhere θ is the survival probability from period t to period t + 1. Note that the replacing
newborns have no assets.
5. Iterate Steps 2–4 until convergence.

55

The transition matrix for employment status is conditioned on (zt , zt+1 ) and denoted byπe (ei0 |ei , zt+1 , zt ) = 

pzt zt+1
1 − qzt zt+1

1 − pzt zt+1 ,
qzt zt+1

where pzt zt+1 is the conditional probability to continue to be unemployed and qzt zt+1 is the conditional probability to continue to be employed. Then the transition matrix for the stationary
equilibrium is given byπe (ei0 |ei ) = 

1−p 
,
1−q
q
p

where
(1 − Πl )u(zh )Πl u(zl )(ρh pzh zh + (1 − ρh )pzh zl ) +((1 − ρl )pzl zh + ρl pzl zl ) ,ussussΠl (1 − u(zl ))(1 − Πl )(1 − u(zh ))(ρh qzh zh + (1 − ρh )qzh zl ) +((1 − ρl )qzl zh + ρl qzl zl ) ,q=
1 − uss
1 − ussp=

and uss = (1 − Πl )u(zh ) + Πl u(zl ). Note that the stationary distribution of employment status
implied by πe (e0 |e) coincides with the employment rate 1 − uss .
In Step 3, given the current period’s distribution µ(xi , kj ) for i = 1, ..., nx and j = 1, ..., ñk , thenext period’s distribution µ(i0 , kj 0 ) for i0 = 1, ..., ne and j 0 = 1, ..., ñk is obtained from transition
probabilities from the current period’s index (i, j) to the next period’s index (i0 , j 0 ). With probability (1 − θ), households die and newly born households with the same properties xi0 = (i0 , yi0 , βi )and no assets replace the old. The transition probability from (i, j) to (i0 , 1) isP(i0 , 1|i, j) =πe (ei0 |ei )πy (yi0 |yi )(1 − θ).
With probability θ, households continue to live in the next period and choose k 0 = ḡ(xi , kj ). Then56

the transition probabilities from (i, j) to (i0 , j 0 ) and (i0 , j 0 + 1) areP(i0 , j 0 |i, j) =πe (ei0 |ei )πy (yi0 |yi )ω ∗ (xi , kj )θ,P(i0 , j 0 + 1|i, j) =πe (ei0 |ei )πy (yi0 |yi ) (1 − ω ∗ (xi , kj )) θ,
where ω ∗ (xi , kj ) =kj 0 +1 −ḡ(xi ,kj )0

kj 0 +1 −kj 00. If j 0 = 1, the transition probability from (i, j) to (i0 , 1) is the sum

0of two probabilities. That is, P(i0 , 1|i, j) = πe (ei0 |ei )πy (yi0 |yi ) (ω ∗ (xi , kj )θ + (1 − θ)).

B.2

Xpa and KS algorithms

Inner loop

Taking as given the forecasting rule Kt+1 = Γ(n−1) (zt , Kt ), the common inner loop

for the two algorithms is as follows.
1. Make an initial guess of the value function v (0) (xi , k; zt , Kt ).
2. Taking v (l−1) (xi , k; zt , Kt ) and Kt+1 = Γ(n−1) (zt , Kt ) as given (l is an index for iteration in
the inner loop),(a) Approximate the expected value function in the next period0h(l−1)
xi zt (k , Kt ) 'XXX
zt+1 e0πz (zt+1 |zt )πe (e0 |ei , zt+1 , zt )πy (y 0 |yi )y0× v (l−1) (e0 , y 0 , βi , k 0 ; zt+1 , Γ(zt , Kt ))for each xi and zt .
(b) Solvev (l) (xi , k; zt , Kt ) =max
0c,k ≥0no
0u(c) + θβi h(l−1)
xi zt (k , Kt )subject to the budget constraintc + k 0 = (1 − τ (zt ; ρ))wt yi ei + (1 + rt )k/θ,57

where wt = (1 − α)zt (Kt /L(zt ))−α and rt = α(Kt /L(zt ))α−1 − δ. Note that L(zt ) =1 − u(zt ) is exogenously given.
3. Iterate Step 2 until v (l) − v (l−1) is sufficiently small.

Outer loop in the KS algorithm

Taking the value function v(xi , k; zt , Kt ) obtained in the

inner loop as given, the outer loop in the KS algorithm is as follows.
1. Make a sequence of {zt } for t = 1, ..., T and an initial distribution µ1 (xi , kj ) approximated by
P P
a histogram at each grid point for i = 1, ..., nx and j = 1, ..., ñk so that i j µ1 (xi , kj ) = 1
holds.
2. In each period t = 1, ..., T , zt and the current period’s distribution µt (xi , kj ) are given.
(a) The aggregate capital is given by Kt =

P P
i

32
j kj µt (xi , kj ).

(b) Taking v(xi , kj ; zt , Kt ) and Kt+1 = Γ(n−1) (zt , Kt ) as given, solve for g(xi , kj ; zt , Kt ) by
the same procedure as in Step 2 for the inner loop. Then aggregate the individual
decision rules with the distribution to obtain

Kt+1 =XXig(xi , kj ; zt , Kt )θµt (xi , kj ).

jUpdate the next period’s distribution µt+1 (xi , kj ) by µ
~ t+1 = P~µt where µ~ t is a stacked
vector of the distribution in period t and P is a transition matrix with the transition
probabilities. With probability (1 − θ), the transition probability from (i, j) to (i0 , 1) isP(i0 , 1|i, j) =πe (ei0 |ei , zt+1 , zt )πy (yi0 |yi )(1 − θ).
32

Note this may not be on the grid points. Therefore we need to approximate g(i , kj ; zt , Kt ) obtained in the inner
loop or redo the optimization.

58

With probability θ, the transition probabilities from (i, j) to (i0 , j 0 ) and (i0 , j 0 + 1) areP(i0 , j 0 |i, j) =πe (e0 |ei , zt+1 , zt )πy (y 0 |yi )ω(xi , kj ; zt , Kt )θ,P(i0 , j 0 + 1|i, j) =πe (e0 |ei , zt+1 , zt )πy (y 0 |yi ) (1 − ω(xi , kj ; zt , Kt )) θ,

where ω(xi , kj ; zt , Kt ) =kj 0 +1 −g(xi ,kj ;zt ,Kt ).
kj 0 +1 −kj 0(n)3. Update the forecasting rule Kt+1 = ΓK (zt , Kt ) by using the sequence of {zt , Kt }Tt=t0 +1 obtained in Step 2 for t = 1, ..., T (the first t0 periods are discarded). Specifically, using thesequence obtained, we construct a subset of the sequence sorted by the realization of zt foreach zi ∈ {z1 , ..., znz } and estimate the forecasting rule by ordinary least squares in the formoflog Kt+1 = bK,0 (zi ) + bK,1 (zi ) log Kt ,for each i = 1, ..., nz .
Outer loop in the Xpa algorithm Taking the value function v(xi , k; z, K) obtained in the
inner loop as given, the outer loop in the Xpa algorithm is as follows.
1. Calculate the bias–correction terms

ζ(xi ) = Kss (xi ) − ḡk (Kss (xi )),where Kss (xi ) =P

j kj µ̄k (kj |xi ) =P

j kj µ̄(xi , kj )/φ(xi ). Also compute ψ(xi ; zm ) = Kss (xi ; zm )/Kssfor the –indexed aggregation where Kss (xi ; zt ) =Pj kj µ̄k (kj |xi ) =Note that the bias–correction terms do not depend on zm .
2. At each grid point (zm , Km ) indexed by m = 1, ..., nz nK ,(a) The aggregate capital Km is given as a grid point.

59

P

j kj µ̄(xi , kj )/φ(xi ; zm ).

0 = Γ(n−1) (z , K ) as given, solve for g(x , k; z , K ) by(b) Taking v(xi , k; zm , Km ) and Km
m
m
i
m
m

the same procedure as in Step 2 for the inner loop. Then evaluate the individual decisionrules at k = K(xi ) ≈ ψ(xi ; zm )Km for each xi and calculate the weighted average toobtain0Km=X[g(xi , ψ(xi ; zm )Km ; zm , Km )θ + ζ(xi )] φ(xi ; zm ).

i

0 obtained in Step 2 at each grid
3. Update the forecasting rules K 0 = Γ(n) (z, K) by using Km

point (zm , Km ).
Numerical setup The individual capital k ranges between [0.0, kmax ] and is spaced into nk = 101
grid points, where kmax = 2000 for the benchmark calibration, kmax = 700 for the +σ(y) calibration,and kmax = 250 for the KS calibration. More grid points are given toward small values considering
concavity of the value function in k. The aggregate capital K ranges between [0.80Kss ,1.15Kss ]around the steady state value of Kss and is evenly spaced into nK = 5 grid points. The distributionof k is approximated by ñk = 2001 grid points. For the labor productivity shock, ny = 7 grid points
are used with Rouwenhorst’s (1995) method to approximate the AR(1) processes. The discountfactor shock takes either of the values in nd = 3 grid points {β̄ − , β̄, β̄ + }. A two–dimensional
cubic spline approximates the expected value function. To solve for k 0 , a variant of Newton’s
method or the golden section search method is used. Convergence criteria are 10−4 for the outerloop, 10−5 for the inner loop, and 10−10 for the other optimizations.

60More numerical results
Figure 11: Unconditional business cycle simulation, KS calibration.

CapitalConsumption

3.68

1.05

3.67

1.04
Log

Log

3.69

3.66
3.65

1.03
1.02

3.64

Xpa
KS

1.01
1

40

80

120 160 200

140

Output120 160 200

0.00
−0.05

1.32

Log

Log

80

Investment

1.34

1.30

−0.10
−0.15

1.28
1

40

80

−0.20

120 160 200

1

40

80

120 160 200

Figure 12: Impulse responses, KS calibration.

TFP

−1
−2
−3
−4

−0.5
−1.0
Xpa
KS

−1.5
0

1

2

3

4

5

0

Output0

Percent

−4

1

2

3

4

5

4

5

Investment

0

−2

−6

Consumption

0.0
Percent

Percent0

Percent

B.3

−5
−10
−15

0

1

2

3

4

−20

5

61

0

1

2

3

Figure 13: Unconditional business cycle simulation, +σ(y) calibration.

Consumption

3.85

1.07

3.84

1.06

Log

Log

Capital

3.83
3.82

1.05
Xpa
KS

1.04

3.81
1

40

80

1.03

120 160 200

1

40

1.40

0.15

1.38

0.10

1.36

0.05

1.34

0.00
1

40

80

120 160 200

Investment

Log

Log

Output80

120 160 200

1

40

80

120 160 200

Figure 14: Impulse responses, +σ(y) calibration.

TFP

−1
−2
−3
−4

−1.0

0

1

2

3

4

5

Xpa
KS
0

Output

−2
−4

1

2

3

4

5

4

5

Investment

0
Percent

Percent

−0.5

−1.5

0

−6

Consumption

0.0
Percent

Percent

0

−5
−10
−15

0

1

2

3

4

5

0

62

1

2

3

Figure 15: den Haan fundamental plots, KMP calibration.

KS:Capital

3.80

Log

3.79
3.78
3.77

Dynamic
StaticActual

3.76
1100

1200

1300

1400

1500

1400

1500

Xpa:Capital

3.80

Log

3.79
3.78
3.77
3.76
1100

1200

1300

Figure 16: den Haan fundamental plots, KS calibration.

KS:Capital

3.69
3.68
Log

3.67
3.66

Dynamic
StaticActual

3.65
3.64
1100

1200

1300

1400

1500

1400

1500

Xpa:Capital

3.69

Log

3.68
3.67
3.66
3.65
3.64
1100

1200

1300

63

Figure 17: den Haan fundamental plots, +σ(y) calibration.

KS:Capital
3.85
Log

3.84
3.83

Dynamic
StaticActual

3.82
3.81
1100

1200

1300

1400

1500

1400

1500

Xpa:Capital
3.85
Log

3.84
3.83
3.82
3.81
1100

1200

1300

64