# Timber Harvesting Model - Cubic Spline Approximation¶

Profit maximizing owner of a commercial tree stand must decide when to clearcut the stand.

In [13]:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
from compecon import DPmodel, BasisSpline, demo, NLP


### Model Parameters¶

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

In [2]:
price = 1.0
kappa = 0.2
smax  = 0.5
gamma = 0.1
delta = 0.9


### State Space¶

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.

In [3]:
n = 200
basis = BasisSpline(n, 0, smax, labels=['biomass'])


### Action Space¶

The action variable is $j\in\{0:\text{'keep'},\quad 1:\text{'clear cut'}\}$.

In [4]:
options = ['keep', 'clear-cut']


### Reward function¶

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.

In [5]:
def reward(s, x, i , j):
return (price * s - kappa) * j


### State Transition Function¶

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)$.

In [6]:
def transition(s, x, i, j, in_, e):
if j:
return np.full_like(s, gamma * smax)
else:
return s + gamma * (smax - s)


### Model Structure¶

The value of the stand, given that it contains $s$ units of biomass at the beginning of the period, satisfies the Bellman equation

$$V(s) = \max\left\{(ps-\kappa) + \delta V(\gamma s_{\max}),\quad \delta V[s+\gamma(s_{\max}-s)] \right\}$$

To solve and simulate this model, use the CompEcon class DPmodel.

In [26]:
model = DPmodel(basis, reward, transition,
discount=delta,
j=options)

S = model.solve()

Solving infinite-horizon model collocation equation by Newton's method
iter change       time
------------------------------
0       3.1e-01    0.0000
1       1.1e-01    0.0186
2       1.1e-02    0.0221
3       4.6e-03    0.0221
4       1.9e-16    0.0378
Elapsed Time =    0.04 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

In [27]:
S.head()

Out[27]:
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

## Analysis¶

To find the biomass level at which it is indifferent to keep or to clear cut, we interpolate as follows:

In [23]:
scrit = NLP(lambda s: model.Value_j(s).dot([1,-1])).broyden(0.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


### Plot the Value Function and Optimal Action¶

In [37]:
fig1 = plt.figure(figsize=[8,5])
gs = gridspec.GridSpec(2, 1,height_ratios=[3.5, 1])

S[['value[keep]', 'value[clear-cut]']].plot(ax= plt.subplot(gs[0]))
plt.title('Action-Contingent Value Functions')
plt.xlabel('')
plt.ylabel('Value of Stand')
plt.xticks([])

ymin,ymax = plt.ylim()
plt.vlines(scrit,ymin, vcrit,linestyle=':')
demo.annotate(scrit, ymin,'$s^*$',ms=0)

S['j*'].cat.codes.plot(ax=plt.subplot(gs[1]))
plt.title('Optimal Action')
plt.ylabel('Action')
plt.ylim([-0.25,1.25])
plt.yticks([0,1],options)
plt.legend([])

Out[37]:
<matplotlib.legend.Legend at 0x23f3f129588>

### Plot Residuals¶

In [34]:
S['resid2'] = 100*S.resid / S.value

fig2 = demo.figure('Bellman Equation Residual','','Percent Residual')
plt.plot(S['resid2'])
plt.hlines(0,0,smax,'k')

Out[34]:
<matplotlib.collections.LineCollection at 0x23f3eeb40f0>

### Simulation¶

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

In [33]:
H = model.simulate(32, 0.0)

fig3 = demo.figure('Timber harvesting simulation','Period','Biomass')
plt.plot(H['biomass'])

Out[33]:
[<matplotlib.lines.Line2D at 0x23f3ee5e5c0>]
In [ ]:
#demo.savefig([fig1, fig2, fig3])