Using dynamic programming with adaptive grid scheme for optimal control problems in economics

Abstract The study of the solutions of dynamic models with optimizing agents has often been limited by a lack of available analytical techniques to explicitly find the global solution paths. On the other hand, the application of numerical techniques such as dynamic programming to find the solution in interesting regions of the state was restricted by the use of fixed grid size techniques. Following Grune (Numer. Math. 75 (3) (1997) 319; University of Bayreuth, submitted, 2003), in this paper an adaptive grid scheme is used for finding the global solutions of discrete time Hamilton–Jacobi–Bellman equations. Local error estimates are established and an adapting iteration for the discretization of the state space is developed. The advantage of the use of adaptive grid scheme is demonstrated by computing the solutions of one- and two-dimensional economic models which exhibit steep curvature, complicated dynamics due to multiple equilibria, thresholds (Skiba sets) separating domains of attraction and periodic solutions. We consider deterministic and stochastic model variants. The studied examples are from economic growth, investment theory, environmental and resource economics.


Introduction
In recent times the lack of closed form solutions of dynamic models with optimizing agents has led to the use of computational methods to solve those models. Closed form solutions are available for special cases such as the linear-quadratic control problem, see Ljungqvist and Sargent (2001, ch. 4) and growth models with log-utility. While those special models suffice for the analytical study of a variety of economic problems in more general cases they are not sufficient for a robust analysis of interesting economic problems. For the more general cases numerical methods have been proposed in the literature. A detailed discussion of a variety of numerical methods and a comparison with our proposed method is provided in Section 3 below.
Our paper is concerned with a family of continuous and discrete time dynamic models with optimizing agents whose solutions are most conveniently approached by the method of dynamic programming. Dynamic programming provides the value function and the control variable in feedback form. Yet, the application of numerical methods such as dynamic programming to find the global dynamics in interesting regions of the state space were restricted by the use of fixed grid size techniques. Following Grüne (1997) in this paper an adaptive grid scheme is used for finding global solutions of models with dynamic optimization. As Santos and Vigo-Aguiar (1998), we use numerical value function iteration but we employ local error estimates based on a flexible grid scheme. Since those numerical methods provide us with approximate solutions only, it is essential to have accuracy estimates for the numerical methods employed.
We consider discounted continuous time and discrete time optimal control problems and flexible grid scheme based on local error estimates. The advantage of the use of flexible grid scheme is demonstrated by computing the value function and the control variable in feedback form of one and two dimensional economic models. In order to study the accuracy of our numerical methods when applied to economic models we first want to test our algorithm by studying a basic growth model of Brock and Mirman (1972) type for which the exact solution is known. This allows us to judge the accuracy of our numerical method and to explore whether the flexible grid scheme is accurately capturing important solution properties, for example, steep curvature of the value function. Our algorithm will be applied to a deterministic as well as to a stochastic version of the Brock and Mirman (1972) growth model. Moreover in the economic literature there exist more complicated dynamic models with optimizing agents which have been a challenge to commonly used numerical techniques. Those models exhibit more complicated dynamics due to the existence of multiple steady state equilibria, 1 thresholds (Skiba-sets) separating domains of attraction and periodic solutions as attractors. Examples of such models can be found in the literature on economic growth and development. 2 Multiple steady states and thresholds can also arise in the 1 If there are local attractors among the equilibria some authors characterize them as indeterminant equilibria. For a survey of models with indeterminacy, see Benhabib and Farmer (1999). 2 In the latter type of models a convex-concave production function arises which leads to thresholds separating paths to low per capita income (poor) countries and high per capita income (rich) countries, see Skiba (1978) and Azariadis and Drazen (1990). In more recent growth models of the Lucas or Romer type multiple steady state growth paths may arise due to externalities and complementarities of inputs, or monopolistic competition, see for example, Benhabib, Perli and Xie (1994) and . dynamic decision problem of the firm, for example due to relative adjustment costs of investment, 3 in resource economics and in ecological management problems. 4 Our paper studies a prototype model from each of those areas and applies the proposed dynamic programming algorithm with adaptive grid scheme to find the global dynamics.
The remainder of the paper is organized as follows. Section 2 describes the basic discretization procedure and its strategy of adaptive grid refinement, while Section 3 compares our approach to other advanced dynamic programming schemes found in the literature. Section 4 solves one dimensional control problems. Here we study the basic growth model for which the exact solution is known and compare the accuracy of our numerical methods with previously used methods. We then study one dimensional dynamic optimization models with multiple domains of attraction and thresholds. Section 5 then presents results for two dimensional examples which also can give rise to thresholds and periodic solutions as attractors. Section 6 concludes the paper.

Numerical Dynamic Programming
In this section we describe a dynamic programming algorithm which enables us to compute optimal value functions as well as optimal trajectories of discounted optimal control problems of the below type (2.1)-(2.2) or (2.3)-(2.4).
After stating the problem formulation we will briefly describe the basic discretization procedure which goes back to Capuzzo Dolcetta (1983) and Falcone (1987) (see also Capuzzo Dolcetta and Falcone (1989) and Bardi and Capuzzo Dolcetta (1997), Appendix A), then discuss the adaptive gridding strategy from Grüne (1997) and, finally, give some hints on implementational aspects.

Problem Formulation
We consider discounted optimal control problems, either given in continuous time or in discrete time t ∈ N 0 given by Our algorithm is easily extended to stochastic optimal control problems, see Section 2.5, below, for details.

Discretization
We start describing the basic discretization technique, which is done using a two step or semi-Lagrangian discretization, first in time and then in space. In the first step, which is only necessary for continuous time problems (2.1)-(2.2), the continuous time optimal control problem is replaced by a first order discrete time approximation given by where β = 1 − δh and x h is defined by the discrete dynamics and h > 0 is the discretization time step. Here U d denotes the set of discrete control sequences u = (u 1 , u 2 , . . .) for u i ∈ U . If the original problem is of type (2.3)-(2.4), then it is already in the form (2.5)-(2.6) with h = 1.
The optimal value function V h is the unique solution of the discrete Hamilton-Jacobi-Bellman equation If we define the dynamic programming operator T h by then V h can be characterized as the unique solution of the fixed point equation The second step of the algorithm now approximates the solution on a grid Γ covering a compact subset Ω of the state space. Note that optimal control problems in economics are usually defined on a non-compact state space and the reduction to a compact domain Ω is a nontrivial task. Typically, we pick a reasonable set Ω and consider only those trajectories which remain in Ω for all future times. Often it is possible to identify sets Ω which are optimally invariant, i.e., which have the property that any optimal trajectory of the unconstrained problem which starts in Ω remains in Ω for all future times. In this case, the restriction to Ω does not change the optimal value functions, otherwise one has to interpret more carefully the final results. In any case, we assume that for any point x ∈ Ω there exists at least one control value u such that x + hf (x, u) ∈ Ω holds.
Denoting the nodes of the grid Γ by x i , i = 1, . . . , P , we are now looking for an approxi- for all nodes x i of the grid, where the value of V Γ h for points x which are not grid points (these are needed for the evaluation of T h ) is determined by interpolation. In our implementation we use rectangular grids with bilinear interpolation in 2d and linear interpolation in 1d.
Note that an approximately optimal control law (in feedback form for the discrete dynamics) can be obtained from this approximation by taking the value u * (x) = u for u realizing the maximum in (2.7), where V h is replaced by V Γ h . This procedure in particular allows the numerical computation of approximately optimal trajectories. For a rigorous convergence analysis of this discretization scheme we refer to Bardi and Capuzzo Dolcetta (1997) and Falcone and Giorgi (1999). Here we only remark that the scheme converges at least with order γ, where γ > 0 depends on the continuity properties of V and V h . In particular, we have γ = 1/2 for continuous time problems and γ = 1 for discrete time problems if V and V h are Lipschitz continuous and γ = 2 for discrete time problems if V h is twice continuously differentiable.

Adaptive Grids
In order to distribute the nodes of the grid efficiently, we make use of a posteriori error estimation. The idea here lies in the fact, that we want to solve the equation V h (x) = T h (V h )(x) for all x ∈ Ω but we do only solve this equation for the nodes x i of the grid. Therefore, for the estimation of the gridding error we estimate the residual of the dynamic programming operator, i.e., the difference between V Γ h (x) and T h (V Γ h )(x) for points which are not nodes of the grid. Thus, for each cell C l of the grid Γ we compute local error estimates η l := max It was shown in Grüne (1997) that the error estimators η l give upper and lower bounds for the real error. More precisely, one can prove the estimates Thus these local error estimates fit into the usual theory for error estimates for PDEs and can be used as indicators for the following local refinement strategy of the grid Γ: (0) Pick an initial grid Γ 0 and set i = 0. Fix a refinement parameter θ ∈ (0, 1) and a tolerance rtol > 0.
(1) Compute the solution V Γ i h on Γ i .
(2) Evaluate the error estimates η l . If η l < rtol for all l then stop.
(3) Refine all cells C j with η j ≥ θ max l η l , set i = i + 1 and go to (1).
As an alternative to specifying rtol one could also stop the refinement iteration when a prescribed maximum number of elements is reached, which is what we did in the numerical experiments in this paper. In all these experiments we used θ = 0.1.
It should be noted that both the upper and the lower bound in (2.11) are necessary for an efficient and reliable adaptation routine: the upper bound ensures that small local errors indeed imply small global errors, which means that refining elements with large local errors -thus reducing the local error -will be the right strategy to reduce the global error (the error estimate is reliable). On the other hand, the lower bound ensures that large local errors do indeed indicate that there is something wrong with the numerical solution, which means that we do only refine if there really are large errors (the estimate is efficient in the sense that we do not make useless refinements). Note that the gap between the upper and lower bound becomes large if β ≈ 1; in this case numerical experience shows that the upper bound is likely to become too conservative in the sense that the estimated errors typically overestimate the real error. However, since the lower bound is independent of β, even in the case β ≈ 1 -which causes problems in many numerical schemes -the efficiency of the adaptive gridding method is maintained, i.e., we do not make unnecessary refinements.
In practice, the error estimates η l cannot be evaluated exactly, because it involves maximization over infinitely many points in the element C l . Instead, we use a number of test points x t 1 , . . . , x t m in each element and evaluateη l = max η(x t i ). Figure 2.1 shows the test points we have used in the 2d examples. Note that this procedure may underestimate the real error estimate on very large elements, thus it is advisable not to choose a too coarse initial grid Γ 0 . This choice of test points naturally leads to the idea of anisotropic refinement of elements, i.e., the refinement of C l not in all but only in selected coordinate directions, thus allowing to use "flat" or "stretched" elements. First note that for each rectangular element in R n we have several possibilities of refinement, as indicated in Figure 2.2 for the 2d case. The simplest choice would be to refine each element in all possible coordinate directions, like in the rightmost case in Figure 2.2. However, for an element C l one can also consider the sets X new,i of potential new nodes which would be added to Γ if the element C l was refined in coordinate direction e i . Figure 2.3 shows these points in 2d. 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. Since the points in X new,i are contained in our selection of test points, anyway, we get these values for free. Now, if the error in one direction is large compared to the overall errorη l , then we only refine the element in the corresponding direction. We have used this strategy in all 2d examples presented in Section 5, but the effect is most apparent in Example 5.1, because here we have a strong curvature in e 1 -direction but no curvature in e 2 -direction.
During the adaptation routine it might happen that the error estimate on a coarser grid causes refinements in regions which after a few more refinement iterations 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 V Γ i onto a "thinned" grid Γ i which is obtained from Γ i by coarsening each element once. Using a specified coarsening tolerance ctol ≥ 0 one can add the following step after Step (2).
(2a) Coarsen all elements C l with η j < θ max l η l and max x∈C l |V 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. We have used this coarsening mechanism in the Examples 5.1 and 5.2, below, with its effect being most apparent in Example 5.2 where the final grid consists of fewer nodes than the initial grid.

Some Implementational Details
In this section we describe some implementational issues for our computations. We only outline the most important facts and refer to the appropriate literature for details. All numerical examples described in this paper were computed using an implementation of our algorithm in C on a Pentium III-PC with 997 MHz and 384 MByte RAM running under a SuSE LINUX system.
For the computation of V Γ h with time step h > 0 on a grid Γ with P nodes x 1 , . . . , x P , we identify the piecewise (bi)linear function For the computation of V one could use the fixed point iteration starting, e.g., from the initial vector V 0 = (0, . . . , 0) T . This procedure is known to converge to the right solution, however, the convergence can be rather slow. A first method to speed up the convergence is obtained if we replace the iteration V j+1 = T h (V j ) by a Gauss-Seidel type iteration V j+1 = S(V j ) which is described in Grüne (1997, Sect. 3) under the name increasing coordinate algorithm.
In our implementation we use a straightforward discretization of U by a finite set U , for which the local maximization can be performed by simply comparing the finitely many values obtained for the finitely many control values from U . Note that other implementations use different strategies for this local optimization, e.g. in Santos and Vigo-Aguiar (1998) and Carlini, Falcone and Ferretti (2002) Brent's Algorithm is used here; both references report good results for this optimization strategy. The advantage of using a finite set of control values U lies in the fact that one can easily couple this value-space iteration with a policy-space iteration as described in Gonzalez and Sagastizabal (1990) or Seek (1997): Once the control values in which the maxima for V i and V i+1 are attained become convergent, we fix these values and solve the corresponding system of linear equations for these fixed control values using a fast iterative solver like the CGS or BICGSTAB method. After convergence of this method we continue with the usual iteration using S until the control values again converge, switch to a linear solver and so on. This combined policy-valuespace iteration turns out to be much more efficient than the plain value space iteration using S, often more than 90 percent faster. We also take advantage of the fact that the solution on the grid Γ i can be interpolated in the new nodes of Γ i+1 and then be used as the initial value for the iteration on Γ i+1 .
The adaptive rectangular grids were stored on a hierarchical tree data structure as described in Grüne, Metscher and Ohlberger (1999). In particular, this allows the compact storage of the complete adaptively refined grid without storing the coordinates of the nodes of the grid; these values are reconstructed when needed while traversing through the data structure. For a given point y ∈ Ω this structure also allows a fast location of the grid element containing y which in turn leads to an efficient computation of the value V Γ h (ϕ(x, u)) which is needed in the evaluation of the operator T h .

Extension to stochastic optimal control problems
Our adaptive approach is easily extended to stochastic discrete time problems of the type and the z t are i.i.d. random variables. Again, this problem could be a priori in discrete time (in which case we have h = 1) or it could be the time discretization of a continuous time stochastic optimal control problem with dynamics governed by an Itô stochastic differential equation, see Camilli and Falcone (1995).
The corresponding dynamic programming operator becomes (2.14) where ϕ(x, u, z) is now a random variable.
If the random variable z is discrete then the evaluation of the expectation E is a simple summation, if z is a continuous random variable then we can compute E via a numerical quadrature formula for the approximation of the integral where p(z) is the probability density of z. Example 5.1 shows the application of our method to such a problem, where z is a truncated Gaussian random variable and the numerical integration was done via the trapezoidal rule.
It should be noted that despite the formal similarity, stochastic optimal control problems have several features different from deterministic ones. First, complicated dynamical behavior like multiple stable steady state equilibria, periodic attractors etc. is less likely because the influence of the stochastic term tends to "smear out" the dynamics in such a way that these phenomena disappear. 5 Furthermore, in stochastic problems the optimal value function typically has more regularity which allows the use of high order approximation techniques, see the following section for details. Finally, stochastic problems can often be formulated in terms of Markov decision problems with continuous transition probability (see Rust (1996) for details), whose structure gives rise to different approximation techniques, in particular allowing to avoid the discretization of the state space.
In these situations, our technique may not be the most efficient approach to these problems, and it has to compete with other efficient techniques, of which some are outlined in the following section. Nevertheless, Example 5.1 shows that adaptive grids are by far more efficient than non-adaptive methods if the same discretization technique (i.e., bilinear approximation in this example) is used for both approaches. It should also be noted that in the smooth case one can obtain estimates for the error in the approximation of the gradient of V h from our error estimates, for details we refer to Grüne (2003).

Comparison with other dynamic programming schemes
In the literature one can find a vast amount of different approaches to dynamic programming problems, many of them using state-of-the art mathematical and numerical techniques for making this appealing computational technique more efficient. In this section we will briefly review some of them and highlight the differences to our approach as described in the preceding section.
Before describing these different approaches, let us recall that the main purpose of our study is to compute the dynamics of the optimally controlled systems, in particular when complicated dynamic behavior like different domains of attractions separated by Skiba-Sets, periodic attractors and other phenomena occur. These phenomena occur in deterministic problems and their occurrence often comes along with a loss of regularity in the related value functions, which are at best Lipschitz but typically nonsmooth.
The majority of numerical dynamic programming approaches, however, do either need sufficient smoothness of the optimal value functions or stochastic problems with a sufficiently rich stochastic influence like in Markov decision problems with at least continuous transition probabilities. 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), Zin (1993, 1997)) or piecewise high-order approximations (Falcone and Ferretti (1998)) can be very efficient, see also the survey of Taylor and Uhlig (1990) where one can find a comparative numerical study of several methods.
Some of these methods (like Spline and piecewise high order approximation) use a grid discretization of the state space similar to our approach and it would be tempting to try whether adaptive discretization ideas based on our local error estimation technique work equally well with these approximation techniques; so far we are not aware of numerical studies in this direction. Again, we would like to recall that this approach needs appropriate smoothness which is not present in most of the examples we are interested in. Smoothness is also needed in other strategies, like in finite difference approximations (Candler (2001)), Gaussian Quadrature discretization (Tauchen and Hussey (1991), Burnside (2001)) and in perturbation techniques (Judd (1996)). While the last should also work if the value function is only piecewise smooth (which is the case in all of our examples), however, this technique computes the optimal value function starting from the stable equilibria and is typically less accurate away from these, which suggests that this may not be the best approach to detect Skiba points or curves, which are located far away from the stable equilibria, furthermore, the approach does not work if there are more complicated stable attractors like, e.g., periodic solutions as in Example 5.3, below.
A crucial aspect of any numerical dynamic programming method approximation is the curse of dimensionality, which states that the numerical effort (induced by the number of grid points) grows exponentially in the dimension of the state space (see Rust (1996) for a comprehensive account on complexity issues). Our adaptive approach will not be able to break the curse of dimensionality, but it helps to considerably reduce the computational cost for any fixed dimension. There are, however, methods that can break the curse of dimensionality by using either randomly distributed grid points (Rust (1997)) or so called low discrepancy grids (Rust (1996), Reiter (1999)). Even though these results are very appealing, these methods are unfortunately not applicable here because they either need Markov Decision Problems with (at least) continuous transition probability like in Rust (1996), which allows to avoid interpolation, or again sufficient smoothness of the optimal value function like in Reiter (1999) where interpolation is replaced by a regression technique. 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.
In this context it should be noted that from the complexity point of view for discretization techniques 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, see Tables 4.1 and 4.2, below.
Summarizing this first part of our comparative discussion, for smooth and "rich" stochastic problems a large amount of sophisticated and efficient techniques have been developed, which do, however, not apply to deterministic or nonsmooth problems. It would be very interesting to see how these methods compare to our adaptive approach in the smooth case, even more, it would be interesting to see whether the ideas of smooth approximation and adaptive discretization can be coupled in order to combine the strong sides of both approaches, but this is beyond the scope of this paper and will be the topic of future research.
As for deterministic and nonsmooth problems, there are still a number of advanced techniques reported in the literature which do all make use of adaptivity in one way or the other. 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 like (2.11) are not given. 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 analysing 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 method, 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. Clearly, selecting a reasonable computational region is a serious problem at its own right, but for problems with multiple steady state equilibria and Skiba sets this approach is not adequate, because these sets have a very low probability (all trajectories move away) but are nevertheless very interesting for understanding the global dynamics of the problem. Furthermore, again this is a heuristic approach which needs to be complemented by a rigorous mathematical foundation. 6 Probably the adaptive approaches with the most solid mathematical background are presented in the papers of Zin (1993, 1997). 7 In 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 upper bound in (2.11), 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 error estimate needs only one step of this trajectory.
Summarizing the second part of our discussion, there are a number of adaptive strategies around which are all reported to show good results, however, they are either heuristic and better suited for other classes of problems (like problems with discontinuities or problems with only one stable equilibrium), or they have nice theoretical features but are practically inconvenient because their implementation is numerically much more expensive than our approach.
Let us finally 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 was used to estimate the quality of the approximating solution, however, to our knowledge it has not been used to control adaptive gridding strategies, and we are not aware of any estimates resembling (2.11), which is a crucial property for an efficient and reliable adaptive gridding scheme.

One Dimensional Examples
In this section we describe the application of our algorithm to a selection of one dimensional optimal control problems. For some one dimensional models the use of adaptive gridding strategies is not really necessary, because these problems can usually be solved with high about the quality of the final solution of a heuristic method. 7 As mentioned above, this approach also uses splines, i.e., a smooth approximation, but the ideas developed in these papers do also work for linear splines which do not require smoothness of the approximated optimal value function.
precision on fixed equidistant grids in a reasonable amount of time and with only low memory consumption. Nevertheless, even for these (numerically) simple problems one can show that adaptive gridding techniques can be much more efficient than algorithms on equidistant grids when there are kinks or steep curvature in the value function. Our first example will illustrate this, because for this problem the exact solution is known, hence the accuracy of the respective methods can be compared directly. For the subsequent examples the exact solution is not known, hence numerical methods are necessary for their analysis.

A basic growth model with explicit solution
We start our numerical investigations with a basic growth model in discrete time, which goes back to Brock and Mirman (1972). This model has been used as a test example for many numerical algorithms, and our method can be regarded as an adaptive version of the numerical analysis performed by Vigo-Aguiar (1995 and1998, Sect. 4). The problem is a discrete time maximization problem of type (2.3)-(2.4) with return function and dynamics given by with constants A > 0 and 0 < α < 1. The exact solution to this problem is known (see Santos and Vigo-Aguiar (1998)) and is given by For our numerical tests we specify A = 5, α = 0.34 and the discount factor to β = 0.95. For these parameters the point x 1 = 2.06734 is the unique optimal equilibrium. The corresponding optimal value function is shown in Figure 4.1.  In Santos and Vigo-Aguiar (1998), this problem was solved with a similar algorithm as proposed here using equidistant grids. In our first computation we have also computed the solution on grids with an equidistant distribution of the nodes. We have computed the solution on the interval [0.1, 10] with control range U = [0.1, 5], whose discretization U consisted of 501 equidistant values. Table 4.1 shows the respective results, the CPU time needed for the computation and also the maximal error (i.e., the error in the L ∞ norm) between the numerical and the exact solution. The achieved accuracy is comparable with the results from Santos and Vigo-Aguiar (1998) (small differences are due to the different local maximization strategy, see Sect. 2.4), the CPU times are, of course, not, because the computations here were performed with a different implementation on a different computer.
# nodes Error CPU time 99 3.2 · 10 −2 1 s 989 6.3 · 10 −4 9 s 9899 5.0 · 10 −5 55 s  Table 4.2 shows the respective results for the adaptive grid algorithm with refinement parameter θ = 0.1. The CPU time in the last row is the overall time for the computation on all levels including all the time needed for the computation of the error estimates and the grid refinements. The estimated error is the error bound obtained from the upper inequality in (2.11), i.e. max η l /(1 − β), in this example 20 · max η l .

A growth model with credit market
Next, we investigate a growth model from Grüne, Semmler and Sieveking (2001) which admits borrowing from the credit market for investment but where the maximal level of debt is given by the present value of an income stream from economic activity. Under the assumption of unrestricted access to capital market, a country, for example, may face a present value borrowing constraint which we want to compute. The interest rate paid on the debt is constant. 8 The model leads to an optimal control problem of type (2.1)-(2.2) with a return function and dynamics given by where x denotes the capital stock, u is the investment and the discount rate δ in (2.1) is the interest rate (for the numerical treatment of non-constant credit cost see Grüne, Semmler and Sieveking (2001)).
We have solved this problem with constants σ = 0.15, A = 0.29, α = 1.1, β = 2, γ = 0.3 and δ = 0.1. The qualitative behavior of the numerical results, which were computed with h = 1/20 and 401 control values in U covering the interval U = [0, 1/4], is shown in the Figures 4.3-4.5. Unlike the previous example, given the parameter constellations of our model, this model exhibits multiple steady state equilibria. We obtain two optimal equilibria x 1 = 0 and x 2 = 0.996, with domains of attraction separated by the Skiba point x * = 0.267. Yet, as it is shown in Grüne, Semmler and Sieveking (2001) there are in fact three steady state equilibria, but the middle one in the vicinity of Skiba point x * = 0.267, is unstable and non-optimal. The kink in the value function, however, is barely visible in Figure 4.3, but the jump in the control variable in Figure 4.4 and also the local mesh refinement (here the final grid consists of 116 nodes) clearly indicate the presence of a Skiba point which represents a threshold below which it is optimal to deplete the capital stock and above which the capital stock and thus per capita income will grow. Figure 4.5 shows the mesh width for the final adaptive grid of our computed example. Also in this example we have evaluated the error bound from inequality (2.11) which is given in Table 4.3.

A 1d investment model with relative adjustment cost
Next, we investigate a continuous time decision problem of the firm where the firm's objective is to maximize, through the choice of an optimal investment strategy, the discounted value of its cash flows subject to some technology constraint. This type of decision problem which also exhibits multiple steady state equilibria may also result, as the example 4.2, in a discontinuous policy function and thus in a discontinuous investment strategy. A model of this type has been introduced, although not numerically solved, by Hartl et al. (2000). The return function and dynamics are given by respectively, where x is the capital stock and ux is the investment (here we have performed the change of variables u = u/x on the model from Hartl et al. (2000) which removes the singularity of the return function g in x = 0).
As the previous example, this model also exhibits multiple stable optimal equilibria for suitable parameters, whose domains of attraction are separated by a Skiba point, see Deissenberg et al. (2001) for a detailed explanation of this phenomenon. Here we have computed the optimal value function for σ = 0.1, a = 1, b = 0.6, c = 0, d = 0.3 and discount rate δ = 0.2. Figure 4.6 shows the corresponding optimal value function computed with time step h = 1/20 and 501 equidistant control values in U discretizing the interval U = [0, 1]. With these parameters, the optimal equilibria are x 1 = 0 and x 2 = 0.565. Yet, the Skiba point in between at x * = 0.010 is barely visible as a kink in the optimal value function in Figure 4.6. The Skiba point is clearly visible as a jump in the optimal feedback control law in Figure 4.7 where the optimal feedback law is shown.
Finally, the Skiba point also shows up in the mesh width of the final adaptive grid, which in this example consists of 399 nodes. The kink in the optimal value function produces large local error, which in turn result in a fine grid around this point. This effect is visible

An ecological management problem
Our final 1d example is an ecological management problem from Brock and Starret (1999). An extensive numerical study of this model was performed in Grüne, Kato and Semmler (2002). We here show a numerical computation with a specific parameter set.
The model is given by a return function and dynamics g(x, u) = au σ − kβx 2 and d dt In the above the return function g(x, u) represents a utility function with two components: the utility of an affector, obtaining a benefit letting phosphorus, u, running into a lake and the disutility of an enjoyer of the lake being negatively affected by the stock of phosphorus, x. The state equation denotes the dynamics of the stock of phosphorus in the lake. We specify the parameters a = 2, σ = β = 1 2 , m = n = 1, ρ = 2, µ = 0.55 and discount rate δ  Table 4.5 again gives the corresponding upper error bounds from (2.11).

Two Dimensional Examples
Next, we describe three two-dimensional problems which we have solved using our algorithm. We first study our above basic growth model of Sect. 4.1, but we employ an extension to a stochastic variant. This gives rise to two state variables. Next we turn to deterministic problems, where the study of the dynamics of two-dimensional problems with multiple steady state equilibria have been quite a challenge for research in economics # nodes estimated Error 101 3.5 · 10 −1 122 3.2 · 10 −1 133 3.4 · 10 −1 135 5.0 · 10 −2 225 1.2 · 10 −1 240 4.2 · 10 −2 since here one expects the separation of domains of attraction given not by threshold points but threshold lines (Skiba lines). Furthermore, two-dimensional models may give rise to periodic solutions as attractors. Examples of both cases are subsequently presented.

A 2d stochastic growth model
Our first 2d problem is a 2d version of the Brock and Mirman (1972) model from Example 4.1. Here the 1d model 4.1 is extended using a second variable modelling a stochastic shock. The model is given by the discrete time equations where A, α and ρ are real constants and the z(t) are i.i.d. random variables with zero mean. The return function in (2.3) is again g(x, u) = ln u.
In our numerical computations we used the parameter values A = 5, α = 0.34, ρ = 0.9 and β = 0.95. As for Example 4.1, the exact solution is known and given by We have computed the solution to this problem on the domain Ω = [0.1, 10] × [−0.32, 0.32]. The integral over the Gaussian variable z was approximated by a trapezoidal rule with 11 discrete values equidistributed in the interval [−0.032, 0.032] which ensures ϕ(x, u, z) ∈ Ω for x ∈ Ω and suitable u ∈ U = [0.5, 10.5]. For evaluating the maximum in T the set U was discretized with 161 points. Table 5.1 shows the results of the resulting adaptive gridding scheme applied with refinement threshold θ = 0.1 and coarsening tolerance ctol = 0.001. Figure 5.1 shows the resulting optimal value function and adapted grid.

A 2d investment model with relative adjustment cost
The following problem from Feichtinger et al. (2000) is a 2d variant of the 1d Example 4.3, where the control variable now is the change of investment rather than investment itself. It is given by Using the parameters σ = 0.25, k 1 = 2, k 2 = 0.0117, c 1 = 0.75, c 2 = 2.5, α = 12 and discount rate δ = 0.04 we have obtained the results shown in the following Figures 5.2-5.3, with time step h = 1/20 and 101 control values U covering U = [−1, 1]. Here we started from a rather fine grid and used the coarsening tolerance ctol = 0.001. As a result, the final grid is coarser than the initial grid but better adapted to the problem and thus provides a smaller error with fewer elements and grid nodes.
The final adaptive grid here consists of 14255 nodes. Note that the Skiba curve is visible in all three figures: as a discontinuity on the optimally controlled vector field ( Figure 5.2, left), as a jump in the optimal control ( Figure 5.2, right), as a kink in the optimal value function ( Figure 5.3, left), and as a line of refined elements in the adaptive grid ( Figure  5.3, right). This last figure in particular shows that the adaptive gridding strategy is very  Note that the estimated error values in Table 5.2 are very conservative here due to two reasons: first, the small denominator 1−β = 0.002 in the right inequality of estimate (2.11) causes a large gap between the upper and lower bound for the error, and our numerical experience shows that the lower bound is often the more reliable one. Second, and more important here, the large error values are concentrated along the Skiba curve, where large errors occur due to the fact that the position of the numerical curve is not precisely at its exact location. Away from this curve, the error estimates are much lower, and since trajectories do only move away from this unstable curve this means that also the error is much lower there. If we evaluate the error bound away from the Skiba curve we obtain an upper bound for the error of about 10 −1 and around the stable equilibria we obtain an even smaller bound of the order of 10 −3 .

A 2d model of interacting resources
As a last example we consider a model for the optimal exploitation of interacting resources from Semmler and Sieveking (1994) and analytically studied in Sieveking and Semmler (1997). It is given by the return function g(x, u) = 1 1 + x 1 u x 1 u + 1 1 + x 2 u x 2 u − u 2 and the dynamics d dt Here x 1 and x 2 are the number of two species (e.g. fish in the ocean) with intraspecific competition, where x 2 is the predator for the prey x 1 . The control variable u models the harvesting effort chosen by the economic agent.
We set the parameters to a 0 = 1.04, a 1 = 0.01, a 2 = 0.07, b 0 = 1.01, b 1 = 0.2, b 2 = 0.01 and the discount rate to δ = 5. We The final adaptive grid for this example consisted of 14208 nodes. The interesting fact in this system is that although the trajectories for constant harvesting effort do always tend to asymptotically stable equilibria, the optimal harvesting strategy here turns out to have an asymptotically stable periodic limit cycle. Figure 5.4 shows the vector field (left) the optimal control (right) and Figure 5.5 (left) the value function and the regions of grid refinement (right). Figure 5.6 shows this optimal periodic trajectory (as a solid line) along with two other optimal trajectories (as dotted lines) which tend to this periodic solution. 1.2 · 10 −2 2399 6.9 · 10 −3 3544 3.3 · 10 −3 6938 1.7 · 10 −3 14208 1.2 · 10 −3

Conclusions
The current paper uses dynamic programming with a flexible gridding strategy for the robust analysis of the dynamics of economic models with optimizing agents. The flexible gridding strategy permits to compute in models with multiple domains of attraction, threshold points (one dimensional examples) and threshold lines (two dimensional examples). The method is also well suited to accurately compute value functions with steep curvature or periodic solutions as attractors. Such solutions are computed in a robust way through local error analysis and grid refinements for prototype economic models of economic growth (deterministic or stochastic versions), dynamic decision problem of the firm and resource and ecological management problems. The global error estimate provided by the method turns out to yield a reliable upper bound for the real error, which, however, can be conservative for discount factor β ≈ 1, in the sense that the real error is overestimated.