Randall Romero Aguilar, PhD
This demo is based on the original Matlab demo accompanying the Computational Economics and Finance 2001 textbook by Mario Miranda and Paul Fackler.
Original (Matlab) CompEcon file: demdp01.m
Running this file requires the Python version of CompEcon. This can be installed with pip by running
!pip install compecon --upgrade
Last updated: 2022-Oct-23
Welfare maximizing social planner must decide how much of a renewable resource to harvest.
s quantity of stock available
q quantity of stock harvested
import numpy as np
import matplotlib.pyplot as plt
from compecon import BasisChebyshev, DPmodel, DPoptions, qnwnorm
α, β, γ, κ, δ = 4.0, 1.0, 0.5, 0.2, 0.9
The steady-state values for this model are
sstar = (α**2 - 1/δ**2)/(2*β) # steady-state stock
qstar = sstar - (δ*α-1)/(δ*β) # steady-state action
print('Steady-State')
print(f"\t{'Stock':12s} = {sstar:8.4f}")
print(f"\t{'Harvest':12s} = {qstar:8.4f}")
Steady-State Stock = 7.3827 Harvest = 4.4938
The state variable is s="Stock", which we restrict to $s\in[6, 9]$.
Here, we represent it with a Chebyshev basis, with $n=8$ nodes.
n, smin, smax = 8, 6, 9
basis = BasisChebyshev(n, smin, smax, labels=['Stock'])
The choice variable q="Harvest" must be nonnegative.
def bounds(s, i=None, j=None):
return np.zeros_like(s), s[:]
The reward function is the utility of harvesting $q$ units.
def reward(s, q, i=None, j=None):
u = (q**(1-γ))/(1-γ)-κ*q
ux= q**(-γ)-κ
uxx = -γ*q**(-γ-1)
return u, ux, uxx
Next period, the stock will be equal that is $s' = \alpha (s-q) - 0.5\beta(s-q)^2$
def transition(s, q, i=None, j=None, in_=None, e=None):
sq = s-q
g = α*sq - 0.5*β*sq**2
gx = -α + β*sq
gxx = -β*np.ones_like(s)
return g, gx, gxx
The value of wealth $s$ satisfies the Bellman equation \begin{equation*} V(s) = \max_q\left\{\frac{q^{1-\gamma}}{1-\gamma} -\kappa q + \delta V(s') \right\} \end{equation*}
To solve and simulate this model,use the CompEcon class DPmodel
model = DPmodel(basis, reward, transition, bounds,
x=['Harvest'],
discount=δ)
Solving the growth model by collocation
S = model.solve()
Solving infinite-horizon model collocation equation by Newton's method iter change time ------------------------------ 0 1.4e+01 0.0049 1 1.9e+01 0.0132 2 3.7e-01 0.0202 3 1.1e-03 0.0257 4 1.5e-08 0.0307 5 1.7e-14 0.0337 Elapsed Time = 0.03 Seconds
DPmodel.solve
returns a pandas DataFrame
with the following data:
S.head()
Stock | value | resid | Harvest | |
---|---|---|---|---|
Stock | ||||
6.000000 | 6.000000 | 32.985613 | 6.795958e-09 | 3.349254 |
6.037975 | 6.037975 | 32.998722 | -1.561702e-09 | 3.379828 |
6.075949 | 6.075949 | 33.011736 | -5.604697e-09 | 3.410456 |
6.113924 | 6.113924 | 33.024658 | -6.686079e-09 | 3.441138 |
6.151899 | 6.151899 | 33.037489 | -5.868451e-09 | 3.471872 |
We are also interested in the shadow price of wealth (the first derivative of the value function).
S['shadow price'] = model.Value(S['Stock'],1)
fig1, ax = plt.subplots()
ax.set(title='Optimal Harvest Policy', xlabel='Stock', ylabel='Harvest')
ax.plot(S['Harvest'])
ax.plot(sstar, qstar,'wo')
ax.annotate(f"$s^*$ = {sstar:.2f}\n$q^*$ = {qstar:.2f}",[sstar, qstar], va='top');
fig2, ax = plt.subplots()
ax.set(title='Value Function', xlabel='Stock', ylabel='Value')
ax.plot(S.value);
fig3, ax = plt.subplots()
ax.set(title='Shadow Price Function', xlabel='Stock', ylabel='Shadow Price')
ax.plot(S['shadow price']);
fig4, ax = plt.subplots()
ax.set(title='Bellman Equation Residuals', xlabel='Stock', ylabel='Residual')
ax.axhline(0,color='w')
ax.plot(S['resid'])
ax.plot(basis.nodes[0], np.zeros_like(basis.nodes[0]), lw=0, marker='.', ms=12)
ax.ticklabel_format(style='sci', axis='y', scilimits=(-1,1))
We simulate 16 periods of the model starting from $s=s_{\min}$
T = 16
data = model.simulate(T, smin)
opts = dict(ha='right', va='top', fontsize='small')
fig5, ax = plt.subplots()
ax.set(title='State and Policy Paths',
xlabel='Period',
ylabel='Stock/Harvest',
xlim=[0, T + 0.5],
xticks=range(0,T+1,4))
ax.plot(data[['Stock', 'Harvest']])
ax.plot(T, sstar, color='C0', marker='*', ms=12)
ax.annotate(f"steady-state stock\n = {sstar:.2f}", [T, sstar-0.01], color='C0', **opts)
ax.plot(T, qstar, color='C1', marker='*', ms=12)
ax.annotate(f"steady-state harvest\n = {qstar:.2f}", [T, qstar-0.01], color='C1', **opts);