Yann Algan
email: yann.algan@ens.fr.

Oliver Allais
email: olivier.allais@ivry.inra.fr.

Wouter J Den Haan
email: wdenhaan@uva.nl.

Des Ponts

E T Chaussées -École

Normale Supérieure

Chris Carroll

Ken Judd

Michel Juillard

Per Krusell

Michael Reiter

Solving heterogeneous-agent models with parameterized cross-sectional distributions

Keywords: Incomplete markets, numerical solutions, projection methods, simulations Incomplete markets, numerical solutions, projection methods, simulations JEL Classification: C63, D52 Résumé

A new algorithm is developed to solve models with heterogeneous agents and aggregate uncertainty that avoids some disadvantages of the prevailing algorithm that strongly relies on simulation techniques and is easier to implement than existing algorithms. A key aspect of the algorithm is a new procedure that parameterizes the cross-sectional distribution, which makes it possible to avoid Monte Carlo integration.

The paper also develops a new simulation procedure that not only avoids cross-sectional sampling variation but is also more than ten times faster than the standard procedure of simulating an economy with a large but finite number of agents. This procedure can help to improve the efficiency of the most popular algorithm in which simulation procedures play a key role.

SOLVING HETEROGENEOUS-AGENT MODELS WITH PARAMETERIZED CROSS-SECTIONAL DISTRIBUTIONS

Yann ALGAN, Oliver ALLAIS, and Wouter J. DEN HAAN * December 19, 2006 1 Introduction Models with heterogeneous agents and aggregate uncertainty are becoming increasingly important. They not only improve the predictions of representative agent models, they also make it possible to study the behavior of sub groups in a general equilibrium framework. Solving such models is difficult, because the set of state variables contains the cross-sectional distribution of agents' characteristics, which is a time-varying infinite dimensional object in the presence of aggregate uncertainty.

The most commonly used algorithm summarizes the cross-sectional distribution with a finite set of moments and calculates the transition law for these state variables using a simulation procedure. [START_REF] Christiano | Algorithms for solving dynamic models with occasionally binding constraints[END_REF] This algorithm is relatively easy to implement, but there are disadvantages to using simulation methods. First, it relies on Monte Carlo integration to calculate the cross-sectional moments, which is known to be an inefficient integration procedure. Second, the observations used to find the transition law are distributed inefficiently. [START_REF] Haan | Solving the stochastic growth model by parameterizing expectations[END_REF] Den Haan (1997) and Reiter (2002) develop alternative algorithms based on traditional approaches in the numerical solutions literature that avoid these disadvantages. These algorithms are, however, quite cumbersome to implement. Krusell and Smith (1998) show that in their heterogeneous-agent model, differences between agents are small in the sense that the marginal propensity to save is very similar among agents. The marginal propensity is only different for those agents that are at or close to the borrowing constraint. There are, however, not many of these agents and they have at most a marginal effect on prices, since their capital holdings are low. Changes in the distribution through wealth redistributions have, therefore, little effect. Consequently, it is possible to solve the model of Krusell and Smith (1998) accurately using only the mean capital stock to characterize the crosssectional distribution. Not only is heterogeneity limited in the model of Krusell and Smith (1998), the law of motion for capital can be accurately described by a linear function. Although still not a trivial exercise, the model is, thus, relatively easy to solve and the disadvantages associated with simulation procedures are unlikely to surface. In fact, we find results that are very similar to those obtained by Krusell and Smith (1998) by using a procedure in which simulation procedures only play a minor role.

Being able to accurately characterize the importance of changes in the crosssectional distribution for agents' decisions with a linear law of motion of the mean is unlikely to remain true in more complex models. [START_REF] Haan | Heterogeneity, aggregate uncertainty and the shortterm interest rate[END_REF] As the heterogeneous-agent models being solved are becoming more complex, the disadvantages of simulation procedures are likely to emerge, just like they do in other numerical problems. In Den Haan and Marcet (1990), one of the authors of this paper proposed a simulation procedure to solve dynamic stochastic models, but personal experience has made clear that there are also quite a few models for which an accurate numerical solution requires alternative techniques. [START_REF] Haan | Solving dynamic models with aggregate shocks and heterogeneous agents[END_REF] It is, thus, important to develop alternative algorithms to solve models with heterogeneous agents in which simulation procedures do not play a key role. Even if the class of heterogeneous-agent models that can be solved with simultationbased algorithms is large, then it is still important that alternative algorithms are developed. The reason is that assessing the accuracy of a numerical solution is hampered by the fact that the true solution is not known and there are many aspects of the solution to evaluate. Consequently, it is not easy to find out which models can be solved accurately with simulation procedures. Alternative solution algorithms are an important tool to determine for which models this is the case and to build (or reduce) confidence in the existing results in the literature that are obtained using simulation procedures.

This paper develops an algorithm that builds on the ideas of Den Haan (1997) and Reiter (2002), but has important efficiency advantages in terms of programming burden and computational efficiency. A key feature of the algorithm is the parameterization of the cross-sectional distribution. This makes it possible to use quadrature instead of Monte Carlo techniques to numerically integrate and to use projection methods to solve for the transition laws. By using reference moments the algorithm can obtain a much more accurate shape of the cross-sectional distribution without using additional state variables, which is an important improvement upon Den Haan (1997). We use a class of polynomials to approximate the cross-sectional distribution. This makes the algorithm easier to implement than the algorithm of Reiter (2002) who uses step functions. More importantly, by using a particular class of approximating polynomials, the problem of finding the coefficients of the approximating function can be reduced to a convex optimization problem, for which reliable algorithms exist. In contrast, for other approximating functions one has to find the coefficients using a nonlinear system of equations, which is a more difficult numerical problem with less reliable convergence properties.

We show how to use the algorithm to solve the model of Krusell and Smith (1997). Papers that describe solution procedures do not make bedtime reading. We hope to increase accessibility of our paper by going through the steps of algorithm using a well known model as an example.

The paper proposes a procedure to simulate an economy with heterogeneous agents that avoids cross-sectional sampling variation and, thus, eliminates an important drawback of simulation procedures. The standard simulation procedures constructs a panel of N T observations and N N agents. By parameterizing the crosssectional distribution and using quadrature integration, however, it is possible to generate an accurate simulation with a continuum of agents. A key aspect of our simulation procedure is to use a continuum instead of a finite number of agents. Note that models with a large number of heterogeneous agents almost always assume a continuum of agents, so that the law of large numbers ensures that idiosyncratic risk is averaged out. In fact, the assumption of a continuum of agents plays a key role, not only in the specification of the state variables and the definition of the equilibrium, but also in the construction of most algorithms. [START_REF] Judd | Projection methods for solving aggregate growth models[END_REF] Thus, by simulating with a continuum instead of a finite number of agents, we stay much closer to the actual model being solved.

Our simulation procedure not only avoids cross-sectional sampling variation, it is also much cheaper. An important role in the reduction of computing time is using our proposed class of approximating functions so that coefficients can be found with a convex optimization procedure. We found that simulating an economy with 10,000 agents for 1000 periods took 13 times as long as simulating the same economy with a continuum of agents. In our own algorithm, the simulation procedure only plays a very minor role. This is, of course, no reason not to use a more accurate simulation procedure. For algorithms that use a simulation procedure to calculate the transition law, however, the improved simulation procedure will have bigger benefits. In the last section, we discuss for what type of model the new simulation procedure is likely to be especially beneficial.

The rest of this paper is organized as follows. The next section describes the production economy of Krusell and Smith (1998). Section 3 briefly discusses existing algorithms and summarizes the contributions of this paper. Section 4 describes the algorithm in detail and Section 5 describes the simulation procedure. Section 6 discusses how to check for accuracy and reports the results. The last section concludes.

The production economy

The economy is a production economy with aggregate shocks in which agents face different employment histories and partially insure themselves through (dis)saving in capital. For more details see Krusell and Smith (1998).

Problem for the individual agent. The economy consists of a unit mass of ex ante identical households. Each period, agents face an idiosyncratic shock ε that determines whether they are employed, ε = 1, or unemployed, ε = 0. An employed agent earns a wage rate of w t . An employed agent earns an after-tax wage rate of (1τ t ) w t and an unemployed agent receives unemployment benefits µw t . Markets are incomplete and the only investment available is capital accumulation. The net rate of return on this investment is equal to r tδ, where r t is the rental rate and δ is the depreciation rate. Agent's i maximization problem is as follows:

max {c i t ,k i t+1 } ∞ t=0 E P ∞ t=0 β t (c i t ) 1-γ -1 1-γ s.t. c i t + k i t+1 = r t k i t + (1 -τ t )w t lε i t + µw t (1 -ε i t ) + (1 -δ)k i t k i t+1 ≥ 0 (1)
Here c i t is the individual level of consumption, k i t is the agent's beginning-of-period capital, and l is the time endowment.

Firm problem. Markets are competitive and the production technology of the firm is characterized by a Cobb-Douglas production function. Consequently, firm heterogeneity is not an issue. Let K t and L t stand for per capita capital and the employment rate, respectively. Per capita output is given by

Y t = a t K α t (lL t ) 1-α (2)
and prices by

w t = (1 -α) a t ( K t lL t ) α (3)
r t = αa t ( K t lL t ) α-1 (4)
Aggregate productivity, a t , is an exogenous stochastic process that can take on two values, 1 -∆ a and 1 + ∆ a .Government

The only role of the government is to tax employed agents and to redistribute funds to the unemployed. We assume that the government's budget is balanced each period. This implies that the tax rate is equal to

τ t = µu t lL t . ( 5
)
where u t = 1 -L t denotes the unemployment rate in period t.

Government The only role of the government is to tax employed agents and to redistribute funds to the unemployed. We assume that the government's budget is balanced each period. This implies that the tax rate is equal to

τ t = µu t lL t . ( 6
)
where u t = 1 -L t denotes the unemployment rate in period t.

Exogenous driving processes. There are two stochastic driving processes. The first is aggregate productivity and the second is the employment status. Both are assumed to be first-order Markov processes. We let π aa 0 εε 0 stand for the probability that a t+1 = a 0 and ε i t+1 = ε 0 when a t = a 0 and ε i t = ε 0 . These transition probabilities are chosen such that the unemployment rate can take on only two values. That is, u t = u b when a t = a b and u t = u g when a t = a g with u b > u g . Equilibrium. Krusell and Smith (1998) consider a recursive equilibria in which the policy functions of the agent depend on his employment status, ε i , his beginningof-period capital holdings, k i , aggregate productivity, a, and the cross-sectional distribution of capital holdings. An equilibrium consists of the following • Individual policy functions that solve the agent's maximization problem.

• A wage and a rental rate that are determined by 3 and 4, respectively.

• A transition law for the cross-sectional distribution of capital, that is consistent with the investment policy function.

Relation to existing algorithms

In this paper we give a brief overview of existing algorithms and then the contributions of this paper.

Existing algorithms

A standard aspect of numerical algorithms that solve models with heterogeneous agents is to summarize the infinite-dimensional cross-sectional distribution of agents' characteristics by a finite set of moments, m. The transition law is then a mapping that generates next period's moment, m 0 , given the values of the moments in the current period, m, and the realization of the aggregate shock. Den Haan (1996), Krusell and Smith (1998), and Rios-Rull (1997) propose to calculate this transition law as follows. First, construct a time series for the cross-sectional moments by simulating an economy with a large but finite number of agents. Second, regress the simulated moments on the set of state variables. The first disadvantage of the simulation procedure is that moments are calculated using Monte Carlo integration, which is known to be an inefficient numerical integration procedure. [START_REF] Judd | Numerical methods in economics[END_REF] The second disadvantage is that the observations are clustered around the mean, since they are taken from a simulated series. An efficient projection procedure, however, requires the explanatory variables to be spread out, for example, by using Chebyshev nodes. [START_REF] Krueger | Computing equilibrium in OLG models with stochastic production[END_REF] Den Haan (1997) parameterizes the cross-sectional distribution with a flexible functional form, P (k; ρ), which makes it possible to use quadrature techniques to do the numerical integration. In addition, his algorithm uses Chebyshev nodes to construct a grid of explanatory variables for the projection step. The coefficients of the approximating density, ρ, are pinned down by the set of moments used, m. The disadvantage of Den Haan (1997) is that the shape of the distribution is completely pinned down by the moments used as state variables and the class of flexible functional forms used. Consequently, a large number of state variables may be needed, not to predict changes in the key moments that matter for agents' behavior, but just to get the shape of the cross-sectional distribution right. Another drawback of Den Haan (1997) is that an inefficient procedure is used to find the coefficients of the approximating density. Reiter (2002) improves upon the algorithm of Den Haan (1997) in an ingenious way by letting the shape of the distribution depend not only on the moments used as state variables, m, but also on a set of reference moments that are obtained by a simulation procedure.

Promising recent alternatives to the standard algorithm have been developed in Preston and Roca (2006) and Reiter (2006). Reiter (2006) first solves a model without aggregate uncertainty using standard projection procedures. Next, by replacing the endogenous variables in the equations of the model with the parameterized numerical solution, he obtains a difference equation in the numerical coefficients. Then he uses perturbation techniques to solve for the sensitivity of the numerical solution to aggregate shocks. This is quite a different approach then the procedure used here, which is good because the more variety among available approaches the better. Preston and Roca (2006) use a "pure" perturbation method to solve the model. [START_REF] Krusell | Income and wealth heterogeneity in the macroeconomy[END_REF] Perturbation methods are likely to work well when the distribution needs to be characterized by many statistics, because dealing with many arguments is the strength of perturbation methods.

The contributions of this paper

The main contributions of this paper are the following.

Calculating the transition law of the cross-sectional distribution. The disadvantage of Reiter (2002) is that the particular implementation of the idea of reference moments is very cumbersome. [START_REF] Miao | Competitive equilibria of economies with a continuum of consumers and aggregate shocks[END_REF] As in Den Haan (1997) and Reiter (2002), this paper develops a procedure to calculate this transition law without relying on simulation procedures to calculate moments and to carry out the projection step. As in Reiter (2002), it uses reference moments, but the modifications introduced make the procedure much more straightforward to implement.

Calculating the approximating density for given moments. This algorithm links a set of moments with a parameterized density. Consequently, an important part of the algorithm is the mapping between the set of moments and the coefficients of the density. One possibility would be to use an equation solver that chooses the set of coefficients so that the moments of the parameterized density are equal to the specified moments. This turns out to be a relatively expensive procedure. By using [START_REF] Krusell | Income and wealth heterogeneity in the macroeconomy[END_REF] Since perturbation methods cannot deal well with the kind of inequality constraint used here, they replace the inequality constraint with a penalty function, which also accomplishes that agents do not have negative capital holdings and try to stay away from low capital stocks. [START_REF] Miao | Competitive equilibria of economies with a continuum of consumers and aggregate shocks[END_REF] Reiter (2002) constructs a reference density G(m), which relates the shape of the distribution to the set of moments that serve as state variables. It is a weighted average of distributions from a simulated economy, where distributions with moments closer to m get more weight. Step functions are used to construct a reference distribution, which has the advantage of being very flexible but has the disadvantage of using a lot of parameters. One problem of the approach in Reiter (2002) is that the moments of the reference density may not be equal to m. This means that one first has to apply operations to ensure that one obtains a new reference function e G(m) for which this is not the case. But even if m contains only first and second moments, then this problem entails more than a linear transformation, since e G(m) has to be a step function that conforms to the specified grid and cannot violate the constraints on the support of the distribution, such as the constraint that k ≥ 0.

a specific class of approximating functions, it is possible to transform this problem into a convex optimization problem, for which much more reliable algorithms exist.

This procedure is likely to be useful outside the literature of numerical solution techniques, since characterizing a cross-sectional distribution with a CDF from a class of flexible functional forms is a common problem in econometrics.

Simulating a panel without cross-sectional sampling variation. This paper develops a procedure to simulate an economy without cross-sectional sampling variation. Standard procedure is to simulate data using a finite number of agents and a finite number of time periods, which means that the outcome depends on the particular random draw used. Sampling variation disappears at a slow rate and could be especially problematic if the number of a particular type of agent is small relative to the total number of agents.

Existing models with a large number of heterogeneous agents typically assume that there is a continuum of agents. [START_REF] Preston | Incomplete markets, heterogeneity and macroeconomic dynamics[END_REF] This implies that conditional on the realization of the aggregate shock there is no cross-sectional sampling variation, a property that plays a key role in the definition of the set of state variables and the definition of the equilibrium. The simulation procedure developed in this paper sticks to the assumption of the model and uses a continuum of agents. By parameterizing the cross-sectional distribution one can use quadrature techniques to obtain a time series for any moment or characteristic of the cross-sectional distribution.

Accuracy tests. It is never trivial to check the accuracy of a numerical solution for dynamic stochastic models, since the true solution is not known. Checking for accuracy is made especially difficult, because there are many aspects to the solution of this type of model. In this paper, we discuss several tests to evaluate the accuracy of the solution for a model with heterogeneous agents.

The algorithm

In this section, we discuss the different steps of the algorithm. In Section 4.1, we start with a discussion of the state variables used, followed by an overview of the algorithm in Section 4.2. The remaining sections describe the steps of the algorithm.

State variables and transition laws

Krusell and Smith (1998) consider a recursive equilibrium in which the policy functions of the agent depend on his employment status, ε i , his beginning-of-period capital holdings, k i , aggregate productivity, a, and the cross-sectional distribution of capital holdings. [START_REF] Reiter | Recursive computation of heterogeneous agent models[END_REF] Let ←f w (k) be the cross-sectional distribution of beginning-ofperiod capital holdings for agents with employment status w ∈ {e, u} and k ≥ 0. The arrow pointing left indicates that the cross-sectional distribution refers to the distribution at the beginning of the period (but after all shocks are observed). Similarly, -→ f w (k) refers to the distribution at the end of the period. The following two steps determine the transition law that links the current-period distribution, ← -

f w (k), with next period's distribution, ← - f w0 (k).
• The end-of-period distribution is determined by a, ←f e , ←f u , and the individual investment function. That is,

-→ f e = -→ Υ e (a, ← - f e , ← - f u ) and -→ f u = -→ Υ u (a, ← - f e , ← - f u ).
• Next period's beginning-of-period distribution, ←f w0 (k), is determined by the end-of-period distribution and the employment-status flows that correspond with the values of a and a 0 . Thus, ← -

f e0 = ← - Υ e (a, a 0 , -→ f e , -→ f u ) and ← - f u0 = ← - Υ u (a, a 0 , -→ f e , -→ f u ). ← - Υ e (•) and ← - Υ u (•)
are simple functions that are determined directly by the transition probabilities. [START_REF] Reiter | Solving heterogeneous agent models by projection and perturbation[END_REF] An alternative to using the cross-sectional distribution of employment and capital holdings is to use all past realizations of the aggregate shocks. [START_REF] Rios-Rull | Computation of equilibria in heterogeneous agent models[END_REF] For the model considered here, we found that a large number of lags is needed. Nevertheless, if one doesn't have a complete description of the cross-sectional distribution, it still may be worthwhile to add some lagged values of a. [START_REF] Veracierto | Plant-level irreversible investment and equilibrium business cycles[END_REF] In our algorithm we, therefore, add the lagged value of a as a state variable. But there is another reason, which will become evident in the remainder of this section.

To deal with the infinite dimension of the cross-sectional distribution, we follow Den Haan (1996, 1997), Krusell and Smith (1997), and Rios-Rull (1997) and describe the cross-sectional distribution with a finite set of moments. The remainder of this section discusses in detail which moments we use.

In this model, agents face a borrowing constraint, k ≥ 0. We, therefore, include the fraction of agents of each type that start the period with zero capital holdings, ← -m e,c and ← -m u,c . Employed agents never choose a zero capital stock. This means that the density -→ f e does not have mass at zero. In contrast, ←f e0 does have mass at zero, because some of the agents that are employed in the current period were unemployed in the last period and chose a zero capital stock. Both ← -m e,c and ← -m u,c can be easily calculated from --→ m u,c -1 and the employment-status flows corresponding to the values of a and a 0 . Thus, instead of using

h a, ← -- m e,c , ← -- m u,c i we also can use h a, a -1 , --→ m u,c -1 i . We prefer to use h a, a -1 , --→ m u,c -1 i
, because a -1 can take on only two values and is, thus, computationally an inexpensive state variable. Moreover, as explained above a -1 could have predictive value that goes beyond the ability to determine ← -m e,c and ← -m u,c . In addition, the algorithm uses centralized moments of the distributions of strictlypositive capital holdings . The set of moments that are used as state variables are stored in the following vector

m = h --→ m u,c -1 , ← -- m e,1 , • • • , ←--- m e,N M , ← -- m u,1 , • • • , ← ---- m u,N M i ,
where ← -m w,j is the j th -order centralized moment for workers with employment status w and strictly-positive capital holdings. The dimension of the vector is

N * M = 2N M + 1.
The aggregate state is thus given by s = [a, a -1 , m]. Since we only use a limited set of moments as state variables, the transition law only needs to specify how this limited set of moments evolve over time. Thus, instead of calculating -→ Υ e (•) and

-→ Υ u (•), we now calculate [ --→ m u,c , --→ m u,1 , • • • , ----→ m u,N M ] = -→ Γ u n (s; ψ Γ u n ) and [ --→ m e,c , --→ m e,1 , • • • , ---→ m e,N M ] = -→ Γ e n (s; ψ Γe n ), where -→ Γ w n (s) is an n th -order polynomial with coefficients ψ Γ w n .
To simplify the notation we will typically write -→ Γ w (s), but one should keep in mind that this is an approximating function with coefficients that are determined by the algorithm. In the implementation of the algorithm, we set-as in Krusell and Smith (1998)-N M equal to 1.

Overview

An important part of this algorithm is to avoid Monte Carlo integration by approximating the densities ←f e and ←f u with a flexible functional forms. To determine this functional form, we use the moments that are used as state variables, m, as well as some additional information that we will refer to as "reference moments". The reference moments are higher-order moments that are helpful in getting the shape of the distribution right.

The algorithm uses the following iterative procedure to solve the model.

• Given transition laws, -→ Γ e (s) and -→ Γ u (s), solve for c(ε, k, s) and k 0 (ε, k, s). This is discussed in Section 4.3.

• Use the solutions for the individual policy functions, c(ε, k, s) and k 0 (ε, k, s)

to obtain information about the "reference moments". This is discussed in Section 4.4.

• Given solutions for the individual policy functions, c(ε, k, s) and k 0 (ε, k, s), solve for -→ Γ e (s) and -→ Γ u (s). This is discussed in Section 4.6. This requires setting up a grid of the aggregate state variables [a, a -1 , m] and a procedure to link the values of the moments m and the reference moments with an explicit cross-sectional density. This procedure is discussed in Section 4.5.

• Iterate until the transitions laws used to solve for the individual policy functions are close to the transition laws implied by the individual policy functions.

Procedure to solve for individual policy functions

The procedure to solve for individual policy functions relies on standard projection methods, except that we modify the standard procedure to deal with the inequality constraint on capital. In this section, we describe how to solve for the individual policy rules taking the aggregate policy rules -→ Γ e (s) and -→ Γ u (s) as given. The firstorder conditions of the agent are given by15

c(ε, k, s) -γ = E [βc(ε 0 , k 0 , s 0 ) -γ (1 + r(s 0 ))] for k 0 > 0, c(ε, k, s) -γ ≥ E [βc(ε 0 , k 0 , s 0 ) -γ (1 + r(s 0 ))] for k 0 = 0, and (7)
c + k 0 = r(s)k + w(s)lε + (1 -δ)k, (8)
In this system w(s) and r(s) only depend on a and the aggregate capital stock and can be solved directly from Equations 3 and 4. The conditional expectation in Equation 7 is a function of the individual and aggregate state variables. To solve the individual problem we approximate this conditional expectation with a flexible functional form. That is,

E £ βc(ε 0 , k 0 , s 0 ) -γ (1 + r(s 0 )) ¤ ≈ Ψ n (k, ε, s; ψ E n ), (9)
where Ψ n (•) is an n th -order polynomial and ψ E n its coefficients. Let k(ε, s) be the capital stock such that

k 0 = 0 if k ≤ k(ε, s). ( 10
)
Then ∂k 0 /∂k = 0 for k < k(ε, s). This implies that E[βc(ε 0 , k 0 , s 0 ) -γ (1 + r(s 0 ))] as a function of k is flat for k < k and nondifferentiable at k = k. When k < k(ε, s) one does not need the approximation Ψ n (•), since k 0 = 0 and one can solve consumption from the budget constraint. To calculate the approximation for the conditional expectation, we only use grid points at which k 0 > k(ε, s). This means that the grid is no longer fixed within the algorithm and we loose some of the optimality properties of using Chebyshev grid points, but we found that with this procedure we can obtain more accurate solutions. Besides this modification, our procedure to solve for the individual policy rules is a standard application of projection methods as discussed in Judd (1992). In particular, we use the following procedure.

• Construct a grid for the values of individual and aggregate state variables.

• Use ψ as the initial value for ψ E n . Given the value ψ, it is straightforward to solve for c(ε, k, s) and k 0 (ε, k, s) from the first-order condition and the budget constraint.

• At each grid point calculate k 0 ,

• For all possible realizations of a 0 and ε 0 calculate βc(ε 0 , k 0 , s 0 ) -γ (1 + r(s 0 )).

This requires calculating m 0 but this is easy since -→ Γ e and -→ Γ u are given. 16 Next, calculate E[βc(ε 0 , k 0 , s 0 ) -γ (1 + r(s 0 ))] by weighting the possibly outcomes with the probabilities.

• Perform a projection to obtain a new estimate for ψ E n , b ψ.

• Use a weighted average of b ψ and ψ as a new initial value for ψ E n .

• Iterate until the coefficients have converged.

As pointed out by Reiter (2006), one doesn't need the law of motion for m 0 to solve the individual problem. Using the parameterized cross-sectional distribution, one could in principle use quadrature methods to directly calculate the values of m 0 when needed. By doing this one could in each iteration not only update the individual policy rules but also the law of motion for the aggregate state variables. This is likely to speed up the algorithm if it is on course towards the fixed point, but the simultaneous updating might make the algorithm less stable.

Procedure to generate reference moments

The reference moments are used to ensure that the functional form of the crosssectional distribution is appropriate without using too many moments explicitly as state variables. Note that an extra state variable increases the dimension of the grid and the set of arguments of each function, whereas an extra reference moment does not. Given the complexity of the system one has to rely on simulations to obtain information about the shape of the distribution. Thus, we also use a simulation procedure to obtain reference moments, but we propose a new simulation procedure that substantially reduces the amount of sampling variation. This new simulation procedure is discussed in Section 5.

The simulation generates a time series with for each observation a set of observations ← -m w,j t for w ∈ {e, u} and j

∈ {N M + 1, • • • , N M }.
The simplest way to proceed would be to use the sample averages as the reference moments, but we let the reference moments depend on a. At each grid point, we know the value of a and a -1 as well as the values of ← -m w,j for w ∈ {e, u} and j ∈ {c, 1, • • • , N M }. We also have a set of higher-order reference 17 Without complicating the algorithm, one could let the higher-order moments depend on all the elements of the aggregate state variables, s, that is

← --- m w,j = Φ w n (s), (11)
where Φ w n is an n th -order flexible functional form. One can obtain the coefficients of Φ w n by a simple regression and at each aggregate grid point it is trivial to use Φ w n to determine the set of reference moments. We have not done so, because Young (2005) points out that higher-order moments do not exhibit a clear relationship on average to lower-order moments. moments ← -m w,j , w ∈ {e, u} and j ∈ {N M + 1, • • • , N M }. Let P (k; ρ w ) be the exponential of an N th M -order polynomial with coefficients ρ w . One way to solve for ρ w is to solve the following system of N M + 1 equations and unknowns.

∞ Z 0 h k - ←-- m w,1 i P (k; ρ w )dk = 0 ∞ Z 0 h (k - ←-- m w,1 ) 2 - ←-- m w,2 i P (k; ρ w )dk = 0 • • • ∞ Z 0 h (k - ←-- m w,1 ) n 2 - ← ---- m w,N M i P (k; ρ w )dk = 0 ∞ Z 0 P (k; ρ w )dk = 1 (12)
It is not clear what for practical applications the convergence properties towards the solution are. In our experience, it is feasible to solve this system but also quite costly. By adopting a different class of approximating polynomials one can reduce this problem to a convex optimization problem for which algorithms with reliable convergence properties exist.

In particular, we parameterize the density with the following N th M -order polynomial:

P (k, ρ w ) = ρ w 0 exp        ρ w 1 h k - ←-- m w,1 i + ρ w 2 • ³ k - ←-- m w,1 ´2 - ←-- m w,2 ¸+ • • • + ρ w N M • ³ k - ←-- m w,1 ´NM - ← ---- m w,N M ¸       . ( 13
)
When the density is written in this particular way, the coefficients, except for ρ w 0 , can be found with the following minimization routine:

min ρ w 1 ,ρ w 2 ,••• ,ρ w N M ∞ Z 0 P (k, ρ w )dk. ( 14
)
The first-order conditions correspond exactly to the first N M equations in (12). ρ w 0 does not appear in these equations, but ρ w 0 is determined by the condition that the density integrates to one.

The (i, j) element of the matrix of second-order derivatives (times ρ w 0 ) is given by

∞ Z 0 • ³ k - ←-- m w,1 ´i - ← -- m w,i ¸•³ k - ←-- m w,1 ´j - ← -- m w,j ¸P (k, ρ w )dk (15)
Thus, the Hessian is a covariance matrix and, thus, positive semi-definite. Consequently, the solution found is a minimum. The fact that this is a convex optimization problem means that one can rely on algorithms with reliable convergence properties. Finding the coefficients of the approximating density by solving the optimization problem 14, turned out to be an order of magnitude faster than finding the coefficients by solving the system of equations given in 12.18

Procedure to solve for aggregate laws of motion

The procedure is characterized by the following steps.

• We construct a grid with values for a, a -1 , --→ m u,c -1 , and ← -m w,j for w ∈ {e, u} and

j ∈ {1, • • • , N M }. Here --→ m u,c
-1 is the fraction of unemployed agents that chose k 0 = 0 last period and ← -m w,j is the j th moment of the distribution of strictlypositive capital holdings for agents with employment status w. Given values for a, a -1 , and --→ m u,c -1 , we can calculate ← -m e,c and ← -m u,c . The grid values for this period's and last period's aggregate state are the two possible realizations and we use Chebyshev nodes to locate the grid points for the other state variables. These are the "x-values".

• Using quadrature methods, we calculate end-of-period moments, --→ m w,j for j ∈ {c, 1, • • • , N M }, at each grid point. These are the "y-values". The parameterization of the cross-sectional distribution discussed in the last section makes it possible to use Simpson quadrature to calculate end-of-period moments.

• Using the y-values and the x-values, we perform a projection step to find the coefficients of the approximating function,

-→ Γ e n (s; ψ Γe n ) and -→ Γ u n (s; ψ Γu n ).

Simulating a panel with a non-random crosssection

An important contribution of this paper is to develop a simulation procedure that eliminates the amount of cross-sectional sampling variation in the simulation of a panel. In the algorithm proposed here simulations play a relatively minor role and are only used to get information on the shape of the cross-sectional distribution. But this simulation procedure could be an important improvement for those algorithms that do rely on simulations to determine the law of motion of the aggregate state variables, such as, the algorithm used in Krusell and Smith (1998). This section discusses the new procedure. Accuracy tests for this procedure are proposed in the next section. The idea of the simulation procedure proposed is to stay close to the idea that there is a continuum of agents in the model. This implies that-conditional on the realization of the aggregate shock-there is no cross-sectional sampling variation with our simulation procedure, just as there is none in the true solution. The standard procedure in the literature uses a large but finite number of agents, N N . Typically a value of 10,000 is used. For the standard calibration this means that the number of (un)employed agents is equal to 9,600 (400) and 9,000 (1,000) in the boom and recession, respectively. Not surprisingly, given the small number of unemployed agents, the sampling variation to determine the moments of the capital holdings of the unemployed is enormous. In this model, the typical capital stock of an unemployed agent is small and it does not matter very much that these moments are estimated inaccurately. In models in which the minority is quantitatively important, however, inaccurate calculation of their average moments is likely to affect aggregate moments. For example, when the minority consists of entrepreneurs, then this sampling variability could turn out to be fatal for the accuracy of the algorithm. If one is interested in a particular subset of the economy then accurate estimation of their moments is important independent of their importance in either number or wealth levels.

The procedure works as follows

• Use a random number generator to draw a time series for the aggregate productivity shock, {a t } N T t=0 . Although the procedures eliminates the cross-sectional sampling variation, there is still sampling variation due to the stochastic nature of the aggregate productivity shock.

• In period 1, the procedure starts with the following information. First, the fractions of employed and unemployed agents with zero capital stock at the beginning of the period, ← -m e,c 1 , and ← -m u,c 1 . Second, N f M centralized moments of the distribution of strictly positive beginning-of-period capital holdings for the unemployed and the employed, ← -m w,j 1 for w ∈ {e, u} and j ∈ {1, • • • , N f M }. Since this procedure is relatively cheap and not part of a complex fixed-point calculation one can set N f M fairly high. We set N f M = N M = 6, but these parameters do not have to be equal to each other.

• The moments ← -m w,j 1 for w ∈ {e, u} and j ∈ {1, • • • , N f M } determine the densities of positive capital holdings for the employed and unemployed, P (k; ρ e 1 ) and P (k; ρ u 1 ). That is, using the procedure discussed in Section 4.5, we find the coefficients of the densities in period 1, ρ e 1 and ρ u 1 , so that the moments of P (k; ρ e 1 ) and P (k; ρ u 1 ) correspond to the specified moments.

• Use P (k; ρ e 1 ), P(k; ρ u 1 ), ← -m e,c 1 , and ← -m u,c 1 , i.e., the distribution of beginning-ofperiod capital holdings together with the individual policy rules to calculate the end-of-period moments, --→ m w,j 1 for w ∈ {e, u} and j ∈ {c, 1,

• • • , N f M }, and --→ m u,c
1 . We use Simpson quadrature to do the integration.

• Use the values of the productivity shocks in period 1 and 2, i.e., a 1 and a 2 , together with the end-of-period moments for period 1 to calculate beginningof-period moments for period 2, ← -m w,j • Use the procedure discussed in 4.5 to find the values for ρ e 2 and ρ u 2 .

• Repeat the procedure for the next period until t = N T .

To ensure that the sample used to obtain information about the cross-sectional distribution has reached (or is at least close to) its ergodic distribution one should disregard an initial set of observations. For the particular model we study in this paper, we found that if the initial distribution is not close to the ergodic set, then one has to disregard a large number of initial observations, since it can take quite a while before the economy has reached the ergodic distribution. After some experimentation, one has a good idea about a reasonable initial distribution and then this is less of a problem.

For some policy functions, it may be the case that some higher-order moments of the cross-sectional distribution do not exist or that higher-order moments are on an explosive path. In our numerical procedure we integrate over a finite range of capital holdings so this problem cannot occur. To make sure that the numerical procedure doesn't hide diverging properties of the true model it is important to check whether the results are robust to changing the upper bound of capital stock.

Results and accuracy

In this section, we discuss the accuracy of the aggregate policy function and the parameterized cross-sectional distribution. We also discuss the accuracy of our simulation approach. Tests to check the accuracy of the individual policy function are standard and these are discussed in A.3. Parameter settings of the numerical procedure, such as the order of the polynomial and the number of grid points, are given in A.2.

Parameter values

In the benchmark specification, parameter values are taken from Krusell and Smith (1998) and are reported in Tables 1 and2. The discount rate, coefficient of relative risk aversion, share of capital in GDP, and the depreciation rate take on standard values. Unemployed people are assumed to earn a fixed fraction of 15% of the wage of the employed. 19 The value of ∆ a is equal to 0.01 so that productivity in a boom, 1+∆ a , is two percent above the value of productivity in a recession, 1-∆ a . Business cycles are symmetric and the expected duration of staying in the same regime is eight quarters. The unemployment rate in a boom, u g , is equal to 4% and the unemployment rate in a recession, u b , is equal to 10%. The average unemployment duration is 2.5 quarters conditional on staying in a recession and equal to 1.5 quarters conditional on staying in a boom. These features correspond with the transition probabilities reported in Table 2.

,ε / s 0 , ε 0 1-∆ a , 0 1-∆ a , 1 1+∆ a , 0 1+∆ a , 1 1-∆ a ,

Aggregate policy function

In this section, we address the accuracy of the aggregate policy function. In Section 6.2.1, we establish the accuracy of the functional form taking the parameterization of the cross-sectional distribution as given. In Section 6.2.2, we establish whether more moments are needed as state variables. In Section 6.2.3, we describe a more demanding accuracy test by taking a multi-period perspective.

Accuracy of functional form of aggregate policy function

The aggregate policy functions,

-→ Γ e n (a, a -1 , m; ψ Γ e n ) and -→ Γ u n (a, a -1 , m; ψ Γ u n
) capture the law of motion of the end-of-period values of the three moments that are used as state variables. The approximation uses a Tensor product polynomial with at most first-order terms for a and a -1 , since a can take on only two values, and up to second-order terms for the elements of m.

Accuracy is evaluated using a grid of the aggregate state variables on which the three variables with continuous support can take on a fine range of values. In particular, { ← -m e,1 } = {35, 35.2, ..., 42.4}, { ← -m u,1 } = {33.5, 33.7, ..., 41.5}, and { --→ m u,c -1 } = {0, 0.05%, ..., 0.2%}. At each grid point, we use the values of m and the reference moments to obtain the corresponding density exactly as they are calculated in the algorithm. Whether this parameterization of the cross-sectional distribution is accurate will be discussed below. Using the parameterized density and the individual policy function, we calculate --→ m u,c , --→ m e,1 , and --→ m u,1 . These explicitly calculated values are compared with those generated by the approximations -→ Γ e (s) and -→ Γ u (s). Table 3 reports for each of the three statistics the average and maximum absolute % error across this fine set of grid points. 20 The errors for the first-order moment of the capital stock of the employed are small. The maximum error is 0.012% and the average error is 0.0059%. The moments for the unemployed are somewhat bigger, but still relatively small. In particular, the maximum error for the first-order 20

Since

--→ m u,c is a small number and already a percentage, we express the error for --→ m u,c in terms of percentage points difference an not as a percentage. moment is equal to 0.92% and the average error is equal to 0.24%. This maximum is attained when the value of ← -m e,1 takes on the highest and ← -m u,1 the lowest grid value, which is an unlikely if not impossible combination to occur. The average and maximum error for --→ m u,c are 0.2 and 0.84 percentage points (pp), respectively.21 There are two reasons why these two numbers are not problematic. First, given that both the actual and the approximation predict very low fractions of agents at the constraint, these errors are of no importance and it wouldn't make sense to spend computing time on improving the part of -→ Γ u (s) that determines --→ m u,c . Second, in Section 6.2.3, we show using a simulation that the economy doesn't get close to points in the state space where such large errors are observed. In fact, using a simulation of 1,000 periods we find an average error of 0.0076 percentage points and a maximum error of 0.071 percentage points.

--→ m e,1 5.9e-3% 1.2e-2% --→ m u,1 2.4e-1% 9.2e-1% --→ m u,c
2.0e-1pp 8.4e-1pp

Number of moments as state variables

The algorithm uses --→ m u,c -1 , ← -m e,1 , and ← -m u,1 as state variables and in this section we analyze whether additional moments should be used as state variables. That is, conditional on staying within the class of cross-sectional distributions pinned down by the reference moments does it make a difference if additional moments are used as state variables. In particular, we check whether changes in the second-order moment matter for the key set of moments the agents predict, i.e., --→ m u,c -1 , ← -m e,1 , and ← -m u,1 . To do this we calculate at each of the aggregate grid points --→ m u,c , --→ m e,1 , and --→ m u,1 in two different ways. First, when ← -m e,2 and ← -m u,2 take on its reference values, i.e., the average observed in the simulated series (conditional on the value of a). Second, when ← -m e,2 and ← -m u,2 take on the maximum values observed in the simulation, but the values of ← -m e,j , and ← -m u,j for j > 2 are still equal to the reference moments.

Table 4 reports for each of the three statistics the average and maximum absolute % change across the grid points when the variance increases. 5.2e-3%

1.1e-1% --→ m u,1 1.9e-1% 6.4e-1% --→ m u,c
1.6e-1pp 5.8e-1pp

The effect of the increase in the variance on --→ m e,1 and --→ m u,1 is small, especially considering that the increase in the variance is enormous. Again, the largest changes occur at unlikely grid points and the changes for --→ m u,c are larger. Of course, it is not surprising that an increase in the variance has an effect on the fraction of agents choosing a zero capital stock, since the increase in the variance increases the fraction of agents close to zero. Given the lack of importance of agents at the constraint, it doesn't make sense to add the second-order moment as a state variable.

Multi-period perspective

A word of caution is warranted in drawing conclusions about accuracy from the type of one-period tests performed in the last two sections. The reason is that small errors can accumulate over time if they do not average out. To investigate this issue we compare values for --→ m e,1 , --→ m u,1 , and --→ m u,c generated by two different procedures. First, we generate these statistics with our simulation procedure that explicitly integrates over the choices made by the agents in the economy. This procedure does not use our approximations -→ Γ e (a, a -1 , m) and -→ Γ u (a, a -1 , m). Second, we generate these statistics using only our approximations -→ Γ e (a, a -1 , m) and -→ Γ u (a, a -1 , m). It is important to point out that this second procedure only uses -→ Γ e (a, a -1 , m) and -→ Γ u (a, a -1 , m), and involves nothing more than basic algebra. That is, the output of our aggregate law of motion will be used as the input in the next period. 22 This comparison, thus, is truly a multi-period accuracy test.

The results for --→ m e,1 , --→ m u,1 , and --→ m u,c are plotted in Figures 1,2, and 3 respectively. The graphs make clear that our approximate aggregate laws of motion do a magnificent job of tracking the movements of --→ m e,1 and --→ m u,1 . In fact, one cannot even distinguish the moments generated by -→ Γ e (s) and -→ Γ u (s) from the corresponding --→ m e,1 generated using the approximation -→ Γ e (s) or the simulation on a continuum of agents moments generated by explicit integration over the individual policy rules. Some differences between the two procedures are visible for --→ m u,c , but our approximate aggregate laws of motion track the changes in --→ m u,c well. As mentioned above, in a sample of 1,000 observations the average and maximum absolute difference are 0.0076 and 0.071 percentage points, respectively.

The parameterized cross-sectional distribution

Parameterization of the cross-sectional distribution with a flexible functional form serves two objectives in our algorithm. First, it enables the algorithm to calculate the aggregate laws of motion, -→ Γ e (s) and -→ Γ u (s), with standard projection techniques, since with a parameterized density (i) next period's moments can be calculated on a prespecified grid and (ii) next period's moments can be calculated with quadrature instead of the less accurate Monte Carlo techniques. Second, it makes it possible to simulate the economy without cross-sectional variation, which improves the procedure to find the reference moments. --→ m u,c generated using the approximation --→ Γ u,c (s) or the simulation on a continuum of agents An accurate representation of the cross-sectional distribution may not be necessary for an accurate solution of the model. What is needed for an accurate solution of the model are accurate aggregate laws of motion, -→ Γ e (s) and -→ Γ u (s), since the agent is only interested in predicting future prices and to determine these one does not need the complete distribution, just the set of statistics that determine prices, that is,

--→ m u,c -1 , ← -- m e,1
, and ← -m u,1 . In this section, we take on the more demanding test to check whether the crosssectional parametrization, P (k; ρ e ) and P (k; ρ u ), are accurate. This is a more demanding test for the following reasons. First, it requires that all N M moments used to pin down the distribution are accurately calculated instead of just the N M moments that are used as state variables. More importantly, because the shape of the cross-sectional distribution is endogenous and time-varying, the functional form used must be flexible enough to capture the unknown and changing shapes.

To check the accuracy of our simulation procedure and, thus, the accuracy of our parameterized densities, we do the following. We start in Section 6.3.1 with a comparison between the simulated time path of moments generated by our parameterized densities, with those generated by a standard simulation using N N agents. The alternative simulation is, of course, subject to sampling variation, but the advantage of the standard simulation procedure is that there is no functional restriction on the cross-sectional distribution at all.

A second accuracy test consists of checking whether the results settle down if N M increases. This is done in Section 6.3.2. The last accuracy test checks whether moments of order higher than N M are calculated precisely.

Comparison between simulation procedures

In this section, we compare the moments generated by our new simulation approach with those generated by the standard simulation procedure. The same individual policy function is used under the two approaches. Note that both approaches are subject to sampling variation of the aggregate shock, a. That seems unavoidable. We plot time paths of generated moments when the value of a alternates between 1 -∆ a and 1 + ∆ a for long periods of time so that the behavior of the economy during a transition between regimes becomes clear. A set of initial observations is discarded so that effects of the initial distribution are no longer present. 23 Figure 4 reports the evolution of the end-of-period first moment of the employed --→ m e,1 . Note that the group of the employed is at least 9,000. Although the sampling variation is still visible in the graph, it is clear that it is small relative to the observed 23 Details about the simulation procedure are given in Table 9. --→ m e,1 generated using a finite and a continuum of agents changes in the calculated moment. Figures 5 and6 show the evolution of the-end-of period first-order moment of the unemployed --→ m u,1 and the fraction of unemployed agents at the borrowing constrained, --→ m u,c . The number of observations in this group is much smaller and the sampling variation is substantial. Sampling variation is especially severe for the fraction of agents at the constraint. In this economy, the number of agents at the constraint is small, however, and the unemployed do not own much capital. Any inaccuracy due to the high sampling variation is, thus, likely to be inconsequential for the properties of any aggregate series. But this is unlikely to carry over to models in which a minority are more influential like, for example, in models with entrepreneurs. Moreover, if one is explicitly interested in the behavior of the group of unemployed then Monte Carlo simulations provide very noisy answers. 24 A good way to document the higher accuracy of the new simulation procedure is to look at the transition from the bad to the good regime. The new simulation procedure clearly shows a sharp increase in the fraction of agents at the constraint when the economy enters the good regime. In contrast, for the standard simulation --→ m u,1 generated using a finite and a continuum of agents 0,00% 0,25% 0,50% 0,75% --→ m u,c generated using a finite and a continuul of agents procedure this increase is not always present and there are many other spikes. The true solution of the model should exhibit such a spike. In this economy, employed agents become unemployed every period. When the economy moves from the highunemployment to the low-unemployment regime, then the flow out of employment into unemployment drops sharply. This means that after the regime change, an unemployed agent is much less likely to have been employed in the recent past. Consequently, after a change to the low-unemployment regime, a larger fraction of unemployment agents will have a zero capital stock.

Increasing N M

In this paper, we use N M moments to characterize the cross-sectional distribution and N M of these are state variables. In Section 6.2.2, we analyzed whether additional moments are necessary as state variable. Here we address the question whether increasing N M makes a difference.

To check the importance of N M , we simulate an economy using different values of N M to parameterize the cross section. We check when the results settle down. The idea of the test is made clear in Figures 7 and8 that plot the second and sixth-order moment of the distribution for the unemployed for different values of N M . The figures make clear that N M has to be sufficiently high. When we increase N M from 2 to 4 then the generated moments change enormously. A further increase from 4 to 5 still causes some changes, but when we increase N M from 5 to 6 then the changes are very minor.

Table 5 reports the results for all moments. It corroborates the results from the figures. Using a value of N M equal to 2 is clearly way too low to generate an accurate set of moments. An increase in N M from 5 to 6, however, only causes minor changes in the generated moments. For example, the average change for ← -m e,2 is only 0.25% and the average change for ← -m u,2 is equally small. The average (maximum) change for ← -m e,6 and ← -m u,6 are equal to 2.33 (3.67) and 2.35 (3.72), respectively. These are, thus, somewhat, higher, but as made clear in the graph, these errors are low relative to the observed variation in the moments.

The shape of the distribution

There are two aspects to our approximation of the cross-sectional distribution. First, the class of functions used and second the value of N M . The coefficients are chosen so that the first N M moments are correct. But higher-order moments are implied by the values of the first N M moments and the class of functions chosen. For example, when one uses a normal distribution, then one can get any mean and variance, but skewness and kurtosis are pinned down. So for any finite value of N M , the class of approximating polynomials used impose certain restrictions on the function form. Here we check those restrictions along a simulated time path by comparing the j th -order moments for j > N M implied by the parameterized cross-section with those calculated by integration of the individual's policy function.

M = 2 N M = 4 N M = 5 N M = 2 N M = 4 N M = 5 ← -- m e,
In particular, we do the following. Draw a long time series for the aggregate productivity shock, a. Let ρ e 1 and ρ w 1 be the parameters of the cross-sectional distribution in the first period and let ← -m e,c 1 and ← -m u,c 1 be the fraction of employed and unemployed agents with zero capital holdings at the beginning of the period. With this information, we calculate the end-of-period values of the first N M moments.25

--→ m w,1 1 = ∞ Z 0 k(ε w , k, s) P (k; ρ w 1 )dk, w ∈ {e, u}. (16)
Table 6: Implied and actual higher-order moments;

N M = 2 N M = 2 Employed Unemployed Error (%) Average Max Average Max ¯← --- = w,3 - ← --- m w,3 ← --- m w,3
5.6% 9.1% 7.0% 15.3%

¯← --- = w,4 - ← --- m w,4 ← --- m w,4 2.8E-1% 3.2E-1% 2.1E-1% 6.4E-1% ¯← --- = w,5 - ← --- m w,5 ← --- m w,5
4.7% 7.5% 5.6% 11.7%

¯← --- = w,6 - ← --- m w,6 ← --- m w,6 5.7E-1% 7.5E-1% 9.3E-1% 2.6% ¯← --- = w,7 - ← --- m w,7 ← --- m w,7
3.7% 6.1% 5.3% 10.3%

¯← --- = w,8 - ← --- m w,8 ← --- m w,8 6.5E-1 1.2% 2.3% 5.7% ¯← --- = w,9 - ← --- m w,9 ← --- m w,9
2.9% 5.0% 5.6% 10.3%

¯← --- = w,10 - ← ---- m w,10 ← ---- m w,10

Concluding comments

In this paper, we have developed a new algorithm to solve models with heterogeneous agents and aggregate uncertainty. Given the rising interest in such models, it is important to have an alternative to the algorithm that is currently predominantly used, especially since it relies on simulation procedures that-given what we know from other numerical problems-are likely to be unsuitable for some models. We have shown how to implement the algorithm for the model developed in Krusell and Smith (1998). This is a relatively simple model and not surprisingly the properties of the model calculated with our solution procedure are very similar to the ones obtained with an algorithm based on a Monte Carlo simulation procedure. We leave it for future research to see whether the same is true for more complex models. We also plan to investigate the importance of replacing the standard Monte Carlo simulation procedure with our simulation procedure in algorithms that use a simulation procedure to calculate aggregate laws of motion. This improvement might turn out to be key in models in which a large fraction of aggregate wealth is held by a small group of the population such as entrepreneurs.

But even if standard Monte Carlo simulation procedures generate accurate aggregate laws of motion, they clearly do not describe well the behavior of agents that are few in number. If one wants to use a model with heterogeneous agents to study the behavior of sub groups, then our more accurate solution should be considered.

A Appendix

A.1 Details on Transition equations

This appendix describes how the change in employment status that occurs at the beginning of each period affects the moments of the cross-sectional distribution. Although, we use centralized moments as state variables, we actually do not use centralized moments here. It is easier to first do the transformation for non-centralized moments and then calculate the centralized moments.

From beginning to end-of-period. Let g a,w be the mass of agents with employment status w when the economy is in regime a. At the beginning of the period we have to following groups of agents:

1. Unemployed with k = 0, whose mass is equal to ← -m u,c g a,u

2. Unemployed with k > 0, whose mass is equal to

³ 1 - ← -- m u,c
´ga,u 3. Employed with k = 0, whose mass is equal to ← -m e,c g a,e

4. Employed with k > 0, whose mass is equal to

³ 1 - ← -- m e,c
´ga,e

Agents in group #1 choose k 0 = 0, while agents in group #2 either set k 0 = 0 or k 0 > 0. Let the fraction of agents that set k 0 = 0 be equal to ζ k 0 =0 u,k>0 . Thus, the fraction of unemployed agents that set k 0 = 0, is equal to

--→ m u,c = ← -- m u,c + ζ k 0 =0 u,k>0 (1 - ← -- m u,c
).

The i th moment of the capital stock chosen by agents in group #2 is equal to

µ k 0 ≥0,i u,k>0 = ζ k 0 =0 u,k>0 × 0 i + (1 -ζ k 0 =0 u,k>0 ) × --→ m u,i .
Thus,

--→ m u,i = µ k 0 ≥0,i u,k>0 (1 -ζ k 0 =0 u,k>0 )
where µ k 0 ≥0,i u,k>0 = R +∞ 0 k(0, k, s) i P (k, ρ u )dk. --→ m e,c = 0, since employed agents never choose a zero capital stock. To calculate --→ m e,i , we need the (mean of the) capital stock chosen by those in group #3, i.e., µ k 0 >0,i e,k=0 = k(1, 0, s) i , and the mean of the capital stock chosen by those in group 4, i.e., From end-of-period to beginning of next period. At the end of the current period, we have to following groups of agents:

1. Unemployed with k 0 = 0, whose mass is equal to --→ m u,c g a,u

2. Unemployed with k 0 > 0, whose mass is equal to (1 ---→ m u,c )g a,u

3. Employed with k 0 > 0, whose mass is equal to g a,e g ww 0 ,aa 0 stands the mass of agents with employment status w that have employment status w 0 in the next period period, conditional on the values of a and a 0 . For all combinations of a and a 0 we have, g uu,aa 0 + g eu,aa 0 + g ue,aa 0 + g ee,aa 0 = 1.

The number of agents with k 0 = 0 as a fraction of all agents that are unemployed in the following period is equal to ← -m u,c0 = g uu,aa 0 g uu,aa 0 + g eu,aa 0 --→ m u,c .

and the fraction of all employed agents at the constraint is equal to ← -m e,c0 = g ue,aa 0 g ue,aa 0 + g ee,aa 0 --→ m u,c

Next period's i th -order moments of the distributions with strictly-positive capital holdings are equal to ← -m u,i0 = g uu,aa 0 (1 ---→ m u,c ) --→ m u,i + g eu,aa 0 --→ m e,i g uu,aa 0 (1 ---→ m u,c ) + g eu,aa 0 and ← -m e,i0 = g ue,aa 0 (1 ---→ m u,c ) --→ m u,i + g ee,aa 0 --→ m e,i g ue,aa 0 (1 ---→ m u,c ) + g ee,aa 0 A.2 Parameters of numerical procedure Tables 8 and9 report the parameter settings used in implementing the numerical procedure.

  17

4. 5

 5 Procedure to find cross-sectional distribution

2 for

 2 w ∈ {e, u} and j ∈ {c, 1, • • • , N f M }. Recall that this simply takes care of the effects of changes in the employment status on the cross-sectional distribution. Details are given in Appendix A.1.

Table 4 :

 4 Effect of increase in variance under reference distribution

  Figure 1:--→ m e,1 generated using the approximation -→ Γ e (s) or the simulation on a continuum of agents

Figure 2 :

 2 Figure 2:--→ m u,1 generated using the approximation -→ Γ u (s) or the simulation on a continuum of agents

  Figure 4:--→ m e,1 generated using a finite and a continuum of agents

Figure 5 :

 5 Figure 5:--→ m u,1 generated using a finite and a continuum of agents

Figure 6 :

 6 Figure 6:--→ m u,c generated using a finite and a continuul of agents

  Figure 7:--→ m u,2 generated using a continuum of agents with different values of N M

  , k, s) i P (k, ρ e )dk.Weighting the two means gives--→ m e,i = ← -m e,c µ k 0 >0,i e,k=0 + (1 -← -m e,c )µ k 0 >0,i e,k>0← -m e,c + (1 -← -m e,c )

M 1 = 5 ,

 15 k = 50, M ---→ m u,c -M ←-- m e,1 = 5, M ← --- m u,1 = 5 M ε = M a = M a -1 = 2

Table 1 :

 1

			Benchmark calibration
	Parameters	β	γ	α	δ	µ	∆ a
	Values	0.99 1 0.36 0.025 0.15 0.01

Table 2 :

 2 Transition probabilities s

Table 3 :

 3 Accuracy of functional form of aggregate policy function Moment approximation error average maximum

Table 5 :

 5 The effect of increasing N M

	Moment	% change relative to N M = 6
		average	maximum
	N

Table 7 :

 7 Implied and actual higher-order moments; N M = 6 N M = 6

	Employed	Unemployed
	Error (%) ¯← ---= w,7 -← ---← m w,7 ---m w,7 ¯← ---= w,8 -← ← ---m w,8 ---m w,8 ¯← ---= w,9 -← ← ---m w,9 ---m w,9 ¯← ---← Average 2.8E-2% 7.3E-1% 1.0E-1% 2.2E-1% Max Average Max 4.3E-2% 1.0E-1% 1.8E-1% 4.3E-1% 9.3E-2% 2.3E-1% 3.8E-1% 8.8E-1% ← ----m w,10 = w,10 -----m w,10

Table 8 :

 8 Construction of the grid

	State variables	k, ε, a -1 , a,	--→ m u,c -1 ,	← --m e,1 ,	← --m u,1
	Number of grid points

Descriptions of this algorithm can be found in Den[START_REF] Haan | Heterogeneity, aggregate uncertainty and the shortterm interest rate[END_REF],[START_REF] Krusell | Income and wealth heterogeneity in the macroeconomy[END_REF], and Rios-Rull (1997).

See Christiano and Fisher (2000) and Judd (1998) for a detailed discussion of the disadvantages of simulation procedures.

In fact,[START_REF] Krueger | Computing equilibrium in OLG models with stochastic production[END_REF] analyze an overlapping generations modle and argue that models in which the cross-sectional distribution can be efficiently described by just the mean are special cases and that in general higher-dimensional characterizations of the state are needed.

See, for example,Christiano and Fischer (2000).

In particular, a crucial property being used is that conditional on realizations of the aggregate shock and this period's cross-sectional distribution, next period's cross-sectional distribution is known with certainty.

See Judd (1998).

Note that in the classic regression problem the standard errors, σ 2 (X 0 X) -1 , are also lower when the x-values are more spread out.

[START_REF] Krueger | Computing equilibrium in OLG models with stochastic production[END_REF] develop and solve an OLG model with a relatively large but finite number of agents and show how Smolyak's algorithm can be used to solve models with 20 to 30 agents.

[START_REF] Miao | Competitive equilibria of economies with a continuum of consumers and aggregate shocks[END_REF] shows the existence of a recursive equilibrium, but also uses expected payoffs as state variables. It is not clear whether a recursive equilibrium exists when the smaller set of state variables is used. For a numerical solution this is less important in the sense that approximation typically entails not using all information.

Details are given in A.1.

This is the approach used in[START_REF] Veracierto | Plant-level irreversible investment and equilibrium business cycles[END_REF]. He solves a model with irreversible investment in which the cross-sectional distribution matters because the investment decision is of the (S,s) variety. Instead of keeping track of the cross-sectional distribution, he keeps track of a history of lower and upper threshold levels.

The numerical cost is fairly low, since a can take on only two values.

We have suppressed the i superscript for notational convenience.

For details see A.1.

We use the BGFS-method as implemented by Peter Spellucci. The upperbound of the Simpson quadrature procedure is equal to 99.

This is the only change relative to Krusell and Smith (1998) who set µ = 0. This has little effect on the properties of the model but avoids the possibility of agents having zero consumption.

This maximum difference for --→ m u,c is also attained at the unlikely combination of a very high value for ← -m e,1 and very low value for ←-m u,1 . At this grid point, the value from our approximation is equal to 0.152% and the recalculated value is 0.996%.

After adjusting, of course, for the employment-status flows.

It may be possible to oversample the unemployed, but this is not trivial in a general equilibrium framework in which the employment status changes. The reason is that oversampling has to be consistent with the appropriate cross-sectional distribution.

By indicating in bold the variable that we are integrating over, we make clear that we are integrating over the argument of the policy function not the outcome of the policy function.

Note that the accuracy measure is actually defined for beginning-of-period moments, but this is simply a transformation of end-of-period values taking into consideration the change in the employment status.

1 ¸j P (k; ρ w 1 )dk, 1 < j ≤ N M , w ∈ {e, u}.

In exactly the same way, we calculate higher-order moments. That is,

1 ¸j P (k; ρ w 1 )dk, j > N M , w ∈ {e, u}.

Now we check whether these higher-order moments (j > N M ) are similar to the moments implied by the parameterized cross-sectional distribution. To do this, we use the first N M end-of-period moments to calculate the coefficients of the corresponding approximating density, -→ ρ e 1 and -→ ρ u 1 . Next, we calculate higher-order moments implied by this parameterized cross-sectional density. That is,

If the shape of the cross-sectional distribution is not too restrictive, then the implied higher-order moments correspond to the explicitly calculated higher-order moments. We perform this exercise using values for N M equal to 2 and 6 and then calculate the average and maximum error observed along the simulation of 2,000 observations with the error term defined as follows: 26 ¯← --

Tables 6 and7 report the errors for N M = 2 and N M = 6, respectively. When N M is equal to 2, then observed error terms are large for the odd-numbered moments. That is, the shape of the distribution implied by our class of approximating functions doesn't capture the correct shape of the distribution when a second-order approximation is used. The results are much better when N M is equal to 6. Now we observe much smaller errors. The largest errors are for the 10-th order moment of the capital stock of the unemployed. For this moment, the maximum error is 1.3%, which is not excessively high.

A.3 Accuracy individual policy function

In this section, we discuss the accuracy of the numerical solution for the individual policy function, taking as given the transition laws for the cross-sectional distribution. We check the accuracy of our policy function using a grid of capital holdings that is much more dense than the one used to obtain the numerical solution and the grid points for the aggregate state variables. In particular, we use k={0, 0.1, 0.2, • • • , 99}. At each grid point, we first calculate the consumption level implied by the numerical approximation, i.e., c app , and next period's capital.

Next, we calculate the conditional expectation, not by using the approximation, but by explicit (numerical) integration of βc 0-γ (r 0 + 1δ) over the possible realizations of a 0 and ε 0 . The numerical approximation for the policy function is used to evaluate next period's consumption values for the different realizations of the idiosyncratic and aggregate random variable. Let c be the current-period value of consumption implied by the explicitly calculated conditional expectation on that grid point. Accuracy is measured using |cc app |/c app . We calculate the average and the maximum across these approximation errors. Table 10 reports the average across aggregate grid points and the maximum. Even the maximum errors are small for most grid points. The average of the maximum errors for the unemployed is equal to 0.67%. For the employed errors are substantially smaller. The maximum percentage error for the unemployed is equal to 5.9%. Such a maximum error is high, but it occurs when consumption takes on a low value, namely 0.690 based on our approximation and 0.651 based on the explicitly calculated conditional expectation. These values are low relative to the average consumption value of the unemployed, which is equal to 2.77.