Asset Pricing with Dynamic Programming ⁄

,


Introduction
In this paper we apply stochastic dynamic programming to intertemporal asset pricing models.In the literature we can find several models.The first generation of intertemporal asset pricing models has employed exogenous dividend stream 1 to match financial statistics of the models to the data.Those models have had difficulties to generate adequate characteristics such the risk-free interest rate, equity premium and the Sharpe-ratio, a measure of the risk-return trade-off.In general in those models the risk-free rate turns out to be too high and the mean equity premium and Sharpe-ratio too low compared to what one finds in time series data.
The next generation of asset pricing models has been based on the stochastic growth model with production originating in Brock and Mirman (1972) and Brock (1979Brock ( , 1982)).The Brock approach extends the asset pricing strategy beyond endowment economies to economies that have endogenous state variables including a capital stock that is used in production.Authors, building on this tradition, 2 have argued that it is crucial how consumption is modelled.In stochastic growth models the randomness occurs in the production function of firms and consumption and dividends are derived endogenously.Yet, the latter type of models have turned out to be even less successful.Given a production shock, consumption can be smoothed through savings and thus asset market facts are even harder to match.Recently numerous extensions of the basic model have been proposed. 3et, one of the most important issues in this endeavor to solve asset pricing models has become the accuracy of the solutions of the asset price model.In this paper we argue that we have to be sufficient confident in the accuracy of the solution method before we can apply it to models that are set up to solve the asset pricing puzzles.The conjecture has been that linearization methods to solve the optimal consumption path in feedback form and to use the growth of marginal utility of consumption as discount factor for pricing assets have been insufficient.Other methods have been developed that use second order approximations (Collard and Juillard, 2001), the perturbation method developed by Judd (1996) and the spline method (Trick and Zin, 1997), to name a few.Our attempt here is to discuss and to apply a solution technique, that builds on the the Hamilton-Jacobi-Bellman equation, the dynamic programming method, to address the accuracy issue and to solve asset pricing models.
The stochastic dynamic programming algorithm, with discretization of the state space and adaptive gridding strategy has already appeared useful in other applications. 4In such a method we do not need to use fixed grids, but adaptive space discretization.In some of the literature, see Munos and Moore (2002) and Trick and Ziu (1997), adaptive methods for dynamic programming have been used, but those methods either do not allow the derivation of rigorous error bounds or are computationally very expensive, see the discussion in Sect. 2. In the method applied here efficient and reliable local error estimation is undertaken and used as a basis for a local refinement of the grid in order to deal with the accuracy problem.This procedure allows for a global dynamic solution method of deterministic as well as stochastic intertemporal decision problems.
In order to test the accuracy of such a stochastic dynamic programming algorithm with flexible grid size as solution method for asset pricing models, we numerically provide the global solution and dynamics of a basic stochastic growth model, as based on Brock and Mirman (1972) and Brock (1978Brock ( , 1982)).The basic model can analytical be solved for the policy function in feedback form.Moreover, the asset prices, the risk-free interest rate, the equity premium and the Sharpe-ratio, can, once the model is solved analytically for the sequence of optimal consumption, easily be solved numerically and those solutions compared to the numerical solutions obtained by our stochastic dynamic programming method.As will be shown the errors, as compared to the analytical solutions, are quite small.Our method uses adaptive space discretization.In the method efficient and reliable local error estimation is undertaken and used as a basis for a local refinement of the grid where we are also able to deal with steep slopes or other non-smooth properties of the value function (such as non-differentiability), obtaining quite accurate solutions even further away from some stationary state 5 .
After the method is tested for its accuracy we can then apply it to more complicated extensions of the model which may include different utility functions, in particular habit formation, adjustment cost of capital, a two sector economy, or heterogenous firms and households. 6In our extensions we focus here on habit formation and adjustment costs of investment.Our new solution method is applied to such an extended asset pricing model.Since, as aforementioned, time separable preferences fail to match financial market characteristics an enormous effort has been invested into models with time non-separable preferences, such as habit formation models, which allow for adjacent complementarity in consumption.Past consumption enters here as a constraint, defined as external habit persistence where the aggregate level of consumption serves as a benchmark level, or internal habit persistence where a household's own past consumption is viewed as a benchmark.
There is a long tradition in economic theory where it is assumed that habits are formed through past consumption. 7In those models high level of consumption in the past depresses current welfare and high current consumption depresses future welfare.This can from it. 5For the former two problems, see Grüne and Semmler (2004), and for the latter, see Becker, Grüne and Semmler (2006). 6For further detailed studies see, for example, Campbell and Cochrane (1999), Jerman (1998), Boldrin, Christiano and Fisher (2001), Cochrane (2001, ch. 21) and Akdezin and Deckert (1997) for extensions along those lines. 7Habit persistence is nowadays used to understand a wide range of issues in growth theory (Carrol et al. 1997, 2000, Alvarez-Cuadrado et al. 2004) macroeconomics (Fuhrer, 2002), and business cycle theory (Boldrin et al, 2001).For a first use of habit persistence in a dynamic decision model see Ryder and Heal (1973).
be written as ratios of current over past consumption (Abel 1990(Abel , 1999) ) or in difference form as (1 − α)c t + α(c t − c t−1 ) with c t current, c t−1 past consumption and α a respective weight.This form of habit formation will be chosen in this paper.This type of habit specification gives rise to time non-separable preferences where risk aversion and intertemporal elasticity of substitution are separated and a time variation of risk aversion will arise.If we define surplus consumption as s t = ct−Xt ct with X t , the habit, and γ, the risk aversion parameter, then the time variation of risk-aversion is γ st : the risk aversion falls with rising surplus consumption and the reverse holds for falling surplus consumption.A high volatility of the surplus consumption will lead to a high volatility of the growth of marginal utility and thus to a high volatility of the stochastic discount factor.
Habit persistence in asset pricing has been introduced by Constantinides (1990) in order to study the equity premium problem.Asset pricing models along this line have been further explored by Campbell and Cochrane (1999), Jerman (1998), and Boldrin et al. (2001).Yet, as mentioned before, introducing habit persistence in stochastic models with production may just produce smoother consumption.As the literature has demonstrated (Jerman 1998, andBoldrin et al. 2001) one also needs adjustment costs of investment to obtain better results on the equity premium.By choosing such a model we will not, following Jerman (1998), allow for elastic labor supply, but rather employ a model with fixed labor supply, since the latter, as shown in Lettau and Uhlig (2000), provides the most favorable case for matching the model with the financial market characteristics.
The paper is organized as follows.Section 2 discusses related literature on solution methods.Section 3 presents the stochastic dynamic programming algorithm.Section 4 introduces the stochastic decision problem in asset pricing.Section 5 applies our method to an analytically solvable basic model of asset pricing and studies the accuracy of our solution method.Section 6 introduces habit formation and adjustment costs of investment and solves for the financial characteristics we want to study.Section 7 interprets our results in the context of previous literature on habit formation models.Section 8 concludes the paper.

Related Literature on Solution Methods
In the literature on solving asset pricing models one can find a vast amount of different approaches. 8We will focus here on those approaches that are employing dynamic programming which is designed to simultaneously approximates both the decision variables and the value function.It is closely related to the Hamilton-Jacobi-Bellman equation.Many of the recent versions of dynamic programming use state-of-the art mathematical and numerical techniques for making this approach more efficient.Here we apply an adaptive gridding algorithm that works for very general Hamilton-Jacobi-Bellman equations, see Section 3 for details.In the present section we briefly review similar approaches and highlight similarities and differences to our approach.
One of the fundamental difficulties with the dynamic programming approach is that the computational load grows exponentially with the dimension of the problem, a phenomenon known as the "curse of dimensionality" (see Rust (1996) for a comprehensive account on complexity issues).In our case, for computing asset pricing in the context of stochastic growth models, starting with Brock and Mirman (1972) as suggested in the literature, the problem to be solved is two dimensional, hence this is not a crucial aspect here.Nevertheless, for the sake of completeness we want to mention approaches like randomly distributed grid points (Rust (1997)) or so called low discrepancy grids (Rust (1996), Reiter (1999)) which are able to break the curse of dimensionality.In principle also Monte-Carlo techniques like in Keane and Wolpin (1994) allow for breaking the curse of dimensionality, but as Rust (1997) points out, the specific algorithm in Keane and Wolpin (1994) uses an interpolation technique which again is subject to exponential growth of the numerical cost in the space dimension.
For low dimensional problems the goal of the numerical strategy is not to avoid the curse of dimensionality but rather to reduce the computational cost for a problem of fixed dimension.For this purpose, two main approaches can be found in the literature, namely higher order approximations and adaptive gridding techniques; the latter will be used in our numerical approach.
The idea of high order approximations lies in exploiting the smoothness of the optimal value function: if the optimal value function turns out to be sufficiently smooth, then methods using approximations by smooth functions, like Chebyshev polynomials (Rust (1996), Judd (1996), Jermann (1998)), Splines (Daniel (1976), Johnson et al. (1993), Trick andZin (1993, 1997)) or piecewise high-order approximations (Falcone and Ferretti (1998)) can be very efficient.Smoothness of the optimal value function is also assumed in the method using second order approximation of the policy function (see Schmidt-Grohe and Uribe 2004 and Collard and Juillard, 2001).Smoothness is also the basis of other high-order strategies, like in finite difference approximations (Candler (2001)), Gaussian Quadrature discretization (Tauchen and Hussey (1991), Burnside (2001)) and in perturbation techniques (Judd (1996)).Yet, the latter should also work if the value function is only piecewise smooth.Some of these methods (like Spline and piecewise high order approximation) use a (fixed) grid discretization of the state space similar to our approach.The combination of adaptive grids with higher order approximation has been investigated for several one dimensional examples in Bauer, Grüne and Semmler (2006).This combination yields good results when the optimal value function is smooth, however, the resulting numerical scheme may become unstable in the presence of non-differentiable optimal value functions.Since, in addition, the practical implementation of spline interpolation on non-uniform grids in more than one space dimension becomes rather complicated and computationally expensive, we have chosen piecewise linear approximations in this paper.
Concerning discretization techniques it should be noted that from the complexity point of view it turns out to be optimal to solve the dynamic programming problem on successively finer grids, using a one-way multigrid strategy (Chow and Tsitsiklis (1991), see also Rust (1996)).In fact, our adaptive gridding algorithm is similar to this approach since the approximation on the previous grid Γ i is always used as the initial value for the computation on the next finer adaptive grid Γ i+1 .This also explains the large reduction in computation time observed for our approach compared to the computation on one fixed equidistant grid.
Let us now turn to our adaptive gridding techniques.Perhaps closest to our approach are the techniques discussed in Munos and Moore (2002).Here a number of heuristic techniques are compared which lead to local and global error indicators which can in turn be used for an adaptive grid generation.Some of the indicators discussed in this paper bear some similarity with our residual based estimator, though rigorous estimates as given below are not given there.In any case, the authors report that these techniques are unsatisfactory and argue for a completely different approach which measures the influence of local errors in certain regions on the global error by analyzing the information flow on the Markov chain related to the discretization of the (deterministic) problem at hand.The reason for this lies in the fact that the model problem treated by Munos and Moore (2002) has a discontinuous optimal value function, which often happens in technical problems with boundary conditions.In fact, also our adaptive scheme performs rather poorly in presence of discontinuities but since our economic problems do always have continuous optimal value functions, Munos' and Moore's conclusions do not apply here.A roughly similar technique is the endogenous oversampling used by Marcet (1994).This is again a heuristic method, which, however, does not lead to adaptive grids but rather selects suitable parts of the state space where the optimally controlled trajectories stay with high probability.
Probably the adaptive approaches with the most solid mathematical background are presented in the papers of Trick andZin (1993, 1997). 9In these papers an alternative approach for the solution of the fully discrete problem is developed using advanced linear programming techniques which are capable of solving huge linear programs with many unknowns and constraints.In Trick and Zin (1993) an adaptive selection of constraints in the linear program is used based on estimating the impact of the missing constraint, a method which is closely related to the chosen solution method but only loosely connected to our adaptive gridding approach.The later paper (Trick and Zin (1997)), however, presents an idea which is very similar to our approach.Due to the structure of their solution they can ensure that the numerical approximation is greater than or equal to the true optimal value function.On the other hand, the induced suboptimal optimal control strategy always produces a value which is lower than the optimal value.Thus, comparing these values for each test point in space one can compute an interval in which the true value must lie, which produces a mathematically concise error estimate that can be used as a refinement criterion.While this approach is certainly a good way to measure errors, which could in particular be less conservative than our measure for an upper bound, we strongly believe that it is less efficient for an adaptive gridding scheme, because (i) the estimated error measured by this procedure is not a local quantity (since it depends on the numerical along the whole suboptimal trajectory), which means that regions may be refined although the real error is large elsewhere, and (ii) compared to our approach it is expensive to evaluate, because for any test point one has to compute the whole suboptimal trajectory, while our residual based error estimate needs only one step of this trajectory.
Let us comment on the idea of a posteriori error estimation.In fact, the idea to evaluate residuals can also be found in the papers of Judd (1996) and Judd and Guu (1997), using, however, not the dynamic programming operator but the associated Euler equation.In these references the resulting residual is used to estimate the quality of the approximating solution, but to our knowledge it has not been used to control adaptive gridding strategies, and we are not aware of any estimates such as ours which is a crucial property for an efficient and reliable adaptive gridding scheme, particularly needed to solve stochastic problems in asset pricing models. 10ummarizing our discussion, there are a number of adaptive strategies around which are all reported to show good results, however, they are either heuristic 11 and better suited for other classes of problems than asset pricing models or they have nice theoretical features but are practically inconvenient because their implementation is numerically much more expensive than our approach.

Stochastic Dynamic Programming
Next we describe the stochastic numerical dynamic programming algorithm that we use to solve the asset pricing characteristics of the stochastic growth model.Details of the adaptive gridding strategy are given in the appendix.
We consider the discrete stochastic dynamic programming equation Here x ∈ Ω ⊂ R 2 , C ⊂ R, Ω and C are compact sets and ε is a random variable with values in R. The mappings ϕ : Ω × C × R → R 2 and g : Ω × C × R → R are supposed to be continuous and Lipschitz continuous in x.Furthermore, we assume that either ϕ(x, c, z) ∈ Ω almost surely for all x ∈ Ω and all c ∈ C, or that suitable boundary values V (x) for x ∈ Ω are specified, such that the right hand side of (3.1) is well defined for all x ∈ Ω.The value β(x, ε) is the (possibly state and ε dependent) discount factor which we assume to be Lipschitz and we assume that there exists β 0 ∈ (0, 1) such that β(x, ε) ∈ (0, β 0 ) holds for all x ∈ Ω.We can relax this condition if no maximization takes place, in this case it suffices that all trajectories end up in a region where β(x, ε) ∈ (0, β 0 ) holds.This is the situation for the asset price problem, see the discussion in Cochrane (2001: 27).
Associated to (3.1) we define the dynamic programming operator given by The solution V of (3.1) is then the unique fixed point of (3.2), i.e., For the numerical solution of (3.3) we use a discretization method that goes back to Falcone (1987) in the deterministic case and was applied to stochastic problems in Santos and Vigo-Aguiar (1998).Here we use unstructured rectangular grids: We assume that Ω ⊂ R n is a rectangular and consider a grid Γ covering Ω with rectangular elements Q l and nodes x j and the space of continuous and piecewise multilinear functions where the e j , j = 1, . . ., n denote the standard basis vectors of the R n , see Grüne (2003) for details of the grid construction.With π Γ : C(Ω, R) → W Γ we denote the projection of an arbitrary continuous function to W Γ , i.e., π Γ (W )(x j ) = W (x j ) for all nodes x j of the grid Γ.
Note that our approach easily carries over to higher order approximations, the use of multilinear approximations is mainly motivated by its ease of implementation, especially for adaptively refined grids. 12Also, the approach can easily be extended to higher dimensions.
We now define the discrete dynamic programming operator by with T from (3.2).Then the discrete fixed point equation has a unique solution V Γ ∈ W Γ which converges to V if the size of the elements Q l tends to zero.The convergence is linear if V is Lipschitz on Ω, see Falcone (1987), and quadratic if V is C 2 , see Santos and Vigo-Aguiar (1998).
For the solution of (3.5) as well as for the computation we need to evaluate the operator T Γ .More precisely, we need to evaluate for all nodes x j of Γ.This first includes the numerical evaluation of the expectation E. If ε is a finite random variable then this is straightforward, if ε is a continuous random variable then the corresponding integral has to be computed, where f is the probability density of ε.In our implementation we approximated this integral by a trapezoidal rule with 10 equidistant intervals.
The second difficulty in the numerical evaluation of T lies in the maximization over c.In our implementation we used Brent's methods with an accuracy of 10 −3 .Here it turned out that due to the errors induced by the spatial discretization on the grid higher accuracies in Brent's method did not yield higher accuracies in the results.
For the solution of the fixed point equation (3.5) we use the Gauss-Seidel type value space iteration which is described in Section 3 of Grüne (1997) (under the name "increasing coordinate algorithm"), where we subsequently compute V i+1 = S Γ (V i ) with S Γ being a Gauss-Seidel type iteration operator (including the maximization over c) obtained from T Γ .This iteration is coupled with a policy space iteration: Once a prescribed percentage of the maximizing u-values in the nodes remains constant from one iteration to another we fix all control values and compute the associated value function by solving a linear system of equations using the CGS or BICGSTAB method (in our examples the CGS method turned out to show more reliable convergence behavior).After convergence of this method we continue with the value space iteration using S Γ until the control values again converge, switch to the linear solver and so on.This combined policy-value space iteration turns out to be much more efficient (often more than 90 percent faster) than the plain Gauss-Seidel value space iteration using S Γ , which in turn is considerably faster than the Banach iteration The gridding algorithm is presented in the appendix.

The Stochastic Decision Problem in Asset Pricing
Our stochastic decision problem arising from the stochastic growth model in the Brock tradition which we want to solve and for which we want to compute certain financial measures is as follows.Before we introduce the basic stochastic growth model that we want to apply our solution technique to, see sect.5, we outline an asset pricing model in a very generic form.The problem we are concerned with is to compute an optimal control, c t , for the dynamic decision problem with X t representing other state variables, for example habit, subject to the dynamics using the constraints c t ≥ 0 and k t ≥ 0 and the initial value Here (k t , z t , X t ) ∈ R 3 is the state and ε t are i.i.d.random variables.We abbreviate This optimal decision problem allows for the computation of c in feedback form, i.e. c t = c(x t ) for some suitable map c : R 2 → R. Based on this c we compute the stochastic (note that m depends on ε t and the derivative u is taken with respect to c t ), which serves as an ingredient for the next step, which consists of solving the asset pricing problem where d(x t ) denotes the dividend at x t and x 0 = x and the dynamics are given by with c from above.
Finally, we use these values to compute the Sharpe-ratio, which represents the ratio of the equity premium to the standard deviation of the equity return.Hereby R f is the risk-free interest rate.
is the risk-free interest rate and is the gross return.
Note that the equality E(m(x)R(x)) = 1 holds, which can serve as a indicator for the accuracy of our numerical solution.
We solve the asset pricing problem in the following three steps: (i) We compute the optimal value function V of the underlying optimal control problem, and compute c from V , (ii) we compute the prices p(x) from c and m, and (iii) we compute the risk-free interest rate, the equitiy premium and the Sharpe-ratio S from c, m and p.
For our baseline stochastic growth model introduced below, which we solve numerically, both c and p are actually available analytically.This allows us to test each single step of our algorithm by replacing the numerically computed c in (ii) and (iii) and/or p in (iii) by their exact values.
For each of the steps we do now sketch our technique for the numerical computation using the algorithm described above in Section 3.
Step (i): Compute the approximate optimal value V Γ )(x) using the method from Section 3 and the adaptive gridding algorithm from the Appendix.
Once V Γ is computed with sufficient accuracy we obtain the optimal control value c(x) in each point by choosing c(x) such that holds.Once c is known, m of equ.(4.2) can be computed from this value.
Step (ii): For computing p(x) we follow the same approach as in Step (i), except that here c(x) is known in advance and hence no maximization needs to be done.
For the computation of p we first solve the dynamic programming equation which is simply a system of linear equations which we solve using the CGS method.This yields a numerical approximation of the function (with the convention 0 s=1 m(x s ) = 1), from which we obtain p by In our numerical computations for the computation of p we have always used the same grid Γ as in the previous computation of V in Step (i).The reason is that it did not seem justified to use a finer grid here, because the accuracy of the entering values c from Step (i) is limited by the resolution of Γ, anyway.However, it might nevertheless be that using a different grid (generated e.g. by an additional adaptation routine) in Step (ii) could increase the numerical accuracy.
Step (iii): The last step is in principle straightforward, because we do now have all the necessary ingredients to compute the risk-free interest rate, the equity premium and Sharpe-ratio S.However, since all the numerical values entering these computations are subject to numerical errors we have to be concerned with the numerical stability of the respective formulas.The first formula for the Sharpe-ratio, where the numerator represents the equity premium, as the spread between the expected equity return and the risk-free interest rate, is This turns out to be considerably less precise than the second formula Since the denominator is the same in both formulas, the difference can only be caused by the different numerators.A further investigation reveals that the numerator of the first formula can be rewritten as while that of the second formula reads Note that in both formulas we have to subtract values which have approximately the same values, which considerably amplifies the numerical errors.As mentioned above, we know that E(m(x)R(x)) = 1, which shows that these formulas are theoretically equivalent.Yet the second formula is more accurate. 14  5 Testing Dynamic Programming for a Basic Model Next we perform numerical computations for the stochastic growth model as present in Brock and Mirman (1972) and Brock (1979Brock ( , 1982)).We use a basic version as suggested and used by Santos nd Vigo-Aguiar (1998) which employs an aggregate capital stock and log utility.The variable X t in equ.(4.1) will therefore be zero.
The dynamics are defined by where A, α and ρ are real constants and the ε t are i.i.d.random variables with zero mean.The return function in (4.1) is u(c) = ln c.
In our numerical computations we used y t = ln z t instead of z t as the second variable.For our numerical experiments we employed the values A = 5, α = 0.34, ρ = 0.9, β = 0.95, 14 The higher accuracy of the second formula can be explained as follows: Assume that we have a small systematic additive (deterministic) numerical error in R(x), e.g., Rnum(x) ≈ R(x) + δ.Such errors are likely to be caused by the interpolation process on the grid.Then, using R f (x) = 1/E(m(x)), in the first formula we obtain while in the second formula we obtain i.e., systematic additive errors cancel out in the second formula.
The following figures show some of the numerically obtained results for this example.Figure 5.1 shows the optimal value function, the resulting final adaptive grid with 2977 nodes, and the corresponding asset price p. Figure 5.2 shows the k-components of several optimal trajectories with different initial values.Observe that all trajectories end up near the "equilibrium" point (k, ln z) = (2, 0) around which they oscillate due to the stochastic influence.For our stochastic growth model the optimal control is known analytically and is given by Note that c here depends linearly on z.
Since u (c) = 1/c, the stochastic discount factor becomes Furthermore, in this model the dividend is given by d(x) = c(x), and from this one easily computes15 Using these values, we can compute the Sharpe-ratio by evaluating the respective expressions numerically.We obtain the value S = 0.007999.
Note that the value is independent of x, but the component parts defining S are not constant at all.In Figure 5.
defining the numerator of the Sharpe-ratio, are shown depending on x, all computed from the exact values.In addition, Figure 5.4 shows a larger plot of the risk free interest rate R f (x) together with a plane at level R f = 1, which allows to identify where the rate is below and where above one.Note that in the equilibrium x = (k, ln z) = (2, 0), around which the optimal trajectories end up in the long run, the value of R f is 1.06.For the test with medium accuracy we used an adaptive grid with 15020 nodes and for the fine accuracy we used 109533 nodes.Table 5.1 shows the respective results with the errors measured in the L ∞ norm.As an alternative, we used the same grids but using the exact c(x) in the grid points with interpolation in between for all our computations (i.e., for p and S).This rules out the errors caused by the maximization routine and leaves us with the pure interpolation errors on the grid.The corresponding accuracies are summarized in Table 5  From the tables it is interesting to see how the original numerical errors in V , c and p propagate through the consecutive numerical approximations.In particular, it turns out that the error in S is much smaller than the error in p and about one order of magnitude smaller than the error in c; in fact in all examples we investigated it is very similar to the error in the underlying optimal value function V from which c is obtained.
We also see that for the exact c even on the coarser grid very good approximations for S are obtained.From our experiments it seems that the crucial point in the numerical procedure is the numerical accuracy for the underlying V and c, while the accuracy of p does not seem to be crucial.Note that the error in the optimal consumption affects the asset price computation in two ways: on the one hand errors in c affect the stochastic discount factor and on the other hand errors in c cause errors in the optimal paths.This combination of two effects seems to be the main cause for the rather large numerical errors in p.On the other hand, the interpolation in the computation of p does not seem to be crucial for the accuracy, as the results in Table 5.2 show.Summarizing, if the consumption strategy c is either known exactly or if it can be determined numerically with a sufficiently high accuracy, then our numerical approach can be expected to be successful.
Note that the errors given above are obtained by maximizing the errors on a 250x250 grid covering the whole computational domain Ω.In a vicinity of the equilibrium the errors are smaller; for instance, the errors in the Sharpe-ratio S are 1.4 • 10 −3 , 9.3 • 10 −4 and 3.7•10 −4 for the fully numerical approach with coarse, medium and fine grid, respectively. 16 Furthermore, we discuss absolute errors here.The reason for this lies in the fact that in our examples the stochastic process determining the magnitude of the Sharpe-ratio does not seem to affect the accuracy, as we will demonstrate in our next numerical test.
We should mention that in this examples for σ = 0.008, for the fully numerical approach the errors for the computation of the Sharpe-ratio S for coarse and medium grids are almost of the same order as the Sharpe-ratio itself from which one might conclude that they do not provide a reasonable approximation.In this case, the relative error for S would be of the order of one, i.e. the numerical results could be considered useless.Yet, we would like to explore whether this undesirable behavior remains true if we modify our model parameters in such a way that the Sharpe-ratio (which is unrealistically small as compared to the empirical Sharpe-ratio 17 ) increases.
To this end we changed our random variables ε i to a Gaussian distributed random variable with standard deviation 18 σ = 0.018, which we now restricted to the interval [−0.072, 0.072].In order not to increase the computational domain Ω we reduced the parameter ρ to ρ = 0.5 which ensures that Ω = [1,4] × [−0.32, 0.32] remains controlled invariant.
The exact value of S for these parameters evaluates to S = 0.017994. 16We want to note here, that this approximation does not rise much neither for the Sharpe-ratio, when we depart for the equilibrium.This is not so in methods that use the second (or higher) order approximations, as for example employed in Schmidt-Grohe and Uribe (2004), where the error rapidly rises with the distance to the equilibrium, for details of a comparison of the solution from that methods and a dynamic programming solution, see Becker, Grüne and Semmler (2006). 17For a report on empirically realistic Sharpe-ratios, see for example, see Section 7. 18 The standard deviation of σ = 0.018 is suggested in the RBC literature, see Boldrin, Christiano and Fisher (2001).
For this set of parameters we have repeated the fully numerical computation on an adaptive grid with 2529 nodes, which is roughly the same amount of nodes as for the coarse accuracy, above.Table 5.3 shows the respective numerical accuracies for σ = 0.018 compared to the accuracies for the model with lower standard deviation σ = 0.008 from Table 5.1.It turns out that the absolute numerical error in the approximation S does not increase when σ increases, implying that the numerical accuracy of our procedure appears to be independent of the standard deviation σ of the underlying random process and thus of the magnitude of the Sharpe-ratio.In other words, when S increases due to larger shocks, then the relative error in S can be expected to decrease. 19  Summing up, our intention so far was to introduce and apply a stochastic dynamic programming algorithm to a basic asset pricing model in order to provide global and reliable solutions.In order to test our algorithm we have applied it to a basic stochastic growth model for which consumption, asset prices and the Sharpe-ratio can analytically be computed.Overall, as we have shown, the errors depend on the discretization step size and grid refinement but are less impacted by the variance of the stochastic component of the model.Our computations show that the optimal consumption, the value function and the Sharpe-ratio can be computed with small absolute errors.We, in particular, could show that the error in the computation of the Sharpe-ratio is independent of the standard deviation of the underlying random process.Our method remains accurate for large shocks and the relative error decreases with increasing Sharpe-ratio. 20Overall the results presented here are very encouraging and our method can be applied to extended versions of the stochastic growth model with some confidence.
where I t = z t Ak α t −c t , where in our numerical computations we used the variable y t = ln z t instead of z t as the second variable.
The utility function in (4.1) is now given by for γ = 1.Since we are working with internal habit, in our case, we have X t = C t−1 .
For our numerical experiments we employed the same parameter values as in section 5.Here too ε t was chosen as a Gaussian distributed random variable with standard deviation σ = 0.008, which we restricted to the interval [−0.032, 0.032].With this choice of parameters it is easily seen that the interval [−0.32, 0.32] is invariant for the second variable y t .
Motivated by our study in Sect. 5 we would like to solve our problem for k t on an interval of the form [a, b]; here we used [0.1, 10].However, the habit persistence implies that for a given habit X t only those value c t are admissible for which c t − bX t > 0 holds, which defines a constraint from below on c t depending on the habit X t .On the other hand, the condition I t ≥ 0 defines a constraint from above on c t depending on k t and y t = ln z t .As a consequence, there exist states x t = (k t , y t , X t ) for which the set of admissible control values c t is empty, i.e., for which the problem is not feasible.On the one hand, we want to exclude these points from our computation, on the other hand we want to have a computational domain which is of a simple shape.A solution for this problem is given by a coordinate transformation which transforms a suitable set Ω of feasible points to the set Ω = [0.1,10] × [−0.32, 0.32] × [0, 7] on which we perform our computation.The coordinate transformation Ψ : R 3 → R 3 we use for this purpose shifts the k-component of a point x = (k, y, X) in such a way that the non-feasible points are mapped to a point xt ∈ Ω.It is given by where s 0 = 0.1 α exp(−0.32)A is chosen such that for y = −0.32 and X = 0 the coordinate change is the identity.This map is built in such a way that for all points x t ∈ Ω = Ψ −1 (Ω) a value c t with c t − bX t ≥ s 0 is admissible.Note that forward invariance of Ω is not automatically guaranteed by this construction and indeed there are parameter combinations for which this property does not hold.
This coordinate transformation allows us to set up our dynamic programming algorithm on a set of feasible points without having to deal with a complicated domain Ω, because numerically we can now work on the simple set Ω using the transformed dynamics instead of (4.2).
In addition to the parameters specified above, for our numerical experiments we have used the following sets of parameters: Note that (a) corresponds to the setting of section 5.
Our first set of numerical results shows the behavior of the k t -component of the optimal trajectories of the optimal control problem (4.1) in Figure 6.1 (a)-(e), the stochastic discount factor along these trajectories in Figure 6.2(a)-(e) and the consumption in Figure 6.3(a)-(e).For all trajectories we used the initial value (k 0 , y 0 , X 0 ) = (2, 0, 4) which is near the point around which the optimal trajectory oscillates.Furthermore, for all trajectories we have used the same sequence of the random variables ε t .The following Table 6.1 shows several characteristic values obtained from our algorithm.The values were obtained by averaging the values along the optimal trajectories described above.As can be observed from the figures 6.2 (b) and (d) the volatility of the stochastic discount factor is increased due to habit formation and time varying risk aversion.It is most volatile for the combination of a high γ and habit persistence, see figure 6.2 (c) and (e).In the latter two cases the curvature of the welfare function affects, as one would expect, the volatility of the discount factor.Yet, as the figures 6.3 (a)-(b) show, the consumption path itself is only very little affected by habit persistence and adjustment costs of capital. 21rom the table 6.1 we can observe that the Sharpe-ratio and equity premium increase strongly with habit persistence and adjustment costs, though not sufficiently to match empirical facts. 22The risk free interest rate is still too high.In particular, as is also well known from asset pricing studies with power utility, the risk free rate rises with higher parameter of relative risk aversion, γ, see table 6.1 (c) and (e).

Evaluating the Results on Habit Formation
It is interesting to compare our numerical results that we have obtained, to previous quantitative studies undertaken for habit formation.We in particular will restrict ourselves to a comparison with the results obtained by Boldrin et  Yet, there are two reasons why the results of those studies differ from ours.First, in comparison to their parameter choice we have chosen parameters that have often been used for the study of the basic stochastic growth models 23 and appear to us less extreme.Table 7.1 reports the parameters and the results.
Both, the study by Jerman (1998) and Boldrin et al. (2001) have chosen a parameter, ϕ = 4.05, in the adjustment costs of investment, a very high value which is at the very upper bound found in the data. 24Since the parameter ϕ smoothes the fluctuation of the capital stock and makes the supply of capital very inelastic, we have rather worked with a ϕ = 0.8 in order to avoid such strong volatility of returns generated by high ϕ.Moreover, both papers use a higher parameter for past consumption, b, than we have chosen.Both papers have also selected a higher standard deviation of the technology shock.Boldrin et al. take σ = 0.018, and Jerman takes a σ = 0.01, whereas we use σ = 0.008 which has been employed in many models. 25Those parameters increase the volatility of the stochastic discount factor, a crucial ingredient to raise the equity premium and the Sharpe-ratio.Jerman, in addition, takes a very high parameter of relative risk aversion, a γ = 5, which also increases the curvature of the welfare function, the volatility of the discount factor and the equity premium when used for the pricing of assets.Jerman also takes a much higher persistence parameter for the technology shocks, a ρ = 0.99, from which one knows that it will make the stochastic discount factor more volatile too.All in all, both studies have chosen parameters which seem to bias the results toward the empirically found financial characteristics.
Second, the studies by Jerman (1998) and Boldrin et al (2001) have employed different procedures to solve the intertemporal decision problem.Boldrin et al. use the Lagrangian multiplier from the corresponding planner's problem to solve for asset prices.The procedure is based on Christiano and Fisher (2000).Jerman uses a log-linear approach to solve the model and the procedure is compared to non-linear approximation.Concerning the accuracy of different solution procedures there are intricate issues involved such as approximation errors due to the order of approximation, the role of parameter choice (steeper curvature creates greater errors), the distance to the equilibrium, absolute or relative errors and so on. 26We also want to note that there is a crucial constraint in habit formation models, namely that the surplus consumption has to remain non-negative when the optimal solution, c t , is computed. 27As we have shown in Section 6 this constraint has to be treated properly in the numerical solution method.
Moreover, we want to note, since the risk aversion, γ st for power utility, falls with the rising consumption surplus ratio, given by s t = ct−Xt ct , the habit persistence model predicts a rising risk aversion (rising Sharpe-ratio) in recessions and falling risk aversion (falling Sharpe-ratio) in booms, for details see Cochrane (2001: 471).Thus, the risk aversion and the Sharpe-ratio move over the business cycle.A state dependent risk aversion and Sharpe-ratio were also observable in our computations.What we have depicted in table 7.1 are averages along the optimal trajectories in the neighborhood of the steady state.
A further remark is needed on the high risk free rate.In table 6.1 it is computed as 1 + R f and for quarterly data.In table 7.1, we have R f for annualized data, in percent terms.As table 7.1 shows the risk free rate is highest for the high γ, and with this the annual risk free rate can go up to 34 percent.(In the basic model with γ = 1 and no habit persistence and to adjustment cost, R f is about 24 percent.)This is much too high as compared to empirical data where the annual risk free rate is usually reported to be about 1 to 2 percent (depending on time periods considered).Yet, as well known in the literature 28 the equity premium and Sharpe-ratio are favorably impacted by a higher γ which, however, produces also a higher risk free rate.Thus, an attempt to increase the equity premium and Sharpe-ratio in the context of the above presumed utility function will make the risk free rate moving in the wrong direction, namely it rises too.We also want to note that the high risk free rate is computed by using a low β.We have chosen β = 0.95, as compared to a β = 0.99 and β = 0.999 in the study by Boldrin et al. and Jerman respectively.As also well known in the literature a higher β will reduce again the risk free rate. 29  Overall, one is inclined to conclude that previous studies, because of the specific parameter 26 For further discussion on those issues, see also Collard and Juillard (2001), in particular on the order of approximation and the role of parameters.choice, have obtained more favorable results in solving the equity premium puzzle.Though also our results for the habit formation model point in the same direction, we show by using a method that has been tested for accuracy, that there are still puzzles remaining for the consumption-based asset pricing models.Finally, we want to note that in our study we have chosen a model variant with no endogenous labor supply, which, as Lettau and Uhlig (2000) show, is the most favorable model for asset pricing in a production economy, since including labor supply as a choice variable, would even reduce the equity premium and the Sharpe-ratio.

Conclusions
In recent times extensive research effort has been devoted to study the asset price characteristics, such as the risk-free interest rate, the equity premium and the Sharpe-ratio, based on the stochastic growth model.Although it includes production this is a consumptionbased asset pricing model.The failure of the model with simple preferences, such as power utility, to match the empirical characteristics of asset prices and returns has given rise to the exploration of numerous extensions of the basic model that allows for different preferences and technology shocks, adjustment cost of capital, two sector economies and heterogenous firms 30 and households. 31  The aim of this paper was two-fold.First, we wanted to apply a solution procedure to solving asset pricing models that has passed an accuracy test.For this purpose we used a basic stochastic growth model and apply our dynamic programming algorithm to that basic model for which consumption, asset prices and the Sharpe-ratio can analytically be computed.We show, the errors of the equity premium and Sharpe-ratio depend on the discretization step size and grid refinement but not on the variance of the shocks.Our computations show that the optimal consumption, the value function and the Sharperatio can be computed with small absolute errors.We also show that the relative errors in computing the Sharpe-ratio can be small.The error in the computation of the Sharperatio is independent of the magnitude of the Sharpe-ratio.In other words the relative error decreases with increasing Sharpe-ratio. 32We also want to note that our method appears to more accurate than the equidistant spline method and the approximation methods using second order approximations of the policy function. 33  Next we study the asset pricing characteristics of a model with more promising extensions than the basic growth model by using our numerical method that appears to solve the asset pricing models with higher accuracy.As extended model we have chosen a model with more complex decision structure, a model with habit persistence, and augmented it, along the line of Boldrin et al (2001) and Jerman (1998), with adjustment costs of investment.Extensive research effort has recently been devoted to such a model variant in order to match the empirical characteristics of asset prices and returns.Our results, based on our dynamic programming algorithm shows that, even if habit persistence is jointly used with adjustment costs of investment, there are still puzzles remaining for consumption based asset pricing models. 34 Appendix: Details of the Gridding Strategy The basic idea of our adaptive gridding algorithm lies in evaluating the residual of the operator T applied to V Γ , as made precise in the following definition.(i) We define the a posteriori error estimate η as a continuous function η ∈ C(Ω, R) by (ii) For any element Q l of the grid Γ we define the elementwise error estimate We define the global error estimate η max by It is shown in Grüne (2003), that for this error estimate the inequalities holds.These inequalities show that the error estimate is reliable and efficient in the sense of numerical error estimator theory, which is extensively used in the numerical solution of partial differential equations.Furthermore, η(x) is continuous and one can show that a similar upper bound holds for the error in the derivative of V and V Γ .
If the size of a grid element tends to zero then also the corresponding error estimate tends to zero, even quadratically in the element size if V Γ satisfies a suitable "discrete C 2 " condition, i.e., a boundedness condition on the second difference quotient.
This observation shows that refining elements carrying large error estimates is a strategy that will eventually reduce the element error and consequently the global error, and thus forms the basis of the adaptive grid generation method which we will describe in the next section.
Clearly, in general the values η l = max x∈Q l η(x) can not be evaluated exactly since the maximization has to be performed over infinitely many points x ∈ Q l .Instead, we approximate η l by ηl = max where X T (Q l ) is a set of test points.In our numerical experiments we have used the test points indicated in Figure 9.1.The adaptive grid itself was implemented on a tree data structure in the programming language C. The adaptive refinement follows the standard practice in numerical schemes and works as follows: (0) Choose an initial grid Γ 0 , set i = 0, fix a refinement threshold θ ∈ (0, 1) (1) Compute V Γ i and the (approximated) error estimates ηl and ηmax .If a desired accuracy or a maximally allowed number of nodes is reached, then stop (2) Refine all elements Q l with ηl ≥ θη max , denote the new grid by Γ i+1 (3) Set i := i + 1 and go to (1) Here for the solution of V Γ i for i ≥ 1 we use the previous solution V Γ i−1 as the initial value for the iteration described in Section 3, which turns out to be very efficient.
During the adaptation routine it might happen that the error estimate causes refinements in regions which later turn out to be very regular.It is therefore advisable to include a coarsening mechanism in the above iteration.This mechanism can, e.g., be controlled by comparing the approximation V Γ i with its projection π Γ i V Γ i onto the grid Γ i which is obtained from Γ i by coarsening each element once.Using a specified coarsening tolerance tol ≥ 0 one can add the following step after Step (2).
(2a) Coarsen all elements Q l with ηl < θη max and This procedure also allows to start from rather fine initial grids Γ 0 , which have the advantage of yielding a good approximation ηl of η l .Unnecessarily fine elements in the initial grids will this way be coarsened afterwards.
In addition, it might be desirable to add additional refinements in order to avoid large differences in size between adjacent elements, e.g., to avoid degeneracies.Such regularization steps could be included as a step (2b) after the error based refinement and coarsening has been performed.In our implementation such a criterion was used; there the difference in refinement levels between two adjacent elements was restricted to at most one.Note that the values in the hanging nodes (these are the nodes appearing at the interface between two elements of different refinement level) have to be determined by interpolation in order to ensure continuity of V Γ .
In addition, our algorithm allows for the anisotropic refinement of elements: consider an element Q of Γ (we drop the indices for notational convenience) and let X new,i be the set of potential new nodes which would be added to Γ if the element Q l was refined in coordinate direction e i , see Figure 9.2.Define the error estimate in these nodes for each coordinate direction e i by η dir,i := max x∈X new,i η(x) and define the overall error measured in these potential new nodes by η dir := max i=1,...,n η dir,i .Note that η dir ≤ η l always holds.If we include all the points in X new := i=1,...,n X new,i in our set of test points X T (Q) (which is reasonable because in order to compute η dir,i we have to evaluate η(x) for x ∈ X new , anyway) then we can also ensure η dir ≤ ηl .Now we refine the element only in those directions for which the corresponding test points yield large values, i.e., if the error estimate η dir,1 is large we refine in x-direction and if the error estimate η dir,2 is large we refine in y-directions (and, of course, we refine in both directions if all test points have large error estimates).
Anisotropic refinement can considerably increase the efficiency of the adaptive gridding strategy, in particular if the solution V has certain anisotropic properties, e.g., if V is linear or almost linear in one coordinate direction.Note that this is the case in our example.On the other hand, a very anisotropic grid Γ can cause degeneracy of the function V Γ like, e.g., large Lipschitz constants or large (discrete) curvature even if V is regular, which might slow down the convergence.However, according to our numerical experience the positive effects of anisotropic grids are usually predominant.

Figure 5 . 4 :
Figure 5.4: The risk free interest rate R f (x) and a plane at level R f = 1

Figure 9 . 1 :
Figure 9.1: Test points X T (Q l ) for a 2d element Q l

Figure 9 . 2 :
Figure 9.2: Potential new nodes X new,1 (left) and X new,2 (right) for a 2d element Q

Table 5 .
1: Approximation errors for the fully numerical approach

Table 5 . 2
: Approximation errors using exact c in the grid points

Table 5 .
3: Approximation for fully numerical approach and different σ al. (2001) and Jerman (1998).Whereas Boldrin et al. use a model with log utility for internal habit, but endogenous labor supply in the household's preferences, Jerman studies the asset price implication of a stochastic growth model, also with internal habit formation but, as in our model, labor effort is not a choice variable.All three papers Boldrin et al. (2001), Jerman (1998) and our variant use adjustment costs of investment in the model with habit formation.Both previous studies claim that habit formation models with adjustment costs come closer to match the financial characteristics of the data.
29We have taken a low value of β in order to make our study comparable to previous studies solving the stochastic growth model with dynamic programming.A too high risk free interest rate can always be reduced by a higher β, seeCochrane (2001, ch.21)andCampbell, Lo and MacKinley (1997, ch.82.).