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-22
Profit maximizing owner of a commercial tree stand must decide when to clearcut the stand.
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
from compecon import DPmodel, BasisSpline, NLP
Assuming that the unit price of biomass is $p=1$, the cost to clearcut-replant is $\kappa=0.2$, the stand carrying capacity $s_{\max}=0.5$, biomass growth factor is $\gamma=10\%$ per period, and the annual discount factor $\delta=0.9$.
price = 1.0
κ = 0.2
smax = 0.5
ɣ = 0.1
δ = 0.9
The state variable is the stand biomass, $s\in [0,s_{\max}]$.
Here, we represent it with a cubic spline basis, with $n=200$ nodes.
n = 200
basis = BasisSpline(n, 0, smax, labels=['biomass'])
The action variable is $j\in\{0:\text{'keep'},\quad 1:\text{'clear cut'}\}$.
options = ['keep', 'clear-cut']
If the farmer clears the stand, the profit is the value of selling the biomass $ps$ minus the cost of clearing and replanting $\kappa$, otherwise the profit is zero.
def reward(s, x, i , j):
return (price * s - κ) * j
If the farmer clears the stand, it begins next period with $\gamma s_{\max}$ units. If he keeps the stand, then it grows to $s + \gamma (s_{\max} - s)$.
def transition(s, x, i, j, in_, e):
if j:
return np.full_like(s, ɣ * smax)
else:
return s + ɣ * (smax - s)
The value of the stand, given that it contains $s$ units of biomass at the beginning of the period, satisfies the Bellman equation
\begin{equation} V(s) = \max\left\{(ps-\kappa) + \delta V(\gamma s_{\max}),\quad \delta V[s+\gamma(s_{\max}-s)] \right\} \end{equation}To solve and simulate this model, use the CompEcon class DPmodel
.
model = DPmodel(basis, reward, transition,
discount=δ,
j=options)
S = model.solve()
Solving infinite-horizon model collocation equation by Newton's method iter change time ------------------------------ 0 3.1e-01 0.0209 1 1.1e-01 0.0299 2 1.1e-02 0.0389 3 4.6e-03 0.0469 4 1.6e-16 0.0559 Elapsed Time = 0.06 Seconds
The solve
method retuns a pandas DataFrame
, which can easily be used to make plots. To see the first 10 entries of S
, we use the head
method
S.head()
biomass | value | resid | j* | value[keep] | value[clear-cut] | |
---|---|---|---|---|---|---|
biomass | ||||||
0.000000 | 0.000000 | 0.067287 | -1.387779e-17 | keep | 0.067287 | -0.132713 |
0.000250 | 0.000250 | 0.067320 | -7.491474e-09 | keep | 0.067320 | -0.132462 |
0.000500 | 0.000500 | 0.067353 | -5.953236e-09 | keep | 0.067353 | -0.132212 |
0.000750 | 0.000750 | 0.067387 | -1.332984e-09 | keep | 0.067387 | -0.131962 |
0.001001 | 0.001001 | 0.067420 | 7.309628e-10 | keep | 0.067420 | -0.131712 |
To find the biomass level at which it is indifferent to keep or to clear cut, we interpolate as follows:
scrit = NLP(lambda s: model.Value_j(s).dot([1,-1])).broyden(0.0)
vcrit = model.Value(scrit)
print(f'When the biomass is {scrit:.5f} its value is {vcrit:.5f} regardless of the action taken by the manager')
When the biomass is 0.30669 its value is 0.17402 regardless of the action taken by the manager
?plt.subplots
Signature: plt.subplots( nrows=1, ncols=1, *, sharex=False, sharey=False, squeeze=True, subplot_kw=None, gridspec_kw=None, **fig_kw, ) Docstring: Create a figure and a set of subplots. This utility wrapper makes it convenient to create common layouts of subplots, including the enclosing figure object, in a single call. Parameters ---------- nrows, ncols : int, default: 1 Number of rows/columns of the subplot grid. sharex, sharey : bool or {'none', 'all', 'row', 'col'}, default: False Controls sharing of properties among x (*sharex*) or y (*sharey*) axes: - True or 'all': x- or y-axis will be shared among all subplots. - False or 'none': each subplot x- or y-axis will be independent. - 'row': each subplot row will share an x- or y-axis. - 'col': each subplot column will share an x- or y-axis. When subplots have a shared x-axis along a column, only the x tick labels of the bottom subplot are created. Similarly, when subplots have a shared y-axis along a row, only the y tick labels of the first column subplot are created. To later turn other subplots' ticklabels on, use `~matplotlib.axes.Axes.tick_params`. When subplots have a shared axis that has units, calling `~matplotlib.axis.Axis.set_units` will update each axis with the new units. squeeze : bool, default: True - If True, extra dimensions are squeezed out from the returned array of `~matplotlib.axes.Axes`: - if only one subplot is constructed (nrows=ncols=1), the resulting single Axes object is returned as a scalar. - for Nx1 or 1xM subplots, the returned object is a 1D numpy object array of Axes objects. - for NxM, subplots with N>1 and M>1 are returned as a 2D array. - If False, no squeezing at all is done: the returned Axes object is always a 2D array containing Axes instances, even if it ends up being 1x1. subplot_kw : dict, optional Dict with keywords passed to the `~matplotlib.figure.Figure.add_subplot` call used to create each subplot. gridspec_kw : dict, optional Dict with keywords passed to the `~matplotlib.gridspec.GridSpec` constructor used to create the grid the subplots are placed on. **fig_kw All additional keyword arguments are passed to the `.pyplot.figure` call. Returns ------- fig : `~.figure.Figure` ax : `.axes.Axes` or array of Axes *ax* can be either a single `~matplotlib.axes.Axes` object or an array of Axes objects if more than one subplot was created. The dimensions of the resulting array can be controlled with the squeeze keyword, see above. Typical idioms for handling the return value are:: # using the variable ax for single a Axes fig, ax = plt.subplots() # using the variable axs for multiple Axes fig, axs = plt.subplots(2, 2) # using tuple unpacking for multiple Axes fig, (ax1, ax2) = plt.subplots(1, 2) fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2) The names ``ax`` and pluralized ``axs`` are preferred over ``axes`` because for the latter it's not clear if it refers to a single `~.axes.Axes` instance or a collection of these. See Also -------- .pyplot.figure .pyplot.subplot .pyplot.axes .Figure.subplots .Figure.add_subplot Examples -------- :: # First create some toy data: x = np.linspace(0, 2*np.pi, 400) y = np.sin(x**2) # Create just a figure and only one subplot fig, ax = plt.subplots() ax.plot(x, y) ax.set_title('Simple plot') # Create two subplots and unpack the output array immediately f, (ax1, ax2) = plt.subplots(1, 2, sharey=True) ax1.plot(x, y) ax1.set_title('Sharing Y axis') ax2.scatter(x, y) # Create four polar axes and access them through the returned array fig, axs = plt.subplots(2, 2, subplot_kw=dict(projection="polar")) axs[0, 0].plot(x, y) axs[1, 1].scatter(x, y) # Share a X axis with each column of subplots plt.subplots(2, 2, sharex='col') # Share a Y axis with each row of subplots plt.subplots(2, 2, sharey='row') # Share both X and Y axes with all subplots plt.subplots(2, 2, sharex='all', sharey='all') # Note that this is the same as plt.subplots(2, 2, sharex=True, sharey=True) # Create figure number 10 with a single subplot # and clears it if it already exists. fig, ax = plt.subplots(num=10, clear=True) File: c:\programdata\anaconda3\lib\site-packages\matplotlib\pyplot.py Type: function
fig1, axs = plt.subplots(2,1,figsize=[8,5], gridspec_kw=dict(height_ratios=[3.5, 1]))
S[['value[keep]', 'value[clear-cut]']].plot(ax=axs[0])
axs[0].set(title='Action-Contingent Value Functions',
xlabel='',
ylabel='Value of Stand',
xticks=[])
ymin,ymax = axs[0].get_ylim()
axs[0].vlines(scrit,ymin, vcrit, linestyle=':')
axs[0].annotate('$s^*$', [scrit, ymin])
S['j*'].cat.codes.plot(ax=axs[1])
axs[1].set(title='Optimal Action',
ylabel='Action',
ylim=[-0.25,1.25],
yticks=[0,1],
yticklabels=options);
S['resid2'] = 100*S.resid / S.value
fig2, ax = plt.subplots()
ax.plot(S['resid2'])
ax.set(title='Bellman Equation Residual',
xlabel='',
ylabel='Percent Residual')
ax.hlines(0,0,smax,'k');
The path followed by the biomass is computed by the simulate()
method. Here we simulate 32 periods starting with a biomass level $s_0 = 0$.
H = model.simulate(32, 0.0)
fig3, ax = plt.subplots()
ax.plot(H['biomass'])
ax.set(title='Timber harvesting simulation',
xlabel='Period',
ylabel='Biomass');