The Federal Reserve Board eagle logo links to home page

Skip to: [Printable Version (PDF)] [Bibliography] [Footnotes]
Finance and Economics Discussion Series: 2012-04 Screen Reader version

Computing DSGE Models with
Recursive Preferences and Stochastic Volatility*

Dario Caldara
Jesus Fernandez-Villaverde
Juan F. Rubio-Ramirez
Wen Yao

Abstract:
This paper compares different solution methods for computing the equilibrium of dynamic stochastic general equilibrium (DSGE) models with recursive preferences such as those in Epstein and Zin (1989 and 1991) and stochastic volatility. Models with these two features have recently become popular, but we know little about the best ways to implement them numerically. To fill this gap, we solve the stochastic neoclassical growth model with recursive preferences and stochastic volatility using four different approaches: second- and third-order perturbation, Chebyshev polynomials, and value function iteration. We document the performance of the methods in terms of computing time, implementation complexity, and accuracy. Our main finding is that perturbations are competitive in terms of accuracy with Chebyshev polynomials and value function iteration while being several orders of magnitude faster to run. Therefore, we conclude that perturbation methods are an attractive approach for computing this class of problems.

Keywords: DSGE models, recursive preferences, perturbation.

JEL Classification: C63, C68, E32.




1 Introduction

This paper compares different solution methods for computing the equilibrium of dynamic stochastic general equilibrium (DSGE) models with recursive preferences and stochastic volatility (SV). Both features have become very popular in finance and in macroeconomics as modelling devices to account for business cycle fluctuations and asset pricing. Recursive preferences, as those first proposed by Kreps and Porteus (1978) and later generalized by Epstein and Zin (1989 and 1991) and Weil (1990), are attractive for two reasons. First, they allow us to separate risk aversion and intertemporal elasticity of substitution (EIS). Second, they offer the intuitive appeal of having preferences for early or later resolution of uncertainty (see the reviews by Backus et al., 2004 and 2007, and Hansen et al. , 2007, for further details and references). SV generates heteroskedastic aggregate fluctuations, a basic property of many time series such as output (see the review by Fernandez-Villaverde and Rubio-Ramirez, 2010), and adds extra flexibility in accounting for asset pricing patterns. In fact, in an influential paper, Bansal and Yaron (2004) have argued that the combination of recursive preferences and SV is the key for their proposed mechanism, long-run risk, to be successful at explaining asset pricing.

But despite the popularity and importance of these issues, nearly nothing is known about the numerical properties of the different solution methods that solve equilibrium models with recursive preferences and SV. For example, we do not know how well value function iteration (VFI) performs or how good local approximations are compared with global ones. Similarly, if we want to estimate the model, we need to assess which solution method is sufficiently reliable yet quick enough to make the exercise feasible. More important, the most common solution algorithm in the DSGE literature, (log-) linearization, cannot be applied, since it makes us miss the whole point of recursive preferences or SV: the resulting (log-) linear decision rules are certainty equivalent and do not depend on risk aversion or volatility. This paper attempts to fill this gap in the literature, and therefore, it complements previous work by Aruoba et al. (2006), in which a similar exercise is performed with the neoclassical growth model with CRRA utility function and constant volatility.

We solve and simulate the model using four main approaches: perturbation (of second and third order), Chebyshev polynomials, and VFI. By doing so, we span most of the relevant methods in the literature. Our results provide a strong guess of how some other methods not covered here, such as finite elements, would work (rather similar to Chebyshev polynomials but more computationally intensive). We report results for a benchmark calibration of the model and for alternative calibrations that change the variance of the productivity shock, the risk aversion, and the intertemporal elasticity of substitution. In that way, we study the performance of the methods both for cases close to the CRRA utility function with constant volatility and for highly non-linear cases far away from the CRRA benchmark. For each method, we compute decision rules, the value function, the ergodic distribution of the economy, business cycle statistics, the welfare costs of aggregate fluctuations, and asset prices. Also, we evaluate the accuracy of the solution by reporting Euler equation errors.

We highlight four main results from our exercise. First, all methods provide a high degree of accuracy. Thus, researchers who stay within our set of solution algorithms can be confident that their quantitative answers are sound.

Second, perturbations deliver a surprisingly high level of accuracy with considerable speed. Both second- and third-order perturbations perform remarkably well in terms of accuracy for the benchmark calibration, being competitive with VFI or Chebyshev polynomials. For this calibration, a second-order perturbation that runs in a fraction of a second does nearly as well in terms of the average Euler equation error as a VFI that takes ten hours to run. Even in the extreme calibration with high risk aversion and high volatility of productivity shocks, perturbation works at a more than acceptable level. Since, in practice, perturbation methods are the only computationally feasible method to solve the medium-scale DSGE models used for policy analysis that have dozens of state variables (as in Smets and Wouters, 2007), this finding has an outmost applicability. Moreover, since implementing second- and third-order perturbations is feasible with off-the-shelf software like Dynare, which requires minimum programming knowledge by the user, our findings may induce many researchers to explore recursive preferences and/or SV in further detail. Two final advantages of perturbation are that, often, the perturbed solution provides insights about the economics of the problem and that it might be an excellent initial guess for VFI or for Chebyshev polynomials.

Third, Chebyshev polynomials provide a terrific level of accuracy with reasonable computational burden. When accuracy is most required and the dimensionality of the state space is not too high, as in our model, they are the obvious choice.

Fourth, we were disappointed by the poor performance of VFI, which, compared with Chebyshev, could not achieve a high accuracy even with a large grid. This suggests that we should relegate VFI to solving those problems where non-differentiabilities complicate the application of the previous methods.

The rest of the paper is organized as follows. In section 2, we present our test model. Section 3 describes the different solution methods used to approximate the decision rules of the model. Section 4 discusses the calibration of the model. Section 5 reports numerical results and section 6 concludes. An appendix provides some additional details.

2 The Stochastic Neoclassical Growth Model with Recursive Preferences and SV

We use the stochastic neoclassical growth model with recursive preferences and SV in the process for technology as our test case. We select this model for three reasons. First, it is the workhorse of modern macroeconomics. Even more complicated New Keynesian models with real and nominal rigidities, such as those in Woodford (2003) or Christiano et al. (2005), are built around the core of the neoclassical growth model. Thus, any lesson learned with it is likely to have a wide applicability. Second, the model is, except for the form of the utility function and the process for SV, the same test case as in Aruoba et al. (2006). This provides us with a set of results to compare to our findings. Three, the introduction of recursive preferences and SV make the model both more non-linear (and hence, a challenge for different solution algorithms) and potentially more relevant for practical use. For example, and as mentioned in the introduction, Bansal and Yaron (2004) have emphasized the importance of the combination of recursive preferences and time-varying volatility to account for asset prices.

The description of the model is straightforward, and we just go through the details required to fix notation. There is a representative household that has preferences over streams of consumption,  c_{t}, and leisure,  1-l_{t} , represented by a recursive function of the form:

\displaystyle U_{t}=\max_{c_{t},l_{t}}\left[ \left( 1-\beta \right) \left( c_{t}^{\upsilon }\left( 1-l_{t}\right) ^{1-\upsilon }\right) ^{\frac{1-\gamma }{\theta } }+\beta \left( {\mathbb{E}}_{t}U_{t+1}^{1-\gamma }\right) ^{\frac{1}{\theta } }\right] ^{\frac{\theta }{1-\gamma }} (1)

The parameters in these preferences include  \beta , the discount factor,  \upsilon , which controls labor supply,  \gamma , which controls risk aversion, and:
\displaystyle \theta =\frac{1-\gamma }{1-\frac{1}{\psi }}    

where  \psi is the EIS. The parameter  \theta is an index of the deviation with respect to the benchmark CRRA utility function (when  \theta =1, we are back in that CRRA case where the inverse of the EIS and risk aversion coincide).

The household's budget constraint is given by:

\displaystyle c_{t}+i_{t}+\frac{b_{t+1}}{R_{t}^{f}}=w_{t}l_{t}+r_{t}k_{t}+b_{t}    

where  i_{t} is investment,  R_{t}^{f} is the risk-free gross interest rate,  b_{t} is the holding of an uncontingent bond that pays 1 unit of consumption good at time  t+1,  w_{t} is the wage,  l_{t} is labor,  r_{t} is the rental rate of capital, and  k_{t} is capital. Asset markets are complete and we could have also included in the budget constraint the whole set of Arrow securities. Since we have a representative household, this is not necessary because the net supply of any security is zero. Households accumulate capital according to the law of motion  k_{t+1}=(1-\delta )k_{t}+i_{t} where  \delta is the depreciation rate.

The final good in the economy is produced by a competitive firm with a Cobb-Douglas technology  y_{t}=e^{z_{t}}k_{t}^{\zeta }l_{t}^{1-\zeta } where  z_{t} is the productivity level that follows:

\displaystyle z_{t}=\lambda z_{t-1}+e^{\sigma _{t}}\varepsilon _{t}{\text{, }}\varepsilon _{t}\sim {\mathcal{N}}\left( 0,1\right) .    

Stationarity is the natural choice for our exercise. If we had a deterministic trend, we would only need to adjust  \beta in our calibration below (and the results would be nearly identical). If we had a stochastic trend, we would need to rescale the variables by the productivity level and solve the transformed problem. However, in this case, it is well known that the economy fluctuates less than when  \lambda <1 , and therefore, all solution methods would be closer, limiting our ability to appreciate differences in their performance.

The innovation  \varepsilon _{t} is scaled by a SV level  \sigma _{t}, which evolves as:

\displaystyle \sigma _{t}=(1-\rho )\overline{\sigma }+\rho \sigma _{t-1}+\eta \omega _{t}   , \displaystyle \omega _{t}\sim {\mathcal{N}}\left( 0,1\right)    

where  \overline{\sigma } is the unconditional mean level of  \sigma _{t},  \rho is the persistence of the processes, and  \eta is the standard deviation of the innovations to  \sigma _{t}. Our specification is parsimonious and it introduces only two new parameters,  \rho and  \eta . At the same time, it captures some important features of the data (see a detailed discussion in Fernandez-Villaverde and Rubio-Ramirez, 2010). The combination of an exponent in the productivity process (  e^{\sigma _{t}}) and a level in the evolution of  \sigma _{t} generates interesting nonlinear dynamics. Another important point is that, with SV, we have two innovations, an innovation to technology,  \varepsilon _{t}, and an innovation to the standard deviation of technology,  \omega _{t}. Finally, the economy must satisfy the aggregate resource constraint  y_{t}=c_{t}+i_{t}.

The definition of equilibrium is standard and we skip it in the interest of space. Also, both welfare theorems hold, a fact that we will exploit by jumping back and forth between the solution of the social planner's problem and the competitive equilibrium. However, this is only to simplify our derivations. It is straightforward to adapt the solution methods described below to solve problems that are not Pareto optimal.

Thus, an alternative way to write this economy is to look at the value function representation of the social planner's problem in terms of its three state variables, capital  k_{t}, productivity  z_{t}, and volatility,  \sigma _{t}:

\displaystyle V\left( k_{t},z_{t},\sigma _{t}\right) =\max_{c_{t},l_{t}}\left[ \left( 1-\beta \right) \left( c_{t}^{\upsilon }\left( 1-l_{t}\right) ^{1-\upsilon }\right) ^{\frac{1-\gamma }{\theta }}+\beta \left( {\mathbb{E}} _{t}V^{1-\gamma }\left( k_{t+1},z_{t+1},\sigma _{t+1}\right) \right) ^{\frac{ 1}{\theta }}\right] ^{\frac{\theta }{1-\gamma }}    
s.t. \displaystyle c_{t}+k_{t+1}=e^{z_{t}}k_{t}^{\zeta }l_{t}^{1-\zeta }+\left( 1-\delta \right) k_{t}    
\displaystyle z_{t}=\lambda z_{t-1}+e^{\sigma _{t}}\varepsilon _{t}{\text{, }}\varepsilon _{t}\sim {\mathcal{N}}\left( 0,1\right)    
\displaystyle \sigma _{t}=(1-\rho )\overline{\sigma }+\rho \sigma _{t-1}+\eta \omega _{t}   , \displaystyle \omega _{t}\sim {\mathcal{N}}\left( 0,1\right) .    

Then, we can find the pricing kernel of the economy

\displaystyle m_{t+1}=\frac{\partial V_{t}/\partial c_{t+1}}{\partial V_{t}/\partial c_{t}} .    

Now, note that:
\displaystyle \frac{\partial V_{t}}{\partial c_{t}}=\left( 1-\beta \right) V_{t}^{1-\frac{ 1-\gamma }{\theta }}\upsilon \frac{(c_{t}^{\upsilon }(1-l_{t})^{1-\upsilon })^{\frac{1-\gamma }{\theta }}}{c_{t}}    

and:
\displaystyle \frac{\partial V_{t}}{\partial c_{t+1}}=\beta V_{t}^{1-\frac{1-\gamma }{ \theta }}({\mathbb{E}}_{t}V_{t+1}^{1-\gamma })^{\frac{1}{\theta }-1}{\mathbb{ E}}_{t}\left( V_{t+1}^{-\gamma }\left( 1-\beta \right) V_{t+1}^{1-\frac{ 1-\gamma }{\theta }}\upsilon \frac{\left( 1-\beta \right) (c_{t+1}^{\upsilon }(1-l_{t+1})^{1-\upsilon })^{\frac{1-\gamma }{\theta }}}{c_{t+1}}\right)    

where in the last step we use the result regarding  \partial V_{t}/\partial c_{t} forwarded by one period. Cancelling redundant terms, we get:
\displaystyle m_{t+1}=\frac{\partial V_{t}/\partial c_{t+1}}{\partial V_{t}/\partial c_{t}} =\beta \left( \frac{c_{t+1}}{c_{t}}\right) ^{\frac{\upsilon (1-\gamma )}{ \theta }-1}\left( \frac{1-l_{t+1}}{1-l_{t}}\right) ^{\frac{(1-\gamma )(1-\upsilon )}{\theta }}\left( \frac{V_{t+1}^{1-\gamma }}{{\mathbb{E}} _{t}V_{t+1}^{1-\gamma }}\right) ^{1-\frac{1}{\theta }}. (2)

This equation shows how the pricing kernel is affected by the presence of recursive preferences. If  \theta =1, the last term,
\displaystyle \left( \frac{V_{t+1}^{1-\gamma }}{{\mathbb{E}}_{t}V_{t+1}^{1-\gamma }} \right) ^{1-\frac{1}{\theta }} (3)

is equal to 1 and we get back the pricing kernel of the standard CRRA case. If  \theta \neq 1, the pricing kernel is twisted by (3).

We identify the net return on equity with the marginal net return on investment:

\displaystyle R_{t+1}^{k}=\zeta e^{z_{t+1}}k_{t+1}^{\zeta -1}l_{t+1}^{1-\zeta }-\delta    

with expected return  {\mathbb{E}}_{t}\left[ R_{t+1}^{k}\right] .

3 Solution Methods

We are interested in comparing different solution methods to approximate the dynamics of the previous model. Since the literature on computational methods is large, it would be cumbersome to review every proposed method. Instead, we select those methods that we find most promising.

Our first method is perturbation (introduced by Judd and Guu, 1992 and 1997 and nicely explained in Schmitt-Grohe and Uribe, 2004). Perturbation algorithms build a Taylor series expansion of the agents' decision rules. Often, perturbation methods are very fast and, despite their local nature, highly accurate in a large range of values of the state variables (Aruoba et al., 2006). This means that, in practice, perturbations are the only method that can handle models with dozens of state variables within any reasonable amount of time. Moreover, perturbation often provides insights into the structure of the solution and on the economics of the model. Finally, linearization and log-linearization, the most common solution methods for DSGE models, are particular cases of a perturbation of first order.

We implement a second- and a third-order perturbation of our model. A first-order perturbation is useless for our investigation: the resulting decision rules are certainty equivalent and, therefore, they depend on  \psi but not on  \gamma or  \sigma _{t}. In other words, the first-order decision rules of the model with recursive preferences coincide with the decision rules of the model with CRRA preferences with the same  \psi and  \overline{\sigma } for any value of  \gamma or  \sigma _{t}. We need to go, at least, to second-order decision rules to have terms that depend on  \gamma or  \sigma _{t} and, hence, allow recursive preferences or SV to play a role. Since the accuracy of second-order decision rules may not be high enough and, in addition, we want to explore time-varying risk premia, we also compute a third-order perturbation. As we will document below, a third-order perturbation provides enough accuracy without unnecessary complications. Thus, we do not need to go to higher orders.

The second method is a projection algorithm with Chebyshev polynomials (Judd, 1992). Projection algorithms build approximated decision rules that minimize a residual function that measures the distance between the left- and right-hand side of the equilibrium conditions of the model. Projection methods are attractive because they offer a global solution over the whole range of the state space. Their main drawback is that they suffer from an acute curse of dimensionality that makes it challenging to extend it to models with many state variables. Among the many different types of projection methods, Aruoba et al. (2006) show that Chebyshev polynomials are particularly efficient. Other projection methods, such as finite elements or parameterized expectations, tend to perform somewhat worse than Chebyshev polynomials, and therefore, in the interest of space, we do not consider them.

Finally, we compute the model using VFI (Epstein and Zin, 1989, show that a version of the contraction mapping theorem holds in the value function of the problem with recursive preferences). VFI is slow and it suffers as well from the curse of dimensionality, but it is safe, reliable, and well understood. Thus, it is a natural default algorithm for the solution of DSGE models.

3.1 Perturbation

We describe now each of the different methods in more detail. We start by explaining how to use a perturbation approach to solve DSGE models using the value function of the household. We are not the first to explore the perturbation of value function problems. Judd (1998) already presents the idea of perturbing the value function instead of the equilibrium conditions of a model. Unfortunately, he does not elaborate much on the topic. Schmitt-Grohe and Uribe (2005) employ a perturbation approach to find a second-order approximation to the value function that allows them to rank different fiscal and monetary policies in terms of welfare. However, we follow Binsbergen et al. (2009) in their emphasis on the generality of the approach.1

To illustrate the procedure, we limit our exposition to deriving the second-order approximation to the value function and the decision rules of the agents. Higher-order terms are derived analogously, but the algebra becomes too cumbersome to be developed explicitly (in our application, the symbolic algebra is undertaken by Mathematica, which automatically generates Fortran 95 code that we can evaluate numerically). Hopefully, our steps will be enough to allow the reader to understand the main thrust of the procedure and obtain higher-order approximations by herself.

First, we rewrite the exogenous processes in terms of a perturbation parameter  \chi ,

\displaystyle z_{t} \displaystyle = \displaystyle \lambda z_{t-1}+\chi e^{\sigma _{t}}\varepsilon _{t}  
\displaystyle \sigma _{t} \displaystyle = \displaystyle (1-\rho )\overline{\sigma }+\rho \sigma _{t-1}+\chi \eta \omega _{t}.  

When  \chi =1, which is just a normalization, we are dealing with the stochastic version of the model. When  \chi =0, we are dealing with the deterministic case with steady state  k_{ss},  z_{ss}=0, and  \sigma _{ss}= \overline{\sigma }. Also, it is convenient for the algebra below to define a vector of states in differences with respect to the steady state,  s_{t}=\left( k_{t}-k_{ss},z_{t-1},\varepsilon _{t},\sigma _{t-1}-\sigma _{ss},\omega _{t},\chi \right) , where  s_{it} is the  i-th component of this vector at time  t for  i\in \left\{ 1,\ldots ,6\right\} . Then, we can write the social planner's value function,  V\left( s_{t}\right) , and the decision rules for consumption,  c\left( s_{t}\right) , investment,  i\left( s_{t}\right) , capital,  k\left( s_{t}\right) , and labor,  l\left( s_{t}\right) , as a function of  s_{t}.

Second, we note that, under differentiability assumptions, the second-order Taylor approximation of the value function around  s_{t}=\mathbf{0} (the vectorial zero) is:

\displaystyle V\left( s_{t}\right) \simeq V_{ss}+V_{i,ss}s_{t}^{i}+\frac{1}{2} V_{ij,ss}s_{t}^{i}s_{t}^{j}    

where:
  1. Each term  V_{...,ss} is a scalar equal to a derivative of the value function evaluated at  \mathbf{0}:  V_{ss}\equiv V\left( \mathbf{0}\right) ,  V_{i,ss}\equiv V_{i}\left( \mathbf{0}\right) for  i\in \left\{ 1,\ldots ,6\right\} , and  V_{ij,ss}\equiv V_{ij}\left( \mathbf{0}\right) for  i,j\in \left\{ 1,\ldots ,6\right\} ,
  2. We use the tensors  V_{i,ss}s_{t}^{i}=\sum_{i=1}^{6}V_{i,ss}s_{i,t} and  V_{ij,ss}s_{t}^{i}s_{t}^{j}=\sum_{i=1}^{6} \sum_{i=1}^{6}V_{ij,ss}s_{i,t}s_{j,t}, which eliminate the symbol  \sum_{i=1}^{6} when no confusion arises.

We can extend this notation to higher-order derivatives of the value function. This expansion could also be performed around a different point of the state space, such as the mode of the ergodic distribution of the state variables. In section 5, we discuss this point further.

Fernandez-Villaverde et al. (2010) show that many of these terms  V_{...,ss} are zero (for instance, those implied by certainty equivalence in the first-order component). More directly related to this paper, Binsbergen et al. (2009) demonstrate that  \gamma does not affect the values of any of the coefficients except  V_{66,ss} and also that  V_{66,ss}\neq 0. This result is intuitive, since the value function of a risk-averse agent is in general affected by uncertainty and we want to have an approximation with terms that capture this effect and allow for the appropriate welfare ranking of decision rules. Indeed,  V_{66,ss} has a straightforward interpretation. At the deterministic steady state with  \chi =1 (that is, even if we are in the stochastic economy, we just happen to be exactly at the steady state values of all the other states), we have:

\displaystyle V\left( 0,0,0,0,0,1\right) \simeq V_{ss}+\frac{1}{2}V_{66,ss}    

Hence  \frac{1}{2}V_{66,ss} is a measure of the welfare cost of the business cycle, that is, of how much utility changes when the variance of the productivity shocks is at steady-state value  \sigma _{ss} instead of zero (note that this quantity is not necessarily negative). This term is an accurate evaluation of the third order of the welfare cost of business cycle fluctuations because all of the third-order terms in the approximation of the value function either have coefficient values of zero or drop when evaluated at the deterministic steady state.

This cost of the business cycle can easily be transformed into consumption equivalent units. We can compute the percentage decrease in consumption  \tau that will make the household indifferent between consuming  \left( 1-\tau \right) c_{ss} units per period with certainty or  c_{t} units with uncertainty. To do so, note that the steady-state value function is just  V_{ss}=c_{ss}^{\upsilon }\left( 1-l_{ss}\right) ^{1-\upsilon }, which implies that:

\displaystyle c_{ss}^{\upsilon }\left( 1-l_{ss}\right) ^{1-\upsilon }+\frac{1}{2} V_{66,ss}=\left( \left( 1-\tau \right) c_{ss}\right) ^{\upsilon }\left( 1-l_{ss}\right) ^{1-\upsilon }    

or:
\displaystyle V_{ss}+\frac{1}{2}V_{66,ss}=\left( 1-\tau \right) ^{\upsilon }V_{ss}    

Then:
\displaystyle \tau =1-\left[ 1+\frac{1}{2}\frac{V_{66,ss}}{V_{ss}}\right] ^{\frac{1}{ \upsilon }}.    

We are perturbing the value function in levels of the variables. However, there is nothing special about levels and we could have done the same in logs (a common practice when linearizing DSGE models) or in any other function of the states. These changes of variables may improve the performance of perturbation (Fernández-Villaverde and Rubio-Ramí rez, 2006). By doing the perturbation in levels, we are picking the most conservative case for perturbation. Since one of the conclusions that we will reach from our numerical results is that perturbation works surprisingly well in terms of accuracy, that conclusion will only be reinforced by an appropriate change of variables.2

The decision rules can be expanded in the same way. For example, the second-order approximation of the decision rule for consumption is, under differentiability assumptions:

\displaystyle c\left( s_{t}\right) \simeq c_{ss}+c_{i,ss}s_{t}^{i}+\frac{1}{2} c_{ij,ss}s_{t}^{i}s_{t}^{j}    

where we have followed the same derivative and tensor notation as before.

As with the approximation of the value function, Binsbergen et al. (2009) show that  \gamma does not affect the values of any of the coefficients except  c_{66,ss}. This term is a constant that captures precautionary behavior caused by risk. This observation tells us two facts. First, a linear approximation to the decision rule does not depend on  \gamma (it is certainty equivalent), and therefore, if we are interested in recursive preferences, we need to go at least to a second-order approximation. Second, given some fixed parameter values, the difference between the second-order approximation to the decision rules of a model with CRRA preferences and a model with recursive preferences is a constant. This constant generates a second, indirect effect, because it changes the ergodic distribution of the state variables and, hence, the points where we evaluate the decision rules along the equilibrium path. These arguments demonstrate how perturbation methods can provide analytic insights beyond computational advantages and help in understanding the numerical results in Tallarini (2000). In the third-order approximation, all of the terms on functions of  \chi ^{2} depend on  \gamma .

Following the same steps, we can derive the decision rules for labor, investment, and capital. In addition we have functions that give us the evolution of other variables of interest, such as the pricing kernel or the risk-free gross interest rate  R_{t}^{f}. All of these functions have the same structure and properties regarding  \gamma as the decision rule for consumption. In the case of functions pricing assets, the second-order approximation generates a constant risk premium, while the third-order approximation creates a time-varying risk premium.

Once we have reached this point, there are two paths we can follow to solve for the coefficients of the perturbation. The first procedure is to write down the equilibrium conditions of the model plus the definition of the value function. Then, we take successive derivatives in this augmented set of equilibrium conditions and solve for the unknown coefficients. This approach, which we call equilibrium conditions perturbation (ECP), gets us, after  n iterations, the  n-th-order approximation to the value function and to the decision rules.

A second procedure is to take derivatives of the value function with respect to states and controls and use those derivatives to find the unknown coefficient. This approach, which we call value function perturbation (VFP), delivers after  \left( n+1\right) steps, the  \left( n+1\right) -th order approximation to the value function and the  n-th order approximation to the decision rules.3 Loosely speaking, ECP undertakes the first step of VFP by hand by forcing the researcher to derive the equilibrium conditions.

The ECP approach is simpler but it relies on our ability to find equilibrium conditions that do not depend on derivatives of the value function. Otherwise, we need to modify the equilibrium conditions to include the definitions of the derivatives of the value function. Even if this is possible to do (and not particularly difficult), it amounts to solving a problem that is equivalent to VFP. This observation is important because it is easy to postulate models that have equilibrium conditions where we cannot get rid of all the derivatives of the value function (for example, in problems of optimal policy design). ECP is also faster from a computational perspective. However, VFP is only trivially more involved because finding the  \left( n+1\right) -th-order approximation to the value function on top of the  n-th order approximation requires nearly no additional effort.

The algorithm presented here is based on the system of equilibrium equations derived using the ECP. In the appendix, we derive a system of equations using the VFP. We take the first-order conditions of the social planner. First, with respect to consumption today:

\displaystyle \frac{\partial V_{t}}{\partial c_{t}}-\mu _{t}=0    

where  \mu _{t} is the Lagrangian multiplier associated with the resource constraint. Second, with respect to capital:
\displaystyle -\mu _{t}+\mathbb{E}_{t}\mu _{t+1}\left( \zeta e^{z_{t+1}}k_{t+1}^{\zeta -1}l_{t+1}^{1-\zeta }+1-\delta \right) =0.    

Third, with respect to labor:
\displaystyle \frac{1-\upsilon }{\upsilon }\frac{c_{t}}{(1-l_{t})}=(1-\zeta )e^{z_{t}}k_{t}^{\zeta }l_{t}^{-\zeta }.    

Then, we have  \mathbb{E}_{t}m_{t+1}\left( \zeta e^{z_{t}}k_{t+1}^{\zeta -1}l_{t+1}^{1-\zeta }+1-\delta \right) =1 where  m_{t+1} was derived in equation (2). Note that, as explained above, the derivatives of the value function in (2) are eliminated.

Once we substitute for the pricing kernel, the augmented equilibrium conditions are:

\displaystyle V_{t}-\left[ \left( 1-\beta \right) \left( c_{t}^{\upsilon }\left( 1-l_{t}\right) ^{1-\upsilon }\right) ^{\frac{1-\gamma }{\theta }}+\beta \left( {\mathbb{E}}_{t}V^{1-\gamma }\left( k_{t+1},z_{t+1}\right) \right) ^{ \frac{1}{\theta }}\right] ^{\frac{\theta }{1-\gamma }}=0    
\displaystyle \mathbb{E}_{t}\left[ \beta \left( \frac{c_{t+1}}{c_{t}}\right) ^{\frac{ 1-\gamma }{\theta }-1}\left( \frac{V_{t+1}^{1-\gamma }}{\mathbb{E} _{t}V_{t+1}^{1-\gamma }}\right) ^{1-\frac{1}{\theta }}\left( \zeta e^{z_{t+1}}k_{t+1}^{\zeta -1}l_{t+1}^{1-\zeta }+1-\delta \right) \right] -1=0    
\displaystyle \frac{1-\upsilon }{\upsilon }\frac{c_{t}}{(1-l_{t})}=(1-\zeta )e^{z_{t}}k_{t}^{\zeta }l_{t}^{-\zeta }=0    
\displaystyle \mathbb{E}_{t}\beta \left( \frac{c_{t+1}}{c_{t}}\right) ^{\frac{1-\gamma }{ \theta }-1}\left( \frac{V_{t+1}^{1-\gamma }}{\mathbb{E}_{t}V_{t+1}^{1-\gamma }}\right) ^{1-\frac{1}{\theta }}R_{t}^{f}-1=0    
\displaystyle c_{t}+i_{t}-e^{z_{t}}k_{t}^{\zeta }l_{t}^{1-\zeta }=0    
\displaystyle k_{t+1}-i_{t}-\left( 1-\delta \right) k_{t}=0    

plus the law of motion for productivity and volatility. Note that all the endogenous variables are functions of the states and that we drop the max operator in front of the value function because we are already evaluating it at the optimum. Thus, a more compact notation for the previous equilibrium conditions as a function of the states is:
\displaystyle F\left( \mathbf{0}\right) =\mathbf{0}    

where  F:\mathbb{R}^{6}\rightarrow \mathbb{R}^{8}.

In steady state,  m_{ss}=\beta and the set of equilibrium conditions simplifies to:

\displaystyle V_{ss}=c_{ss}^{\upsilon }\left( 1-l_{ss}\right) ^{1-\upsilon }    
\displaystyle \left( \zeta k_{ss}^{\zeta -1}l_{ss}^{1-\zeta }+1-\delta \right) =1/\beta    
\displaystyle \frac{1-\upsilon }{\upsilon }\frac{c_{ss}}{(1-l_{ss})}=(1-\zeta )k_{ss}^{\zeta }l_{ss}^{-\zeta }    
\displaystyle R_{ss}^{f}=1/\beta    
\displaystyle c_{ss}+i_{ss}=k_{ss}^{\zeta }l_{ss}^{1-\zeta }    
\displaystyle i_{ss}=\delta k_{ss}    

a system of 6 equations on 6 unknowns,  V_{ss},  c_{ss},  k_{ss},  i_{ss} ,  l_{ss}, and  R_{ss}^{f} that can be easily solved (see the appendix for the derivations). This steady state is identical to the steady state of the real business cycle model with a standard CRRA utility function and no SV.

To find the first-order approximation to the value function and the decision rules, we take first derivatives of the function  F with respect to the states  s_{t} and evaluate them at  \mathbf{0}:

\displaystyle F_{i}\left( \mathbf{0}\right) =\mathbf{0} for \displaystyle i\in \left\{ 1,\ldots ,6\right\} .    

This step gives us 48 different first derivatives (8 equilibrium conditions times the 6 variables of  F). Since each dimension of  F is equal to zero for all possible values of  s_{t}, their derivatives must also be equal to zero. Therefore, once we substitute the steady-state values and forget about the exogenous processes (which we do not need to solve for), we have a quadratic system of 36 equations on 36 unknowns:  V_{i,ss},  c_{i,ss},  i_{i,ss},  k_{i,ss},  l_{i,ss}, and  R_{i,ss}^{f} for  i\in \left\{ 1,\ldots ,6\right\} . One of the solutions is an unstable root of the system that violates the transversality condition of the problem and we eliminate it. Thus, we keep the solution that implies stability.

To find the second-order approximation, we take derivatives on the first derivatives of the function  F, again with respect to the states and the perturbation parameter:

\displaystyle F_{ij}\left( \mathbf{0}\right) =\mathbf{0} for \displaystyle i,j\in \left\{ 1,\ldots ,6\right\} .    

This step gives us a new system of equations. Then, we plug in the terms that we already know from the steady state and from the first-order approximation and we get that the only unknowns left are the second-order terms of the value function and other functions of interest. Quite conveniently, this system of equations is linear and it can be solved quickly. Repeating these steps (taking higher-order derivatives, plugging in the terms already known, and solving for the remaining unknowns), we can get any arbitrary order approximation. For simplicity, and since we checked that we were already obtaining a high accuracy, we decided to stop at a third-order approximation (we are particularly interested in applying the perturbation for estimation purposes and we want to document how a third-order approximation is accurate enough for many problems without spending too much time deriving higher-order terms).

3.2 Projection

Projection methods take basis functions to build an approximated value function and decision rules that minimize a residual function defined by the augmented equilibrium conditions of the model. There are two popular methods for choosing basis functions: finite elements and the spectral method. We will present only the spectral method for several reasons: first, in the neoclassical growth model the decision rules and value function are smooth and spectral methods provide an excellent approximation. Second, spectral methods allow us to use a large number of basis functions, with the consequent high accuracy. Third, spectral methods are easier to implement. Their main drawback is that since they approximate the solution with a spectral basis, if the decision rules display a rapidly changing local behavior or kinks, it may be difficult for this scheme to capture those local properties.

Our target is to solve the decision rule for labor and the value function  \{l_{t},V_{t}\} from:

\displaystyle \mathcal{H}(l_{t},V_{t})=\left[ \begin{array}{c} u_{c,t}-\beta \left( \mathbb{E}_{t}V_{t+1}^{1-\gamma }\right) ^{\frac{1}{ \theta }-1}\mathbb{E}_{t}\left[ V_{t+1}^{\frac{\left( 1-\gamma \right) \left( \theta -1\right) }{\theta }}u_{c,t+1}\left( \zeta e^{z_{t+1}}k_{t+1}^{\zeta -1}l_{t+1}^{1-\zeta }+1-\delta \right) \right] \\ V_{t}-\left[ (1-\beta )(c_{t}^{\upsilon }(1-l_{t}^{\upsilon }))^{\frac{ 1-\gamma }{\theta }}+\beta \mathbb{E}_{t}(V_{t+1}^{1-\gamma })^{\frac{1}{ \theta }}\right] ^{\frac{\theta }{1-\gamma }} \end{array} \right] =\mathbf{0}    

where, to save on notation, we define  V_{t}=V\left( k_{t},z_{t},\sigma _{t}\right) and:
\displaystyle u_{c,t}=\frac{1-\gamma }{\theta }\upsilon \frac{\left( c_{t}^{\upsilon }\left( 1-l_{t}\right) ^{1-\upsilon }\right) ^{\frac{1-\gamma }{\theta }}}{ c_{t}}    

Then, from the static condition
\displaystyle c_{t}=\frac{\upsilon }{1-\upsilon }(1-\zeta )e^{z_{t}}k_{t}^{\zeta }l_{t}^{-\zeta }(1-l_{t})    

and the resource constraint, we can find  c_{t} and  k_{t+1}.

Spectral methods solve this problem by specifying the decision rule for labor and the value function  \{l_{t},V_{t}\} as linear combinations of weighted basis functions:

\displaystyle l(k_{t},z_{j},\sigma _{m};\rho ) \displaystyle = \displaystyle \Sigma _{i}\rho _{ijm}^{l}\psi _{i}(k_{t})  
\displaystyle V(k_{t},z_{j},\sigma _{m};\rho ) \displaystyle = \displaystyle \Sigma _{i}\rho _{ijm}^{V}\psi _{i}(k_{t})  

where  \{\psi _{i}(k)\}_{i=1,...,n_{k}} are the  n_{k} basis functions that we will use for our approximation along the capital dimension and  \rho =\{\rho _{ijm}^{l},\rho _{ijm}^{V}\}_{i=1,...,n_{k};j=1,...,J;m=1,..,M} are unknown coefficients to be determined. In this expression, we have discretized the stochastic processes  \sigma _{t} for volatility and  z_{t} for productivity using Tauchen's (1986) method as follows. First, we have a grid of  M points  G_{\sigma }=\left\{ e^{\sigma _{1}},e^{\sigma 2},...,e^{\sigma _{M}}\right\} for  \sigma _{t} and a transition matrix  \Pi ^{M} with generic element  \pi _{i,j}^{M}=Prob\left( e^{\sigma _{t+1}}=e^{\sigma _{j}}\vert e^{\sigma _{t}}=e^{\sigma _{i}}\right) . The grid covers 3 standard deviations of the process in each direction. Then, for each  M point, we find a grid with  J points  G_{z}^{m}=\left\{ z_{1}^{m},z_{2}^{m},...,z_{J}^{m}\right\} for  z_{t} and transition matrixes  \Pi ^{J,m} with generic element  \pi _{i,j}^{J,m}=Prob\left( z_{t+1}^{m}=z_{j}^{m}\vert z_{t}^{m}=z_{i}^{m}\right) . Again, and conditional on  e^{\sigma _{m}}, the grid covers 3 standard deviations in each direction. Values for the decision rule outside the grids  G_{\sigma } and  G_{z}^{m} can be approximated by interpolation. We make the grids for  z_{t} depend on the level of volatility  m to adapt the accuracy of Tauchen's procedure to each conditional variance (although this forces us to interpolate when we switch variances). Since we set  J=25 and  M=5, the approximation is quite accurate along the productivity axis (we explored other choices of  J and  M to be sure that our choice was sensible).

A common choice for the basis functions are Chebyshev polynomials because of their flexibility and convenience. Since their domain is [-1,1], we need to bound capital to the set [  \underline{k},\overline{k}], where  \underline{k} (  \overline{k}) is chosen sufficiently low (high) to bind only with extremely low probability, and define a linear map from those bounds into [-1,1]. Then, we set  \psi _{i}(k_{t})=\widetilde{\psi }_{i}(\phi _{k}(k_{t})) where  \widetilde{\psi }_{i}(\cdot ) are Chebyshev polynomials and  \phi _{k}(k_{t}) is the linear mapping from [  \underline{k},\overline{k}] to [-1,1].

By plugging  l(k_{t},z_{j},\sigma _{m};\rho ) and  V(k_{t},z_{j},\sigma _{m};\rho ) into  \mathcal{H}(l_{t},V_{t}), we find the residual function:

\displaystyle \mathcal{R}(k_{t},z_{j},\sigma _{m};\rho )=\mathcal{H}(l(k_{t},z_{j},\sigma _{m};\rho ),V(k_{t},z_{j},\sigma _{m};\rho ))    

We determine the coefficients  \rho to get the residual function as close to  \mathbf{0} as possible. However, to do so, we need to choose a weight of the residual function over the space  \left( k_{t},z_{j},\sigma _{m}\right) . A collocation point criterion delivers the best trade-off between speed and accuracy (Fornberg, 1998) by making the residual function exactly equal to zero in  \{k_{i}\}_{i=1}^{n_{k}} roots of the  n_{k}-th order Chebyshev polynomial and in the Tauchen points (also, by the Chebyshev interpolation theorem, if an approximating function is exact at the roots of the  n_{k}-th order Chebyshev polynomial, then, as  n_{k}\rightarrow \infty , the approximation error becomes arbitrarily small). Therefore, we just need to solve the following system of  n_{k}\times J\times M\times 2 equations:
\displaystyle \mathcal{R}(k_{i},z_{j},\sigma _{m};\rho )=\mathbf{0} for any \displaystyle i,j,m    collocation points    

on  n_{k}\times J\times M\times 2 unknowns  \rho . We solve this system with a Newton method and an iteration based on the increment of the number of basis functions. First, we solve a system with only three collocation points for capital and 125 (  125=25\times 5) points for technology. Then, we use that solution as a guess for a system with more collocation points for capital (with the new coefficients being guessed to be equal to zero) and iterate on the procedure. We stop the iteration when we have 11 polynomials in the capital dimension (therefore, in the last step we solve for  2,750=11\times 25\times 5\times 2 coefficients). The iteration is needed because otherwise the residual function is too cumbersome to allow for direct solution of the 2,750 final coefficients.

3.3 Value Function Iteration

Our final solution method is VFI. Since the dynamic algorithm is well known, our presentation is most brief. Consider the following Bellman operator:

\displaystyle TV\left( k_{t},z_{t},\sigma _{t}\right) =\max_{c_{t},l_{t}}\left[ \left( 1-\beta \right) \left( c_{t}^{\upsilon }\left( 1-l_{t}\right) ^{1-\upsilon }\right) ^{\frac{1-\gamma }{\theta }}+\beta \left( {\mathbb{E}} _{t}V^{1-\gamma }\left( k_{t+1},z_{t+1},\sigma _{t+1}\right) \right) ^{\frac{ 1}{\theta }}\right] ^{\frac{\theta }{1-\gamma }}    
s.t. \displaystyle c_{t}+k_{t+1}=e^{z_{t}}k_{t}^{\zeta }l_{t}^{1-\zeta }+\left( 1-\delta \right) k_{t}    
\displaystyle z_{t}=\lambda z_{t-1}+e^{\sigma _{t}}\varepsilon _{t}{\text{, }}\varepsilon _{t}\sim {\mathcal{N}}\left( 0,1\right)    
\displaystyle \sigma _{t}=(1-\rho )\overline{\sigma }+\rho \sigma _{t-1}+\eta \omega _{t}   , \displaystyle \omega _{t}\sim {\mathcal{N}}\left( 0,1\right) .    

To solve for this Bellman operator, we define a grid on capital,  G_{k}=\left\{ k_{1},k_{2},...,k_{M}\right\} , a grid on volatility and on the productivity level. The grid on capital is just a uniform distribution of points over the capital dimension. As we did for projection, we set a grid  G_{\sigma }=\left\{ e^{\sigma _{1}},e^{\sigma 2},...,e^{\sigma _{M}}\right\} for  \sigma _{t} and a transition matrix  \Pi ^{M} for volatility and  M grids  G_{z}^{m}=\left\{ z_{1}^{m},z_{2}^{m},...,z_{J}^{m}\right\} for  z_{t} and transition matrixes  \Pi ^{J,m} using Tauchen's (1986) procedure. The algorithm to iterate on the value function for this grid is:
  1. Set  n=0 and  V^{0}(k_{t},z_{t},\sigma _{t})=c_{ss}^{\upsilon }\left( 1-l_{ss}\right) ^{1-\upsilon } for all  k_{t}\in G_{k} and all  z_{t}\in G_{z}.
  2. For every  \{k_{t},z_{t},\sigma _{t}\}, use the Newton method to find  c_{t}^{\ast },  l_{t}^{\ast },  k_{t+1}^{\ast } that solve:
    \displaystyle c_{t}=\frac{\upsilon }{1-\upsilon }(1-\zeta )e^{z_{t}}k_{t}^{\zeta }l_{t}^{-\zeta }(1-l_{t})    
    \displaystyle \left( 1-\beta \right) \upsilon \frac{\left( c_{t}^{\upsilon }\left( 1-l_{t}\right) ^{1-\upsilon }\right) ^{\frac{1-\gamma }{\theta }}}{c_{t}} =\beta \left( \mathbb{E}_{t}\left( V_{t+1}^{n}\right) ^{1-\gamma }\right) ^{ \frac{1}{\theta }-1}\mathbb{E}_{t}\left[ \left( V_{t+1}^{n}\right) ^{-\gamma }V_{1,t+1}^{n}\right]    
    \displaystyle c_{t}+k_{t+1}=e^{z_{t}}k_{t}^{\zeta }l_{t}^{1-\zeta }+(1-\delta )k_{t}    

  3. Construct  V^{n+1} from the Bellman equation:
    \displaystyle V^{n+1}=((1-\beta )(c_{t}^{\ast \upsilon }(1-l_{t}^{\ast })^{1-\upsilon })^{ \frac{1-\gamma }{\theta }}+\beta (\mathbb{E}_{t}(V(k_{t+1}^{\ast },z_{t+1},\sigma _{t+1})^{1-\gamma }))^{\frac{1}{\theta }})^{\frac{\theta }{ 1-\gamma }}    

  4. If  \frac{\vert V^{n+1}-V^{n}\vert}{\vert V^{n}\vert}\geq 1.0e^{-7}, then  n=n+1 and go to 2. Otherwise, stop.

To accelerate convergence and give VFI a fair chance, we implement a multigrid scheme as described by Chow and Tsitsiklis (1991). We start by iterating on a small grid. Then, after convergence, we add more points to the grid and recompute the Bellman operator using the previously found value function as an initial guess (with linear interpolation to fill the unknown values in the new grid points). Since the previous value function is an excellent grid, we quickly converge in the new grid. Repeating these steps several times, we move from an initial 23,000-point grid into a final one with 375,000 points (3,000 points for capital, 25 for productivity, and 5 for volatility).

4 Calibration

We now select a benchmark calibration for our numerical computations. We follow the literature as closely as possible and select parameter values to match, in the steady state, some basic observations of the U.S. economy (as we will see below, for the benchmark calibration, the means of the ergodic distribution and the steady-state values are nearly identical). We set the discount factor  \beta =0.991 to generate an annual interest rate of around 3.6 percent. We set the parameter that governs labor supply,  \theta  =0.357 , to get the representative household to work one-third of its time. The elasticity of output to capital,  \zeta  =0.3, matches the labor share of national income. A value of the depreciation rate  \delta =0.0196 matches the ratio of investment-output. Finally,  \lambda =0.95 and  \log \overline{\sigma }=0.007 are standard values for the stochastic properties of the Solow residual. For the SV process, we pick  \rho =0.9 and  \eta =0.06, to match the persistence and standard deviation of the heteroskedastic component of the Solow residual during the last 5 decades.

[Table 1 here]

Since we want to explore the dynamics of the model for a range of values that encompasses all the estimates from the literature, we select four values for the parameter that controls risk aversion,  \gamma , 2, 5, 10, and 40, and two values for EIS  \psi , 0.5, and 1.5, which bracket most of the values used in the literature (although many authors prefer smaller values for  \psi , we found that the simulation results for smaller  \psi 's do not change much from the case when  \psi  =0.5). We then compute the model for all eight combinations of values of  \gamma and  \psi , that is  \left\{ 2,0.5\right\} ,  \left\{ 5,0.5\right\} ,  \left\{ 10,0.5\right\} , and so on. When  \psi  =0.5 and  \gamma =2, we are back in the standard CRRA case. However, in the interest of space, we will report only a limited subset of results that we find are the most interesting ones.

We pick as the benchmark case the calibration  \{\gamma ,\psi ,\log \overline{\sigma },\eta \}=\{5,0.5,0.007,0.06\}. These values reflect an EIS centered around the median of the estimates in the literature, a reasonably high level of risk aversion, and the observed volatility of productivity shocks. To check robustness, we increase, in the extreme case, the risk aversion, the average standard deviation of the productivity shock, and the standard deviation of the innovations to volatility to  \{\gamma ,\psi ,\log \overline{\sigma },\eta \}=\{40,0.5,0.021,0.1\}. This case combines levels of risk aversion that are in the upper bound of all estimates in the literature with huge productivity shocks. Therefore, it pushes all solution methods to their limits, in particular, making life hard for perturbation since the interaction of the large precautionary behavior induced by  \gamma and large shocks will move the economy far away from the deterministic steady state. We leave the discussion of the effects of  \psi =1.5 for the robustness analysis at the end of the next section.

5 Numerical Results

In this section we report our numerical findings. First, we present and discuss the computed decision rules. Second, we show the results of simulating the model. Third, we report the Euler equation errors. Fourth, we discuss the effects of changing the EIS and the perturbation point. Finally, we discuss implementation and computing time.

5.1 Decision Rules

Figure 1 plots the decision rules of the household for consumption and labor in the benchmark case. In the top two panels, we plot the decision rules along the capital axis when  z_{t}=0 and  \sigma _{t}=\overline{\sigma }. The capital interval is centered on the steady-state level of capital of 9.54 with a width of  \pm 40\%, [5.72,13.36]. This interval is big enough to encompass all the simulations. We also plot, with two vertical lines, the 10 and 90 percentile of the ergodic distribution of capital.4 Since all methods provide nearly indistinguishable answers, we observe only one line in both panels. It is possible to appreciate miniscule differences in labor supply between second-order perturbation and the other methods only when capital is far from its steady-state level. Monotonicity of the decision rules is preserved by all methods.5 We must be cautious, however, mapping differences in choices into differences in utility. The Euler error function below provides a better view of the welfare consequences of different approximations.

In the bottom two panels, we plot the decision rules as a function of  \sigma _{t} when  k_{t}=k_{ss} and  z_{t}=0. Both VFI and Chebyshev polynomials capture well precautionary behavior: the household consumes less and work more when  \sigma _{t} is higher. In comparison, it would seem that perturbations do a much poorer job at capturing this behavior. However, this is not the right interpretation. In a second-order perturbation, the effect of  \sigma _{t} only appears in a term where  \sigma _{t} interacts with  z_{t}. Since here we are plotting a cut of the decision rule while keeping  z_{t}=0, that term drops out and the decision rule that comes from a second-order approximation is a flat line. In the third-order perturbation,  \sigma _{t} appears by itself, but raised to its cube. Hence, the slope of the decision rule is negligible. As we will see below, the interaction effects of  \sigma _{t} with other variables, which are hard to visualize in two-dimensional graph, are all that we need to deliver a satisfactory overall performance.

[Figure 1 here]

We see bigger differences in the decision rules as we increase the risk aversion and the variance of innovations. Figure 2 plots the decision rules for the extreme calibration following the same convention than before. In this figure, we change the capital interval where we compute the top decision rules to [3,32] (roughly 1/3 and 3 times the steady-state capital) because, owing to the high variance of the calibration, the equilibrium paths fluctuate through much wider ranges of capital.

[Figure 2 here]

We highlight several results. First, all the methods deliver similar results in our original capital interval for the benchmark calibration. Second, as we go far away from the steady state, VFI and the Chebyshev polynomial still overlap with each other (and, as shown by our Euler error computations below, we can roughly take them as the "exact" solution), but second- and third-order approximations start to deviate. Third, the decision rule for consumption approximated by the third-order perturbation changes from concavity into convexity for values of capital bigger than 15. This phenomenon (also documented in Aruoba et al. 2006) is due to the poor performance of local approximation when we move too far away from the expansion point and the polynomials begin to behave wildly. Numerically, this issue is irrelevant because the problematic region is visited with nearly zero probability.

5.2 Simulations

Applied economists often characterize the behavior of the model through statistics from simulated paths of the economy. We simulate the model, starting from the deterministic steady state, for 10,000 periods, using the decision rules for each of the eight combinations of risk aversion and EIS discussed above. To make the comparison meaningful, the shocks are common across all paths. We discard the first 1,000 periods as a burn-in to eliminate the transition from the deterministic steady state of the model to the middle regions of the ergodic distribution of capital. This is usually achieved in many fewer periods than the ones in our burn-in, but we want to be conservative in our results. The remaining observations constitute a sample from the ergodic distribution of the economy.

For the benchmark calibration, the simulations from all of the solution methods generate almost identical equilibrium paths (and therefore we do not report them). We focus instead on the densities of the endogenous variables as shown in figure 3 (remember that volatility and productivity are identical across the different solution methods). Given the low risk aversion and SV of the productivity shocks, all densities are roughly centered around the deterministic steady-state value of the variable. For example, the mean of the distribution of capital is only 0.2 percent higher than the deterministic value. Also, capital is nearly always between 8.5 and 10.5. This range will be important below to judge the accuracy of our approximations.

[Figure 3 here]

Table 2 reports business cycle statistics and, because DSGE models with recursive preferences and SV are often used for asset pricing, the average and variance of the (quarterly) risk-free rate and return on capital. Again, we see that nearly all values are the same, a simple consequence of the similarity of the decision rules.

[Table 2 here]

The welfare cost of the business cycle is reported in table 3 in consumption equivalent terms. The computed costs are actually negative. Besides the Jensen's effect on average productivity, this is also due to the fact that when we have leisure in the utility function, the indirect utility function may be convex in input prices (agents change their behavior over time by a large amount to take advantage of changing productivity). Cho and Cooley (2000) present a similar example. Welfare costs are comparable across methods. Remember that the welfare cost of the business cycle for the second- and third-order perturbations is the same because the third-order terms all drop or are zero when evaluated at the steady state.

[Table 3 here]

When we move to the extreme calibration, we see more differences. Figure 4 plots the histograms of the simulated series for each solution method. Looking at quantities, the histograms of consumption, output, and labor are the same across all of the methods. The ergodic distribution of capital puts nearly all the mass between values of 6 and 15. This considerable move to the right in comparison with figure 3 is due to the effect of precautionary behavior in the presence of high risk aversion, large productivity shocks, and high SV. Capital also visits low values of capital more than in the benchmark calibration because of large, persistent productivity shocks. In any case, the translation is more pronounced to the right than to the left.

[Figure 4 here]

Table 4 reports business cycle statistics. Differences across methods are minor in terms of means (note that the mean of the risk-free rate is lower than in the benchmark calibration because of the extra accumulation of capital induced by precautionary behavior). In terms of variances, the second-order perturbation produces less volatility than all other methods. This suggests that a second-order perturbation may not be good enough if we face high variance of the shocks and/or high risk aversion. A third-order perturbation, in comparison, eliminates most of the differences and delivers nearly the same implications as Chebyshev polynomials or VFI.

[Table 4 here]

Finally, table 5 presents the welfare cost of the business cycle. Now, in comparison with the benchmark calibration, the welfare cost of the business cycle is positive and significant, slightly above 1.1 percent. This is not a surprise, since we have both a large risk aversion and productivity shocks with an average standard deviation three times as big as the observed one. All methods deliver numbers that are close.

[Table 5 here]

5.3 Euler Equation Errors

While the plots of the decision rules and the computation of densities and business cycle statistics that we presented in the previous subsection are highly informative, it is also important to evaluate the accuracy of each of the procedures. Euler equation errors, introduced by Judd (1992), have become a common tool for determining the quality of the solution method. The idea is to observe that, in our model, the intertemporal condition:

\displaystyle u_{c,t}=\beta (\mathbb{E}_{t}V_{t+1}^{1-\gamma })^{\frac{1}{\theta }-1} \mathbb{E}_{t}\left( V_{t+1}^{\frac{(\gamma -1)(1-\theta )}{\theta } }u_{c,t+1}R\left( k_{t},z_{t},\sigma _{t};z_{t+1},\sigma _{t+1}\right) \right) (4)

where  R\left( k_{t},z_{t},\sigma _{t};z_{t+1},\sigma _{t+1}\right) =1+\zeta e^{z_{t+1}}k_{t+1}^{\zeta -1}l_{t+1}^{1-\zeta }-\delta is the gross return of capital given states  k_{t},  z_{t},  \sigma _{t}, and realizations  z_{t+1} and  \sigma _{t+1} should hold exactly for any given  k_{t}, and  z_{t}. However, since the solution methods we use are only approximations, there will be an error in (4) when we plug in the computed decision rules. This Euler equation error function  EE^{i}\left( k_{t},z_{t},\sigma _{t}\right) is defined, in consumption terms:
\displaystyle EE^{i}\left( k_{t},z_{t},\sigma _{t}\right) =1-\frac{\left[ \frac{\beta ({ \mathbb{E}}_{t}\left( V_{t+1}^{i}\right) ^{1-\gamma })^{\frac{1}{\theta }-1}{ \mathbb{E}}_{t}\left( \left( V_{t+1}^{i}\right) ^{\frac{(\gamma -1)(1-\theta )}{\theta }}u_{c,t+1}^{i}R\left( k_{t},z_{t},\sigma _{t};z_{t+1},\sigma _{t+1}\right) \right) }{\frac{1-\gamma }{\theta }\upsilon (1-l_{t}^{i})^{\left( 1-\upsilon \right) \frac{1-\gamma }{\theta }}}\right] ^{\frac{1}{\upsilon \frac{1-\gamma }{\theta }-1}}}{c_{t}^{i}}    

This function determines the (unit free) error in the Euler equation as a fraction of the consumption given the current states and solution method  i . Following Judd and Guu (1997), we can interpret this error as the optimization error incurred by the use of the approximated decision rule and we report the absolute errors in base 10 logarithms to ease interpretation. Thus, a value of -3 means a $1 mistake for each $1000 spent, a value of -4 a $1 mistake for each $10,000 spent, and so on.
[Figure 5 here]

Figure 5 displays a transversal cut of the errors for the benchmark calibration when  z=0 and  \sigma _{t}=\overline{\sigma }. Other transversal cuts at different technology and volatility levels reveal similar patterns. The first lesson from figure 5 is that all methods deliver high accuracy. We know from figure 3 that capital is nearly always between 8.5 and 10.5. In that range, the (log10) Euler equation errors are at most -5, and most of the time they are even smaller. For instance, the second- and third-order perturbations have an Euler equation error of around -7 in the neighborhood of the deterministic steady state, VFI of around -6.5, and Chebyshev an impressive -11/-13. The second lesson from figure 5 is that, as expected, global methods (Chebyshev and VFI) perform very well in the whole range of capital values, while perturbations deteriorate as we move away from the steady state. For second-order perturbation, the Euler error in the steady state is almost four orders of magnitude smaller than on the boundaries. Third-order perturbation is around half an order of magnitude more accurate than second-order perturbation over the whole range of values (except in a small region close to the deterministic steady state).

There are two complementary ways to summarize the information from Euler equation error functions. First, we report the maximum error in our interval (capital between 60 percent and 140 percent of the steady state and the grids for productivity and volatility). The maximum Euler error is useful because it bounds the mistake owing to the approximation. The second procedure for summarizing Euler equation errors is to integrate the function with respect to the ergodic distribution of capital and productivity to find the average error. We can think of this exercise as a generalization of the Den Haan-Marcet test (Den Haan and Marcet, 1994). The top- left panel in Figure 6 reports the maximum Euler error (darker bars) and the integral of the Euler error for the benchmark calibration. Both perturbations have a maximum Euler error of around -2.7, VFI of -3.1, and Chebyshev, an impressive -9.8. We read this result as indicating that all methods perform adequately. Both perturbations have roughly the same integral of the Euler error (around -5.3), VFI a slightly better -6.4, while Chebyshev polynomials do fantastically well at -10.4 (the average loss of welfare is $1 for each $500 billion). But even an approximation with an average error of $1 for each $200,000, such as the one implied by third-order perturbation, must suffice for most relevant applications.

[Figure 6 here]

We repeat our exercise for the extreme calibration. As we did when we computed the decision rules of the agents, we have changed the capital interval to [3,32]. The top-right panel in Figure 6 reports maximum Euler equation errors and their integrals. The maximum Euler equation error is large for perturbation methods while it is rather small using Chebyshev polynomials. However, given the very large range of capital used in the computation, this maximum Euler error provides a too negative view of accuracy. We find the integral of the Euler equation error to be more instructive. With a second-order perturbation, we have -4.02 and with a third-order perturbation we have -4.12. To evaluate this number, remember that we have extremely high risk aversion and large productivity shocks. Even in this challenging environment, perturbations deliver a high degree of accuracy. VFI does not display a big loss of precision compared to the benchmark case. On the other hand, Chebyshev polynomials deteriorate somewhat, but the accuracy it delivers it is still of $1 out of each $1 million spent.

5.4 Robustness: Changing the EIS and Changing the Perturbation Point

In the results we reported above, we kept the EIS equal to 0.5, a conventional value in the literature, while we modified the risk aversion and the volatility of productivity shocks. However, since some researchers prefer higher values of the EIS (see, for instance, Bansal and Yaron, 2004, a paper that we have used to motivate our investigation), we also computed our model with  \psi  =1.5. Basically our results were unchanged. To save on space, we concentrate only on the Euler equation errors (decision rules and simulation paths are available upon request). In the bottom-left panel in Figure 6, we report the maxima of the Euler equation errors and their integrals with respect to the ergodic distribution. The relative size and values of the entries in this table are quite similar to the values reported for the benchmark calibration (except, partially, VFI that performs a bit better). The bottom-right panel in Figure 6 repeats the same exercise for the extreme calibration. Again, the entries in the table are very close to the ones in the extreme calibration (and now, VFI does not perform better than when  \psi  =0.5).

As a final robustness test, we computed the perturbations not around the deterministic steady state (as we did in the main text), but around a point close to the mode of the ergodic distribution of capital. This strategy, if perhaps difficult to implement because of the need to compute the mode of the ergodic distribution,6 could deliver better accuracy because we approximate the value function and decision rules in a region where the model spends more time. As we suspected, we found only trivial improvements in terms of accuracy. Moreover, expanding at a point different from the deterministic steady state has the disadvantage that the theorems that ensure the convergence of the Taylor approximation might fail (see theorem 6 in Jin and Judd, 2002).

5.5 Implementation and Computing Time

We briefly discuss implementation and computing time. For the benchmark calibration, second-order perturbation and third- order perturbation algorithms take only 0.02 second and 0.05 second, respectively, in a 3.3GHz Intel PC with Windows 7 (the reference computer for all times below), and it is simple to implement: 664 lines of code in Fortran 95 for second order and 1133 lines of code for third order, plus in both cases, the analytical derivatives of the equilibrium conditions that Fortran 95 borrows from a code written in Mathematica 6.7 The code that generates the analytic derivatives has between 150 to 210 lines, although Mathematica is much less verbose. While the number of lines doubles in the third order, the complexity in terms of coding does not increase much: the extra lines are mainly from declaring external functions and reading and assigning values to the perturbation coefficients. An interesting observation is that we only need to take the analytic derivatives once, since they are expressed in terms of parameters and not in terms of parameter values. This allows Fortran to evaluate the analytic derivatives extremely fast for new combinations of parameter values. This advantage of perturbation is particularly relevant when we need to solve the model repeatedly for many different parameter values, for example, when we are estimating the model. For completeness, the second-order perturbation was also run in Dynare (although we had to use version 4.0, which computes analytic derivatives, instead of previous versions, which use numerical derivatives that are not accurate enough for perturbation). This run was a double-check of the code and a test of the feasibility of using off-the-shelf software to solve DSGE models with recursive preferences.

The projection algorithm takes around 300 seconds, but it requires a good initial guess for the solution of the system of equations. Finding the initial guess for some combination of parameter values proved to be challenging. The code is 652 lines of Fortran 95. Finally, the VFI code is 707 lines of Fortran 95, but it takes about ten hours to run.

6 Conclusions

In this paper, we have compared different solution methods for DSGE models with recursive preferences and SV. We evaluated the different algorithms based on accuracy, speed, and programming burden. We learned that all of the most promising methods (perturbation, projection, and VFI) do a fair job in terms of accuracy. We were surprised by how well simple second-order and third-order perturbations perform even for fairly non-linear problems. We were impressed by how accurate Chebyshev polynomials can be. However, their computational cost was higher and we are concerned about the curse of dimensionality. In any case, it seems clear to us that, when accuracy is the key consideration, Chebyshev polynomials are the way to go. Finally, we were disappointed by VFI since even with 125,000 points in the grid, it only did marginally better than perturbation and it performed much worse than Chebyshev polynomials in our benchmark calibration. This suggests that unless there are compelling reasons such as non-differentiabilities or non-convexities in the model, we better avoid VFI.

A theme we have not developed in this paper is the possibility of interplay among different solution methods. For instance, we can compute extremely easily a second-order approximation to the value function and use it as an initial guess for VFI. This second-order approximation is such a good guess that VFI will converge in few iterations. We verified this idea in non-reported experiments, where VFI took one-tenth of the time to converge once we used the second-order approximation to the value function as the initial guess. This approach may even work when the true value function is not differentiable at some points or has jumps, since the only goal of perturbation is to provide a good starting point, not a theoretically sound approximation. This algorithm may be particularly useful in problems with many state variables. More research on this type of hybrid method is a natural extension of our work.

We close the paper by pointing out that recursive preferences are only one example of a large class of non-standard preferences that have received much attention by theorists and applied researchers over the last several years (see Backus, Routledge, and Zin, 2004). Having fast and reliable solution methods for this class of new preferences will help researchers to sort out which of these preferences deserve further attention and to derive empirical implications. Thus, this paper is a first step in the task of learning how to compute DSGE models with non-standard preferences.

Table 1: Calibrated Parameters
Parameter  \beta  \upsilon  \zeta  \delta  \lambda  \log \overline{\sigma }  \rho  \eta
Value 0.991 0.357 0.3 0.0196 0.95 0.007 0.9 0.06

Table 2a: Business Cycle Statistics(Mean) - Benchmark Calibration
Mean  c  y  i  R^{f}(\%)  R^{k}(\%)
Second-Order Perturbation 0.7253 0.9128 0.1873 0.9070 0.9078
Third-Order Perturbation 0.7257 0.9133 0.1875 0.9062 0.9069
Chebyshev Polynomial 0.7256 0.9130 0.1875 0.9063 0.9066
Value Function Iteration 0.7256 0.9130 0.1875 0.9063 0.9066

Table 2b:Business Cycle Statistics(Variance (%)) - Benchmark Calibration
Variance  c  y  i  R^{f}(\%)  R^{k}(\%)
Second-Order Perturbation 0.0331 0.1084 0.0293 0.0001 0.0001
Third-Order Perturbation 0.0330 0.1079 0.0288 0.0001 0.0001
Chebyshev Polynomial 0.0347 0.1117 0.0313 0.0001 0.0001
Value Function Iteration 0.0347 0.1117 0.0313 0.0001 0.0001

Table 3: Welfare Costs of Business Cycle - Benchmark Calibration
2nd-Order Pert. 3rd-Order Pert. Chebyshev Value Function
-2.0864e(-5) -2.0864e(-5) -3.2849e(-5) -3.2849e(-5)

Table 4a: Business Cycle Statistics(Mean) - Extreme Calibration
Mean  c  y  i  R^{f}(\%)  R^{k}(\%)
Second-Order Perturbation 0.7338 0.9297 0.1950 0.8432 0.8562
Third-Order Perturbation 0.7344 0.9311 0.1955 0.8416 0.8529
Chebyshev Polynomial 0.7359 0.9329 0.1970 0.8331 0.8402
Value Function Iteration 0.7359 0.9329 0.1970 0.8352 0.8403

Table 4b: Business Cycle Statistics(Variance (%)) - Extreme Calibration
Variance  c  y  i  R^{f}(\%)  R^{k}(\%)
Second-Order Perturbation 0.2956 1.0575 0.2718 0.0004 0.0004
Third-Order Perturbation 0.3634 1.2178 0.3113 0.0004 0.0005
Chebyshev Polynomial 0.3413 1.1523 0.3425 0.0005 0.0006
Value Function Iteration 0.3414 1.1528 0.3427 0.0005 0.0006

Table 5: Welfare Costs of Business Cycle - Extreme Calibration
2nd-Order Pert. 3rd-Order Pert. Chebyshev Value Function
1.1278e-2 1.1278e-2 1.2855e-2 1.2838e-2
Figure 1: Decision Rules - Benchmark Calibration.
Figure 1: Decision rules of the household for benchmark calibration. Four panels. The figure plots the decision rules for households computed using four different solution algorithms: second and third order perturbation, Chebyshev polynomials, and value function iteration. Top-left panel: Decision rule for consumption as function of the capital stock, for steady state level of technology and variance of technology. Data plotted as curves. X axis displays capital stock, Y axis consumption level. This panel shows that consumption is increasing in the capital stock, and that all solution methods deliver identical decision rules. Top-right panel: Decision rule for labor as function of the capital stock, for steady state level of technology and variance of technology. Data plotted as curves. X axis displays capital stock, Y axis hours worked as fraction of total time. This panel shows that hours worked are decreasing in the capital stock, and that all solution methods deliver identical decision rules. Bottom-left panel: Decision rule for consumption as function of the volatility of technology, for steady state level of technology and capital stock. Data plotted as curves. X axis displays volatility levels, Y axis consumption level. This panel shows that consumption decreases with the level of volatility, although for the low levels associated to the benchmark calibration the slope is nearly zero. Second and third order perturbation produce a decision rule for consumption with a higher slope than Chebyshev polynomials and value function iteration, but numerical differences are negligible. Bottom-right panel: Decision rule for labor as function of the volatility of technology, for steady state level of technology and capital stock. Data plotted as curves. X axis displays volatility levels, Y axis hours worked as share of total time. This panel shows that labor increases with the level of volatility, although for the low levels associated to the benchmark calibration the slope is nearly zero. Second and third order perturbation produce a decision rule for labor with a lower slope than Chebyshev polynomials and value function iteration, but numerical differences are negligible.
Figure 2: Decision Rules - Extreme Calibration.
Figure 2: Decision rules of the household for extreme calibration. Four panels. The figure plots the decision rules for households computed using four different solution algorithms: second and third order perturbation, Chebyshev polynomials, and value function iteration. Top-left panel: Decision rule for consumption as function of the capital stock, for steady state level of technology and variance of technology. Data plotted as curves. X axis displays capital stock, Y axis consumption level. This panel shows that consumption is increasing in the capital stock. Theoretically, consumption is concave in the capital stock. Concave functions are produced by the Chebyshev polynomials, the value function iteration, and the second-order perturbation, although the latter method displays higher concavity (and it is less precise). Concavity is not preserved by the third order perturbation method for values of capital far from the approximation point. Yet, for values of capital within the 10th and 90th percentile of the ergodic distribution, plotted as vertical lines, all solution methods deliver nearly identical consumption decision rules. Top-right panel: Decision rule for labor as function of the capital stock, for steady state level of technology and variance of technology. Data plotted as curves. X axis displays capital stock, Y axis hours worked as fraction of total time. This panel shows that hours worked are decreasing in the capital stock. Theoretically, labor is convex in the capital stock. Convex functions are produced by the Chebyshev polynomials, the value function iteration, and the second-order perturbation, although the latter method displays higher convexity (and it is less precise). Convexity is not preserved by the third order perturbation method for values of capital far from the approximation point. Yet, for values of capital within the 10th and 90th percentile of the ergodic distribution, plotted as vertical lines, all solution methods deliver nearly identical labor decision rules. Bottom-left panel: Decision rule for consumption as function of the volatility of technology, for steady state level of technology and capital stock. Data plotted as curves. X axis displays volatility levels, Y axis consumption level. This panel shows that consumption decreases with the level of volatility, although for the low levels associated to the benchmark calibration the slope is nearly zero. Second and third order perturbation produce a decision rule for consumption with a higher slope than Chebyshev polynomials and value function iteration, but numerical differences are negligible. Bottom-right panel: Decision rule for labor as function of the volatility of technology, for steady state level of technology and capital stock. Data plotted as curves. X axis displays volatility levels, Y axis hours worked as share of total time. This panel shows that labor increases with the level of volatility, although for the low levels associated to the benchmark calibration the slope is nearly zero. Second and third order perturbation produce a decision rule for labor with a lower slope than Chebyshev polynomials and value function iteration, but numerical differences are negligible.
Figure 3: Densities - Benchmark Calibration.
Figure 3: Density functions estimated from simulation of the model under benchmark calibration. Six panels. The figure plots the density functions computed using four different solution algorithms: second and third order perturbation, Chebyshev polynomials, and value function iteration. Data plotted as curves. Top-left panel: Density function for consumption. Data plotted as curves. X axis displays consumption level, Y axis frequency. This panel shows that consumption is normally distributed around the mean of its ergodic distribution, and that all solution methods deliver nearly identical results. Top-right panel: Density function for hours worked. Data plotted as curves. X axis displays hours worked as fraction of total time, Y axis frequency. This panel shows that hours worked are normally distributed around the mean of its ergodic distribution, and that all solution methods deliver nearly identical results.  Middle-left panel: Density function for capital. Data plotted as curves. X axis displays capital level, Y axis frequency. This panel shows that capital is normally distributed around the mean of its ergodic distribution, and that all solution methods deliver nearly identical results. Middle-right panel: Density function for output. Data plotted as curves. X axis displays output level, Y axis frequency. This panel shows that output is normally distributed around the mean of its ergodic distribution, and that all solution methods deliver nearly identical results. Bottom-left panel: Density function for return on risk-free bond. Data plotted as curves. X axis displays capital level, Y axis frequency. This panel shows that the return on the risk-free bond is normally distributed around the mean of its ergodic distribution, and that all solution methods deliver nearly identical results. Bottom-right panel: Density function for the return on equity. Data plotted as curves. X axis displays output level, Y axis frequency. This panel shows that the return on equity is normally distributed around the mean of its ergodic distribution, and that all solution methods deliver nearly identical results.
Figure 4: Densities - Extreme Calibration.
Figure 4: Density functions estimated from simulation of the model under the extreme. Six panels. The figure plots the density functions computed using four different solution algorithms: second and third order perturbation, Chebyshev polynomials, and value function iteration. Data plotted as curves. Top-left panel: Density function for consumption. Data plotted as curves. X axis displays consumption level, Y axis frequency. This panel shows the distribution of consumption around the mean of its ergodic distribution. Small departure from normality are due to the presence of stochastic volatility, with levels of volatility being sufficiently different from each other. Chebyshev polynomials and value function iteration deliver identical results, while second and third order perturbations deliver slightly different distributions, although differences are small. Top-right panel: Density function for hours worked. Data plotted as curves. X axis displays hours worked, Y axis frequency. This panel shows the distribution of hours worked around the mean of its ergodic distribution. Chebyshev polynomials and value function iteration deliver identical results, while second and third order perturbations deliver slightly different distributions, although differences are small. Middle-left panel: Density function for capital. Data plotted as curves. X axis displays capital level, Y axis frequency. This panel shows the distribution of capital around the mean of its ergodic distribution. Departures from normality are due to the presence of stochastic volatility, with levels of volatility being sufficiently different from each other. Chebyshev polynomials and value function iteration deliver identical results, while second and third order perturbations deliver different distributions, both around the mode and at the tails. Middle-right panel: Density function for output. Data plotted as curves. X axis displays output level, Y axis frequency. This panel shows the distribution of output around the mean of its ergodic distribution. All solution methods provide nearly identical distributions. Bottom-left panel: Density function for the return on the risk-free bond. Data plotted as curves. X axis displays return on the risk-free bond, Y axis frequency. This panel shows the distribution of the return on the risk-free bond around the mean of its ergodic distribution. Small departure from normality are due to the presence of stochastic volatility, with levels of volatility being sufficiently different from each other. Chebyshev polynomials and value function iteration deliver identical results, while second and third order perturbations deliver slightly different distributions, although differences are small.. Bottom-right panel: Density function for the return on equity. Data plotted as curves. X axis displays the return on equity, Y axis frequency. This panel shows the distribution of the equity return around the mean of its ergodic distribution. All solution methods provide nearly identical distributions.
Figure 5: Euler Equation Error - Benchmark Calibration.
Figure 5. Euler Equation Error plotted for the benchmark calibration. The figure plots the log10 of the euler equation error computed for different algorithms: second and third order perturbation, Chebyshev polynomials, and value function iteration. Data plotted as curves. X axis displays capital stock, Y axis log10 of the Euler Equation Error.  Figure 5 displays a transversal cut of the errors for the benchmark calibration when both technology and its volatility take steady state values. The first lesson from Figure 5 is that all methods deliver high accuracy. We know from figure 3 that capital is nearly always between 8.5 and 10.5. In that range, the log10 Euler equation errors are at most -5, and most of the time they are even smaller. For instance, the second- and third-order perturbations have an Euler equation error of around -7 in the neighborhood of the deterministic steady state, VFI of around -6.5, and  Chebyshev an impressive -11/-13. The second lesson from Figure 5 is that, as expected, global methods (Chebyshev and VFI) perform very well in the whole range of capital values, while perturbations deteriorate as we move away from the steady state. For second-order perturbation, the Euler error in the steady state is almost four orders of magnitude smaller than on the boundaries. Third-order perturbation is around half an order of magnitude more accurate than second-order perturbation over the whole range of values (except in a small region close to the deterministic steady state).
Figure 6: Maximum Euler Equation Errors.
   Figure 6. Maximum Euler Equation Errors. Data plotted as bars. X axis lists the four different solution alghoritms: second and third order perturbation, Chebyshev polynomials, and value function iteration. Y axis log10 of the Euler Equation Errors. Top- left panel:  Maximum Euler error (darker bars) and the integral of the Euler error (light bars) for the benchmark calibration. Both perturbations have a maximum Euler error of around -2.7, VFI of -3.1, and Chebyshev, an impressive -9.8. We read this result as indicating that all methods perform adequately. Both perturbations have roughly the same integral of the Euler error (around -5.3), value function iteration a slightly better -6.4, while Chebyshev polynomials do fantastically well at -10.4 (the average loss of welfare is $1 for each $500 billion). But even an approximation with an average error of $1 for each $200,000, such as the one implied by third-order perturbation, must suffice for most relevant applications. Top-right panel:  Maximum Euler error (darker bars) and the integral of the Euler error (light bars) for the extreme calibration. The maximum Euler equation error is -1.8 perturbation methods while it is -4.4 using Chebyshev polynomials. However, given the very large range of capital used in the computation, this maximum Euler error provides a too negative view of accuracy. We find the integral of the Euler equation error to be more instructive. With a second-order perturbation, we have -4.02 and with a third-order perturbation we have -4.12. Value function iteration does not display a big loss of precision compared to the benchmark case. On the other hand, Chebyshev polynomials deteriorate somewhat, but the accuracy it delivers it is still of $1 out of each $1 million spent. Bottom left panel: Maximum Euler error (darker bars) and the integral of the Euler error (light bars) for the benchmark calibration assuming an intertemporal elasticity of substitution equal to 1.5 instead of 0.5. The relative size and values of the entries in this table are quite similar to the values reported for the benchmark calibration in the top-left panel. Maximum Euler error (darker bars) and the integral of the Euler error (light bars) for the extreme calibration assuming an intertemporal elasticity of substitution equal to 1.5 instead of 0.5. The relative size and values of the entries in this table are quite similar to the values reported for the extreme calibration in the top-right panel.

7 Appendix

In this appendix, we present the steady state of the model and the alternative perturbation approach, the value function perturbation (VFP).

7.1 Steady State of the Model

To solve the system:

\displaystyle V_{ss}=c_{ss}^{\upsilon }\left( 1-l_{ss}\right) ^{1-\upsilon }    
\displaystyle \left( \zeta k_{ss}^{\zeta -1}l_{ss}^{1-\zeta }+1-\delta \right) =1/\beta    
\displaystyle \frac{1-\nu }{\nu }\frac{c_{ss}}{(1-l_{ss})}=(1-\zeta )k_{ss}^{\zeta }l_{ss}^{-\zeta }    
\displaystyle m_{ss}R_{ss}^{f}=1/\beta    
\displaystyle c_{ss}+i_{ss}=k_{ss}^{\zeta }l_{ss}^{1-\zeta }    
\displaystyle i_{ss}=\delta k_{ss}    

note first that:
\displaystyle \frac{k_{ss}}{l_{ss}}=\left( \frac{1}{\zeta }\left( \frac{1}{\beta } -1+\delta \right) \right) ^{\frac{1}{\zeta -1}}=\Omega    

Now, from the leisure-consumption condition:

\displaystyle \frac{c_{ss}}{1-l_{ss}}=\frac{\upsilon }{1-\upsilon }\left( 1-\zeta \right) \Omega ^{\zeta }=\Phi \Rightarrow c_{ss}=\Phi \left( 1-l_{ss}\right)    

Then:
\displaystyle c_{ss}+\delta k_{ss}=k_{ss}^{\zeta }l_{ss}^{1-\zeta }=\Omega ^{\zeta }l_{ss}\Rightarrow c_{ss}=\left( \Omega ^{\zeta }-\delta \Omega \right) l_{ss}    

and:
\displaystyle \Phi \left( 1-l_{ss}\right) =\left( \Omega ^{\zeta }-\delta \Omega \right) l_{ss}\Rightarrow    
\displaystyle l_{ss}=\frac{\Phi }{\Omega ^{\zeta }-\delta \Omega +\Phi }    
\displaystyle k_{ss}=\frac{\Phi \Omega }{\Omega ^{\zeta }-\delta \Omega +\Phi }    

from which we can find  V_{ss} and  i_{ss}.

7.2 Value Function Perturbation (VFP)

We mentioned in the main text that instead of perturbing the equilibrium conditions of the model, we could directly perturb the value function in what we called value function perturbation (VFP). To undertake the VFP, we write the value function as:

\displaystyle V\left( k_{t},z_{t},\sigma _{t};\chi \right) =\max_{c_{t},l_{t}}\left[ \left( 1-\beta \right) \left( c_{t}^{\upsilon }\left( 1-l_{t}\right) ^{1-\upsilon }\right) ^{\frac{1-\gamma }{\theta }}+\beta {\mathbb{E}} _{t}V^{1-\gamma }\left( k_{t+1},z_{t+1},\sigma _{t+1};\chi \right) ^{\frac{1 }{\theta }}\right] ^{\frac{\theta }{1-\gamma }}    

To find a second-order approximation to the value function, we take derivatives of the value function with respect to controls (  c_{t},l_{t}), states (  k_{t},z_{t},\sigma _{t}), and the perturbation parameter  \chi . We collect these 6 equations, together with the resource constraint, the value function itself, and the exogenous processes in a system:

\displaystyle \widetilde{F}\left( k_{t},z_{t},\chi \right) =\mathbf{0}    

where the hat over  F emphasizes that now we are dealing with a slightly different set of equations than the  F in the main text.

After solving for the steady state of this system, we take derivatives of the function  \widetilde{F} with respect to  k_{t},  z_{t},  \sigma _{t} , and  \chi :

\displaystyle \widetilde{F}_{i}\left( k_{ss},0,\sigma _{ss};0\right) =\mathbf{0} for  \displaystyle i=\left\{ 1,2,3,4\right\}    

and we solve for the unknown coefficients. This solution will give us a second-order approximation of the value function but only a first-order approximation of the decision rules. By repeating these steps  n times, we can obtain the  n+1-order approximation of the value function and the  n -order approximation of the decision rules. It is straightforward to check that the coefficients obtained by ECP and VFP are the same. Thus, the choice of one approach or the other should be dictated by expediency.

Bibliography

Anderson, E.W., L.P. Hansen, E.R. McGrattan, and T.J. Sargent, (1996).
"On the Mechanics of Forming and Estimating Dynamic Linear Economies," in H.M. Amman et al. (eds.) Handbook of Computational Economics, Elsevier.
Aruoba S.B., J. Fernandez-Villaverde, and J. Rubio-Ramirez, (2006).
"Comparing Solution Methods for Dynamic Equilibrium Economies." Journal of Economic Dynamics and Control 30, 2477-2508.
Backus D.K., B.R. Routledge, and S.E. Zin, (2004).
"Exotic Preferences for Macroeconomists." NBER Macroeconomics Annual 2004, 319-390.
Backus D.K., B.R. Routledge, and S.E. Zin, (2007).
"Asset Pricing in Business Cycle Analysis." Mimeo, New York University.
Bansal R. and A. Yaron, (2004).
"Risks for the Long Run: A Potential Resolution of Asset Pricing Puzzles." Journal of Finance 59, 1481-1509.
Benigno P. and M. Woodford, (2006).
"Linear-Quadratic Approximation of Optimal Policy Problems." Mimeo, Columbia University.
Binsbergen J.H. van, J. Fernandez-Villaverde, R.S.J. Koijen, and J. Rubio-Ramirez, (2009).
"Likelihood Estimation of DSGE Models with Epstein-Zin Preferences." Mimeo , University of Pennsylvania.
Cho J-O. and T.C. Cooley, (2000).
"Business Cycle Uncertainty and Economic Welfare." Mimeo, New York University.
Chow C.S., and J.N. Tsitsiklis, (1991).
"An Optimal One-Way Multigrid Algorithm for Discrete-Time Stochastic Control." IEEE Transaction on Automatic Control 36, 898-914.
Christiano L., M. Eichenbaum, and C.L. Evans, (2005).
"Nominal Rigidities and the Dynamic Effects of a Shock to Monetary Policy." Journal of Political Economy 113, 1-45.
Den Haan W. J., and A. Marcet, (1994).
"Accuracy in Simulations." Review of Economic Studies 61, 3-17.
Epstein L., and S.E. Zin, (1989).
"Substitution, Risk Aversion, and the Temporal Behavior of Consumption and Asset Returns: A Theoretical Framework." Econometrica 57, 937-969.
Epstein L., and S.E. Zin, (1991).
"Substitution, Risk Aversion, and the Temporal Behavior of Consumption and Asset Returns: An Empirical Analysis." Journal of Political Economy 99, 263-286.
Fernandez-Villaverde J. and J. Rubio-Ramirez, (2006).
"Solving DSGE Models with Perturbation Methods and a Change of Variables." Journal of Economic Dynamics and Control 30, 2509-2531.
Fernandez-Villaverde J. and J. Rubio-Ramirez, (2010).
"Macroeconomics and Volatility: Data, Models, and Estimation" NBER Working Paper 16618.
Fernandez-Villaverde J., P. Guerron-Quintana, and J. Rubio-Ramirez, (2010).
"Fortune or Virtue: Time-Variant Volatilities Versus Parameter Drifting in U.S. Data." NBER Working Paper 15928.
Fornberg B., (1998).
A Practical Guide to Pseudospectral Methods. Cambridge University Press, Cambridge.
Hansen L.P., and T.J. Sargent, (1995).
"Discounted Linear Exponential Quadratic Gaussian Control." IEEE Transactions on Automatic Control 40, 968-971.
Hansen L.P., J. Heaton, J. Lee, and N. Roussanov, (2007).
"Intertemporal Substitution and Risk Aversion." In James J. Heckman and Edward E. Leamer (eds.), Handbook of Econometrics 6, 3967-4056.
Jin H. and K.L. Judd, (2002).
"Perturbation Methods for General Dynamic Stochastic Models." Mimeo, Hoover Institution.
Judd K.L., (1992).
"Projection Methods for Solving Aggregate Growth Models." Journal of Economic Theory 58, 410-452.
Judd K.L., (1998).
Numerical Methods in Economics. MIT Press.
Judd K.L. and S.M. Guu, (1992).
"Perturbation Solution Methods for Economic Growth Models." In H. Varian (ed.), Economic and Financial Modelling with Mathematica. Springer Verlag.
Judd K.L. and S.M. Guu, (1997).
"Asymptotic Methods for Aggregate Growth Models." Journal of Economic Dynamics and Control 21, 1025-1042.
Kreps D.M. and E.L. Porteus, (1978).
"Temporal Resolution of Uncertainty and Dynamic Choice Theory." Econometrica 46, 185-200.
Krueger D. and F. Kubler, (2005).
"Pareto Improving Social Security When Financial Markets Are Incomplete." American Economic Review 96, 737-755.
Levine P., J. Pearlman, and R. Pierse, (2007).
"Linear-Quadratic Approximation, External Habit, and Targeting Rules." Mimeo, University of Surrey.
Schmitt-Grohe S. and M. Uribe, (2004).
"Solving Dynamic General Equilibrium Models Using a Second-Order Approximation to the Policy Function." Journal of Economic Dynamics and Control 28, 755-775.
Schmitt-Grohe S. and M. Uribe, (2005).
"Optimal Fiscal and Monetary Policy in a Medium Scale Macroeconomic Model." NBER Macroeconomic Annual 2005, 382-425.
Smets F. and R. Wouters, (2007).
"Shocks and Frictions in US Business Cycles: A Bayesian DSGE Approach." American Economic Review 97, 586-606.
Tallarini. T.D., (2000).
"Risk-Sensitive Real Business Cycles." Journal of Monetary Economics 45, 507-532.
Tauchen G., (1986).
"Finite State Markov-Chain Approximations to Univariate and Vector Autoregressions." Economics Letters 20, 177-181.
Weil P., (1990).
"Nonexpected Utility in Macroeconomics." Quarterly Journal of Economics 105, 29-42.
Woodford M.D., (2003).
Interest and Prices. Princeton University Press.



Footnotes

* We thank Michel Juillard for his help with computational issues and Larry Christiano, Dirk Krueger, Pawel Zabczyk, and participants at several seminars for comments. Beyond the usual disclaimer, we must note that any views expressed herein are those of the authors and not necessarily those of the Board of Governors of the Federal Reserve System or the Federal Reserve Bank of Atlanta. Finally, we also thank the NSF for financial support. Return to Text
† Federal Reserve Board, Washington, DC 20551, [email protected]. Return to Text
‡ University of Pennsylvania, NBER, CEPR, and FEDEA. Address: 160 McNeil, 3718 Locust Walk, Philadelphia, PA 19104, and [email protected]. Return to Text
♠ Duke University, Federal Reserve Bank of Atlanta, and FEDEA. Address: 213 Social Sciences Building, Box 90097, Durham, NC 27708-0097, and [email protected]. Return to Text
♦ University of Pennsylvania. Address: 160 McNeil, 3718 Locust Walk, Philadelphia, PA 19104, and [email protected]. Return to Text
1. The perturbation method is related to Benigno and Woodford (2006) and Hansen and Sargent (1995). Benigno and Woodford present a linear-quadratic (LQ) approximation to solve optimal policy problems when the constraints of the problem are non-linear (see also Levine et al., 2007). This allows them to find the correct local welfare ranking of different policies. Our perturbation can also deal with non-linear constraints and obtains the correct local approximation to welfare and policies, but with the advantage that it is easily generalizable to higher-order approximations. Hansen and Sargent (1995) modify the LQ problem to adjust for risk. In that way, they can handle some versions of recursive utilities. Hansen and Sargent's method, however, requires imposing a tight functional form for future utility and to surrender the assumption that risk-adjusted utility is separable across states of the world. Perturbation does not suffer from those limitations. Return to Text
2. This comment begets the question, nevertheless, of why we did not perform a perturbation in logs, since many economists will find it more natural than using levels. Our experience with the CRRA utility case (Aruoba et al., 2006) and some computations with recursive preferences not included in the paper suggest that a perturbation in logs does slightly worse than a perturbation in levels. Return to Text
3. The classical strategy of finding a quadratic approximation of the utility function to derive a linear decision rule is a second-order example of VFP (Anderson et al.,1996). A standard linearization of the equilibrium conditions of a DSGE model when we add the value function to those equilibrium conditions is a simple case of ECP. This is done, for instance, in Schmitt-Grohe and Uribe (2005). Return to Text
4. There is the technical consideration of which ergodic distribution to use for this task, since this is an object that can only be found by simulation. All over the paper, we use the ergodic simulation generated by VFI. We checked that the results are robust to using the ergodic distributions coming from the other methods. Return to Text
5. Similar figures could be plotted for other values of  z_{t} and  \sigma _{t}. We omit them because of space considerations. Return to Text
6. For example, the algorithm of finding a perturbation around the steady state, simulate from it, find a second perturbation around the model of the implied ergodic simulation, and so on until convergence, may not settle in any fixed point. In our exercise, we avoid this problem because we have the ergodic distribution implied by VFI. This is an unfair advantage for perturbations at the mode of the ergodic distribution but it makes our point below about the lack of improvement in accuracy even stronger. Return to Text
7. We use lines of code as a proxy for the complexity of implementation. We do not count comment lines. Return to Text

This version is optimized for use by screen readers. Descriptions for all mathematical expressions are provided in LaTex format. A printable pdf version is available. Return to Text