Computational Methods for Production-Based Asset Pricing Models with Recursive Utility

2.1 Preferences

Our economy admits a representative agent whose utility function follows Epstein and Zin (1989) and Weil (1990):

where 0 < β < 1 is the subjective discount factor, 𝔼_t is the conditional expectations operator, C_t denotes aggregate consumption, C¯t=(Ct,Ct+1,…), γ denotes the agent’s coefficient of relative risk aversion, ψ denotes the agent’s inter-temporal elasticity of substitution (IES), and θ=1−γ1−1/ψ. A particularly desirable feature of this utility function is that it separates the IES and risk-aversion parameters, as opposed to the standard constant relative risk aversion (CRRA) utility, where IES and risk aversion are inversely related. In theory, it is not clear that there should be a tight link between these two parameters, as risk aversion is atemporal while IES is temporal.

As shown in Epstein and Zin (1989) and Weil (1989), the log of the stochastic discount factor, m_t+1, for these preferences is

where Δc_t+1 denotes log consumption growth and r_a,t+1 denotes the log gross return on the aggregate wealth portfolio. Of particular importance is the presence of r_a,t+1 in the specification of the discount factor, which makes innovations to expected consumption growth a priced risk factor; in the standard model with CRRA utility, θ = 1 and the last term disappears. It is this feature of Epstein-Zin-Weil utility that makes agents concerned about shocks to expected future consumption growth and that allows us to amplify the equity premium (Mehra and Prescott, 1985) while evading the risk-free rate puzzle (Weil (1989)).

2.2 Technology

A single firm owns the capital stock and produces a consumption good via Cobb-Douglas technology, using labor and capital as inputs:

where Z_t is the labor-augmenting stochastic technology level, H_t denotes the number of hours worked, K_t represents capital and α is the share of capital in the production function. Under the assumption that H¯ is the agent’s total leisure endowment, it is clear that utility is going to be maximized when Ht=H¯, ∀t, since H_t does not appear in the utility function. Normalizing H¯=1, the production function simplifies to

The log technology process, zt=ln⁡(Zt), evolves exogenously according to

We limit ourselves to the special case of φ = 1, since this allows us to retain only one state variable; our computational results are similar for the case of persistent, yet trend stationary z_t (when |φ| < 1). Hence, the log of TFP is a random walk with drift parameter μ, and in this special case shocks to technology are permanent.

2.3 Capital accumulation

Following Jermann (1998), we allow capital adjustment costs in the accumulation equation

where

is an increasing, concave function which induces large changes in the capital stock to be more costly than successive small changes. The parameter ξ governs the degree of concavity and has the desirable feature that as ξ → ∞, ϕ(x) becomes the identity function (with an appropriate specification of α₁ and α₂); that is, capital adjustment costs disappear. At the other extreme, as ξ → 0, I_t → 0, ∀t, and ϕ(x)≡exp⁡(μ)−1+δ, allowing us to obtain an endowment economy where all output is consumed each period and the capital stock grows deterministically at rate exp(μ). Hence, for intermediate values, the adjustment cost parameter, ξ, allows us flexibility in matching the relative volatilities of consumption and output. As mentioned above, the remaining parameters are defined so as to eliminate adjustment costs in the deterministic steady state: α1=(exp⁡(μ)−1+δ)1/ξ and α2=11−ξ(exp⁡(μ)−1+δ) (see B for a derivation).

2.4 Equilibrium

In equilibrium, the aggregate resource constraint is binding:

In this basic environment, the welfare theorems are satisfied and the solution to the social planner’s problem yields the same allocations as a competitive equilibrium. The planner’s problem is

subject to

By solving the first-order conditions, we obtain the intertemporal Euler equation

where

is an alternative expression for the Epstein-Zin-Weil stochastic discount factor, equivalent to exp⁡(mt+1) in Equation (2). As shown in C, the return on equity is

which is the term scaling the stochastic discount factor in Equation (13). Hence, the Euler equation can be written compactly as Et[Mt+1Rt+1E]=1. This expression for the return to equity will be useful when evaluating the quality of model solutions.

To preserve stationarity in the economy, we normalize all variables by the level of the contemporaneous technology process:

Since the preferences represented by Equation (1) are homothetic (Epstein and Zin, 1989), the normalized system of equilibrium conditions can then be expressed as

where

Augmenting Equations (17)–(20) with appropriately scaled versions of Equations (14), (3), and (9) yields the system of equations that we solve with the methods of the following section.

3.1 Perturbation

Perturbation methods, suggested for economic models by Judd and Guu (1997) and Judd (1998), and widely popularized by Schmitt-Grohé and Uribe (2004), build an asymptotically valid polynomial approximation of a function around a point where the solution is known. In general notation, perturbation seeks a local approximation to a function, F, where

and where F(x(0), 0) is known. The typical specialization in economics is for F to represent a system of nonlinear stochastic difference equations,

where the deterministic steady state, Et[f(xt+1(0),xt(0),0)]=f(xss,xss,0) is known. The canonical economic example is the neoclassical growth model, where f is a system of equations including the inter-temporal Euler equation and constraints, and where the polynomial approximation to f is a Taylor expansion. However, there is no a priori reason to restrict our attention to the inter-temporal Euler equation; since we are interested in computing financial moments and since bond prices in a recursive utility model depend on the value function, it is natural for us to approximate the value function directly. Judd and Guu (1997) and Judd (1998) are early examples of using the value function to generate perturbation conditions. Our particular solution method utilizes both the value function and the intertemporal Euler equation.

We use system (17)–(19) to build approximations of the value and policy functions:

where

and

To obtain these approximations, we take successive derivatives of Equations (17) and (18) with respect to K^_t and σ_z, and evaluate the resulting systems of equations at the deterministic steady state to obtain closed form solutions for the coefficients in Equations (25) and (26).

For example, evaluating Equations (17)–(20) at the deterministic steady state is sufficient to determine Z^ss, K^ss, V^ss(0,0) and C^ss(0,0). Taking first derivatives (with respect to K^_t and σ_z) of Equations (17) and (18) and again evaluating at the deterministic steady state allows us to solve for V^ss(1,0), V^ss(0,1), C^ss(1,0) and C^ss(0,1). Continuing in this fashion leads to the approximations in Equations (23) and (24), where the order of approximation is equivalent to the number of times we have differentiated Equations (17) and (18). We emphasize that the full system (17)–(20) is only used to determine the deterministic steady state, and the polynomial approximations are constructed by taking derivatives of only Equations (17) and (18). As mentioned in Aruoba, Fernández-Villaverde, and Rubio-Ramírez (2006), the first order solution involves a quadratic matrix equation, but each order of approximation thereafter only necessitates the solution of a linear system. Hence, higher order solutions only require a matrix inversion, albeit of rapidly increasing size.

With approximations V~pert(K^,σz) and C~pert(K^,σz) in hand, we can compute any other variable in the economy, where the accuracy of the approximation of those variables will depend on the underlying accuracy of our approximations for V^_t and C^_t. However, we can also approximate other variables of interest by augmenting system (17)–(19) with additional equilibrium conditions. In our case, we are interested in approximating both the risk-free rate and log⁡(V^t/C^t)_ (which we will use in computing welfare costs). The requisite equilibrium conditions are

Adding Equations (27) and (28) to system (17)–(19) allows us to obtain approximations

as outlined above.

3.2 Projection

Similar to perturbation, projection methods seek a polynomial approximation to a function, F, as in Equation (21), or more commonly to the special case, f, as in Equation (22). However, rather than using the known solution at ε = 0 to construct a local approximation, we specify a polynomial expansion, x^, with coefficients chosen to minimize f(x^) globally, over the domain of x. As before, for the neoclassical growth model, f would be comprised of the inter-temporal Euler equation and constraints, and a projection solution would specify a polynomial expansion of the consumption policy that would minimize the Euler equation error.

The analogous approach to our problem would be to use Equations (17) and (18) to obtain approximations,

and

where M is the order of approximation and φ_j, j=1,2,…, represent a set of linearly independent polynomial basis functions. That is, given an order of approximation M, we could specify a grid of N ≥ M + 1 points for K^ and evaluate Equations (17) and (18) [coupled with the constraint (19)] at those points to obtain a system of 2N equations in 2(M + 1) unknowns. We could then use a nonlinear solution method to find the coefficients a_j and b_j for j = 1, …, M, in Equations (31) and (32), that best satisfy (17) and (18).

As there is no theorem to guarantee convergence of the preceding approach, we follow an alternative methodology, suggested by Campanale, Castro, and Clementi (2010) and Kaltenbrunner and Lochstoer (2008), which is to couple polynomial approximations of the value and policy functions with value function iteration. Specifically, we seek a polynomial approximation to the value function as in Equation (31). Letting N_k ≥ M, we specify a (not necessarily equally spaced) grid for K^, spanning the values (0.1K^ss,1.9K^ss). Additionally, we set Nε=⌈M+12⌉ and confine ε to the order N_ε Gauss-Hermite abscissae. To ease notation, we suppress time subscripts, collect the N_k values of K^ in the vector K^ and group the basis functions evaluated at each value of K^ in the matrix Φ(K^), where Φ(K^)ij=φj(K^i). Using this notation, V~proj(K^)=Φ(K^)a. As this solution method is not standard practice, we provide pseudocode in A.

4.1 Calibration

We fix several parameters of our model and report them in Table 1. These values are a widely accepted parameterization of the US economy in the literature; in particular, the depreciation rate and share of capital, δ and α, respectively, are identical to those of Jermann (1998). The quarterly growth rate, μ, implies annual growth of 1.6 percent and the intertemporal elasticity of substitution (IES) parameter, ψ, is set in the middle of the range (1,2] advocated by Bansal, Kiku, and Yaron (2007). In reality, we considered alternative values of μ and ψ but do not report the corresponding solutions and simulations as they do not alter the qualitative nature of our results. Finally, the adjustment cost parameter, ξ, was chosen so that the ratio of volatilities of log consumption growth to log output growth matches empirical estimates (in the vicinity of 0.5), which depend on the time period and frequency of the data (see discussion below).

To understand our parameterization of the TFP volatility, it is instructive to consider the data moments reported in Table 2. The table contains means and volatilities for GDP, aggregate consumption and the 90 T-bill, both at annual and quarterly frequencies, for several sample periods. We highlight two important features of the data. First, the volatility of log output growth is markedly different between pre-war and post-war samples, the former being roughly 2.5 to 3 times as great as the latter. Second, the mean of the risk-free rate increases and its volatility decreases as the time horizon is curtailed to include fewer years. In the case of the mean, values in later samples are up to twice as large as that of the pre-war sample.

As a result of the variance in sample moments across sub-periods, we observe a wide range of calibrated values for the TFP volatility, σ_z, and the discount factor, β, in the literature. Since Std (Δy)≈(1−α)σz, models that calibrate to post-war, quarterly data often specify much smaller values of σ_z than models which calibrate to pre-war, annual data. Hence, we allow σz∈{0.01,0.02,0.03,0.04}, which correspond to Std (Δy)∈{0.0064,0.0128,0.0192,0.0256}, a range that encompasses the moments reported in Table 2.

For the remaining parameters, the subjective discount factor, β, and coefficient of relative risk aversion, γ, we entertain β ∈ [0.980, 0.998] and γ ∈ {2, 5, 10}. We choose these values because they not only encompass accepted values in the literature, but they allow a broad enough range of parameterizations to investigate their effect on the sensitivity of the value function to σ_z. In general, we are primarily concerned with β > 0.99, as these higher values are requisite for matching moments of the risk-free asset.

4.2 Model implications

4.3 Graphical evidence

We now provide graphical evidence to support the model solutions of Section 4.2. Figure 1 shows policy function approximations for 5th order Chebyshev projection and perturbation of orders 1, 2 and 3, all for the case of σ_z = 0.04 and β = 0.998 (for C^_t and V^_t, the “Pert” and “NPert” solutions are identical). Figure 2 and Figure 3 depict similar approximations for the value and bond price functions.

Figure 1:

Consumption policy approximations for Chebyshev projection (5th order) and perturbation of orders 1, 2 and 3, where σ_z = 0.04 and β = 0.998.

Figure 2:

Value function approximations for Chebyshev projection (5th order) and perturbation of orders 1, 2 and 3, where σ_z = 0.04 and β = 0.998.

Note the widely different scales on the vertical axes.

Figure 3:

Bond price function approximations for Chebyshev projection (5th order) and perturbation of orders 1, 2 and 3, where σ_z = 0.04 and β = 0.998.

Note the widely different scales on the vertical axes.

These plots clarify the results reported in Table 4: both projection and perturbation produce policy functions that are in close agreement, while there is a wide discrepancy in their solutions for the value and bond price functions, even among the different order perturbations. Since the quantity dynamics of the solutions are not sensitive to the value function, it is not surprising that the simulated macroeconomic moments of the two methods do not differ by a great amount. However, the risk-free rate and log⁡(V^t/C^t) both depend directly on the level and shape of the value function, and hence are quite different across methods. Figure 4–Figure 6 depict the same approximations for the case of σ_z = 0.01, and demonstrate that when the TFP volatility is low, the solution methods are more similar, as expected from the simulation output in Table 3. In this case, the consumption policy approximations overlap to an even greater extent, the value function approximations are separated by only a (relatively) small level shift and the bond price functions are shifted closer together and also overlap more. The discrepancies in the value and bond price functions further diminish as we shrink σ_z toward zero.

We note that in our model it is possible to analytically solve for the return on equity in terms of aggregate variables (see C), and hence it is unaffected by poor local approximations of the value function. It follows that local approximations of the equity premium are only affected by the value function via the risk-free rate. This result may not extrapolate to more general models where analytical expressions of the return on equity are not available.

Figure 4:

Consumption policy approximations for Chebyshev projection and perturbation of orders 1, 2 and 3, where σ_z = 0.01 and β = 0.998.

Figure 5:

Value function approximations for Chebyshev projection and perturbation of orders 2 and 3, where σ_z = 0.01 and β = 0.998.

Figure 6:

Bond price function approximations for Chebyshev projection and perturbation of orders 2 and 3, where σ_z = 0.01 and β = 0.998.

To understand why the value functions for the two solution methods diverge for large σ_z, it is useful to think of the value function and consumption policy as functions of both the state variable, K^, and the TFP volatility parameter, σ_z. Figure 7–Figure 9 show policy, value and bond price function approximations for 5th order Chebyshev projection and 3rd order perturbation, when σ_z ∈ [0, 0.04]. Thus, the approximations in Figure 1–Figure 6 are simply cross sections (fixing σ_z) of the functions depicted in Figure 7–Figure 9; e.g. the value function in panel (a) of Figure 2 corresponds to the cross-section of Figure 8 at σ_z = 0.04. It becomes apparent from inspecting these surfaces that the value function exhibits a high degree of curvature in the direction of σ_z, with the amount of curvature increasing as σ_z approaches zero, whereas the consumption policy is quite flat. As a result, similar to the square root function considered in 53, we anticipate that the radius of convergence of a Taylor polynomial approximation of the value function will diminish for approximations centered at points very close to σ_z = 0. This is a exactly what a perturbation solution is: a local Taylor approximation atσ_z = 0. It follows that, for the particular case of the value function, we have no guarantee of convergence for values of σ_z far away from zero, and in fact we should not be surprised to see divergent behavior, as suggested by the previous plots. This result is congruent with Den Haan and de Wind (2009) who find that nonlinearities in DSGE models can render high-order perturbation solutions that are explosive. On the other hand, the radius of convergence for the consumption policy is likely to be quite large, and we expect local Taylor approximations to converge for a wide range of σ_z – this is corroborated by Figure 7, where we are unable to distinguish the two surfaces. We are careful to note, however, that due to the nonlinearity and numerical difficulty of the problem, it is not possible to derive the true radius of convergence, and hence we cannot be certain that it is the cause of the large numerical discrepancies in the local approximation; rather, we simply use the nonlinearities of the approximated surfaces as a guide, suggesting the radius of convergence might be responsible.

Figure 7:

Consumption policy approximations for 5th order Chebyshev projection and 3rd order perturbation, for σ_z ∈ [0, 0.04].

Figure 8:

Value function approximations for 5th order Chebyshev projection and 3rd order perturbation, for σ_z ∈ [0, 0.04].

The lower surface corresponds to perturbation.

Figure 9:

Bond price function approximations for 5th order Chebyshev projection and 3rd order perturbation, for σ_z ∈ [0, 0.04].

The lower surface corresponds to perturbation.

The upshot of the foregoing results is that a global projection method is more robust to value function curvature, since it seeks to minimize an error equation expressed in terms of the true (unknown) value function, rather than approximating the truth at a distant focal point.

4.4 Solution evaluation and sensitivity analysis

In the preceding analysis we have merely shown some conditions under which the two solution methods we consider are different; we have not formally investigated their relative accuracy. We now undertake the important task of determining which of the solutions is a closer approximation to the unknown truth and do so for a variety of model parameter values. While our primary evaluation criterion will be Euler equation errors, we will conclude the section with a discussion of the Den-Haan-Marcet statistic (Den Haan and Marcet, 1994).

4.4.1 Pricing errors

The fundamental asset pricing equation is 1=Et[Mt+1Rt+1], where M_t+1 is the time t stochastic discount factor, and R_t+1 is the return for any asset between t and t + 1. Thus, from Equation (18) we have

(33)1=Et[Mt+1Rt+1e]=Et[Mt+1 ϕ′(I^tK^t)((α−1)Y^t+1+C^t+1K^t+1+ϕ(I^t+1K^t+1)+1−δϕ′(I^t+1K^t+1))].

Since the stochastic discount factor M_t+1 incorporates current consumption, C^_t, in its denominator [see Equation (14)], we can interpret pricing errors as a fraction of contemporaneous consumption. As suggested by Judd and Guu (1997), Aruoba, Fernández-Villaverde, and Rubio-Ramírez (2006), and Caldara et al. (2012), base 10 logarithms of pricing errors in Equation (33) can be interpreted in the following manner: a value of −1 corresponds to a 10% consumption error, a value of −2 corresponds to a 1% consumption error, a value of −3 corresponds to a 0.1% consumption error, etc.. Combining Equation (33) with the long simulations of ε_t [see Equation (20)] used in Section 4.2, we can compute the mean of the pricing errors implied by the model, for each solution method. The expectation is approximated by a Gauss-Hermite quadrature rule, with the order chosen so as to exactly compute the integral for the finite polynomial solutions. The pricing errors are reported graphically in the upper rows of Figure 10–Figure 13. The individual plots depict how mean Euler equation errors vary with β, where β ∈ [0.980, 0.998] – in general, perturbation solution quality degrades as β rises. Moving across the upper rows, from left to right, we are then able to observe the effect of increasing risk aversion, γ, and moving between the four figures we observe the effect of increasing TFP volatility, σ_z – as with β, the quality of the perturbation solution degrades as each of these parameters increases. At the lower extreme, in Figure 10, when σ_z = 0.01 and γ = 2, a 3rd order perturbation dominates a 5th order projection for virtually all values of β that we consider. However, both solutions produce errors that most would consider economically insignificant (less than 0.01% of consumption). Holding σ_z fixed and increasing γ, the perturbation errors rise to levels as high as 1% of consumption, for high values of β. These qualitative results become more pronounced in Figure 11–Figure 13, where at the upper extreme (σ_z = 0.04 and γ = 10), perturbation errors exceed 10% of consumption, for high values of β. It is this final case that deserves particular attention: models that calibrate to annual, pre-war data generally require high values of σ_z (on the order of 0.04) and β (on the order of 0.998 or above) in order to match output volatility and the level of the risk-free rate. We see that these calibrations, matched with moderate levels of risk aversion (above 5) can lead to poor local approximations. On the other hand, models that calibrate to quarterly, post-war data typically obtain much smaller values of σ_z (on the order of 0.01 or below), and do not suffer from poor local approximations.

Figure 10:

Mean log₁₀ Euler equation errors (first row), mean risk-free rate (second row) and mean log ratio of value function to consumption policy, plotted as functions of β for different values of γ.

σ_z = 0.01 in all cases.

Figure 11:

Mean log₁₀ Euler equation errors (first row), mean risk-free rate (second row) and mean log ratio of value function to consumption policy, plotted as functions of β for different values of γ.

σ_z = 0.02 in all cases.

Figure 12:

Mean log₁₀ Euler equation errors (first row), mean risk-free rate (second row) and mean log ratio of value function to consumption policy, plotted as functions of β for different values of γ.

σ_z = 0.03 in all cases.

Figure 13:

Mean log₁₀ Euler equation errors (first row), mean risk-free rate (second row) and mean log ratio of value function to consumption policy, plotted as functions of β for different values of γ.

σ_z = 0.04 in all cases.

The fundamental characteristic driving these results is the curvature of the value function with respect to TFP volatility: as shown in Figure 8, the value function can exhibit a high degree of curvature in the direction of σ_z. In cases where the curvature is extreme, a local method such as perturbation will have difficulty approximating the function at points far from the deterministic steady state (the locus of approximation), a result which is corroborated by Figure 10–Figure 13. This is especially relevant for models that require a high calibrated value of σ_z. The parameters γ and β have an effect insofar as they increase the sensitivity of the value function to changes in σ_z; i.e. the sensitivity increases with each of these parameters.

The second and third rows of Figure 10–Figure 13 depict the mean of Rtf and the mean of log⁡(V^t/C^t), respectively, across simulations. As in Section 4.2, we compute these values both nonlinearly, via Equations (14), (27) and (28), and directly, via Equations (29) and (30). As with the Euler equation errors, the nonlinear perturbation risk-free rate deviates dramatically from that of projection as β, γ and σ_z rise. This discrepancy is most pronounced for σ_z = 0.04, γ ≥ 5 and β ≥ 0.99. On the other hand, the direct perturbation computation appears to be quite similar to projection for all parameter values. In truth, the local method exhibits a small amount of divergent behavior as well, but the graphical evidence is washed out by the scale of the nonlinear deviation. The moments in Table 3 and Table 4 give an idea of the magnitude of divergence.

The reason for the discrepancy in the risk-free rate computations is the same as for the Euler equation errors: for high values of β, γ and σ_z, perturbation provides a poor approximation to the value function. Since the risk-free rate depends directly on the value function in models with recursive utility [see Equations (27) and (14)], it is likewise poorly approximated by perturbation, insofar as the value function approximation is poor. This effect is most severe when we compute the risk-free fate nonlinearly. Viewed from the opposite perspective, if our approximations for C^_t and V^_t were highly accurate, a nonlinear computation of Rtf would likewise be highly accurate. Thus, the large deviations in Figure 11–Figure 13 are just further evidence of the poor local approximation of V^_t for high β, γ and σ_z. The interesting aspect of our results, though, is that we can ameliorate the effect of the value function by directly computing the risk-free rate via a Taylor expansion. This latter method anchors the risk-free at its deterministic steady state and weakens the computational relationship between Rtf and V^_t.

The final rows of Figure 10–Figure 13 lend more insight to the foregoing results. As mentioned in Section 4.2, we include the log of V^t/C^t in our analysis, since this variable is instrumental in welfare evaluations. Once again, the approximations deteriorate as β, γ and σ_z increase, which we attribute to the poor local approximation of the value function. However, in the case of log⁡(V^t/C^t) the deviations are more severe (for the direct Taylor expansion method) than for the risk-free. The reason for this is that the risk-free rate depends on the value function in both numerator and denominator [see Equations (14) and (27)], which mitigates the error propagation of the value function approximation. The same could be true of welfare computations log⁡(V^t∗/C^t∗)−log⁡(V^t∗∗/C^t∗∗) when V^t∗ and V^t∗∗ are computed with the same σ_z (i.e. welfare effects are evaluated for a variable other than σ_z). However, since the approximations of V^_t are likely to be very different for different values σ_z (it will be well approximated for low volatility and badly approximated for high volatility), we expect that welfare computations attributed to TFP volatility will be very poorly approximated.

As a final note, we repeated all of the previous analysis for μ = 0 and ψ = 0.5; these changes only caused slight shifts, leaving the results above qualitatively the same. For space considerations, we do not include them. We recognize, however, the particular interplay between β and μ: lowering the growth rate increases the model’s ability to tolerate high values of β before local approximations become very poor. That is, we can think of the model as depending on the single growth adjusted subjective discount factor β∗=βexp⁡(μ) – lowering μ or β effectively decreases β^*, and hence the model’s sensitivity to σ_z.