# Timber Harvesting -- using cubic splines¶

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

In [2]:
price = 1.0  # price of biomass
kappa = 0.2  # clearcut-replant cost
smax  = 0.5  # stand carrying capacity
gamma = 0.1  # biomass growth parameter
delta = 0.9  # discount factor


### Code the growth function¶

In [3]:
def h(s): return np.array(s + gamma*(smax - s))


## SOLUTION¶

### Code the approximant and the residual¶

In [4]:
ns = 200
vhat = BasisSpline(ns,0,smax,k=3)

In [5]:
def vhat1(s): return price*s - kappa + delta * vhat(h(0))
def vhat0(s): return delta * vhat(h(s))

In [6]:
def resid(c,s=vhat.nodes):
vhat.c = c
return vhat(s) - np.maximum(vhat0(s), vhat1(s))


### Solve collocation equation¶

In [7]:
cc = NLP(resid).broyden(vhat.c)


### Compute critical biomass¶

In [8]:
scrit = NLP(lambda s: vhat0(s)-vhat1(s)).broyden(0.0)[0]


## ANALYSIS¶

### Compute refined state grid¶

In [9]:
ss = np.linspace(0,smax,1000)


### Plot Conditional Value Functions¶

In [10]:
fig1 =demo.figure('Conditional Value Functions','Biomass','Value of Stand')
plt.plot(ss,vhat0(ss),label='Grow')
plt.plot(ss,vhat1(ss),label='Clear-Cut')
plt.legend()

vcrit = vhat(scrit)
ymin = plt.ylim()[0]
plt.vlines(scrit, ymin,vcrit,'grey',linestyles='--')
demo.annotate(scrit,ymin,'$s^*$',ms=10)
demo.bullet(scrit,vcrit)
print(f'Optimal Biomass Harvest Level = {scrit:.4f}')

Optimal Biomass Harvest Level = 0.3067


### Plot Value Function Residual¶

In [11]:
fig2 = demo.figure('Bellman Equation Residual', 'Biomass', 'Percent Residual')
plt.plot(ss, 100*resid(cc,ss) / vhat(ss))
plt.hlines(0,0,smax,linestyles='--')

Out[11]:
<matplotlib.collections.LineCollection at 0x1c59c8d9348>

### Compute ergodic mean annual harvest¶

In [12]:
s = h(0)
for n in range(100):
if s > scrit: break
s = h(s)

print(f'Ergodic Mean Annual Harvest = {s/n:.4f} after {n+1} iterations')

Ergodic Mean Annual Harvest = 0.0362 after 10 iterations

In [13]:
#demo.savefig([fig1,fig2])