Technical Appendix: The Optimum Quantity of Debt^{†}

S. Rao Aiyagari and Ellen R. McGrattan Federal Reserve Bank of Minneapolis November 1995

† The views expressed herein are those of the authors and not necessarily those of the Federal Reserve Bank of Minneapolis or the Federal Reserve System.

Technical Appendix: The Optimum Quantity of Debt

In this appendix, we describe the numerical methods used to compute an equilibrium in the economy with an inelastic labor supply and in the economy with an elastic labor supply (i.e., our benchmark economy). Although the economy with inelastically supplied labor is a special case of the benchmark economy, the equilibrium in the inelastic labor supply case is much easier to compute and is therefore treated separately. In each case, we start with the consumer’s problem, assuming the consumer takes prices as given. We then show how the equilibrium prices are determined. To verify that the methods work well for our problem, we apply them to some related test problems that have known solutions.

1. The case with inelastically supplied labor

1.1. Computing the consumer’s decision functions

The consumer chooses sequences of consumption and asset holdings to maximize expected utility; i.e.,^{1}

max_{{c̃_{t},ã_{t+1}}}E[∑_{t=0}^{∞}(β(1 + g)^{1−ν})^{t}c̃_{t}^{1−ν}/(1 − ν)|ã_{0}, e_{0}]  (1)
subject to c̃_{t} + (1 + g)ã_{t+1} ≤ (1 + r̄)ã_{t} + w̄e_{t} + χ,  (2)
ã_{t} ≥ 0,  (3)

with the after-tax interest rate r̄, the after-tax wage rate w̄, and a lump-sum transfer χ taken as given. The productivity level e_{t} is assumed to follow a Markov chain with Ee_{t} = 1. The probability of transiting from state i to state j is given by π_{i,j}, i, j = 1, …, m.

To incorporate the inequality constraints in Eq. (3), we replace the objective in Eq. (1) with

E[∑_{t=0}^{∞}(β(1 + g)^{1−ν})^{t}\{\frac{c̃_{t}^{1−ν}}{1 − ν} + \frac{1}{3}ζ min(ã_{t}, 0)^{3}\}|ã_{0}, e_{0}],  (4)

1 As in the paper, we normalize variables by dividing through by output, and we assume that ã_{t} ≥ 0

is sufficiently restrictive for the parameters that we consider.

where ζ is some positive parameter. Note that if ã_{t} is negative, then we subtract −ζã_{t}^{3} from the consumer’s value function. With this respecification, the optimization problem is as follows: given ζ, we choose sequences of consumption and assets that are optimal for Eq. (4) subject to Eq. (2). The parameter ζ is set so that Eq. (3) is approximately satisfied. In practice, the optimization is done either by iterating over a sequence 1, 10, 10^{2}, etc. for ζ, until the constraints are satisfied to within the tolerance, or by starting with a reasonably large value of ζ.^{2}

Because the decision rules are stationary, we can compute the functions c(x, i) and α(x, i) that satisfy the following first-order conditions for i = 1, …, m and x ∈ [0, x_{max}]:

(1 + g)c(x, i)^{−ν} = β(1 + g)^{1−ν}{∑_{j} π_{i,j} (1 + r̄)c(x′, j)^{−ν} + ζ min(α(x, i), 0)^{2}}, (5)
c(x, i) + (1 + g)α(x, i) = (1 + r̄)x + w̄e(i) + χ,  (6)

where x′ = α(x, i) and e(i) is the productivity level in state i. Note that these conditions assume an interior solution for c(x, i). If we substitute the expression for c(x, i) in (6) into Eq. (5), we have a functional equation in α. We denote this expression by R(x, i; α); i.e.,

R(x, i; α) = (1 + g){(1 + r̄)x + w̄e(i) + χ − (1 + g)α(x, i)}^{−ν}

− β(1 + g)^{1−ν}{∑_{j} π_{i,j} {(1 + r̄)(1 + r̄)α(x, i) + w̄e(j)

+ χ − (1 + g)α(α(x, i), j)}^{−ν} + ζ min(α(x, i), 0)^{2}}.  (7)

The computational task is, therefore, to find an approximation for α(x, i) – say α^{h}(x, i)

– that implies R(x, i; α^{h}) ≈ 0 for all x and i. We accomplish this task by applying a finite element method. In particular, we do the following. We choose some discretization of the domain of our functions. Since only x is continuous, we need to specify some partition of [0,x_{max}] where x_{max} is such that no x > x_{max} would be chosen by the consumer. We refer to each subinterval of x as an element. On each element, we choose a set of basis functions for approximating α; that is, we assume α^{h} can be represented as a weighted sum of basis functions, where the weights and basis functions may be different for each element. In our case, we choose linear basis functions for all elements; e.g.

2 See R. Fletcher (1987), Practical Methods of Optimization (New York: Wiley) for a discussion of the

relationship of ζ to the Lagrange multipliers of the constraints in Eq. (3).

α^{h}(x, i) = ψ_{e}^{i}N_{e}(x) + ψ_{e+1}^{i}N_{e+1}(x), N_{e}(x) = \frac{x_{e+1} − x}{x_{e+1} − x_{e}} , N_{e+1}(x) = \frac{x − x_{e}}{x_{e+1} − x_{e}} ,

on the element [x_{e}, x_{e+1}]. This choice is motivated by our test problem (see below). Notice that α^{h}(x_{e}, i) = ψ_{e}^{i} and α^{h}(x_{e+1}, i) = ψ_{e+1}^{i}. If we consider the approximation globally, we need to compute values for ψ_{e}^{i} for all nodes e (assume there are n) and for all levels of productivity i. We choose these values for ψ_{e}^{i} by setting the weighted residual to zero; i.e.,

∫ R(x, i; α^{h})N_{e}(x)dx = 0,  i = 1, …, m, e = 1, … n.  (8)

In effect, we solve a problem of the following form: find ψ⃗ such that h(ψ⃗) = 0, where ψ⃗ is the vector of coefficients that we are searching over and h is the system of equations in Eq. (8).

There are several practical points worth noting. The first point is that, for parameterizations in which the no-borrowing constraint binds, the penalty function only serves to get α(x, i) > −ϵ, where ϵ is small but positive. Thus, if the function is truly equal to zero for low values of x and i, then the algorithm will not yield a mass point at x = 0. To deal with this problem, we use a two-step procedure. At the first step, we apply the penalty function method and choose a sufficiently fine mesh to resolve the kink. The kink is defined as the grid point at which the second derivative is maximized. At the second step, we use the candidate solution to determine the boundary conditions α^{h}(x, i) = 0, x < x^{∗}, where x^{∗} is the grid point at which the second derivative of α^{h}(x, i) is highest. These boundary constraints are imposed on the solution prior to the final run. The second point concerns the procedure for solving h(ψ⃗) = 0. We solve this system of equations by applying a Newton method. Therefore, we need to compute the derivative of h(ψ⃗) with respect to ψ⃗, and we need to invert it. We calculate analytical derivatives because it saves computing time. With respect to inverting the Jacobian, we can rely on the fact that it is very sparse. The sparseness is due to the fact that the approximation is done element by element; thus, the basis functions are nonzero on relatively small subdomains. In practice, however, we only use sparse matrix routines if the dimension of the matrix to be inverted is bigger than 2000 × 2000.

1.2. Computing the distribution of assets

To compute equilibrium prices, we must first construct the joint distribution over assets and productivities. This distribution can be derived from the decision rule for asset holdings; i.e.,

x_{t+1} = α(x_{t}, e_{t}),

where x_{t} is asset holdings in t and e_{t} is the earnings shock in t. We want to compute the invariant distribution for asset holdings; namely, H(x, i) = P r(x_{t} < x | e_{t} = e(i)). To compute the distribution, we solve the following functional equation:

H(x, i) = ∑_{j=1}^{m} π_{j,i}H(α^{−1}(x, j), j)I(x ≥ α(0, j)),  (9)

where π is the transition matrix for the Markov chain governing earnings and I is an indicator function (i.e., I(x > y) is equal to one if x > y and is equal to zero otherwise). The form of the functional equation in (9) is motivated below.

Suppose X_{1} is a random variable with density function f_{1}. Assume that X_{2} is a second random variable that is given by X_{2} = h(X_{1}) for some function h. In this case, the density function for X_{2} is

f_{2}(x_{2}) = ∫_{−∞}^{∞} f_{1}(x_{1})δ(x_{2} − h(x_{1}))dx_{1},  (10)

where δ is the Dirac delta function. If we integrate the left-hand side of Eq. (10), we get

F_{2}(x_{2}) = ∫_{−∞}^{x_{2}} f_{2}(s)ds = ∫_{−∞}^{x_{2}} ∫_{−∞}^{∞} f_{1}(x_{1})δ(s − h(x_{1}))dx_{1}ds
= ∫_{−∞}^{∞}[∫_{−∞}^{x_{2}} δ(s − h(x_{1}))ds]f_{1}(x_{1})dx_{1}ds
= ∫_{−∞}^{∞} I(h(x_{1}) ≤ x_{2})f_{1}(x_{1})dx_{1}
= ∫_{−∞}^{∞} I(h(x_{1}) ≤ x_{2})dF_{1}(x_{1}).  (11)

Now assume that h is strictly increasing and, hence, invertible. In this case, we can rewrite (11) as follows:

F_{2}(x_{2}) = ∫_{−∞}^{∞} I(h^{−1}(h(x_{1})) ≤ h^{−1}(x_{2}))dF_{1}(x_{1})
= ∫_{−∞}^{∞} I(x_{1} ≤ h^{−1}(x_{2}))dF_{1}(x_{1})
= ∫_{−∞}^{h^{−1}(x_{2})} dF_{1}(x_{1})
= F_{1}(h^{−1}(x_{2})).

In the example with assets, we have α increasing but not strictly. To motivate the functional equation in Eq. (9), we consider cases in which the borrowing constraints do not bind and cases in which they do and conclude that both cases imply the same functional equation for H. First, consider states of the world where the borrowing constraint does not bind (e.g., values of i such that α(x, i) > 0 for all x ≥ 0). Using (11), we know that the equation to be solved is

H_{n+1}(x′, i) = ∑_{j=1}^{m} π_{j,i} ∫_{−∞}^{∞} I(α(x, j) ≤ x′)dH_{n}(x, j),  (12)

where H_{n} is the cumulative distribution function for today’s assets, H_{n+1} is the cumulative distribution function for tomorrow’s assets, and both H_{n} and H_{n+1} are functions of the earnings state. If an invariant distribution exists, it is given by the fixed point of Eq. (12); i.e., H = H_{n} = H_{n+1}. In this case, I(α(x, j) ≤ x′) is equal to one for x′ ≥ α(0, j) and is equal to zero for x′ < α(0, j). Therefore,

∫_{−∞}^{∞} I(α(x, j) ≤ x′)dH_{n}(x, j) = \{ \begin{matrix}H_{n}(α^{−1}(x′, j), j) & if x′ ≥ α(0, j) \\ 0 & if x′ < α(0, j), \end{matrix}

and (12) can be rewritten as follows:

H_{n+1}(x′, i) = ∑_{j=1}^{m} π_{j,i}H_{n}(α^{−1}(x′, j), j)I(x′ ≥ α(0, j)).  (13)

Second, consider states of the world where the borrowing constraint binds (e.g., values of i such that α(x, i) = 0 for x ∈ [0, x^{∗}] where x^{∗} > 0). In these cases, if we assume that α^{−1}(0, j) = x^{∗} where x^{∗} is the point at which we see a kink, then we again have the mapping in (13). If we assume that H_{n} = H_{n+1} at a fixed point, the functional equation for the invariant distribution is given by Eq. (9).

To compute H, we again apply the finite element method with linear basis functions. In this case, the residual is the difference between the right- and left-hand sides of Eq. (9). For this problem, we do not have to worry about inequality constraints directly, but we do have to deal with them indirectly. If inequality constraints bind in the consumer’s problem, then the decision functions for low productivity levels will be set to zero for some interval [0, x^{∗}]. Thus, there will be a mass point at x = 0 and throughout the distribution at points traversed prior to reaching the zero-asset position. The mass points in the distribution imply that the solution to (9) has discontinuities, possibly at a countably infinite number of points. Thus, we need to know if the finite element method, as we implement it, will yield a good approximation to the distribution H.

1.3. Computing an equilibrium

Assume that we have in place numerical algorithms for computing α(x, i) and H(x, i). We now describe an algorithm for computing the equilibrium after-tax interest rate. All other prices and quantities can be determined from the interest rate.

We start with bounds on the interest rate r; that is, we assume that the equilibrium interest rate is in the interval [r_{l}, r_{u}]. Our inputs are the share of capital in production θ, the government consumption-to-GDP ratio γ, the debt-to-GDP ratio b, the discount factor β, the rate of capital depreciation δ, the utility parameters η and μ, the growth rate g, the initial bounds on r, the productivity levels e(i), i = 1, …, m, and the transition probability matrix π. Our initial guess for the equilibrium r is \frac{1}{2}(r_{l} + r_{u}). If lump-sum taxes are assumed, then we set the taxes equal to γ + χ + (r − g)b. If proportional taxes are assumed, then we back out the tax rate from the government budget constraint; i.e.,

τ_{y} = \frac{γ + χ + rb − gb}{1 + rb − δθ/(r + δ)}.  (14)

With τ_{y} and r, we have the after-tax interest rate r̄ = (1 − τ_{y})r. We can also compute the after-tax wage rate

w̄ = (1 − τ_{y})(1 − θ).  (15)

We now have the inputs needed for computing finite element approximations of α and H. The next step is to calculate the mean asset holdings Eã_{t}. Since our approximation for H is piecewise linear, an estimate of the mean asset holdings is simply

\frac{1}{2} ∑_{j} ∑_{i}(Ĥ(i + 1, j) − Ĥ(i, j))(x(i + 1) + x(i)),

where x(i) is the ith grid point over asset holdings and Ĥ(i, j) is the finite element approximation of H(x(i), j). Given this estimate of Eã_{t}, we check to see if markets are cleared; i.e., if

Eã_{t} = \frac{θ}{r + δ} + b.  (16)

If the right-hand side of (16) is larger than the left-hand side, then we set r_{l} = r and repeat the above steps. Otherwise, we set r_{u} = r and repeat. Iterations are made until Eq. (16) is satisfied to within some tolerance level.

2. The case with elastically supplied labor

2.1. Computing the consumer’s decision functions

The consumer chooses sequences of consumption, asset holdings, and leisure to maximize expected utility; i.e.,

max_{{c̃_{t},ã_{t+1},ℓ_{t}}}E[∑_{t=0}^{∞}(β(1 + g)^{η(1−μ)})^{t}(c̃_{t}^{η}ℓ_{t}^{1−η})^{1−μ}/(1 − μ)|ã_{0}, e_{0}]  (17)
subject to c̃_{t} + (1 + g)ã_{t+1} ≤ (1 + r̄)ã_{t} + w̄e_{t}(1 − ℓ_{t}) + χ,  (18)
ã_{t} ≥ 0,  (19)
ℓ_{t} ≤ 1,  (20)

with the after-tax interest rate r̄, the after-tax wage rate w̄, and a lump-sum transfer χ taken as given. In specifying the optimization problem, we have left out two inequality constraints; namely, c_{t} ≥ 0 and ℓ_{t} ≥ 0. When we compute the decision functions, these inequality constraints are ignored, but we check the solution to make sure that they are satisfied. The productivity level e_{t} is assumed to follow a Markov chain with Ee_{t} = 1. The probability of transiting from state i to state j is given by π_{i,j}, i, j = 1, …, m.

To incorporate the inequality constraints in equations (19) and (20), we replace (17) with

E[∑_{t=0}^{∞}(β(1 + g)^{η(1−μ)})^{t}\{\frac{(c̃_{t}^{η}ℓ_{t}^{1−η})^{1−μ}}{1 − μ} + \frac{1}{3}ζ(min(ã_{t}, 0)^{3} + min(1 − ℓ_{t}, 0)^{3})\}|ã_{0}, e_{0}], (21)

where ζ is some positive parameter. Note that if ã_{t} is less than zero, then we subtract −ζã_{t}^{3} from the consumer’s value function. If ℓ_{t} is greater than one, then we subtract −ζ(1 − ℓ_{t})^{3}. With this respecification, the optimization problem is as follows: given ζ, we choose sequences of consumption, asset holdings, and leisure that are optimal for Eq. (21) subject to Eq. (18). The parameter ζ is set so that Eq. (19) is approximately satisfied. In practice, the optimization is done either by iterating over a sequence 1, 10, 10^{2}, etc. for ζ, until the constraints are satisfied to within the tolerance, or by starting with a reasonably large value of ζ.

Because the decision rules are stationary, we can compute the functions c(x, i), α(x, i), and ℓ(x, i) that satisfy the following first-order conditions:

η(1 + g)c(x, i)^{η(1−μ)−1}ℓ(x, i)^{(1−η)(1−μ)} = β(1 + g)^{η(1−μ)}{∑_{j} π_{i,j} η(1 + r̄)

· c(x′, j)^{η(1−μ)−1}ℓ(x′, j)^{(1−η)(1−μ)} + ζ min(α(x, i), 0)^{2}},  (22)

(1 − η)c(x, i)^{η(1−μ)}ℓ(x, i)^{(1−η)(1−μ)−1}

= ηw̄e(i)c(x, i)^{η(1−μ)−1}ℓ(x, i)^{(1−η)(1−μ)}
+ ζ min(1 − ℓ(x, i), 0)^{2},  (23)
c(x, i) + (1 + g)α(x, i) = (1 + r̄)x + w̄e(i)(1 − ℓ(x, i)) + χ,  (24)
x′ = α(x, i),  (25)

for i = 1, …, m and x ∈ [0, x_{max}], where e(i) is the productivity level in state i. Notice that c(x, i) is a function of x, i, e(i), the parameters, and the function α(·). Therefore, with a candidate solution for α, we can back out c using (24). This is not the case for leisure, however, since (23) is an implicit function of ℓ(x, i) (if we assume that x, i, e(i), the parameters, α(·), and c(·) are known). However, we can construct a numerical solution for ℓ(x, i) by applying a Newton method with given values for η, μ, ζ, r̄, w̄, χ, g, x, i, e(i), and α(x, i). We simply iterate on

ℓ^{k+1} = ℓ^{k} − f(ℓ^{k})/J(ℓ^{k}), k = 0, 1, …, (26)

where

f(ℓ) = (1 − η)c(ℓ)^{η(1−μ)}ℓ^{(1−η)(1−μ)−1} − ζ min(1 − ℓ, 0)^{2}
− w̄e(i)ηc(ℓ)^{η(1−μ)−1}ℓ^{(1−η)(1−μ)},
c(ℓ) = (1 + r̄)x + w̄e(i)(1 − ℓ) + χ − (1 + g)α(x, i),  (27)
J(ℓ) = −2w̄e(i)(1 − η)η(1 − μ)c(ℓ)^{η(1−μ)−1}ℓ^{(1−η)(1−μ)−1}
+ (1 − η)((1 − η)(1 − μ) − 1)c(ℓ)^{η(1−μ)}ℓ^{(1−η)(1−μ)−2}
+ w̄^{2}e(i)^{2}η(η(1 − μ) − 1)c(ℓ)^{η(1−μ)−2}ℓ^{(1−η)(1−μ)}
+ 2ζ min(1 − ℓ, 0).  (28)

The function f is the Euler equation in (23), the function c is consumption derived from the budget constraint, and the function J is the derivative of f with respect to ℓ. We start the iterations in Eq. (26) with an initial guess ℓ^{0}, and we stop when |ℓ^{k+1} − ℓ^{k}| is less than some tolerance parameter.

Let ℓ^{∗}(x, i; α) be the solution to Eq. (26). Then we can write the first-order conditions in Eqs. (22)-(25) in terms of one residual; i.e.,

R(x, i; α) = η(1 + g)c(ℓ^{∗}(x, i; α))^{η(1−μ)−1}ℓ^{∗}(x, i; α)^{(1−η)(1−μ)}

− β(1 + g)^{η(1−μ)}{∑_{j} π_{i,j} η(1 + r̄)c(ℓ^{∗}(α(x, i), j; α))^{η(1−μ)−1}

· ℓ^{∗}(α(x, i), j; α)^{(1−η)(1−μ)} + ζ min(α(x, i), 0)^{2}},

where c(·) is defined in (27). If we ignore the constraint on leisure (i.e., ℓ_{t} ≤ 1), the residual is given by

R(x, i; α) = η(1 + g)\{\frac{1 − η}{ηw̄e(i)}\}^{(1−η)(1−μ)}\{(w̄e(i) + (1 + r̄)x − (1 + g)α(x, i) + χ)η\}^{−μ}

− β(1 + g)^{η(1−μ)}\{∑_{j} π_{i,j} η(1 + r̄)\{\frac{1 − η}{ηw̄e(j)}\}^{(1−η)(1−μ)}\{(w̄e(j)

+ (1 + r̄)α(x, i) − (1 + g)α(α(x, i), j) + χ)η}^{−μ} + ζ min(α(x, i), 0)^{2}}.(29)

We show later that other calculations are also simplified if an interior solution for leisure is assumed. However, for some parameterizations, we need to impose the constraint on leisure. Therefore, we describe how to compute the equilibrium assuming the constraint is violated in some regions of the state space.

As in the inelastic labor case, the computational task is to find an approximation for α(x, i), say α^{h}(x, i), that implies R(x, i; α^{h}) ≈ 0 for all x and i. We follow the same procedure described for the inelastic labor case (i.e., we apply the finite element method).

2.2. Computing the distribution of assets

Computation of the cumulative distribution is the same for the elastic and inelastic labor cases. Therefore, the description in section 1.2 carries over to this case.

2.3. Computing an equilibrium

Assume that we have in place numerical algorithms for computing α(x, i) and H(x, i). We now describe an algorithm for computing the equilibrium interest rate and aggregate hours. All other prices and quantities can be determined from these two.

We start with an initial guess for N. With N fixed, we apply a bisection method to calculate the interest rate that clears the asset market. Assume that the equilibrium r is in the interval [r_{l}, r_{u}] and that our initial guess for the equilibrium r is \frac{1}{2}(r_{l}+r_{u}). If lump-sum taxes are assumed, then we set the taxes equal to γ + χ + (r − g)b. If proportional taxes are assumed, then we back out the tax rate from the government budget constraint. Given the tax rate and the guess for N, we can determine the after-tax interest rate (1 − τ_{y})r and the after-tax wage rate (1 − τ_{y})(1 − θ)/N. These two rates are used to compute the finite element approximation of α. To check to see that the asset market has cleared, we also compute a finite element approximation of H. H is used to calculate the mean asset holdings Eã_{t}. Given this estimate of Eã_{t}, we check to see if markets are cleared; i.e., if

Eã_{t} = \frac{θ}{r + δ} + b.  (30)

If the right-hand side of (30) is larger than the left-hand side, then we set r_{l} = r; otherwise, we set r_{u} = r. Iterations are made until Eq. (30) is satisfied to within some tolerance level. Once Eq. (30) is satisfied, we update our guess for N. The updating scheme is simply Newton-Raphson:

N^{k+1} = N^{k} − J(N^{k})^{−1}f(N^{k}),

where

f(N) = ∑_{i} ∫ e(i)(1 − ℓ^{∗}(x, i; α))dH(x, i) − N,  (31)

and J is the derivative of f with respect to N. Note that the first term on the right-hand side of (31) is a function of N because the optimal decision rules depend on it.

3. Test problems

3.1. Computing the consumer’s decision functions – a test case Assume that the household solves

max_{{c_{t},a_{t+1}}} ∑_{t=0}^{∞} β^{t}u(c_{t})
subject to c_{t} + a_{t+1} = (1 + r)a_{t} + w.

This specification assumes that there is no uncertainty (e_{t} = 1); therefore, wages are constant.

The dynamic program for this example involves the following form for Bellman’s equation:

v(x) = max_{0≤y≤Rx+w}{u(Rx + w − y) + βv(y)},  (32)

where y_{t} = x_{t+1} is the asset position next period and R = 1 + r is the gross return. A conjectured solution is as follows:

\{ y = 0  if x ∈ [0, m_{1}], ; y = \frac{−m_{1}^{2}}{m_{2}−m_{1}} + \frac{m_{1}}{m_{2}−m_{1}} x  if x ∈ (m_{1}, m_{2}], ; y = \frac{(m_{3}m_{1}−m_{2}^{2})}{m_{3}−m_{2}} + \frac{m_{2}−m_{1}}{m_{3}−m_{2}} x if x ∈ (m_{2}, m_{3}], ; ⋮  (33)

where m_{j}, j = 1, 2, …, will be calculated below. Note that the solution assumes that if x = m_{j+1}, then y = m_{j}.

The Lagrangian for the maximization in the right-hand side of Eq. (32) is given by

L = u(xR + w − y) + βv(y) + p(Rx + w − y) + qy,

where p and q are multipliers. The first-order conditions for this problem are

− u′(Rx + w − y) + βv′(y) − p + q = 0,
p ≥ 0, Rx + w − y ≥ 0, p(Rx + w − y) = 0,
q ≥ 0, y ≥ 0, qy = 0.  (34)

If we assume that the conjecture above is correct, then when x ∈ [0, m_{1}) we have the y ≥ 0 constraint binding. Therefore, if we assume that Rx + w > 0, then it must be true that y < Rx + w, p = 0, and

v′(0) = \frac{1}{β} u′(Rx + w) − \frac{q}{β} ≤ \frac{1}{β} u′(Rx + w).

Furthermore, from Bellman’s equation, we get

v(x) = u(Rx + w) + βv(0), for x ∈ (0, m_{1}) and v(0) = \frac{β}{1 − β} u(w),

and taking derivatives, we get

v′(x) = Ru′(Rx + w) < \frac{1}{β} u′(Rx + w),  (35)

since βR < 1.

Consider next the interval (m_{1}, m_{2}]. The conjectured solution is such that in this interval, the constraint y ≥ 0 is not binding. If we assume that Rx + w > y, then the first-order conditions imply

v′(y) = \frac{1}{β} u′(Rx + w − y).  (36)

If y = 0 at x = m_{1}, then

v′(0) = \frac{1}{β} u′(Rm_{1} + w).  (37)

Using Eq. (35) evaluated at x = 0 and Eq. (37), we get

u′(Rm_{1} + w) = βRu′(w),

which gives us an equation for m_{1}. For example, if u(c) = c^{1−μ}/(1 − μ), then

m_{1} = \frac{w(1 − (βR)^{\frac{1}{μ}})}{(βR)^{\frac{1}{μ}} R},

and y = 0 in the interval [0, m_{1}].

Now we want to compute the asset function for the next interval (m_{1}, m_{2}]. If the conjecture in Eq. (33) is correct, then Eq. (36) holds, as does

v′(x) = Ru′(Rx + w − g(x)),  (38)

which is the derivative of the value function once y is replaced by the optimal policy y = g(x). The conjecture assumes that y = m_{1} when x = m_{2}, and by Eq. (36) and Eq. (38), we have,

u′(Rm_{2} + w − m_{1}) = βRu′(Rm_{1} + w).

Note that this equation can be solved for m_{2}. If we follow the same logic for the remaining m’s, we find that, in general,

u′(Rm_{j+1} + w − m_{j}) = βRu′(Rm_{j} + w − m_{j−1}), j = 1, …, and m_{0} = 0.  (39)

Thus, given m_{1} and m_{2}, we can compute m_{3} and so on.

What we have done is conjectured a solution and derived the functions analytically. It is easy to show that the solution is, in fact, piecewise linear and that the conjecture is correct.

Now we consider the finite element approximation. Let β = 0.95, w = 1.0, r = 0.02, u(c) = c^{1−μ}/(1 − μ), and μ = 3. We can use the formula in Eq. (39) to derive the exact solution. In Figure 1a, we plot the true solution and the finite element approximation assuming the grid is given by m_{j}, j = 0, …, 14. Notice that the two functions are in very good agreement. In Figure 1b, we look at exactly the same picture, but we focus on the corner, where the constraints bind. The finite element approximation in this case does not assume any boundary constraints are imposed. Notice that there is a slight deviation of the approximate and exact solutions because the penalty function ensures only approximate satisfaction of the inequality constraints. The approximate solution in Figure 1c imposes that the asset holdings are zero at the first two grid points. Notice that the fit is exact. In Figures 2a-2c, we show the results of a more realistic experiment. In this case, we start with a grid that is “stretched” exponentially: the intervals get exponentially larger the farther they are from the origin. To obtain the results of Figure 2a, we use 15 grid points (n = 15). Notice in Figure 2a that the grid points of the approximate solution do not line up exactly with the kinks in the exact solution. Notice also that boundary conditions are not imposed. In Figure 2b, we add two grid points around the value of x where the decision function becomes positive. The picture is similar to Figure 2a. In Figure 2c, we impose that asset holdings are zero for the first three grid points. Notice how close we get to the exact solution even though we do not use the best possible discretization. Finally, in Figure 3, we plot the decision functions for three different sized grids: n = 5, n = 9, and n = 17. In all three cases, we impose some boundary conditions. If we plot the solutions for all x, the picture looks like Figure 1a. In Figure 3, we focus on only the corner. Notice how close the approximate function is to the exact function, even with n = 5.

3.2. Computing the distribution of assets – a test case

Suppose the productivities take on two possible values and the decision functions are given by

α(x, i) = \{ \begin{matrix}max(0, −0.25 + x), & if i = 1 \\ 0.5 + 0.5x, & if i = 2, \end{matrix}

with π_{1,1} = π_{2,2} = 0.8. In this case,

H(0, 1) = 0.8H(0.25, 1),
H(0, 2) = 0.2H(0.25, 1),
H(0.25, 1) = 0.8H(0.5, 1),
H(0.25, 2) = 0.2H(0.5, 1),
H(0.5, 1) = 0.8H(0.75, 1) + 0.2H(0, 2),
H(0.5, 2) = 0.2H(0.75, 1) + 0.8H(0, 2),
H(0.75, 1) = 0.8H(1.0, 1) + 0.2H(0.5, 2),
H(0.75, 2) = 0.2H(1.0, 1) + 0.8H(0.5, 2),
H(0.875, 1) = 0.8H(1.125, 1) + 0.2H(0.75, 2),
H(0.875, 2) = 0.2H(1.125, 1) + 0.8H(0.75, 2),
H(1.0, 1) = 0.5,
H(1.0, 2) = 0.5.

If we assume that H(x, i) = 0.5 for x > 1, then the above expressions can easily be solved. We can first determine H(0, j), H(0.25, j), H(0.5, j), and H(0.75, j) for j =1,2 by solving Ax = b, where

A = [ \begin{matrix}1 & −.8 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & −.8 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & −.8 & −.2 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & −.2 & 0 \\ 0 & −.2 & 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & −.2 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & −.2 & −.8 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & −.8 & 1 \end{matrix} ] b = [ \begin{matrix}0 \\ 0 \\ 0 \\ .4 \\ 0 \\ 0 \\ 0 \\ .1 \end{matrix} ] .

The solution is H(0, 1) = 0.225, H(0.25, 1) = 0.282, H(0.5, 1) = 0.352, H(0.75, 1) = 0.426, H(0, 2) = 0.056, H(0.25, 2) = 0.070, H(0.5, 2) = 0.130, and H(0.75, 2) = 0.204. Note that we can back out the other points from these solutions by applying the formula in Eq. (9). In Figure 4a, we plot the finite element approximation and the exact solution for a relatively coarse grid (n = 13) with evenly spaced partitions. Notice that a grid this coarse will not resolve the discontinuities. In Figure 4b, we refine the grid by doubling the number of partitions (i.e., n = 25); in this case, they are still evenly spaced. In Figures 4c and 4d, we recompute with an evenly spaced grid but with finer partitions (i.e., n = 49 and n = 97). In both of these cases, it is clear where the discontinuities occur. In Figure 4e, we plot the finite element approximation for a grid adapted to better resolve the discontinuities. Note that we get very good agreement with the above solution.

In Figures 4f and 4g, we focus on regions where jumps occur to illustrate how the method performs. Although the solution is relatively accurate, there are slight deviations at the discontinuities due to the fact that the method is essentially “centrally differencing.” That is, to determine the solution at x = .5, information is used from above and below this point. The result is nonmonotonicity near discontinuities. To fix this nonmonotonicity, we simply adjust the final solution as follows. In regions just to the left of large gradients (i.e., discontinuities), if we see Ĥ(i, j) < Ĥ(i − 1, j), then we set Ĥ(i, j) = Ĥ(i − 1, j). In regions to the right of discontinuities, if we see Ĥ(i, j) < Ĥ(i − 1, j), then we look at points k = i + 1, i + 2, …, until we find Ĥ(k, j) ≥ Ĥ(k − 1, j). Once we find such a point, we set Ĥ(l, j) = Ĥ(k − 1, j), l = i − 1, …, k − 1. The solution shown in Figure 4h applies this “trick.” If we compare this solution to Figure 4f, we see that the nonmonotonicity is eliminated and that the solution is closer to the exact solution.

Finally, we should note that for this example, if values of x greater than one are not included in the grid, then the approximation is very inaccurate. This inaccuracy is due to the fact that values of H(x, i), x < 1 are used to construct an estimate of H(x, i), x > 1. For this example, this extrapolation will lead to the inaccurate approximation shown in Figure 5. Notice that the approximation in Figure 5 lies above the true solution and violates the terminal condition that H(x, i) = 0.5.

Figure 1a. Decision functions for the test case with correct discretization, no boundary conditions imposed, and n = 15.

Figure 1b. Decision functions for the test case with correct discretization, no boundary conditions imposed, and n = 15. (Only x ∈ [0, .07] shown.)

Figure 1c. Decision functions for the test case with correct discretization, boundary conditions imposed, and n = 15. (Only x ∈ [0, .07] shown.)

Figure 2a. Decision functions for the test case without correct discretization, without boundary conditions imposed, and n = 15. (Only x ∈ [0, .07] shown.)

Figure 2b. Decision functions for the test case without correct discretization, without boundary conditions imposed, and n = 17. (Only x ∈ [0, .07] shown.)

Figure 2c. Decision functions for the test case without correct discretization, with boundary conditions imposed, and n = 17. (Only x ∈ [0, .07] shown.)

Figure 3. Decision functions for the test case without correct discretization, with boundary conditions imposed, and various n. (Only x ∈ [0, .09] shown.)

Figure 4a. Distribution for the test case with evenly spaced mesh and n = 13. Figure 4b. Distribution for the test case with evenly spaced mesh and n = 25.

Figure 4c. Distribution for the test case with evenly spaced mesh and n = 49. Figure 4d. Distribution for the test case with evenly spaced mesh and n = 97.

Figure 4e. Distribution for the test case with an adapted grid and n = 73. 24 Figure 4f. Distribution for the test case with an adapted grid. (Only x ∈ [0.475, 0.525] shown.)

Figure 4g. Distribution for the test case with an adapted grid. (Only x ∈ [0.84, 1.02] shown.)

Figure 4h. Distribution for the test case with an adapted grid and nonmonotonicities eliminated. (Only x ∈ [0.475, 0.525] shown.)

Figure 5. Inaccurate solution for the distribution.