How to solve dynamic stochastic models computing expectations just once
Abstract
EconStor is a publication server for scholarly economic literature, provided as a non-commercial public service by the ZBW.
Full text
Judd, Kenneth L.; Maliar, Lilia; Maliar, Serguei; Tsener, Inna Article How to solve dynamic stochastic models computing expectations just once Quantitative Economics Provided in Cooperation with: The Econometric Society Suggested Citation: Judd, Kenneth L.; Maliar, Lilia; Maliar, Serguei; Tsener, Inna (2017) : How to solve dynamic stochastic models computing expectations just once, Quantitative Economics, ISSN 1759-7331, The Econometric Society, New Haven, CT, Vol. 8, Iss. 3, pp. 851-893, https://doi.org/10.3982/QE329 This Version is available at: https://hdl.handle.net/10419/195556 Standard-Nutzungsbedingungen: Die Dokumente auf EconStor dürfen zu eigenen wissenschaftlichen Zwecken und zum Privatgebrauch gespeichert und kopiert werden. Sie dürfen die Dokumente nicht für öffentliche oder kommerzielle Zwecke vervielfältigen, öffentlich ausstellen, öffentlich zugänglich machen, vertreiben oder anderweitig nutzen. Sofern die Verfasser die Dokumente unter Open-Content-Lizenzen (insbesondere CC-Lizenzen) zur Verfügung gestellt haben sollten, gelten abweichend von diesen Nutzungsbedingungen die in der dort genannten Lizenz gewährten Nutzungsrechte. Terms of use: Documents in EconStor may be saved and copied for your personal and scholarly purposes. You are not to copy documents for public or commercial purposes, to exhibit the documents publicly, to make them publicly available on the internet, or to distribute or otherwise use the documents in public. If the documents have been made available under an Open Content Licence (especially Creative Commons Licences), you may exercise further usage rights as specified in the indicated licence. https://creativecommons.org/licenses/by-nc/4.0/
Quantitative Economics 8 (2017), 851–893 1759-7331/20170851 How to solve dynamic stochastic models computing expectations just once Kenneth L. Judd Hoover Institution, Stanford University Lilia Maliar Department of Economics, Stanford University Serguei Maliar Department of Economics, Santa Clara University Inna Tsener Department of Applied Economics, University of the Balearic Islands We introduce a computational technique—precomputation of integrals—that makes it possible to construct conditional expectation functions in dynamic stochastic models in the initial stage of a solution procedure. This technique is very general: it works for a broad class of approximating functions, including piecewise polynomials; it can be applied to both Bellman and Euler equations; and it is compatible with both continuous-state and discrete-state shocks. In the case of normally distributed shocks, the integrals can be constructed in a closed form. After the integrals are precomputed, we can solve stochastic models as if they were deterministic. We illustrate this technique using oneand multi-agent growth models with continuous-state shocks (and up to 60 state variables), as well as Aiyagari’s (1994) model with discrete-state shocks. Precomputation of integrals saves programming efforts, reduces computational burden, and increases the accuracy of solutions. It is of special value in computationally intense applications. MATLAB codes are provided. Keywords. Dynamic model, precomputation, numerical integration, dynamic programming, value function iteration, Bellman equation, Euler equation, envelope condition method, endogenous grid method, Aiyagari model. JEL classification. C61, C63, C68. Kenneth L. Judd: [email protected] Lilia Maliar: [email protected] Serguei Maliar: [email protected] Inna Tsener: [email protected] We are indebted to the editor and three anonymous referees for many thoughtful comments and suggestions. Errors are ours. Support from the Hoover Institution and Department of Economics at Stanford University, University of Alicante, University of the Balearic Islands, Santa Clara University, and the MINECO/FEDER Grant ECO2015-70540-P is gratefully acknowledged. Copyright ©2017 The Authors. Quantitative Economics. The Econometric Society. Licensed under the Creative Commons Attribution-NonCommercial License 4.0. Available at http://www.qeconomics.org. DOI: 10.3982/QE329
852 Judd, Maliar, Maliar, and Tsener Quantitative Economics 8 (2017) 1. Introduction Existing global methods for solving dynamic stochastic models compute conditional expectation functions in their iterative cycles.1Recomputing expectation functions in each iteration is costly and the cost grows rapidly (i) when the number of random variables increases (because the dimensionality of integrals increases), (ii) when more accurate integration methods are used (because the number of integration nodes increases), and (iii) when models become more complex (because numerical solvers are used more intensively, and this involves additional evaluations of integrals). In this paper, we introduce a computational technique that makes it possible to construct conditional expectation functions in the initial stage of the solution procedure; we refer to this technique as precomputation of expectation functions or precomputation of integrals. The idea is simple and can be seen from the following example: Assume that the value function of a stylized growth model is approximated using a linear polynomial function V(kz)≈b0+b1k+b2z,whereb0,b1,andb2are polynomial coefficients, kis capital, and zis productivity that follows a first-order autoregressive process, z=zρexp(ε),withρ∈(−11)and εbeing a random shock. The key step of our precomputation analysis is to notice that expectation of an ordinary polynomial function can be derived in a closed form as E[b0+b1k+b2z]=b0+b1k+b2zρI,where I≡E[exp(ε)]. Given a distribution function of ε, the integral Ican be constructed either analytically or numerically. In particular, it can be constructed analytically in the case of normally distributed shocks, used by a vast majority of economic models, namely, for ε∼N(0σ2),wehaveI=exp(σ2 2). Importantly, we need to construct Ijust once, in the stage of initialization. Within the main iterative cycle on bs, expectation functions can be evaluated by using a closed-form expression that contains no random variables, that is, E[V(k z)]≈b0+b1k+b2zρI. In effect, precomputation of integrals allows us to solve a stochastic problem as if it was a deterministic problem. In our example, the advantage of using precomputation of integrals is twofold: First, we attain higher accuracy of numerical solutions because we construct integrals exactly, whereas the related literature constructs integrals approximately, by using some numerical integration method (e.g., Monte Carlo, quasi-Monte Carlo, quadrature, monomials).2Second, we are able to reduce the cost of constructing numerical solutions because in the main iterative cycle we evaluate future value function in just one composite future state, whereas the existing global solution methods approximate future value function as a weighted average across a possibly large number of future states. Of course, our example is very special. However, it turns out that integrals can be precomputed in a variety of other contexts: First, precomputation of integrals can be implemented for any set of equations that contain conditional expectation functions, including the Bellman and Euler equations. Second, integrals can be precomputed not only 1For reviews of methods for solving dynamic economic models, see Taylor and Uhlig (1990), Rust (1996, 2008), Gaspar and Judd (1997), Judd (1998), Marimon and Scott (1999), Santos (1999), Christiano and Fisher (2000), Miranda and Fackler (2002), Aruoba, Fernández-Villaverde, and Rubio-Ramírez (2006), Stachursky (2009), Den Haan (2010), Kollmann, Maliar, Malin, and Pichler (2011), and Maliar and Maliar (2014). 2See Judd, Maliar, and Maliar (2017) for a discussion of alternative accuracy measures of numerical solutions.
Quantitative Economics 8 (2017) How to solve dynamic stochastic models 853 for ordinary polynomial functions but also for any other approximating families whose bases are separable in endogenous and exogenous state variables, including orthogonal polynomial families such as Chebyshev, Smolyak, Hermite, and piecewise polynomial functions, as well as many nonpolynomial families.3Third, precomputation of integrals can be combined with other computational techniques used by existing global solution methods, including a variety of solution domains, integration rules, fitting methods, and iterative schemes for finding unknown parameters of approximating functions. Fourth, precomputation of expectation functions is also possible for models with a discrete set of shocks and a discrete set of controls. Fifth, integrals can be computed analytically not only for univariate, but also for multivariate normally distributed shocks, including the case when shocks are correlated. Finally, in those cases when integrals cannot be constructed analytically, we can construct them numerically using very accurate methods since they should be constructed just once (i.e., this is a one time fixed cost). We emphasize that precomputation of integrals is not a new solution method but an analytical and numerical manipulation of the model’s equations that simplifies the construction of conditional expectation functions. Moreover, precomputation of integrals is not related to any specific solution method: any numerical solution method that can be applied to solve the original model’s equations can also be applied to solve the model’s equations, obtained after precomputing the integrals. In particular, we show that the precomputation technique can enhance the performance of five existing solution methods: conventional value function iteration, the endogenous grid method of Carroll (2006), the envelope condition method of Maliar and Maliar (2013), and two versions of the Euler equation methods. In the context of the stylized one-agent neoclassical growth model, we show that precomputation of integrals reduces the running time of the value-iterative and Euler equation methods up to three and five times, respectively, depending on the degree of the polynomial approximation and the specific solution method considered. Furthermore, we show that precomputation of integrals can be implemented in more complex models such as a growth model with elastic labor supply. MATLAB codes are provided in a supplementary file on the journal website, http://qeconomics.org/supp/329/code_and_data.zip. It is noteworthy that precomputation of integrals leads to much larger gains in terms of accuracy and speed in models with multiple shocks than in models with one shock. We solve a multicountry growth model with up to 30 heterogeneous countries (60 state variables, including 30 correlated shocks) using a generalized stochastic simulation algorithm (GSSA) in line with Judd, Maliar, and Maliar (2011). In the latter paper, integrals are approximated numerically using deterministic integration methods such as a Gauss Hermite product rule and two monomial integration methods; for these methods, the number of integration nodes grows with the number of shocks exponentially, quadratically, and linearly, respectively. In contrast, precomputation of integrals always means just one integration node. We find that precomputation of integrals can reduce the running time by many orders of magnitude with multivariate shocks compared to numerical approximations 3See Judd (1998) for a survey of polynomial approximating functions. Also, see Krueger and Kubler (2004), and Judd, Maliar, Maliar, and Valero (2014) for a discussion of Smolyak approximating functions.
854 Judd, Maliar, Maliar, and Tsener Quantitative Economics 8 (2017) of integrals. In particular, in the model with 30 countries, only the solution method using precomputation of integrals was able to deliver accurate second-order numerical solutions, while similar solution methods, which use numerical approximations of integrals, were too expensive. Finally, we show that expectation functions can be also precomputed in models with a discrete set of shocks. As an illustration, we use Aiyagari’s (1994) model in which decision functions are parameterized by piecewise linear polynomial functions. Here, the reduction in running time depends critically on the number of states in the associated Markov chain for exogenous shocks. In our baseline case of a seven-state Markov chain, the running time is reduced by around 50%, but far larger reductions in cost are possible when the number of states increases. Our analysis suggests that the gains from precomputation of expectation functions will be especially large in models with multiple shocks in which Markov chains are characterized by a very large number of states. In the case of value-iterative methods, precomputation of integrals is always beneficial: it does not affect the way in which value function iteration is implemented; it just makes it faster. In the case of the Euler equation methods, precomputation has benefits but may also have costs. To precompute expectation functions, we must reparameterize the Euler equation in a particular way and we must construct a policy function for a specific variable that is the integrand of the expectation function in the Euler equation.4In some applications, it could happen that this new policy function is harder to approximate accurately than conventional policy functions for capital or consumption. We observe this effect in the context of Aiyagari’s (1994) model: an algorithm iterating on the reparameterized Euler equation delivers slightly less accurate solutions than a similar algorithm iterating on the conventional Euler equation because the integrand has a spike in the area of the kink. However, such an accuracy reduction is very minor, even in this special case. The rest of the paper is as follows. In Section 2, we introduce the technique of precomputation of integrals in the context of an optimal growth model with one continuous-state shock. In Section 3, we show the precomputation results for models with multivariate continuous-state shocks. In Section 4, we show how to precompute expectation functions in models with discrete-state shocks. In Section 5,wedescribe possible generalizations of the precomputation technique, and we discuss some of its limitations. Finally, in Section 6, we conclude. Appendices are available in a supplementary file on the journal website, http://qeconomics.org/supp/329/supplement.pdf. 2. Univariate continuous-state expectations We show the technique for precomputing univariate (one-dimensional) continuousstate integrals in the context of the standard one-agent neoclassical stochastic growth model. 4It turns out that the reparameterization of the Euler equation that is necessary for precomputing integrals leads to the same system of equations that does the version of the envelope condition method that iterates on the derivative of the value function (ECM-DVF); see Maliar and Maliar (2013) and Arellano, Maliar, Maliar, and Tsyrennikov (2016). Thus, our analysis suggests that the ECM method produces systems of equations that are suitable for precomputation of integrals by construction.
Quantitative Economics 8 (2017) How to solve dynamic stochastic models 855 2.1 Neoclassical stochastic growth model The representative agent solves max {kt+1ct}∞ t=0 E0 ∞ t=0 βtu(ct)(1) s.t. ct+kt+1=(1−δ)kt+ztf(kt) (2) lnzt+1=ρlnzt+εt+1(3) where (k0z0)is given, Etis an operator of conditional expectation, ct,kt,andztare consumption, capital, and productivity level, respectively, β∈(01),δ∈(01],and ρ∈(−11)are parameters, εt+1∼N(0σ2)is a productivity shock, and uand fare the utility and production functions, respectively, both of which are strictly increasing, continuously differentiable, and concave. Bellman equation We can characterize a solution to the model (1)–(3)byusingadynamic programming approach. The Bellman equation is V(kz)=max kc u(c) +βEVkz (4) s.t. k=(1−δ)k +zf (k) −c (5) lnz=ρln z+ε(6) where ε∼N(0σ2)and E[V(k z)]≡E[V(k z)|kz]is expectation conditional on state (kz); here and later in the text, the primes on variables denote their next-period values. We solve for value function V(kz)that satisfies (4)–(6). Euler equation We can also characterize a solution to (1)–(3) by a set of the first-order conditions. The Euler equation is u(c) =βEuc1−δ+zfk(7) We solve for the policy functions such as c=C(kz) and k=K(kz) that satisfy (5), (6), and (7). Numerical approximation of expectation functions (integrals) To implement iteration on Bellman and Euler equations, we must construct and evaluate expectation functions, E[V(k z)]and E[u(c)(1−δ+zf(k))], respectively. The related literature approximates expectation functions using numerical integration methods, such as Monte Carlo, quasi-Monte Carlo, quadrature, and monomial. All such methods approximate the value of an integral of a given function Gby a weighted sum of the integrand function evaluated in a finite number of nodes, EGε=+∞ −∞ Gεωεdε≈ J j=1 wjGε j(8)
856 Judd, Maliar, Maliar, and Tsener Quantitative Economics 8 (2017) where ω(·)is a probability density function of ε,and{ε j}and {wj}are the integration nodes and weights, respectively, j=1J. However, integration methods differ in the number and placement of the integration nodes and in the choice of integration weights. Typically, there is a trade-off between the accuracy and cost: integration formulas with more nodes (and thus, with a higher evaluation cost) lead to more accurate approximations. Importantly, all the existing solution methods recompute expectation functions using (8) each time when they are evaluated in the solution procedure, which involves a substantial cost. We will show that this cost can be significantly reduced, while attaining the highest possible accuracy in approximation of integrals. 2.2 Precomputation of univariate integrals under ordinary polynomial approximations In this section, we show that for a class of ordinary polynomial functions, conditional expectation functions can be constructed prior to solving the model, that is, precomputed. Our precomputation analysis builds on the following fact: If a state-contingent function is parameterized by an ordinary polynomial function, then its conditional expectation function can be analytically characterized for any vector of the polynomial coefficients and any state of the world. (We had shown this result for a linear polynomial function in the Introduction, and we now generalize this result for polynomials of higher degrees.) Let a function P(k z) be approximated with a complete ordinary polynomial function, that is, P(kz;b) =b0+b1k+b2z+b3k2+b4kz +b5z2+···+bnzL(9) where b≡(b0b1bn)∈Rn+1is a vector of polynomial coefficients and Lis a degree of polynomial. Taking into account that k=K(kz) is known at present and that future productivity depends on a random draw z=zρexp(ε), we can represent conditional expectation of P(kz;b) as EPkz;b =Eb0+b1k+b2zρexpε+b3k2+b4kzρexpε+···+bnzLρ expLε =b0+b1k+b2zρEexpε+b3k2+b4kzρEexpε+··· +bnzLρEexpLε =b0I0+b1I1k+b2I2zρ+b3I3k2+b4I4kzρ+···+bnInzLρ =b 0+b 1k+b 2zρ+b 3k2+b 4kzρ+···+b nzLρ ≡Pkzρ;b (10) where b≡(b 0b 1b n)∈Rn+1. The coefficients b iand bi,fori=01n, are related by b i=biIi(11)
Quantitative Economics 8 (2017) How to solve dynamic stochastic models 857 where Ii=Eexpliε(12) with libeing a power on zin the ith monomial term. For example, zdoes not enter in the monomial terms i=0136,sol0=l1=l3=l6=0and the corresponding Is are I0=I1=I3=I6=1;zdoes enter linearly in monomial terms i=247,sothatwe have l2=l4=l7=1and I2=I4=I7=E[exp(ε)]; similarly, we obtain l5=l8=2and I5=I8=E[exp(2ε)],...,andIn=E[exp(Lε)]. The integrals I0Independ only on the properties of exogenous shock εand can be computed up-front without specifying the values of the coefficients b(which are unknown before the model is solved), that is, the integrals can be precomputed. Once Is are constructed, evaluation of conditional expectation becomes trivial in the iterative procedure. Namely, conditional expectation of a polynomial function is given by the same polynomial function but evaluated at a different coefficients’ vector, that is, E[P(kz;b)]=P(kzρ;b), where a relation between band bis determined by (11) and (12). In particular, for the case of normally distributed innovations ε∼N(0σ2),the integrals (12) can be computed exactly (in a closed form); we do this in Section 2.4.For those probability distributions for which analytical characterizations are not possible, integrals can be precomputed numerically very accurately, since we need to construct them just once at the beginning of the iterative cycle (i.e., this is a one time fixed cost). Remark 1. The precomputation result (11)and(12) also holds for piecewise polynomial approximations. For example, let us parameterize P(kz) on a rectangular grid [k1kM]×[z1zN]by using a collection of piecewise linear polynomial functions. Then, in each local area [kmkm+1]×[znzn+1], we have a local polynomial function P[kmkm+1]×[znzn+1](kz) =β(mn) 0+β(mn) 1k+β(mn) 2z where n∈{1N}and m∈{1M}. The expectation function in each local area is given by EP[kmkm+1]×[znzn+1]kz=β(mn) 0+β(mn) 1k+β(mn) 2zρEexpε i.e., (11)and(12) hold in each local area. The precomputation analysis can be generalized to higher-order piecewise polynomial approximations in a similar way. 2.3 Characterizing the solution under precomputation of integrals In this section, we show how to precompute expectations in the Bellman and Euler equations. 2.3.1 Bellman equation with precomputation of integrals Precomputation of integrals is straightforward in the Bellman equation. We parameterize the true value function V(kz)with a flexible functional form V (kz;b) given by polynomial (9). Using (10),
858 Judd, Maliar, Maliar, and Tsener Quantitative Economics 8 (2017) we rewrite the Bellman equation (4)–(6)as V (kz;b) =max kc u(c) +β Vkzρ;b(13) s.t. k=(1−δ)k +zf (k) −c (14) b i=biIii=01n (15) where =indicates that the equality is satisfied approximately, b≡(b0b1bn)∈Rn+1 and b≡(b 0b 1b n)∈Rn+1,andtheIisin(15) are determined by (12). In the transformed Bellman equation (13)–(15), the effect of uncertainty on the solution is summarized by expected values in (12) that determine a relation between band b. To solve the transformed Bellman equation, we proceed in two steps. First, construct {I0In}using (12); second, find bthat solves (13)–(15). Apart from the fact that band bdiffer, the transformed Bellman equation (13)–(15) is standard and can be solved using a variety of solution algorithms available in the literature; we show some algorithms in Section 2.5. The only difference is that we are able to construct expectation functions faster and/or more accurately. 2.3.2 Euler equation with precomputation of integrals Precomputation of integrals is trickier in the Euler equation. The precomputation technique introduced in Section 2.2 assumes that a function we parameterize is the same as a function for which we need to compute expectation. This is true for the Bellman equation approach when we parameterize V(kz)and compute E[V(k z)]. However, this is not true for the existing Euler equation approaches that parameterize such policy functions as the expectation function (Den Haan and Marcet (1990)), the capital function (Judd (1992)), and the labor function (Maliar and Maliar (2005b)), because these functions do not coincide with the integrand of the expectation function in the Euler equation E[u(c)(1−δ+zf(k))]. None of these parameterizations allows us to precompute the integrals in the Euler equation in the way we did in the Bellman equation.5 To adapt the precomputation technique to the Euler equation methods, we reparameterize the Euler equation, namely, we introduce a new variable qthat represents the integrand of the expectation function in the Euler equation (7): q≡u(c)1−δ+zf (k)(16) In terms of qand q, the Euler equation (7)is q 1−δ+zf (k) =βEq(17) We now have the same function on the left side q=Q(kz) as the one inside the expectation on the right side E[Q(kz)]and we can use the precomputation result. By 5In particular, the parameterized expectation algorithm of Den Haan and Marcet (1990) parameterizes the expectation function but to precompute integrals, we need to parameterize the integrand of the expectation function.
Quantitative Economics 8 (2017) How to solve dynamic stochastic models 865 Envelope condition method The ECM of Maliar and Maliar (2013)reliesonaconventional exogenous grid (kz) and simplifies the root-finding using a different mechanism. Namely, to construct policy functions, ECM uses the envelope condition instead of the first-order conditions (FOCs) used by conventional VFI and EGM. For problem (4)–(6), the envelope condition provides a convenient closed-form expression for consumption policy function: c=u−1V1(kz) 1−δ+zf (k)(26) Given c, we compute kfrom budget constraint (5) and, finally, evaluate E[V(k z)]on the right side of the Bellman equation (4)toconstructV(kz)on the left side. In this model, ECM is simpler than Carroll’s (2006) EGM since all the policy functions can be constructed analytically and a numerical solver need not be used at all (not even once). We compare the accuracy and speed of the ECM with and without precomputation of integrals in Table 5. The gains from precomputation of integrals are higher for ECM than for previously considered VFI and EGM. Since ECM avoids finding the roots of nonlinear equations, the cost of approximating integrals constitutes a larger fraction of total costs. As a result, savings in the cost of integration are also more substantial in percentage terms for ECM than for the other methods. Remark 3. Our results should not be interpreted as a comparison of ECM and EGM. First, as we said earlier, the cost of EGM in this specific model can be reduced by using the change of variables proposed in Carroll (2006); second, in more complicated models (e.g., the one with valued leisure), neither ECM nor EGM can avoid root-finding completely, although they both can simplify it.7Maliar and Maliar (2013,2014) show that in Table 5. Accuracy and cost of the ECM method with and without precomputation. γ=1/3γ=3 No Precomputation Precomputation No Precomputation Precomputation Polynomial Degree L1L∞CPU L1L∞CPU L1L∞CPU L1L∞CPU 2nd −402 −352 046 −402 −352 014 −343 −243 061 −343 −243 018 3rd −538 −464 062 −538 −464 019 −438 −311 103 −438 −311 032 4th −665 −577 083 −665 −577 026 −527 −382 148 −527 −382 046 5th −797 −685 087 −797 −685 027 −605 −445 183 −605 −445 059 Note: The statistics L1and L∞are, respectively, the average and maximum of absolute approximation errors across optimality condition and test points (in log10 units) on a stochastic simulation of 10,000 observations; CPU is the time necessary for computing a solution (in seconds); γis the coefficient of risk aversion. 7For the studied model, Carroll’s (2006) technique requires us to introduce a new variable Y≡zf(k)+ (1−δ)kand to rewrite the Bellman equation as V(Yz)=max{u(Y −k)+βE[V(Yz)]}. Using an endogenous grid on k(instead of k) makes it possible to compute E[V(Yz)]with just one evaluation per iteration. However, in a similar model with labor, Y=zf(k)+(1−d)kdepends on labor , and E[V(Yz)] cannot be computed without solving for labor policy function; see Barillas and Fernandez-Villaverde (2007) for a method that recomputes the endogenous grid by iterating on the labor policy function, and see Maliar and Maliar (2013) for a method that computes labor using a numerical solver.
866 Judd, Maliar, Maliar, and Tsener Quantitative Economics 8 (2017) a model with valued leisure, the performance of ECM and EGM is very similar in terms of both accuracy and computational expense. Our point is that for all such methods, precomputation of integrals can help save on costs and increase the accuracy of solutions. 2.5.3 Numerical results for Euler equation methods We next demonstrate that precomputation of integrals can also reduce the cost of Euler equation methods. Our benchmark solution method is similar to projection methods in Judd (1992)inthatitusesa rectangular grid and deterministic (quadrature) integration methods, but it also uses fixed-point iteration with damping as in Den Haan and Marcet (1990)andJudd, Maliar, and Maliar (2011). Parameterizing Qfunction Our first Euler equation method parameterizes the new variable qand uses a change of variables that eliminates consumption. In terms of q, the Euler equation is given by (17), and the budget constraint is k=(1−δ)k +zf (k) −u−1q 1−δ+zf (k)(27) The results of our numerical experiments for the Euler equation algorithm considered are provided in Table 6. The main finding in the table is that savings in costs from precomputation of integrals are much larger for the Euler equation method than for the previously considered VFI and EGM methods, and are similar to those under the ECM method. Here, the running time is typically reduced by a factor of 3. Such a large reduction in costs is due to the fact that the considered Euler equation method avoids finding the roots of nonlinear equations and spends the largest fraction of total costs on integration. Parameterizing capital function jointly with Qfunction In some applications, it might be preferable to approximate policy functions other than qbecause such policy functions are better behaved or more convenient for some computational purposes.8Also, it Table 6. Accuracy and cost of the Euler equation algorithm, parameterizing Qfunction, with and without precomputation. γ=1/3γ=3 No Precomputation Precomputation No Precomputation Precomputation Polynomial Degree L1L∞CPU L1L∞CPU L1L∞CPU L1L∞CPU 2nd −402 −352 053 −402 −352 016 −344 −246 076 −344 −246 023 3rd −538 −464 073 −538 −464 022 −438 −311 137 −438 −311 041 4th −665 −577 098 −665 −577 030 −526 −382 192 −526 −382 059 5th −797 −685 106 −797 −685 032 −605 −445 234 −605 −445 076 Note: The statistics L1and L∞are, respectively, the average and maximum of absolute residuals across optimality condition and test points (in log 10 units) on a stochastic simulation of 10,000 observations; CPU is the time necessary for computing a solution (in seconds); γis the coefficient of risk aversion. 8For example, a capital policy function may have less curvature than consumption and labor policy functions and, hence, it can be approximated more accurately. Moreover, knowing the capital policy function
Quantitative Economics 8 (2017) How to solve dynamic stochastic models 867 could happen that recovering the decision variables such as cand kfrom qrequires us to solve some nonlinear equations numerically. In this case, the advantages of precomputing the integrals might be offset by disadvantages of numerical root-finding. While precomputation of integrals is not directly suitable for approximating functions other than q, it is still possible to approximate other policy functions jointly with q(and thus, sidestep the root-finding step). As an example, we show how to use the precomputation technique for approximating the capital policy function k=K(kz). Premultiplying both sides of Euler equation (17)bykand rearranging the terms, we obtain an equivalent representation of the Euler equation (as long as k>0), k=βEq q1−δ+zf (k)k = K(kz;v) (28) where K(·; v) is a flexible functional form and vis the coefficients vector. Thus, we have expressed kin two ways: one is as today’s choice variable kparameterized by K(·; v); the other is as a combination of the model’s variables that involves conditional expectation of a random variable E[q]. We can compute the capital policy function as a fixed point of (28) by iterating on vuntil convergence. However, to perform such iterations, we also need to construct q,definedby(16), and E[q].Weagaincomparetwootherwise identical numerical solution methods: one that recomputes E[q]in each step of the iterative procedure using a five-node Gauss Hermite quadrature method; the other that precomputes E[q]in the beginning. In Table 7, we provide the results for the Euler equation method approximating the capital policy function Kjointly with Q. The accuracy of this version of the Euler equation algorithm is similar to that of the Euler equation algorithm that parameterizes only Q; however, the running time is somewhat larger. This is because such a method is less numerically stable than the other methods and we had to use damping to enhance numerical stability, namely, we update the policy function only by 15% from one iteration to another. Table 7. Accuracy and speed of the Euler equation algorithm, jointly parameterizing functions Qand K, with and without precomputation. γ=1/3γ=3 No Precomputation Precomputation No Precomputation Precomputation Polynomial Degree L1L∞CPU L1L∞CPU L1L∞CPU L1L∞CPU 2nd −402 −352 318 −402 −352 090 −344 −246 486 −344 −246 087 3rd −538 −464 369 −538 −464 105 −438 −311 726 −438 −311 156 4th −665 −577 436 −665 −577 126 −526 −382 944 −526 −382 224 5th −742 −653 446 −742 −653 133 −605 −445 1077 −605 −445 270 Note: The statistics L1and L∞are, respectively, the average and maximum of absolute residuals across optimality condition and test points (in log 10 units) on a stochastic simulation of 10,000 observations; CPU is the time necessary for computing a solution (in seconds); γis the coefficient of risk aversion. makes it possible to compute endogenous state variables in all grid points without solving for control variables, which is convenient for parallel computation.
868 Judd, Maliar, Maliar, and Tsener Quantitative Economics 8 (2017) The gains from precomputation of integrals are even larger for this version of the Euler equation method than for the one parameterizing the Qfunction; in particular, under high risk aversion coefficient, the running time is reduced by up to a factor of 5.Combining precomputation of integrals with approximations of several functions simultaneously will be useful in more complex models in which there are several Euler equations and several conditional expectation functions that need to be precomputed; this case is studied in Section 3. Conventional Euler equation method parameterizing the capital function In the case of the Euler equation methods, precomputation of integrals has benefits, but may also have costs. Specifically, to be able to precompute expectation functions, we must reparameterize the Euler equation in a particular way and we must compute a decision function for a specific variable, which is the integrand of the expectation function. It could happen that in some applications, this alternative implementation of the Euler equation iteration produces less accurate solutions than the conventional Euler equation method. To see whether or not this is the case in our model, we solve the model by using the conventional Euler equation method that solves for the capital policy function k= K(kz) satisfying k=(1−δ)k +zf (k) −u−1βEuckz1−δ+zfk(29) where the above Euler equation follows by combining (5)and(7). The results are shown in Table 8. As we observe, the accuracy of the conventional Euler equation algorithm is similar to that of other methods considered in this section. The running time for this method is slightly smaller than that of the Euler equation method that parameterizes capital policy function jointly with Q, however, is significantly larger than that of the Euler equation method that parameterizes only function Q. Thus, in this particular case, a reparameterization of the Euler equation in terms of Qhas no negative side effects on either accuracy or cost. Table 8. Accuracy and cost of a conventional Euler equation algorithm parameterizing Kwithout precomputation. γ=1/3γ=3 Polynomial Degree L1L∞CPU L1L∞CPU 2nd −402 −353 054 −344 −246 079 3rd −538 −464 073 −438 −311 137 4th −665 −577 097 −526 −382 187 5th −794 −683 102 −605 −445 223 Note: The statistics L1and L∞are, respectively, the average and maximum of absolute residuals across optimality condition and test points (in log10 units) on a stochastic simulation of 10,000 observations; CPU is the time necessary for computing a solution (in seconds); γis the coefficient of risk aversion.
Quantitative Economics 8 (2017) How to solve dynamic stochastic models 869 Model with elastic labor supply We finally show that precomputation of the expectation functions can simplify the construction of numerical solutions to more complex models such as the model with elastic labor supply. In such a model, the representative agent solves V(kz)=max kclu(cl) +βEVkz (30) s.t. k=(1−δ)k +zf (k l) −c (31) lnz=ρln z+ε(32) where ldenotes labor, uis a utility function that is strictly increasing in consumption and strictly decreasing in labor, fis a production functions that is strictly increasing in both arguments, and both utility and production functions, are continuously differentiable and concave. Our goal is to solve for value function V(kz)and policy functions c=C(kz),l=L(k z),andk=K(kz). As far as value-iterative methods are concerned, value function can be precomputed in the model with elastic labor supply in the same way as in the model with inelastic labor supply studied in Section 2.5.2. Let us show how the precomputation analysis can be generalized to the Euler equation methods. An interior solution to (30)–(32)satisfies first-order conditions u1(cl) =βEu1cl1−δ+zf1kl(33) u2(cl) =−u1(cl)zf2(kl) (34) Again, we rewrite the Euler equation in the form that is suitable for precomputation, q 1−δ+zf1(k l) =βEq(35) where qdenotes the integrand of the Euler equation (33), q≡u1(cl)1−δ+zf1(kl)(36) Since we must parameterize qin terms of the state variables q=Q(kz), the construction of the intratemporal choice requires us to solve numerically two equations (34)and (36) with respect to two unknowns cand lgiven (k z). How does this problem compare to the one we face under the conventional Euler equation methods? The answer depends on what function is parameterized in terms of state variables. If we parameterize the capital function k=K(kz), we must also solve numerically two equations (31)and (34) with respect to two unknowns cand lgiven (kz), so the cost will be comparable; the same is true if we parameterize the consumption function c=C(kz). However, if we parameterize the labor function l=L(k z), the intratemporal choice cand ksatisfying (31)and(34) can be derived in a closed form under conventional utility functions such as Cobb–Douglas or addilog, and the cost will be considerably lower; see Maliar and Maliar (2005b). Finally, it is possible to reduce the cost of iteration on either the conventional or transformed Euler equations by constructing the intratemporal choice
870 Judd, Maliar, Maliar, and Tsener Quantitative Economics 8 (2017) Table 9. Accuracy and cost of the ECM and Euler equation algorithms in the model with elastic labor supply with and without precomputation. γ=1/3γ=3 No Precomputation Precomputation No Precomputation Precomputation Polynomial Degree L1L∞CPU L1L∞CPU L1L∞CPU L1L∞CPU 2nd −342 −243 573 −342 −243 373 −343 −247 413 −343 −247 339 3rd −444 −319 803 −444 −319 573 −444 −320 767 −444 −320 641 4th −541 −395 1001 −541 −395 735 −541 −394 1440 −541 −394 945 5th −628 −462 1150 −628 −462 861 −628 −462 1802 −628 −462 1257 Note: The statistics L1and L∞are, respectively, the average and maximum of absolute approximation errors across optimality condition and test points (in log 10 units) on a stochastic simulation of 10,000 observations; CPU is the time necessary for computing a solution (in seconds). functions outside the main iterative cycle on some prespecified grid of points; see Maliar and Maliar (2005b,2014) for this different kind of precomputation results. As in the case of inelastic labor supply, the transformed Euler equation (35) is identical to that under ECM-DVF (up to notation), and, therefore, we can use the code accompanying the paper of Maliar and Maliar (2013) as a basis for our precomputation analysis. Following that paper, we parameterize the model (30)–(32) by an additively separable utility function u(cl) =c1−γ−1 1−γ+B(1−l)1−μ−1 1−μfor which the two equations (34)and(36) can be combined into one equation relating qand l: q=B(1−l)−μ zf2(k l) 1−δ+zf1(kl)(37) We assume that the production function is Cobb–Douglas f(kl) =Akαl1−αwith α=033 and A=1/β−(1−δ) αl1−α(this value normalizes the deterministic steady state of capital to 1). We fix γ=5and μ=5, and we calibrate the parameters β,δ,andBto reproduce the steady state labor, capital–output, and consumption–output ratios of l=1/3, πk=10,andπc=3/4, respectively, and we assume ρ=095 and σ=001.Theimplementation of the solution methods with elastic labor supply is similar to that with inelastic labor supply; see Appendix B for details. To solve for labor choices, we solve (37) numerically in each grid point by using a Newton solver. In Table 9, we compare the results for the value-iterative and Euler equation methods with and without precomputation of integrals. Again, the introduction of precomputation of integrals has no effect on accuracy, but reduces the running time by about 30% for both value-iterative and Euler equation methods. Here, the reduction in cost is not as large as in the model with inelastic labor supply because a large fraction of the total cost comes from applying a numerical solver to multiple grid points. 3. Multivariate continuous-state expectations Although precomputation of integrals can visibly speed up computations even in models with just one shock, the real case of interest is models with multivariate shocks, in
Quantitative Economics 8 (2017) How to solve dynamic stochastic models 871 which the cost of integrals evaluation can increase rapidly with the dimensionality of the problem. As an example, we describe a multicountry model with Nstochastic shocks analyzed in Judd, Maliar, and Maliar (2011). This model provides a convenient framework for testing how the accuracy and cost of numerical solution methods change with the number of state variables, in particular, with the number of exogenous shocks. 3.1 Heterogeneous neoclassical stochastic growth model A social planner maximizes a weighted sum of the expected lifetime utilities of Nagents (interpreted as countries) subject to the aggregate resource constraint, that is, max {ch tkh t+1}h=1N t=0∞ E0 N h=1 τh∞ t=0 βtuch t(38) subject to N h=1 ch t= N h=1zh tfkh t+(1−δ)kh t−kh t+1(39) where ch t,kh t,zh t,u,f,andτhare consumption, capital, productivity level, utility function, production function, and welfare weight of a country h∈{1N}, respectively, βis the discount factor, and δis the depreciation rate. The initial condition (k0z0)is given, where kt≡(k1 tkN t)and zt≡(z1 tzN t). The process for productivity shocks in country his given by lnzh t=ρlnzh t−1+εh t(40) where εt=(ε1 tεN t)∼N(0NΣ) is generated by a multivariate normal distribution with zero mean 0Nand variance–covariance matrix Σ. Here, we assume that all the countries have identical utility and production functions; however, it is straightforward to extend our analysis for the case when the utility and production functions differ across countries; see Maliar, Maliar, and Judd (2011). Bellman equation We can represent problem (38)–(40) in a dynamic programming form, V(kz)=max {(kh)ch}h=1N N h=1 τhuch+βEVkz(41) s.t. N h=1 ch= N h=1zhfkh+(1−δ)kh−kh(42) lnzh=ρlnzh+εh(43) where Vis an optimal value function, ε=((ε1)(εN))∼N(0NΣ);k≡(k1 kN),andz≡(z1zN).
872 Judd, Maliar, Maliar, and Tsener Quantitative Economics 8 (2017) Euler equation An interior solution to the planner’s problem (38)–(40)satisfiesasetof Nfirst-order conditions uch=βEuch1−δ+zhfkh(44) where h=1N. Under the Euler equation approach, we solve for policy functions ch=Ch(kz)and (kh)=Kh(kz),h=1N, that satisfy (39), (40), and (44). 3.2 Precomputation of multivariate integrals under polynomial approximations It turns out that multivariate integrals can be precomputed in essentially the same way as univariate integrals in Section 2.2.Letk=(k1kN)and z=(z1zN)be the economy’s state variables and let a policy function P(kz)be approximated with ordinary polynomial function P(kz;b) =b0+b1k1+···+bNkN+bN+1z1+··· +b2NzN+b11k12+···+bNNkN2 + ij∈{1N}i=j bij kikj+ i∈{1N}j∈{N+12N} bij kizj + ij∈{N+12N} bij zizj+··· +bN+1N+1z1L+···+b2N2NzNL (45) where b=(b0b1b2Nb11b2N2Nb11b2N2N)∈Rn+1is a vector of polynomial coefficients and Lis a polynomial degree. Again, taking into account that k≡((k1)(kN))is known and that (zh)=(zh)ρexp((εh))has a known distribution function, we can represent a conditional expectation of P(kz;b) as EPkz;b=Eb0+b1k1+···+bNkN+bN+1z1+··· +b2NzN+b11k12+···+bNNkN2 + ij∈{1N} bij kikj+ i∈{1N}j∈{N+12N} bij kizj + ij∈{N+12N} bij zizj+··· +bN+1N+1z1L+···+b2N2NzNL =Pkzρ;b (46) where b≡(b 0b 1b 2Nb 11b 2N2Nb 11b 2N2N)∈Rn+1, and the coefficients b iand biare related by b i=biIii=01n (47)
Quantitative Economics 8 (2017) How to solve dynamic stochastic models 873 where Ii=Eexpl iε=Eexpl1 iε1+l2 iε2+···+lN iεN(48) with (l1 ilN i)≡libeing the powers on (z1)(zN)in the ith monomial term, respectively. As in the one-dimensional case, {I0In}depend only on the properties of exogenous shock εand can be computed up-front (i.e., precomputed). Once Isareconstructed, we again obtain that conditional expectation of a polynomial function is given by the same polynomial function but evaluated at a different coefficients vector, that is, E[P(kz;b)]=P(kzρ;b); this result is parallel to what we had in the univariate case in Section 2.2. 3.3 Characterizing the solution under precomputation of integrals We now show how to precompute expectations in the Bellman and Euler equations in thecaseofmultivariateshocks. 3.3.1 Bellman equation with precomputation of integrals Using the precomputation result (46), we rewrite the Bellman equation (41)as V(kz;b) =max {(kh)ch}h=1N N h=1 τhuch+β Vkzρ;b(49) s.t. (42), (43), and (46), where a relation between band bis summarized by (47)and(48). To solve the transformed Bellman equation, we proceed in two steps. First, we construct {I0In}using (48); second, we find bthat solves the Bellman equation (49). 3.3.2 Euler equation with precomputation of integrals Let us introduce a new variable qhthat represents the integrand of the expectation function in (44), qh≡uch1−δ+zhfkh(50) where h=1N.Intermsofqhand (qh), the Euler equation (44)is qh 1−δ+zhfkh=βEqh(51) We parameterize the function qh=Qh(kz)with a flexible functional form Qh(kz;bh), h=1N. Using the precomputation result (46), we formulate a version of the Euler equation (51)intermsof Qh, Qhkz;bh 1−δ+zhfkh =β Qhkzρ;bh(52) where bhand (bh)are related by (47)and(48) for all h=1N. To solve the transformed Euler equation, we proceed in two steps. First, we construct {I0In}using (48); second, we find bthat solves (42), (43), (47), (48), and (52).
874 Judd, Maliar, Maliar, and Tsener Quantitative Economics 8 (2017) 3.4 Analytical versus numerical construction of integrals In this section, we compare the accuracy and computational expense of analytical and numerical approximation of integrals. Analytical construction of integrals under precomputation For an important special case of a multivariate normal distribution ε=((ε1)(εN))∼N(0NΣ),wecan compute Iianalytically. As in the unidimensional case, we complete the square of the expression inside of the exponential function, Ii=Eexpl1 iε1+···+lN iεN=Eexpl iε =1 (2π)N/2|Σ|1/2+∞ −∞ ···+∞ −∞ expl iεexp−1 2εΣ−1εdε =1 (2π)N/2|Σ|1/2+∞ −∞ ···+∞ −∞ exp−1 2εΣ−1ε−2l iεdε =1 (2π)N/2|Σ|1/2+∞ −∞ ···+∞ −∞ exp−1 2ε−ΣliΣ−1ε−Σli−l iΣlidε =expl iΣli 2 (53) where in the last equality, we use the fact that +∞ −∞ ···+∞ −∞ f(x)dx=1for a density function f(x)of a normally distributed variable xwith mean Σliand variance–covariance matrix Σ. Numerical construction of integrals Similarly to the univariate case, multivariate numerical integration methods approximate integrals using a weighted sum (8), where the integration nodes and weights are constructed in the multidimensional space. Gauss Hermite quadrature can be extended to the multivariate case using a product rule; however, product rules are prohibitively costly even for a moderately large number of shocks N; in particular, the cost of GH(2),GH(5),andGH(10)grows exponentially as 2N,5N, and 10N, respectively. A tractable and sufficiently accurate alternative to product rules is nonproduct monomial integration methods; see Stroud (1971) for a collection of monomial formulas, and see Judd, Maliar, and Maliar (2011) for a description of two monomial integration methods with 2Nand 2N2+1nodes, referred to as M1 and M2, respectively. A quasi-Monte Carlo method is another class of numerical integration methods that is tractable in problems with high dimensionality; see Rust (1996)andGeweke (1996)fora discussion (we do not study such methods). Comparison of the accuracy Constructing the expectation functions of polynomial approximations (46) includes evaluation of expectations of exponents of different linear combination of shocks in (48). Our numerical experiments show that for the same sum of the powers on (z1)(zN)in the ith monomial term, l1 i+···+lN i, numerical integration methods produce the largest approximation errors whenever all lis, except one, are zeros. For example, E[exp(2ε1)]is approximated less accurately than E[exp(ε1+ε2)];
Quantitative Economics 8 (2017) How to solve dynamic stochastic models 881 Bellman equation In the model (57)–(60)withJexogenous states, the optimal value function can be represented by Jstate-contingent value functions Vj(k) that satisfy the Bellman equation V(k) =max kc u(c) +β J j=1 πjVjk(61) s.t. k=Wz +(1+R)k −c (62) k≥−φ (63) where j ∈{1J}and πj satisfies (59). Euler equation The policy functions for the problem (57)–(60) can also be represented as Jstate-contingent functions. A consumption policy function Cj(k) satisfies the Euler equation uC(k)−η=β J j=1 πjuCjk(1+R)(64) where j ∈{1J}and η≥0is a Lagrange multiplier, which is associated with the borrowing constraint (63), that satisfies a Kuhn–Tucker condition η(k+φ) =0. 4.2 Precomputation of expectation functions Let us show how to precompute conditional expectation functions under the assumption of ordinary polynomial approximation. Each state-contingent value or policy function is approximated by a polynomial function Pk;b=b 0+b 1k+b 2k2+···+b nkL(65) where =1J,b≡(b 0b 1b n)∈Rn+1is a vector of polynomial coefficients, and Lis the degree of polynomial. Given today’s state , conditional expectation of (65)is EPk;bj|kz= J j=1 πjPk;bj=π1Pk;b1+···+πJPk;bJ =π1b1 0+b1 1k+···+b1 nkL+··· +πJbJ 0+bJ 1k+···+bJ nkL =π1b1 0+···+πJbJ 0+π1b1 1+···+πJbJ 1k+··· +π1b1 n+···+πJbJ nkL =Pk;b (66)
882 Judd, Maliar, Maliar, and Tsener Quantitative Economics 8 (2017) where (b)≡((b 0)(b 1)(b n))∈Rn+1is given by b i= J j=1 πjbj ii=01n (67) Condition (67) provides a basis for precomputation of conditional expectation functions in the discrete-state case. As in the continuous-state case, conditional expectation of a polynomial function is given by the same polynomial function, but evaluated at a different coefficients vector, that is, E[P(k;bj)|kz]=P(k;(b)),where(b)is determined by (67). Remark 4. Similar to the continuous-state case, the precomputation result (66)and (67) also holds for piecewise polynomial approximations. As an example, consider again a collection of piecewise linear polynomial functions on a grid [k1kM]such that in each state =1J, we have a local polynomial function P[kmkm+1]k;b(m)=b(m) 0+b(m) 1k where m∈{1M}. Then the expectation function in each interval [kmkm+1]is given by EP[kmkm+1]k;b(m)= J j=1 πjb(jm) 0+b(jm) 1k=b(m) 0+b(m) 1k where (b(m) 0)=[π1b(1m) 0+···+πJb(Jm) 0]and (b(m) 1)=[π1b(1m) 1+···+πJb(Jm) 1], i.e., (66)and(67) hold in each local area. The precomputation results for higher-order piecewise polynomial approximations can be shown in a similar manner. 4.3 Characterizing the solution under precomputation of expectation functions In this section, we show how to precompute expectations in the Bellman and Euler equations in the discrete-shock case. 4.3.1 Bellman equation with precomputation of expectation functions Using precomputation result (66), we rewrite the Bellman equation (61)–(63)as Vk;b =max kc u(c) +β Vk;b (68) s.t. (62), (63), (66), and (67), where ∈{1J}. In the transformed Bellman equation (68), the effect of uncertainty on the solution is captured by a relation between bjand (b), which is described by (66) and (67).
Quantitative Economics 8 (2017) How to solve dynamic stochastic models 883 4.3.2 Euler equation with precomputation of expectation functions To make the Euler equation (64) suitable for precomputation of expectation functions, we use the same change of variables as in the continuous-state case: q≡u(c)[1+R](69) We next rewrite the Euler equation (64) and budget constraint (62) by eliminating cand by expressing them in terms of qapproximated using a polynomial function Q(k;bj), Qk;b 1+R−η=β J j=1 πjQk;bj(70) where j ∈{1J}and Lagrange multiplier ηsatisfies the same conditions as in (64). Finally, under polynomial approximation (65), we can use the precomputation result (66) to rewrite the Euler equation (70)as Qk;b 1+R−η =β Qk;b(71) where ∈{1J}and (b)is determined by (67). Again, after expectation functions are precomputed, all the effect of uncertainty on the solution is compressed into a mapping between the vectors band (b), which is described by (66)and(67). 4.4 Numerical assessment of gains from precomputation of expectation functions in Aiyagari’s (1994)model We assess gains from precomputation of expectation functions in Aiyagari’s (1994) model. With discrete-state shocks, expectation functions are computed exactly both with and without precomputation analysis. Thus, all the gains will come in terms of a cost reduction, unlike in the case of continuous-state shocks, in which we may have both a lower cost and higher accuracy. We consider a version of the Euler equation solution method described in Maliar and Maliar (2006), and we show how to precompute expectation functions in the MATLAB code accompanying that paper. A detailed description of the solution algorithm is provided in Appendix D. 4.4.1 Methodology and implementation details Following the literature, we study a stationary equilibrium in which aggregate variables are constant over time. A stationary equilibrium is defined as a stationary probability measure x, optimal consumption and capital policy functions Cj(k) and Kj(k), respectively, j=1J,andpositivereal number (R W ) such that (i) xsatisfies x=K×ZP(kzB)dx for all B∈B, (ii) Cj(k) and Kj(k) solve (57)–(60)forgiven(RW ), (iii) R=f1(KN) −δand W=f2(KN), (iv) N=Zztdx and K=K×ZK(kz)dx,
884 Judd, Maliar, Maliar, and Tsener Quantitative Economics 8 (2017) where Z≡{z1zJ}is a set of all possible productivity levels, K≡[−φk]is an interval of possible asset holdings, Bis a collection of Borel subsets of all possible individual states K×Z,andP(kzB) is conditional probability that an agent with today’s state (kz) will be in a set B∈Bin the next period. To approximate the capital decision functions Kj(k),j=1J, we use piecewise linear polynomial approximation. We use an unevenly spaced grid of 1000 points: we placed a dense grid of 900 grid point in the area of the borrowing limit in which the solution is highly nonlinear, and we placed the remaining 100 grid points to cover the rest of the solution domain. Precomputation of expectation functions can be easily implemented in MATLAB in the one-dimensional case by using a routine interp1,whichproduces linear, cubic, and spline approximations, and which returns both the knots and coefficients. We take these coefficients as input, and we construct the precomputation coefficients as implied by (66)and(67). In our calibration procedure, we follow Aiyagari (1994). The model’s period is 1 year. We assume u(ct)=c1−γ t−1 1−γwith γ∈{ 1 33}and f(KN) =KαN1−αwith α=036 (we use a normalization N=1for convenience). We use β=096 and δ=008.Weset the debt limit at φ=0and the upper bound on capital at k=δ1/(α−1).AsinAiyagari (1994), we assume that idiosyncratic shocks follow an AR(1)process: logzt+1=ρlog zt+ σ(1−ρ2)1/2εt+1,whereρ∈{0609},σ∈{0204},andεt+1∼N(01),andwediscretize this process into a seven-state Markov chain using Tauchen’s (1986)procedure;seealso Tauchen and Hussey (1991) for a discussion of this method.9We solve for the equilibrium interest rate by using stochastic simulation and a bisection method as in Aiyagari (1994), and we construct stationary probability distribution as in Rios-Rull (1997). Finally, as a measure of accuracy, we report the mean and maximum of absolute unit-free residuals in the transformed Euler equation (71) on a stochastic simulation of 10,000 observations. 4.4.2 Numerical results We illustrate a numerical solution to Aiyagari’s (1994)modelin Figure 1. The properties of the solution to Aiyagari’s (1994) model are well known. The capital function is close to linear except in the kink area, and the consumption function has strong nonlinearities in the area of the kink. The constructed probability distribution indicates that our upper bound is chosen large enough so that probability of reaching this bound is low (in our simulations, it was actually never reached). But the limit on borrowing is important and is occasionally reached. In Table 13, we compare the accuracy and cost for two versions of the algorithm with and without precomputation of expectation functions. As we can see, the residuals are again practically identical for the methods with and without precomputation. The solutions are very accurate everywhere in the solution domain (average error is of order 10−6percent), except for a very close neighborhood of the 9We discretize the process for shocks to provide a numerical illustration of the precomputation technique for discrete-state shocks. Alternatively, we could have solved Aiyagari’s (1994) model by applying the precomputation technique to the original AR(1)process with continuous-state shock.
Quantitative Economics 8 (2017) How to solve dynamic stochastic models 885 Figure 1. Numerical solution to Aiyagari’s (1994)model. borrowing limit, where the maximum residuals are much larger and sensitive to small changes in the construction of the numerical solutions. To gain insight into the accuracy deterioration near the kink, we note that qhas a spike in the lowest asset and productivity states in Figure 1. As follows from the ECM analysis in Section 2.3.3,qis the marginal value function, and it is largest when consumption is lowest. It is typically better to approximate numerically functions that do not have much curvature, for example, it is usually better to approximate consumption rather than the marginal utility of consumption (especially with high risk aversion). This reasoning suggests that a highly nonlinear variable qmay not be the best candidate for approximating in the context of Aiyagari’s (1994) model. To check on this conjecture, namely, to see how the choice of a function to parameterize affects the accuracy of solutions, we also solved the model by using the stanTable 13. Accuracy and cost of the Euler equation algorithm with and without precomputation. γ=1/3γ=3 No Precomputation Precomputation No Precomputation Precomputation Parameterization L1L∞CPU L1L∞CPU L1L∞CPU L1L∞CPU ρ=06,σ=02−634 −299 1503−634 −299 726−493 −242 1547−493 −242 728 ρ=06,σ=04−618 −358 1406−618 −358 700−467 −187 1433−467 −187 681 ρ=09,σ=02−678 −314 1510−678 −314 688−536 −271 1556−536 −271 712 ρ=09,σ=04−662 −317 1357−662 −317 659−504 −196 1344−504 −196 651 Note: The statistics L1and L∞are, respectively, the average and maximum of absolute residuals across optimality condition and test points (in log 10 units) on a stochastic simulation of 10,000 observations; CPU is the time necessary for computing a solution (in seconds); γis the coefficient of risk aversion.
886 Judd, Maliar, Maliar, and Tsener Quantitative Economics 8 (2017) Table 14. Accuracy and cost of the Euler equation algorithm parameterizing K. γ=1/3γ=3 Parameterization L1L∞CPU L1L∞CPU ρ=06,σ=02−623 −304 1115−475 −302 1154 ρ=06,σ=04−605 −367 1042−446 −188 1066 ρ=09,σ=02−667 −316 1093−515 −293 1153 ρ=09,σ=04−647 −307 1009−482 −210 1000 Note: The statistics L1and L∞are, respectively, the average and maximum of absolute residuals across optimality condition and test points (in log10 units) on a stochastic simulation of 10,000 observations; CPU is the time necessary for computing a solution (in seconds); γis the coefficient of risk aversion. dard Euler equation algorithm that parameterizes only the capital policy function and that sidesteps precomputation. In Table 14, we report the average and maximum Euler equation errors produced by this standard method. The maximum residuals produced by the method parameterizing kare somewhat smaller than those produced by the algorithm parameterizing q, although the average residuals are very similar. These results suggest that both methods are very accurate away from the kink, but their accuracy deterioration near the kink is somewhat larger for the method parameterizing qthan for the method parameterizing k. More generally, these result suggest that in some problems, it might be preferable to approximate decision functions others than qand to compute expectations in the conventional way even if this implies a higher cost. Our main result is that precomputation of expectation functions reduces the running time by about 50%. This reduction in running times corresponds to a Markov chain with seven states that we used following Aiyagari (1994); see also Maliar, Maliar, and Valli (2010). If we use a Markov chain with more states, the savings in costs from precomputation of expectation functions will be much larger. We conjecture that the gains from precomputation of expectations functions will be especially sizable in models with multivariate shocks in which Markov chains may have a very large number of states. 5. Approximating functions consistent with precomputation of expectation functions All our precomputation results are obtained for ordinary polynomial functions. In this section, we establish other families of approximating functions for which expectation functions can be precomputed. The cases of continuous-state and discrete-state exogenous shocks are analyzed in Sections 5.1 and 5.2, respectively. 5.1 Continuous-state problems Let x∈Rnxand z∈Rnzbe vectors of endogenous and exogenous state variables, respectively. We make the following three assumptions.
Quantitative Economics 8 (2017) How to solve dynamic stochastic models 887 Assumption 1. Each basis function of an approximating polynomial is multiplicatively separable in xand z,that is, P(xz;b) = n i=0 biψi(x)ϕi(z) (72) where b≡(b0bn)∈Rn+1,ψi(x)ϕi(z) is the ith basis function,by convention, ψ0(x)ϕ0(z) is equal to 1,and ψi:Rnx→Rand ϕi:Rnz→Rfor i=1n. Assumption 2. Next-period endogenous state variables xare t-measurable. Assumption 3. Exogenous state variables follow a stochastic process z=Z(zε), where ε∈Ω⊆Rnεis a vector of disturbances with a density function w:Rnε→R+ that is known at tand that satisfies ε∈Ωw(ε)dε=1.Moreover,all integrals of type ε∈Ωϕi(Z(zε))w(ε)dεexist and are finite,up to the order of approximation in (72). Proposition 1. Under Assumptions 1–3,we have EPxz;b|xz=PxZ(zμ);b(z)(73) where μ≡E[ε]and b(z) ≡(b 0(z) b n(z)),and bare related by b i(z) =biε∈Ω ϕiZzε ϕiZ(zμ)wεdε(74) Proof. This result is obtained as EPxz;b|xz=ε∈Ω Pxz;bwεdε Assumption 1 =ε∈Ω n i=0 biψixϕizwεdε Assumption 2 = n i=0 biψixε∈Ω ϕiZzεwεdε Assumption 3 = n i=0 biψixϕiZ(zμ)ε∈Ω ϕiZzε ϕiZ(zμ)wεdε (74) = n i=0 b i(z)ψixϕiZ(zμ)(72) =PxZ(zμ);b(z) A few comments about our assumptions are in order. Assumption 1holds not only for ordinary polynomials (9), but also for orthogonal polynomial families (like Chebyshev, Hermite, and Legendre families), as well as for many nonpolynomial families (for example, ψi(x) and ϕi(z) can be trigonometric functions). Assumption 1also holds for piecewise approximations such as piecewise linear functions and cubic splines, as well
888 Judd, Maliar, Maliar, and Tsener Quantitative Economics 8 (2017) as mixed approximations that combine higher-order polynomials for some variables with piecewise functions for others. Furthermore, Assumption 1is also consistent with anisotropic approximating functions that allow use for different degrees of polynomial approximation in different model’s variables; see, for example, an anisotropic Smolyak method introduced in Judd et al. (2014). Anisotropic families may be especially useful in high-dimensional problems if there is more curvature in some directions than in others or if some of the state variables are exogenous so that higher-order smoothness in those variables is not really required. Assumption 2means that next-period endogenous state variables such as kare known with certainty at present. There are some interesting economic models in which this assumption is not satisfied, including asset-pricing models, dynamic games, input– output models with private shocks, and research and development models. It would be of interest to extend precomputation technique to the case when all or some endogenous state variables are random. Finally, Assumption 3means that endogenous state variables have finite moments, which is normally satisfied in economic models. The assumption that the random shocks have a density function wis used for expositional convenience: it is sufficient to assume instead that random shocks have a distribution function that produces finite integrals in all basis functions considered. To construct an expectation function, we choose the next-period reference point Z(z0). This choice is a matter of convenience. In our precomputation examples (as well as in many macroeconomic models), the stochastic process for an exogenous state variable is multiplicatively separable in the current state term and future disturbances, which makes the integral (74) independent of the current economy’s state z,thatis, b i(z) =b ifor all z.Forexample,in(22), the process is Z(zε)=zρexp(ε)and, hence, we have b i(z) =biε∈Ω zρexp(ε) zρexp(0)w(ε)dε≡b ifor all z. However, in general, the precomputed integral b i(z) will depend on the current economy’s state z, and we must construct a mapping b i(z), either analytically or numerically. 5.2 Discrete-state problems Let x∈Rnxand z∈{z1zJ}∈Rnzbe vectors of endogenous and exogenous state variables, respectively. We formulate a set of Assumptions 1–3for the discrete-shock case that are parallel to Assumptions 1–3for the continuous-shock case. Assumption 1.For each state ∈{1J},an approximating function is given by Px;b= n i=0 b iψi(x) (75) where b≡(b 0b n)∈Rn+1;ψi(x) is the ith basis function,by convention,ψ0(x) is equal to 1,and ψi:Rnx→Rfor i=1n. Assumption 2.Next-period endogenous state variables xare t-measurable.
Quantitative Economics 8 (2017) How to solve dynamic stochastic models 889 Assumption 3.Exogenous state variables follow a stochastic process that has a countable number of states z∈[z1zJ]and transition probabilities are given by πj = Pr(z=zj|z=z),where ∈{1J}. Proposition 2. Under Assumptions 1–3,we have EPx;bj|xz=Px;b(76) where (b)≡((b 0)(b 1)(b n))∈Rn+1and the coefficients (b)and bjare related by b i= J j=1 πjbj ii=01n (77) The proof follows the derivations in (67). In the case of discrete-state shocks, our generalizations go in two dimensions. First, we allow for any additively separable approximating family, in addition to the ordinary polynomial functions used in (65). Second, we allow for multiple endogenous and exogenous state variables. This means that if each shock i∈{1N}has Jistates, then we need to approximate J=J1×··· ×JNdecision functions of exogenous state variables.10 In this case, the number of exogenous states grows exponentially with the number of shocks. To alleviate the curse of dimensionality, one can use some nonproduct rules for selecting a smaller set of states that can approximate the process for multivariate shocks sufficiently well; see, for example, a cluster grid and epsilon distinguishable set techniques introduced in Maliar and Maliar (2015). Finally, integrals can also be precomputed in models with discrete control variables or a mixture of discrete and continuous control variables. Indeed, precomputation results (10), (46), and (66) hold independently of whether control variables are continuous or discrete. Discrete controls are used in the recent literature on structural estimation in which dynamic stochastic models must be solved repeatedly a large number of times, and savings on cost from precomputation of expectation functions can be of value in this computationally intense class of problems. 6. Conclusion A vast majority of the existing solution methods in the economics literature use approximating families of functions studied in the present paper. For such methods, we can precompute integrals in the stage of initialization and in effect, we can transform a stochastic problem into a deterministic problem. The technique of precomputation of integrals is very general and can be applied to essentially any set of equations that contains expectation functions. It works for both continuousand discrete-state shocks, and it can be combined with other computational techniques used by the existing solution 10In the multivariate case, it is also possible to start from a continuous-state multivariate Markov process and to discretize this process into a Markov chain with a countable number of states by using the analysis of Tauchen (1986) and Tauchen and Hussey (1991).
890 Judd, Maliar, Maliar, and Tsener Quantitative Economics 8 (2017) methods, including a variety of solution domains, integration rules, fitting methods, and iterative schemes for finding unknown parameters of the approximating functions. Integrals can be constructed in a closed form for models with unior multivariate normally distributed shocks, which is the case of special interest to economics. In those cases in which integrals cannot be constructed analytically, we can precompute them numerically by applying very accurate computational methods since this is a one time fixed cost. In addition, precomputation of integrals is very simple to implement. For small problems with few shocks that must be solved just once, precomputation of integrals is useful but not critical. Nevertheless, some interesting economic models in the recent literature may have dozens of exogenous shocks, including largescale new Keynesian models used by central banks for projection and policy analysis, large-scale overlapping generation models, heterogeneous agents models, and climate change models; see Maliar and Maliar (2014) for a discussion and further examples of large-scale applications.11 Furthermore, the literature on structural estimation must recompute a solution to dynamic economic models a large number of times under different parameter vectors. For these and other computationally intense applications, precomputation of integrals can be the only tractable alternative. References Aiyagari, R. (1994), “Uninsured idiosyncratic risk and aggregate saving.” Quarterly Journal of Economics, 109 (3), 659–684. [851,854,860,880,883,884,885,886] Arellano, C., L. Maliar, S. Maliar, and V. Tsyrennikov (2016), “Envelope condition method with an application to default risk model.” Journal of Economic Dynamics and Control, 69, 436–459. [854,859,862] Aruoba, S., J. Fernández-Villaverde, and J. Rubio-Ramírez (2006), “Comparing solution methods for dynamic equilibrium economies.” Journal of Economic Dynamics and Control, 30, 2477–2508. [852,863] Barillas, F. and J. Fernandez-Villaverde (2007), “A generalization of the endogenous grid method.” Journal of Economic Dynamics & Control, 31, 2698–2712. [865] Bewley, T. (1977), “The permanent income hypothesis: A theoretical formulation.” Journal of Economic Theory, 16 (2), 252–292. [880] Carroll, C. D. (1992), “The buffer-stock theory of saving: Some macroeconomic evidence.” Brookings Papers on Economic Activity, 1992 (2), 61–156. [880] Carroll, K. (2006), “The method of endogenous grid points for solving dynamic stochastic optimal problems.” Economics Letters, 91, 312–320. [853,862,863,864,865] 11In particular, precomputation of integrals can significantly reduce the computational expense in the analysis of Lepetyuk, Maliar, and Maliar (2017), who construct a global nonlinear solution to a version of the Canadian Central Bank ToTEM model with 21 state variables in the presence of zero lower bound on the nominal interest rate.