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
Public authority must decide how much water to release from a reservoir so as to maximize benefits derived from agricultural and recreational uses.
s reservoiur level at beginning of summer
x quantity of water released for irrigation
import numpy as np
import matplotlib.pyplot as plt
from compecon import BasisChebyshev, DPmodel, DPoptions, qnwlogn
import seaborn as sns
import pandas as pd
a0, a1, b0, b1 = 1, -2, 2, -3
μ, σ, δ = 1.0, 0.2, 0.9
The deterministic steady-state values for this model are
xstar = 1.0 # action
sstar = 1.0 + (a0*(1-δ)/b0)**(1/b1) # stock
The state variable is s="Reservoir Level", which we restrict to $s\in[2, 8]$.
Here, we represent it with a Chebyshev basis, with $n=15$ nodes.
n, smin, smax = 15, 2, 8
basis = BasisChebyshev(n, smin, smax, labels=['Reservoir'])
m = 3 #number of rainfall shocks
e, w = qnwlogn(m, np.log(μ)-σ**2/2,σ**2) # rainfall shocks and proabilities
The choice variable x="Irrigation" must be nonnegative.
def bounds(s, i=None, j=None):
return np.zeros_like(s), 1.0*s
The reward function is
def reward(s, x, i=None, j=None):
sx = s-x
u = (a0/(1+a1))*x**(1+a1) + (b0/(1+b1))*sx**(1+b1)
ux = a0*x**a1 - b0*sx**b1
uxx = a0*a1*x**(a1-1) + b0*b1*sx**(b1-1)
return u, ux, uxx
Next period, reservoir level wealth will be equal to current level minus irrigation plus random rainfall:
def transition(s, x, i=None, j=None, in_=None, e=None):
g = s - x + e
gx = -np.ones_like(s)
gxx = np.zeros_like(s)
return g, gx, gxx
The value of wealth $s$ satisfies the Bellman equation \begin{equation*} V(s) = \max_x\left\{\left(\frac{a_0}{1+a_1}\right)x^{1+a1} + \left(\frac{b_0}{1+b_1}\right)(s-x)^{1+b1}+ \delta V(s-x+e) \right\} \end{equation*}
To solve and simulate this model,use the CompEcon class DPmodel
water_model = DPmodel(basis, reward, transition, bounds,e=e,w=w,
x=['Irrigation'],
discount=δ )
The DPmodel.lqapprox
solves the linear-quadratic approximation, in this case arround the steady-state. It returns a LQmodel which works similar to the DPmodel object. We compute the solution at the basis nodes to use it as initial guess for the Newton's methods.
water_lq = water_model.lqapprox(sstar,xstar).solution(basis.nodes)
Solving the growth model by collocation. We take the values for the Value and Policy functions at the basis nodes obtained from the linear-quadratic approximation as initial guess values for the Newton's method.
S = water_model.solve(v=water_lq['value'], x=water_lq['Irrigation'])
S.head()
Solving infinite-horizon model collocation equation by Newton's method iter change time ------------------------------ 0 7.1e-01 0.0160 1 1.2e-01 0.0302 2 1.2e-02 0.0420 3 1.3e-04 0.0510 4 1.9e-08 0.0596 5 4.6e-15 0.0656 Elapsed Time = 0.07 Seconds
Reservoir | value | resid | Irrigation | |
---|---|---|---|---|
Reservoir | ||||
2.000000 | 2.000000 | -14.086357 | 8.445569e-07 | 0.628126 |
2.040268 | 2.040268 | -13.985981 | -6.417778e-07 | 0.638659 |
2.080537 | 2.080537 | -13.888843 | -7.616363e-07 | 0.649069 |
2.120805 | 2.120805 | -13.794753 | -3.379168e-07 | 0.659357 |
2.161074 | 2.161074 | -13.703538 | 1.731421e-07 | 0.669527 |
DPmodel.solve
returns a pandas DataFrame
with the following data:
We are also interested in the shadow price of wealth (the first derivative of the value function).
S['shadow price'] = water_model.Value(S['Reservoir'],1)
S.head()
Reservoir | value | resid | Irrigation | shadow price | |
---|---|---|---|---|---|
Reservoir | |||||
2.000000 | 2.000000 | -14.086357 | 8.445569e-07 | 0.628126 | 2.534519 |
2.040268 | 2.040268 | -13.985981 | -6.417778e-07 | 0.638659 | 2.451652 |
2.080537 | 2.080537 | -13.888843 | -7.616363e-07 | 0.649069 | 2.373665 |
2.120805 | 2.120805 | -13.794753 | -3.379168e-07 | 0.659357 | 2.300174 |
2.161074 | 2.161074 | -13.703538 | 1.731421e-07 | 0.669527 | 2.230828 |
fig1, ax = plt.subplots()
ax.set(title='Optimal Irrigation Policy', xlabel='Reservoir Level', ylabel='Irrigation')
ax.plot(S['Irrigation'])
ax.plot(sstar, xstar, 'o', color='C5');
ax.annotate(f'$s^*$ = {sstar:.2f}\n$x^*$ = {xstar:.2f}', [sstar, xstar], va='top', color='C5');
fig2, ax = plt.subplots()
ax.set(title='Value Function', xlabel='Reservoir Level', ylabel='Value')
ax.plot(S['value']);
fig3, ax = plt.subplots()
ax.set(title='Shadow Price Function', xlabel='Reservoir Level', ylabel='Shadow Price')
ax.plot(S['shadow price']);
fig4, ax = plt.subplots()
ax.set(title='Bellman Equation Residual', xlabel='Reservoir Level', 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='.', markersize=9)
ax.ticklabel_format(style='sci', axis='y', scilimits=(-1,1))
We simulate 21 periods of the model starting from $s=s_{\min}$
T = 31
nrep = 100_000
data = water_model.simulate(T, np.tile(smin,nrep))
subdata = data.query('_rep < 100')
fig6, ax = plt.subplots()
ax.set(title='Simulated and Expected Reservoir Level',
xlabel='Year',
ylabel='Reservoir Level',
xlim=[0, T + 0.5])
ax.plot(data[['time','Reservoir']].groupby('time').mean())
ax.plot(subdata.pivot('time','_rep','Reservoir'),lw=1, color='C0', alpha=0.2)
ax.annotate(f'steady-state reservoir\n = {sstar:.2f}',[T, sstar], color='C0', ha='right', va='top')
Text(31, 3.7144176165949068, 'steady-state reservoir\n = 3.71')
fig7, ax = plt.subplots()
ax.set(title='Simulated and Expected Irrigation',
xlabel='Year',
ylabel='Irrigation',
xlim=[0, T + 0.5])
ax.plot(subdata.pivot('time','_rep','Irrigation'),lw=1, color='C2', alpha=0.2)
ax.plot(data[['time','Irrigation']].groupby('time').mean(), color='C2', lw=4)
ax.annotate(f'steady-state irrigation\n = {xstar:.2f}',[T, xstar], color='C2', ha='right', va='top')
Text(31, 1.0, 'steady-state irrigation\n = 1.00')
subdata = data[data['time']==T][['Reservoir','Irrigation']]
stats =pd.DataFrame({'Deterministic Steady-State': [sstar, xstar],
'Ergodic Means': subdata.mean(),
'Ergodic Standard Deviations': subdata.std()}).T
stats
Reservoir | Irrigation | |
---|---|---|
Deterministic Steady-State | 3.714418 | 1.000000 |
Ergodic Means | 3.764014 | 0.999627 |
Ergodic Standard Deviations | 0.358049 | 0.062222 |
fig8, ax = plt.subplots()
ax.set(title='Ergodic Reservoir and Irrigation Distribution',xlabel='Wealth',ylabel='Probability')
sns.kdeplot(subdata['Reservoir'], ax=ax, label='Reservoir')
sns.kdeplot(subdata['Irrigation'], ax=ax, label='Irrigation')
ax.legend();