_{1}

^{*}

The optimality conditions for macroeconomic problems with limited commitment often contain partial derivatives of the optimal value function, corresponding to the outside option. This paper contributes to the literature on recursive contracts by proposing an algorithm for approximating the gradient of the value function using simulation-based methods. Our method combines numerical solution and simulation of the model, Monte-Carlo integration and numerical differentiation. It does not suffer from the curse of dimensionality and is therefore convenient for models involving many state variables. The algorithm inherits the speed and accuracy limitations of the numerical solution method it relies on. Our accuracy analysis is limited to a few classical examples from macroeconomic literature.

The purpose of this paper is to propose a simple algorithm for computing partial derivatives of the optimal value function. Macroeconomic problems involving incentive compatibility constraints have received wide attention in the literature due to recent advances in dynamic optimization techniques (see [

This issue has been previously considered in [

In order to be able to use finite differences to approximate the gradient at a given point, one would need to know the values of the optimal value function at a certain set of points. Our algorithm obtains approximations of these values with arbitrary precision. Moreover, achieving this accuracy is feasible for all points in the state-space which have economic relevance.

The initial step of our algorithm involves obtaining numerical solution to a problem using a procedure which satisfies three criteria. First, it approximates some unknown function with flexible functional forms of finite elements. Second, it can deliver an accurate solution as the number of the finite elements in the function goes to infinity. Third and last, the resulting numerical solution must be such that it can be formulated as a set of policy functions approximated with flexible functional forms. The next step involves using Monte-Carlo integration in order to evaluate the conditional expectation of the discounted sum of future instantaneous utilities. The final step involves applying the method of finite differences to approximate the values of the partial derivatives of the value function.

The attractive features of the algorithm include its rather wide scope of applicability and simplicity of implementation. It can be used to study the questions of risk sharing under imperfect enforcement of contracts, as well as partnerships with limited commitment when several state variables appear in the model corresponding to the outside option. Such models may include habit formation preferences, several types of capital, or reputational co-state variables. The suggested method is computationally inexpensive. It does not suffer from the curse of dimensionality and therefore it is particularly convenient for models involving many state variables.

The rest of the paper is organized as follows. Section 2 discusses an example, where our algorithm proves to be useful. Section 3 sketches the idea behind the algorithm. Section 4 deals with implementation of the algorithm, while Section 5 compares it with some available alternatives. Section 6 concludes.

To fix ideas, we start with an example of a macroeconomic model where our computational algorithm proves to be useful. The key feature of this example is that solving it boils down to designing an optimal social contract which takes into account not only technological but also incentive and legal constraints. Our example illustrates the need for computing the gradient of the value function and its practical implementation. Furthermore, it shows that our algorithm is applicable to some widely used models, to which the method in [

Consider a model of international risk sharing, which distinguishes itself from the canonical model [

The planner’s problem is to choose the sequences of consumption { c i t } and investment { i i t } to maximize a weighted sum of utilities

max { c i t , i i t } E 0 ∑ i = 1 I λ i ∑ t = 0 ∞ β t u ( c i t , h i t ) (1)

subject to an aggregate feasibility constraint

∑ i = 1 I c i t + ∑ i = 1 I i i t = ∑ i = 1 I f ( k i t , θ i t ) , (2)

individual participation constraints for each i = 1 , ⋯ , I ,

E t ∑ j = 0 ∞ β j u ( c i t + j , h i t + j ) ≥ V i a ( k i t , h i t , θ i t ) , (3)

the equations of motion for the capital,

k i t + 1 = ( 1 − δ ) k i t + i i t , (4)

the laws of motion for habits,

h i t + 1 = h i t + λ ( c i t − h i t ) , (5)

and non-negativity constraints c i t , i i t ≥ 0 . We assume that productivity shocks θ t follow a first order stationary vector autoregressive process, and that the initial values for the state variables k i 0 , h i 0 , θ i 0 , and the initial non-negative weights, λ i , are given. In addition, the usual restrictions apply to the discount factor, β ∈ ( 0,1 ) , the capital depreciation rate, δ ∈ ( 0,1 ) , and the persistence of habits, λ ∈ ( 0,1 ) .

The outside option V i a ( k i t , h i t , θ i t ) in the participation constraint (3) represents the optimal value function corresponding to the autarkic environment.

max { c i t , i i t } t = 0 ∞ E 0 ∑ t = 0 ∞ β t u ( c i t , h i t ) (6)

subject to

c i t + i i t = f ( k i t , θ i t ) ,

k i t + 1 = f ( k i t , θ i t ) − c i t + ( 1 − δ ) k i t ,

h i t + 1 = h i t + λ ( c i t − h i t ) ,

with the initial values being equal to the values of the state variables k i t , h i t , θ i t at the moment of deviation from the optimal plan.

In addition to Equations (2)-(5), the optimal allocations, for all i , s = 1 , ⋯ , I , must satisfy the risk sharing condition,

Λ i , t Λ s , t = ξ s t ξ i t ,

where

Λ i , t = u c ( i , t ) + λ β E t ∑ j = 0 ∞ β j ( 1 − λ ) j [ ξ i t + 1 + j ξ i t u h ( i , t + 1 + j ) − μ i t + 1 + j ξ i t ∂ V i a ∂ h i t + 1 + j ( i , t + 1 + j ) ] , (7)

the intertemporal condition,

Λ i , t = β E t [ ξ i t + 1 + j ξ i t Λ i , t + 1 ( f k ( k i t + 1 , θ i t + 1 ) + 1 − δ ) − μ i t + 1 ξ i t ∂ V i a ∂ k i t + 1 ( i , t + 1 ) ] , (8)

the complementary slackness condition,

μ i t [ E t ∑ j = 0 ∞ β j u ( c i t + j , h i t + j ) − V i a ( k i t , h i t , θ i t ) ] = 0 ,

and the law of motion for the co-state variables M i t + 1 = M i t + μ i t , where ξ i t = λ i + M i t + 1 , μ i t ≥ 0 , and M i 0 = 0 . In the equations above u c ( i , t ) denotes

∂ u ( c i t , h i t ) ∂ c i t , andsimilar abbreviations apply to other terms1.

The gradient of the optimal value function V i a enters the intertemporal condition (8) and the risk-sharing condition (7). Approximation of this gradient is the purpose of the algorithm proposed in [

algorithms can approximate ∂ V i a ∂ k i t ( k i t , h i t , θ i t ) , we will use that fact to compare their accuracy in Section 1. Our approach can also approximate ∂ V i a ∂ h i t ( k i t , h i t , θ i t ) ,

whereas [

Typically, in the models with participation constraints the reservation value is the value function of the outside alternative, evaluated at the current values of the endogenous state variables, x ¯ , and exogenous shocks, s ¯ . Consider the optimal value function at a point ( x ¯ , s ¯ ) as an outcome of a standard optimization problem for the outside alternative:

V ( x ¯ , s ¯ ) = max { a t } E 0 ∑ t = 0 ∞ β t r ( x t , a t , s t ) (9)

subject to

x t + 1 = l ( x t , a t , s t ) , a t ∈ A ( x t , s t ) , (10)

x 0 = x ¯ , s 0 = s ¯ ,

where r is an instantaneous utility function, β ∈ ( 0,1 ) the discount factor, { s t } an exogenous Markov stochastic process, x t a vector of endogenous state variables, a t a vector of control variables, A a feasibility correspondence and l the law of motion for the endogenous state variables. The functional equation to this problem can be derived using the standard dynamic programming techniques. It yields a time invariant policy function f such that the optimal allocations satisfy a t = f ( x t , s t ) .

The purpose of our algorithm is to find a pointwise approximation to the partial derivative ∂ ∂ x i V ( x ¯ , s ¯ ) of the value function with respect to its i-th argument. The algorithm takes the following three steps:

Step I (Numerical Solution) Solve the model in (9) with a spectral method and formulate the solution in terms of approximated policy functions a t = f ^ ( ω ; x t , s t ) , which depend on the state variables and some coefficients, ω .

Step II (Monte Carlo Integration) Simulate N sequences of the realizations of the stochastic process { s t n } t = 1 T of size T with a starting value s 0 n = s ¯ , for all n = 1 , ⋯ , N . For a each sequence { s t n } t = 0 T , simulate the series of the endogenous variables { x t n , a t n } t = 0 T using approximated policy functions f ^ , the equations for motion for the state variables (10), and the initial values x 0 n = x ¯ . Using the simulated series calculate the discounted sums of the instantaneous returns and average over N:

V ( x ¯ , s ¯ ) ≃ 1 N ∑ n = 1 N ∑ t = 0 T β t r ( x t n , a t n , s t n ) .

Step III (Numerical Differentiation) Repeat Step II to obtain approximations of the value function at two points, for instance V ( x ¯ + ϵ ι i , s ¯ ) and V ( x ¯ − ϵ ι i , s ¯ ) , where ι i denotes a conformable vector of zeros with one on its i-th coordinate, and ϵ is a small positive number. Calculate the value of the partial derivative using, for example, Stirling’s finite difference formula:

∂ ∂ x i V ( x ¯ , s ¯ ) ≃ V ( x ¯ + ϵ ι i , s ¯ ) − V ( x ¯ − ϵ ι i , s ¯ ) 2 ϵ .

The optimal choice of the method for calculating the derivatives in Step III is problem specific and its accuracy depends on the smoothness of the value function. The approaches available include a variety of difference formulas, Richardson Extrapolation, or curve fitting with cubic splines. These are described at length in the standard numerical methods texts such as [

A brief note should be made at this point on the accuracy of the algorithm. In principle, arbitrary accuracy of the approximation can be achieved, by simultaneously increasing the dimension of the approximating family of functions in Step I, increasing the size of Monte Carlo iterations in Steps I and II, and decreasing the denominator ϵ in Step III. In practical applications, however, there are several sources of the approximation errors. First, in order to obtain the values of the optimal value function at a point, one relies on the approximations of the policy functions implied by the numerical solution to the model. Second, because we consider stochastic models, there is an additional error stemming from the evaluation of the integral in the computation of expected discounted returns. Finally, numerical differentiation introduces two more sources of error: the truncation error and the roundoff error. The truncation error comes from omitting higher order terms in the Taylor series expansion. The roundoff error is associated with storing real numbers in computer’s floating-point format. Section 0 discusses some practical accuracy issues in the context of an example.

This section describes a practical computational strategy for implementing the algorithm using the example from Section 2. The optimality conditions include partial derivatives of the value function corresponding to the dynamic programming formulation of the agents outside option, i.e. autarky. The functional equation for the autarkic problem is:

V ( k , h , θ ) = max ( c , i ) ∈ A ( k , θ ) { u ( c , h ) + β E [ V ( k ′ , h ′ , θ ′ ) | ( k , h , θ ) ] }

h ′ = h + λ ( c − h ) ,

k ′ = ( 1 − δ ) k + i ,

A ( k , θ ) = { ( c , i ) ∈ ℝ + 2 : c + i = f ( k , θ ) } .

The objective of the algorithm is to find the values of the partials V h ( ⋅ ) and V k ( ⋅ ) at a point ( k ¯ , h ¯ , θ ¯ ) , which is likely to happen in equilibrium. Since the analytical expression for these derivatives is in general unavailable, we have no choice but to rely on numerical differentiation. Another complication which arises here is that the closed form solution to the optimal value function is generally unavailable too. Hence, one needs to approximate value function at two points, e.g. V ( k ¯ + ε , h ¯ , θ ¯ ) and V ( k ¯ − ε , h ¯ , θ ¯ ) with arbitrary accuracy in order to be able to use the finite differencing approach.

The first step of the algorithm involves solving the model with a spectral method which can approximate the policy functions with arbitrary accuracy. This example relies on a version of stochastic simulation algorithm, which formulates the solution in terms of approximated policy functions. The Euler equation for the problem is given by:

Λ t = β E t [ Λ t + 1 ( f k ( k t + 1 , θ t + 1 ) + 1 − δ ) ] ,

where marginal utility of consumption is

Λ t = u c ( c t , h t ) + β λ E t [ ∑ j = 0 ∞ β j ( 1 − λ ) j u h ( c t + j + 1 , h t + j + 1 ) ] .

To simplify the exposition we will consider the case of non-persistent habits, which corresponds to λ = 1 . We assume the functional forms standard in the growth literature. The instantaneous utility function is given by

u ( c t , h t ) = ( c t − b h t ) 1 − σ 1 − σ , where b ∈ ( 0,1 ) and σ > 0 . The production function

is Cobb-Douglas and is given by f ( k t , θ t ) = θ t k t α . The stochastic process for productivity is log θ t = ρ log θ t − 1 + ε t , where ρ ∈ ( 0,1 ) , and { ε t } are independent normally distributed random variables with zero mean and variance σ ε 2 . In this example, we restrict attention to one particular set of the parameters which are summarized in

The algorithm follows the three steps: 1) numerical solution; 2) Monte Carlo integration; 3) numerical differentiation.

Step I (Numerical Solution) The sequences of optimal allocations { c t , h t + 1 , k t + 1 } t = 0 ∞ must satisfy the following system of stochastic difference equations:

( c t − b h t ) − σ = β E t [ b ( c t + 1 − b h t + 1 ) − σ × ( 1 + ( α θ t + 1 k t + 1 α − 1 + 1 − δ ) ( 1 b − β ( c t + 2 − b h t + 2 c t + 1 − b h t + 1 ) − σ ) ) ] , (11)

k t + 1 = θ t k t α + ( 1 − δ ) k t − c t , (12)

h t + 1 = c t . (13)

For the expositional purpose, we solve the model with a version of a stochastic simulation algorithm, which is easiest to implement (see e.g. [

1) Fix the initial conditions and draw a series of { θ t } t = 1 T that obeys the law of motion for the exogenous shocks with T sufficiently large. To ensure sufficient accuracy of solution we chose T = 50000 for all the numerical examples considered. The computational burden of this is still rather low since the model needs to be solved only once.

2) Substitute the conditional expectations in (11) with the flexible functional forms that depend on the state variables k t , h t , θ t and some coefficients, ω , to yield

Preferences | Technology | |||||
---|---|---|---|---|---|---|

σ | β | b | α | δ | ρ | σ ε |

3 | 0.95 | 0.5 | 0.36 | 0.06 | 0.95 | 0.007 |

( c t ( ω ) − b h t ( ω ) ) − σ = β ψ ( ω ; k t ( ω ) , h t ( ω ) , θ t ) ,

where ψ ( ω ; k t ( ω ) , h t ( ω ) , θ t ) = exp ( P n ( ω ; log k t ( ω ) , log h t ( ω ) , log θ t ) ) , and P n denotes polynomial of degree n. By using the exponent of the logarithmic polynomial expansion we guarantee that the left hand side of (11) remains positive. Given c t ( ω ) , the next period values for the capital and habit stocks follow directly from the laws of motion (12) and (13).

3) Using the realizations of { θ t } t = 0 T , repeat the previous step in order to obtain recursively a series of the endogenous variables { c t ( ω ) , k t + 1 ( ω ) , h t + 1 ( ω ) } t = 0 T , for this particular parameterization ω .

4) Run the following non-linear regression

Y t ( ω ) = exp ( P n ( ξ ; log k t ( ω ) , log h t ( ω ) , log θ t ) ) + η t ,

where the role of the dependent variable Y t ( ω ) is performed by the expression inside the conditional expectation in Equation (11).

5) Letting S ( ω ) be the result of the regression in the previous step, use an iterative procedure to find the fixed point of S, and the set of coefficients ω f = S ( ω f ) . This would provide the solution for the endogenous variables { c t ( ω f ) , k t + 1 ( ω f ) , h t + 1 ( ω f ) } t = 0 T for this particular realization of the stochastic process { θ t } t = 1 T along with the approximated policy functions:

c t ( k t , h t , θ t ) = b h t + [ β ψ ( ω f ; k t , h t , θ t ) ] − 1 σ ,

k t + 1 ( k t , h t , θ t ) = θ t k t α + ( 1 − δ ) k t − b h t − [ β ψ ( ω f ; k t , h t , θ t ) ] t − 1 σ ,

h t + 1 ( k t , h t , θ t ) = b h t + [ β ψ ( ω f ; k t , h t , θ t ) ] − 1 σ .

Step II (Monte Carlo Integration) Our objective is to find approximations of partials at a range of points. Supposing that the point of interest is ( k ¯ , h ¯ , θ ¯ ) the algorithm proceeds as follows:

Simulate N sequences of the realizations of the stochastic process { θ n } t = 0 T ¯ of size T ¯ with a starting value θ 0 n = θ ¯ , for all n = 1 , ⋯ , N . For a each sequence { θ n } t = 0 T simulate the series of the endogenous variables { k t n , h t n , c t n } t = 0 T ¯ using approximated policy functions, the laws of motion (12)-(13), and the corresponding initial values k 0 n = k ¯ , h 0 n = h ¯ . Using the simulated series calculate the discounted sums of the instantaneous utilities and average over N,

V ( k ¯ , h ¯ , θ ¯ ) ≃ 1 N ∑ n = 1 N ∑ t = 0 T ¯ β t ( c t n − b h t n ) 1 − σ 1 − σ .

Step III (Numerical Differentiation) To obtain V k ( k ¯ , h ¯ , θ ¯ ) get approximations of the optimal value function at V ( k + ϵ , h , θ ) and V ( k − ϵ , h , θ ) , where ϵ is a small positive number. Calculate the approximated value of the partial derivative using Stirling’s finite difference formula:

∂ V ( k ¯ , h ¯ , θ ¯ ) ∂ k ≃ V ( k ¯ + ϵ , h ¯ , θ ¯ ) − V ( k ¯ − ϵ , h ¯ , θ ¯ ) 2 ϵ . (14)

The partial with respect to the habit stock is obtained in a similar way. The length of the simulated series T ¯ can be very moderate due to discounting of the future utilities. The optimal value of ϵ is both computer and problem specific.

This section considers the accuracy of the algorithm in the context of our example. First, we compare performance of our algorithm with the approach in [

Consider the optimality conditions for the autarkic problem, written in the sequence form:

u c ( c t , h t ) + β λ E t [ V h ( k t + 1 , h t + 1 , θ t + 1 ) ] = β E t [ V k ( k t + 1 , h t + 1 , θ t + 1 ) ] , (15)

V k ( k t , h t , θ t ) = β E t [ V k ( k t + 1 , h t + 1 , θ t + 1 ) ] ( f k ( k t , θ t ) + 1 − δ ) , (16)

V h ( k t , h t , θ t ) = u h ( c t , h t ) + β ( 1 − λ ) E t [ V h ( k t + 1 , h t + 1 , θ t + 1 ) ] . (17)

Condition (17) can be used to compare our algorithm with the method in [

V h ( k t , h t , θ t ) = u h ( c t , h t ) + E t [ ∑ j = 1 ∞ β j ( 1 − λ ) j u h ( c t + j , h t + j ) ] .

An approximation of this derivative can be obtained by parameterizing the right hand side with flexible functional forms in the state variables ( k t , h t , θ t ) and running a non-linear regression using the simulated series from the numerical solution of the model.

The framework we have chosen for a worked out example embeds several well known special cases. For instance, if λ = 0 , it reduces to the Brock-Mirman stochastic growth model. In this case, the analytical form of the one-period return function r, which maps the graph A of the feasibility correspondence Γ into the real numbers is known. The correspondence describing the feasibility constraints is given by

Γ ( k t , θ t ) = [ f ( 1 − δ ) k t , f ( k t , θ t ) + ( 1 − δ ) k t ] ,

and the instantaneous return function becomes

r : A → ℝ given by r ( k t , k t + 1 , θ t ) = u ( f ( k t , θ t ) + ( 1 − δ ) k t − k t + 1 ) ,

where A = { ( k t , k t + 1 , θ t ) ∈ ℝ + 3 : k t + 1 ∈ Γ ( k t , θ t ) } . Hence, by virtue of the Benveniste-Scheinkman theorem the derivative of interest can be expressed as

V k ( k t , θ t ) = u ′ ( f ( k t , θ t ) + ( 1 − δ ) k t − g ( k t , θ t ) ) [ f k ( k t , θ t ) + 1 − δ ] ,

where g is the optimal policy function for capital stock. This special case allows us to compare the simulation from our algorithm with the example where the only source of approximation errors is the approximation of the policy function g. This will allow us to isolate the contribution of the approximation errors in evaluation of the integrals and numerical differentiation to the overall approximation error of the algorithm. As shown in

Our final special case compares the approximation of the derivative with its known exact solution. It is well known that for the functional forms f ( k t , θ t ) = θ t k t α , u ( c t ) = log c t , and δ = 1 , the optimal policy function is defined by the simple law of motion k t + 1 = α β θ t k t α . Moreover, the derivative of the value function has the following analytical solution:

V k ( k t , θ t ) = α θ t k t α − 1 ( 1 − α β ) θ t k t α = α 1 − α β 1 k t .

By replacing the approximated policy function with the known closed form solution, we can isolate the effect of the errors stemming from Monte Carlo integration and numerical differentiation on the accuracy of the approximation. ^{−9} of the value of the derivative. This suggests that obtaining accurate approximation of the policy functions is crucial for the accuracy of the whole algorithm.

This paper contributes to the literature on recursive contracts by proposing an algorithm for computing the gradient of the optimal value function using simulation-based techniques. The proposed procedure is conceptually straightforward, computationally inexpensive, and simple to implement. It allows researchers to extend existing risk-sharing models with limited commitment by including additional endogenous state variables. For example, one may extend a multi-country single-good model in [

In terms of accuracy, the algorithm demonstrates performance comparable with Marcet and Marimon’s method [

While our algorithm has wide applicability, it inherits its speed and accuracy trade-offs from the underlying numerical solution method. Our experiments suggest that obtaining accurate approximation of the policy functions is crucial for the accuracy of the whole algorithm.

An additional limitation on the algorithm’s computational speed is imposed by model’s simulation and Monte-Carlo integration. However, both steps can be parallelized along the lines proposed in [

The author declares no conflicts of interest regarding the publication of this paper.

Dmitriev, A. (2019) On Approximating the Gradient of the Value Function. Theoretical Economics Letters, 9, 126-138. https://doi.org/10.4236/tel.2019.91011