On Approximating the Gradient of the Value Function

Alexandre Dmitriev^{1,2}

Show more

1. Introduction

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 [1] [2] and references therein). Often the optimality conditions for this class of problems involve partial derivatives with respect to endogenous state variables of the optimal value function corresponding to the dynamic programming formulation of an outside option. Although many numerical methods can provide an approximation for the value function, there is no reason to believe that a derivative of this approximation will be close in any sense to the actual value of the derivative. In this note we suggest an algorithm for accurately computing these partial derivatives by simulation.

This issue has been previously considered in [3] , in the context of a stochastic growth model with capital accumulation under one-sided lack of commitment. To circumvent the problem of finding the values of the derivatives in [3] , the authors proposed a method based on the ideas of Benveniste and Scheinkman [4] . Unfortunately, their method has limited applicability since it depends on the availability of an analytical solution for the derivatives as conditional expectations of the known functions of the model solution. This paper proposes a simple algorithm to fill this gap in the literature.

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.

2. Applicability of the Algorithm: An Example

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 [3] cannot be applied.

Consider a model of international risk sharing, which distinguishes itself from the canonical model [5] in two respects. First, as in [6] , we introduce a friction in the credit markets. We assume that the international loans are feasible only to the extent to which they can be enforced by the threat of exclusion from participation in any other international risk sharing arrangement. Second, we incorporate habit formation preferences into the model. The motivation for doing this is threefold. Habits help us illustrate the features of the algorithm by expanding the set of endogenous states in the model. Habit formation preferences tend to improve performance of the international business cycle models [7] [8] . Finally, empirical studies suggest that habit formation is consistent with the observed consumption behavior [9] . To simplify the exposition, we assume inelastic labor supply.

The planner’s problem is to choose the sequences of consumption $\left\{{c}_{it}\right\}$ and investment $\left\{{i}_{it}\right\}$ to maximize a weighted sum of utilities

$\underset{\left\{{c}_{it},{i}_{it}\right\}}{\mathrm{max}}{E}_{0}{\displaystyle \underset{i=1}{\overset{I}{\sum}}}\text{\hspace{0.05em}}{\lambda}_{i}{\displaystyle \underset{t=0}{\overset{\infty}{\sum}}}\text{\hspace{0.05em}}{\beta}^{t}u\left({c}_{it},{h}_{it}\right)$ (1)

subject to an aggregate feasibility constraint

$\underset{i=1}{\overset{I}{\sum}}}\text{\hspace{0.05em}}{c}_{it}+{\displaystyle \underset{i=1}{\overset{I}{\sum}}}\text{\hspace{0.05em}}{i}_{it}={\displaystyle \underset{i=1}{\overset{I}{\sum}}}\text{\hspace{0.05em}}f\left({k}_{it},{\theta}_{it}\right),$ (2)

individual participation constraints for each $i=1,\cdots ,I$,

${E}_{t}{\displaystyle \underset{j=0}{\overset{\infty}{\sum}}}\text{\hspace{0.05em}}{\beta}^{j}u\left({c}_{it+j},{h}_{it+j}\right)\ge {V}_{i}^{a}\left({k}_{it},{h}_{it},{\theta}_{it}\right),$ (3)

the equations of motion for the capital,

${k}_{it+1}=\left(1-\delta \right){k}_{it}+{i}_{it},$ (4)

the laws of motion for habits,

${h}_{it+1}={h}_{it}+\lambda \left({c}_{it}-{h}_{it}\right),$ (5)

and non-negativity constraints ${c}_{it}\mathrm{,}{i}_{it}\ge 0$. We assume that productivity shocks ${\theta}_{t}$ follow a first order stationary vector autoregressive process, and that the initial values for the state variables ${k}_{i0}\mathrm{,}{h}_{i0}\mathrm{,}{\theta}_{i0}$, and the initial non-negative weights, ${\lambda}_{i}$, are given. In addition, the usual restrictions apply to the discount factor, $\beta \in \left(\mathrm{0,1}\right)$, the capital depreciation rate, $\delta \in \left(\mathrm{0,1}\right)$, and the persistence of habits, $\lambda \in \left(\mathrm{0,1}\right)$.

The outside option ${V}_{i}^{a}\left({k}_{it}\mathrm{,}{h}_{it}\mathrm{,}{\theta}_{it}\right)$ in the participation constraint (3) represents the optimal value function corresponding to the autarkic environment.

$\underset{{\left\{{c}_{it},{i}_{it}\right\}}_{t=0}^{\infty}}{\mathrm{max}}{E}_{0}{\displaystyle \underset{t=0}{\overset{\infty}{\sum}}}\text{\hspace{0.05em}}{\beta}^{t}u\left({c}_{it},{h}_{it}\right)$ (6)

subject to

${c}_{it}+{i}_{it}=f\left({k}_{it},{\theta}_{it}\right),$

${k}_{it+1}=f\left({k}_{it},{\theta}_{it}\right)-{c}_{it}+\left(1-\delta \right){k}_{it},$

${h}_{it+1}={h}_{it}+\lambda \left({c}_{it}-{h}_{it}\right),$

with the initial values being equal to the values of the state variables ${k}_{it}\mathrm{,}{h}_{it}\mathrm{,}{\theta}_{it}$ at the moment of deviation from the optimal plan.

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

$\frac{{\Lambda}_{i,t}}{{\Lambda}_{s,t}}=\frac{{\xi}_{st}}{{\xi}_{it}},$

where

$\begin{array}{c}{\Lambda}_{i,t}={u}_{c}\left(i,t\right)+\lambda \beta {E}_{t}{\displaystyle \underset{j=0}{\overset{\infty}{\sum}}}\text{\hspace{0.05em}}{\beta}^{j}{\left(1-\lambda \right)}^{j}[\frac{{\xi}_{it+1+j}}{{\xi}_{it}}{u}_{h}\left(i,t+1+j\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-\frac{{\mu}_{it+1+j}}{{\xi}_{it}}\frac{\partial {V}_{i}^{a}}{\partial {h}_{it+1+j}}\left(i,t+1+j\right)],\end{array}$ (7)

the intertemporal condition,

${\Lambda}_{i,t}=\beta {E}_{t}\left[\frac{{\xi}_{it+1+j}}{{\xi}_{it}}{\Lambda}_{i,t+1}\left({f}_{k}\left({k}_{it+1},{\theta}_{it+1}\right)+1-\delta \right)-\frac{{\mu}_{it+1}}{{\xi}_{it}}\frac{\partial {V}_{i}^{a}}{\partial {k}_{it+1}}\left(i,t+1\right)\right],$ (8)

the complementary slackness condition,

${\mu}_{it}\left[{E}_{t}{\displaystyle \underset{j=0}{\overset{\infty}{\sum}}}\text{\hspace{0.05em}}{\beta}^{j}u\left({c}_{it+j},{h}_{it+j}\right)-{V}_{i}^{a}\left({k}_{it},{h}_{it},{\theta}_{it}\right)\right]=0,$

and the law of motion for the co-state variables ${M}_{it+1}={M}_{it}+{\mu}_{it}$, where ${\xi}_{it}={\lambda}_{i}+{M}_{it+1}$, ${\mu}_{it}\ge 0$, and ${M}_{i0}=0$. In the equations above ${u}_{c}\left(i\mathrm{,}t\right)$ denotes

$\frac{\partial u\left({c}_{it}\mathrm{,}{h}_{it}\right)}{\partial {c}_{it}}$ , 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 [3] and in this paper. Because both

algorithms can approximate $\frac{\partial {V}_{i}^{a}}{\partial {k}_{it}}\left({k}_{it}\mathrm{,}{h}_{it}\mathrm{,}{\theta}_{it}\right)$, we will use that fact to compare their accuracy in Section 1. Our approach can also approximate $\frac{\partial {V}_{i}^{a}}{\partial {h}_{it}}\left({k}_{it}\mathrm{,}{h}_{it}\mathrm{,}{\theta}_{it}\right)$,

whereas [3] cannot, because the analytical solution for this partial derivative as an expectation of the known functions of the model’s solution is not available. This will be further discussed in Section 1.

3. The Algorithm

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, $\stackrel{\xaf}{x}$, and exogenous shocks, $\stackrel{\xaf}{s}$. Consider the optimal value function at a point $\left(\stackrel{\xaf}{x},\stackrel{\xaf}{s}\right)$ as an outcome of a standard optimization problem for the outside alternative:

$V\left(\stackrel{\xaf}{x},\stackrel{\xaf}{s}\right)=\underset{\left\{{a}_{t}\right\}}{\mathrm{max}}{E}_{0}{\displaystyle \underset{t=0}{\overset{\infty}{\sum}}}\text{\hspace{0.05em}}{\beta}^{t}r\left({x}_{t},{a}_{t},{s}_{t}\right)$ (9)

subject to

${x}_{t+1}=l\left({x}_{t},{a}_{t},{s}_{t}\right),\text{\hspace{0.17em}}{a}_{t}\in A\left({x}_{t},{s}_{t}\right),$ (10)

${x}_{0}=\stackrel{\xaf}{x},{s}_{0}=\stackrel{\xaf}{s},$

where r is an instantaneous utility function, $\beta \in \left(\mathrm{0,1}\right)$ the discount factor, $\left\{{s}_{t}\right\}$ 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\left({x}_{t},{s}_{t}\right)$.

The purpose of our algorithm is to find a pointwise approximation to the partial derivative $\frac{\partial}{\partial {x}_{i}}V\left(\stackrel{\xaf}{x},\stackrel{\xaf}{s}\right)$ 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}=\stackrel{^}{f}\left(\omega ;{x}_{t},{s}_{t}\right)$, which depend on the state variables and some coefficients, $\omega $.

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

$V\left(\stackrel{\xaf}{x},\stackrel{\xaf}{s}\right)\simeq \frac{1}{N}{\displaystyle \underset{n=1}{\overset{N}{\sum}}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}{\displaystyle \underset{t=0}{\overset{T}{\sum}}}\text{\hspace{0.05em}}{\beta}^{t}r\left({x}_{t}^{n}\mathrm{,}{a}_{t}^{n}\mathrm{,}{s}_{t}^{n}\right)\mathrm{.}$

Step III (Numerical Differentiation) Repeat Step II to obtain approximations of the value function at two points, for instance $V\left(\stackrel{\xaf}{x}+\u03f5{\iota}_{i}\mathrm{,}\stackrel{\xaf}{s}\right)$ and $V\left(\stackrel{\xaf}{x}-\u03f5{\iota}_{i}\mathrm{,}\stackrel{\xaf}{s}\right)$, where ${\iota}_{i}$ denotes a conformable vector of zeros with one on its i-th coordinate, and $\u03f5$ is a small positive number. Calculate the value of the partial derivative using, for example, Stirling’s finite difference formula:

$\frac{\partial}{\partial {x}_{i}}V\left(\stackrel{\xaf}{x},\stackrel{\xaf}{s}\right)\simeq \frac{V\left(\stackrel{\xaf}{x}+\u03f5{\iota}_{i}\mathrm{,}\stackrel{\xaf}{s}\right)-V\left(\stackrel{\xaf}{x}-\u03f5{\iota}_{i}\mathrm{,}\stackrel{\xaf}{s}\right)}{2\u03f5}\mathrm{.}$

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 [10] [11] [12] .

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 $\u03f5$ 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.

4. Implementation of the Algorithm

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\left(k\mathrm{,}h\mathrm{,}\theta \right)=\underset{\left(c\mathrm{,}i\right)\in A\left(k\mathrm{,}\theta \right)}{\mathrm{max}}\left\{u\left(c\mathrm{,}h\right)+\beta E\left[V\left({k}^{\prime}\mathrm{,}{h}^{\prime}\mathrm{,}{\theta}^{\prime}\right)|\left(k\mathrm{,}h\mathrm{,}\theta \right)\right]\right\}$

${h}^{\prime}=h+\lambda \left(c-h\right),$

${k}^{\prime}=\left(1-\delta \right)k+i\mathrm{,}$

$A\left(k,\theta \right)=\left\{\left(c,i\right)\in {\mathbb{R}}_{+}^{2}:c+i=f\left(k,\theta \right)\right\}.$

The objective of the algorithm is to find the values of the partials ${V}_{h}(\cdot )$ and ${V}_{k}(\cdot )$ at a point $\left(\stackrel{\xaf}{k}\mathrm{,}\stackrel{\xaf}{h}\mathrm{,}\stackrel{\xaf}{\theta}\right)$, 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\left(\stackrel{\xaf}{k}+\epsilon \mathrm{,}\stackrel{\xaf}{h}\mathrm{,}\stackrel{\xaf}{\theta}\right)$ and $V\left(\stackrel{\xaf}{k}-\epsilon \mathrm{,}\stackrel{\xaf}{h}\mathrm{,}\stackrel{\xaf}{\theta}\right)$ 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:

${\Lambda}_{t}=\beta {E}_{t}\left[{\Lambda}_{t+1}\left({f}_{k}\left({k}_{t+1},{\theta}_{t+1}\right)+1-\delta \right)\right],$

where marginal utility of consumption is

${\Lambda}_{t}={u}_{c}\left({c}_{t},{h}_{t}\right)+\beta \lambda {E}_{t}\left[{\displaystyle \underset{j=0}{\overset{\infty}{\sum}}}\text{\hspace{0.05em}}{\beta}^{j}{\left(1-\lambda \right)}^{j}{u}_{h}\left({c}_{t+j+1},{h}_{t+j+1}\right)\right].$

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

$u\left({c}_{t},{h}_{t}\right)=\frac{{\left({c}_{t}-b{h}_{t}\right)}^{1-\sigma}}{1-\sigma}$, where $b\in \left(\mathrm{0,1}\right)$ and $\sigma >0$. The production function

is Cobb-Douglas and is given by $f\left({k}_{t},{\theta}_{t}\right)={\theta}_{t}{k}_{t}^{\alpha}$. The stochastic process for productivity is $\mathrm{log}{\theta}_{t}=\rho \mathrm{log}{\theta}_{t-1}+{\epsilon}_{t}$, where $\rho \in \left(\mathrm{0,1}\right)$ , and $\left\{{\epsilon}_{t}\right\}$ are independent normally distributed random variables with zero mean and variance ${\sigma}_{\epsilon}^{2}$. In this example, we restrict attention to one particular set of the parameters which are summarized in Table 1.

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 ${\left\{{c}_{t},{h}_{t+1},{k}_{t+1}\right\}}_{t=0}^{\infty}$ must satisfy the following system of stochastic difference equations:

$\begin{array}{l}{\left({c}_{t}-b{h}_{t}\right)}^{-\sigma}\\ =\beta {E}_{t}\left[b{\left({c}_{t+1}-b{h}_{t+1}\right)}^{-\sigma}\times \left(1+\left(\alpha {\theta}_{t+1}{k}_{t+1}^{\alpha -1}+1-\delta \right)\left(\frac{1}{b}-\beta {\left(\frac{{c}_{t+2}-b{h}_{t+2}}{{c}_{t+1}-b{h}_{t+1}}\right)}^{-\sigma}\right)\right)\right],\end{array}$ (11)

${k}_{t+1}={\theta}_{t}{k}_{t}^{\alpha}+\left(1-\delta \right){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. [13] ). It takes the following steps:

1) Fix the initial conditions and draw a series of ${\left\{{\theta}_{t}\right\}}_{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}\mathrm{,}{h}_{t}\mathrm{,}{\theta}_{t}$ and some coefficients, $\omega $, to yield

Table 1. Parameterization of the model.

${\left({c}_{t}\left(\omega \right)-b{h}_{t}\left(\omega \right)\right)}^{-\sigma}=\beta \psi \left(\omega \mathrm{;}{k}_{t}\left(\omega \right)\mathrm{,}{h}_{t}\left(\omega \right)\mathrm{,}{\theta}_{t}\right)\mathrm{,}$

where $\psi \left(\omega \mathrm{;}{k}_{t}\left(\omega \right)\mathrm{,}{h}_{t}\left(\omega \right)\mathrm{,}{\theta}_{t}\right)=\mathrm{exp}\left({P}_{n}\left(\omega \mathrm{;}\mathrm{log}{k}_{t}\left(\omega \right)\mathrm{,}\mathrm{log}{h}_{t}\left(\omega \right)\mathrm{,}\mathrm{log}{\theta}_{t}\right)\right)$, 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}\left(\omega \right)$, 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 ${\left\{{\theta}_{t}\right\}}_{t=0}^{T}$, repeat the previous step in order to obtain recursively a series of the endogenous variables ${\left\{{c}_{t}\left(\omega \right),{k}_{t+1}\left(\omega \right),{h}_{t+1}\left(\omega \right)\right\}}_{t=0}^{T}$, for this particular parameterization $\omega $.

4) Run the following non-linear regression

${Y}_{t}\left(\omega \right)=\mathrm{exp}\left({P}_{n}\left(\xi ;\mathrm{log}{k}_{t}\left(\omega \right),\mathrm{log}{h}_{t}\left(\omega \right),\mathrm{log}{\theta}_{t}\right)\right)+{\eta}_{t},$

where the role of the dependent variable ${Y}_{t}\left(\omega \right)$ is performed by the expression inside the conditional expectation in Equation (11).

5) Letting $S\left(\omega \right)$ 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 ${\omega}_{f}=S\left({\omega}_{f}\right)$. This would provide the solution for the endogenous variables ${\left\{{c}_{t}\left({\omega}_{f}\right),{k}_{t+1}\left({\omega}_{f}\right),{h}_{t+1}\left({\omega}_{f}\right)\right\}}_{t=0}^{T}$ for this particular realization of the stochastic process ${\left\{{\theta}_{t}\right\}}_{t=1}^{T}$ along with the approximated policy functions:

${c}_{t}\left({k}_{t},{h}_{t},{\theta}_{t}\right)=b{h}_{t}+{\left[\beta \psi \left({\omega}_{f};{k}_{t},{h}_{t},{\theta}_{t}\right)\right]}^{-\frac{1}{\sigma}},$

${k}_{t+1}\left({k}_{t},{h}_{t},{\theta}_{t}\right)={\theta}_{t}{k}_{t}^{\alpha}+\left(1-\delta \right){k}_{t}-b{h}_{t}-{\left[\beta \psi \left({\omega}_{f};{k}_{t},{h}_{t},{\theta}_{t}\right)\right]}_{t}^{-\frac{1}{\sigma}},$

${h}_{t+1}\left({k}_{t},{h}_{t},{\theta}_{t}\right)=b{h}_{t}+{\left[\beta \psi \left({\omega}_{f};{k}_{t},{h}_{t},{\theta}_{t}\right)\right]}^{-\frac{1}{\sigma}}.$

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 $\left(\stackrel{\xaf}{k}\mathrm{,}\stackrel{\xaf}{h}\mathrm{,}\stackrel{\xaf}{\theta}\right)$ the algorithm proceeds as follows:

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

$V\left(\stackrel{\xaf}{k}\mathrm{,}\stackrel{\xaf}{h}\mathrm{,}\stackrel{\xaf}{\theta}\right)\simeq \frac{1}{N}{\displaystyle \underset{n=1}{\overset{N}{\sum}}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}{\displaystyle \underset{t=0}{\overset{\stackrel{\xaf}{T}}{\sum}}}\text{\hspace{0.05em}}{\beta}^{t}\frac{{\left({c}_{t}^{n}-b{h}_{t}^{n}\right)}^{1-\sigma}}{1-\sigma}\mathrm{.}$

Step III (Numerical Differentiation) To obtain ${V}_{k}\left(\stackrel{\xaf}{k}\mathrm{,}\stackrel{\xaf}{h}\mathrm{,}\stackrel{\xaf}{\theta}\right)$ get approximations of the optimal value function at $V\left(k+\u03f5\mathrm{,}h\mathrm{,}\theta \right)$ and $V\left(k-\u03f5\mathrm{,}h\mathrm{,}\theta \right)$, where $\u03f5$ is a small positive number. Calculate the approximated value of the partial derivative using Stirling’s finite difference formula:

$\frac{\partial V\left(\stackrel{\xaf}{k}\mathrm{,}\stackrel{\xaf}{h}\mathrm{,}\stackrel{\xaf}{\theta}\right)}{\partial k}\simeq \frac{V\left(\stackrel{\xaf}{k}+\u03f5\mathrm{,}\stackrel{\xaf}{h}\mathrm{,}\stackrel{\xaf}{\theta}\right)-V\left(\stackrel{\xaf}{k}-\u03f5\mathrm{,}\stackrel{\xaf}{h}\mathrm{,}\stackrel{\xaf}{\theta}\right)}{2\u03f5}\mathrm{.}$ (14)

The partial with respect to the habit stock is obtained in a similar way. The length of the simulated series $\stackrel{\xaf}{T}$ can be very moderate due to discounting of the future utilities. The optimal value of $\u03f5$ is both computer and problem specific.

5. Numerical Accuracy: A Comparison

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 [3] when such comparison is feasible. Next, we present several special cases, which isolate the contributions of different sources to the overall approximation error.

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

${u}_{c}\left({c}_{t}\mathrm{,}{h}_{t}\right)+\beta \lambda {E}_{t}\left[{V}_{h}\left({k}_{t+1}\mathrm{,}{h}_{t+1}\mathrm{,}{\theta}_{t+1}\right)\right]=\beta {E}_{t}\left[{V}_{k}\left({k}_{t+1}\mathrm{,}{h}_{t+1}\mathrm{,}{\theta}_{t+1}\right)\right]\mathrm{,}$ (15)

${V}_{k}\left({k}_{t}\mathrm{,}{h}_{t}\mathrm{,}{\theta}_{t}\right)=\beta {E}_{t}\left[{V}_{k}\left({k}_{t+1}\mathrm{,}{h}_{t+1}\mathrm{,}{\theta}_{t+1}\right)\right]\left({f}_{k}\left({k}_{t}\mathrm{,}{\theta}_{t}\right)+1-\delta \right)\mathrm{,}$ (16)

${V}_{h}\left({k}_{t}\mathrm{,}{h}_{t}\mathrm{,}{\theta}_{t}\right)={u}_{h}\left({c}_{t}\mathrm{,}{h}_{t}\right)+\beta \left(1-\lambda \right){E}_{t}\left[{V}_{h}\left({k}_{t+1}\mathrm{,}{h}_{t+1}\mathrm{,}{\theta}_{t+1}\right)\right]\mathrm{.}$ (17)

Condition (17) can be used to compare our algorithm with the method in [3] . The latter requires solving the model numerically and expressing the derivatives of interest in terms of conditional expectations and functions of the equilibrium path of the model. Applying recursive substitution and the law of iterated expectations to (17) yields:

${V}_{h}\left({k}_{t},{h}_{t},{\theta}_{t}\right)={u}_{h}\left({c}_{t},{h}_{t}\right)+{E}_{t}\left[{\displaystyle \underset{j=1}{\overset{\infty}{\sum}}}\text{\hspace{0.05em}}{\beta}^{j}{\left(1-\lambda \right)}^{j}{u}_{h}\left({c}_{t+j},{h}_{t+j}\right)\right].$

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

Figure 1 shows the approximations of the derivative obtained using our algorithm and the approach in [3] . The approximated values of ${V}_{h}\left({k}_{t}\mathrm{,}{h}_{t}\mathrm{,}{\theta}_{t}\right)$ are plotted for a range of a state variable while keeping the remaining states fixed at their deterministic steady state values. The histograms plot the sample distributions of capital and habit stock. A few observations can be made based on Figure 1. First, the two algorithms produce indistinguishable results when the state variables take the values which often occur in equilibrium. Second, for the points which are unlikely to occur in equilibrium, the approximations differ significantly. To see this feature, consider the range of values of capital stock in excess of 6.5. The plots of the approximate derivatives reported in the upper panel of Figure 1 do not coincide. Moreover, the upper tail of the histogram suggests that such values of ${k}_{t}$ are not unlikely to happen in equilibrium. Notice, that while considering a relatively high value of ${k}_{t}$ we kept the remaining arguments of ${V}_{h}\left({k}_{t}\mathrm{,}{h}_{t}\mathrm{,}{\theta}_{t}\right)$ at their deterministic steady state values. However, the points where capital is very high while consumption (and hence the habit stock) are at the steady state level are rather unusual. This is an expected result, since Monte Carlo integration delivers good approximations only in the region of the state space which is frequently visited by the model in equilibrium.

The framework we have chosen for a worked out example embeds several well known special cases. For instance, if $\lambda =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 $\Gamma $ into the real numbers is known. The correspondence describing the feasibility constraints is given by

$\Gamma \left({k}_{t},{\theta}_{t}\right)=\left[f\left(1-\delta \right){k}_{t},f\left({k}_{t},{\theta}_{t}\right)+\left(1-\delta \right){k}_{t}\right],$

and the instantaneous return function becomes

$r\mathrm{:}A\to \mathbb{R}\text{\hspace{0.17em}}\text{given}\text{\hspace{0.17em}}\text{by}\text{\hspace{0.17em}}r\left({k}_{t}\mathrm{,}{k}_{t+1}\mathrm{,}{\theta}_{t}\right)=u\left(f\left({k}_{t}\mathrm{,}{\theta}_{t}\right)+\left(1-\delta \right){k}_{t}-{k}_{t+1}\right)\mathrm{,}$

where $A=\left\{\left({k}_{t}\mathrm{,}{k}_{t+1}\mathrm{,}{\theta}_{t}\right)\in {\mathbb{R}}_{+}^{3}\mathrm{:}{k}_{t+1}\in \Gamma \left({k}_{t}\mathrm{,}{\theta}_{t}\right)\right\}$. Hence, by virtue of the Benveniste-Scheinkman theorem the derivative of interest can be expressed as

${V}_{k}\left({k}_{t}\mathrm{,}{\theta}_{t}\right)={u}^{\prime}\left(f\left({k}_{t}\mathrm{,}{\theta}_{t}\right)+\left(1-\delta \right){k}_{t}-g\left({k}_{t}\mathrm{,}{\theta}_{t}\right)\right)\left[{f}_{k}\left({k}_{t}\mathrm{,}{\theta}_{t}\right)+1-\delta \right]\mathrm{,}$

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 Figure 1, the approximation delivered by our algorithm is very close to the approximation which relies on the Benveniste-Scheinkman theorem. Once again, in the region of the state space which is often visited by the model in equilibrium the two approximations are virtually identical. This allows us to tentatively suggest that the main contribution to the approximation error of the algorithm comes from the approximation of the policy functions.

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\left({k}_{t},{\theta}_{t}\right)={\theta}_{t}{k}_{t}^{\alpha}$, $u\left({c}_{t}\right)=\mathrm{log}{c}_{t}$, and $\delta =1$, the optimal policy function is defined by the simple law of motion ${k}_{t+1}=\alpha \beta {\theta}_{t}{k}_{t}^{\alpha}$. Moreover, the derivative of the value function has the following analytical solution:

${V}_{k}\left({k}_{t},{\theta}_{t}\right)=\frac{\alpha {\theta}_{t}{k}_{t}^{\alpha -1}}{\left(1-\alpha \beta \right){\theta}_{t}{k}_{t}^{\alpha}}=\frac{\alpha}{1-\alpha \beta}\frac{1}{{k}_{t}}.$

Figure 1. Comparison with the algorithm of Marcet and Marimon [3] .

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. Figure 2 compares the approximated derivative obtained using the exact policy function for
${k}_{t}$ with its analytical counterpart. The reported graphs are visually indistinguishable for the range of six standard deviations of
${k}_{t}$ around its deterministic steady state value. The approximation errors stemming from Monte Carlo integration and numerical differentiation are of an order of 10^{−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.

6. Concluding Remarks

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 [6] by including different types of capital or time-non-separable preferences. Alternatively, a multi-country multi-good model in [14] can be extended to include capital accumulation and cross-country trade in investment goods.

Figure 2. Comparison with a Closed Form Solution.

In terms of accuracy, the algorithm demonstrates performance comparable with Marcet and Marimon’s method [3] in our benchmark example. In contrast to [3] , our method is flexible enough to handle dynamic models with large numbers of state variables even when derivatives of interest cannot be expressed in terms of conditional expectations and functions of the equilibrium path of the model.

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 [15] in order to reduce the computational time burden. Exploring the costs and benefits of a parallel implementation of the algorithm is left for the future research.

NOTES

^{1}Derivation of the optimality condition is available upon request.

References

[1] Cole, H. and Kubler, F. (2012) Recursive Contracts, Lotteries and Weakly Concave Pareto Sets. Review of Economic Dynamics, 15, 479-500.

[2] Golosov, M., Tsyvinski, A. and Werquin, N. (2016) Recursive Contracts and Endogenously Incomplete Markets. In: Handbook of Macroeconomics, Vol. 2, Elsevier, Amsterdam, 725-841.

https://doi.org/10.1016/bs.hesmac.2016.03.007

[3] Marcet, A. and Marimon, R. (1992) Communication, Commitment, and Growth. Journal of Economic Theory, 58, 219-249.

https://doi.org/10.1016/0022-0531(92)90054-L

[4] Benveniste, L.M. and Scheinkman, J.A. (1979) On the Differentiability of the Value Function in Dynamic Models of Economics. Econometrica, 47, 727-732.

https://doi.org/10.2307/1910417

[5] Backus, D.K., Kehoe, P.J. and Kydland, F.E. (1992) International Real Business Cycles. Journal of Political Economy, 100, 745-775.

https://doi.org/10.1086/261838

[6] Kehoe, P.J. and Perri, F. (2002) International Business Cycles with Endogenous Incomplete Markets. Econometrica, 70, 907-928.

https://doi.org/10.1111/1468-0262.00314

[7] Dmitriev, A. and Roberts, I. (2012) International Business Cycles with Complete Markets. Journal of Economic Dynamics and Control, 36, 862-875.

https://doi.org/10.1016/j.jedc.2011.12.006

[8] Dmitriev, A. and Krznar, I. (2012) Habit Persistence and International Comovements. Macroeconomic Dynamics, 16, 312-330.

https://doi.org/10.1017/S1365100510000957

[9] Chen, X. and Ludvigson, S.C. (2009) Land of Addicts? An Empirical Investigation of Habit-Based Asset Pricing Models. Journal of Applied Econometrics, 24, 1057-1093.

[10] Judd, K.L. (1998) Numerical Methods in Economics. The MIT Press, Cambridge, MA.

[11] Mathews, J.H. and Fink, K.K. (2004) Numerical Methods Using Matlab. 4th Edition, Prentice-Hall, Upper Saddle River, NJ.

[12] Press, W.H., Teukolsky, S.A., Vetterling, W.T. and Flannery, B.P. (1992) Numerical Recipes in C: The Art of Scientific Computing. 2nd Edition, Cambridge University Press, Cambridge.

[13] den Haan, W. and Marcet, A. (1990) Solving the Stochastic Growth Model by Parameterizing Expectations. Journal of Business & Economic Statistics, 8, 31-34.

[14] Bodenstein, M. (2008) International Asset Markets and Real Exchange Rate Volatility. Review of Economic Dynamics, 11, 688-705.

https://doi.org/10.1016/j.red.2007.12.003

[15] Creel, M. (2008) Using Parallelization to Solve a Macroeconomic Model: A Parallel Parameterized Expectations Algorithm. Computational Economics, 32, 343-352.

https://doi.org/10.1007/s10614-008-9142-6