Timber Harvesting – using cubic splines
Contents
Timber Harvesting – using cubic splines¶
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: demdp01b.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
About¶
Profit maximizing owner of a commercial tree stand must decide when to clearcut the stand.
import numpy as np
from compecon import NLP, BasisSpline
import matplotlib.pyplot as plt
price = 1.0 # price of biomass
κ = 0.2 # clearcut-replant cost
smax = 0.5 # stand carrying capacity
ɣ = 0.1 # biomass growth parameter
δ = 0.9 # discount factor
Code the growth function¶
def h(s): return np.array(s + ɣ*(smax - s))
SOLUTION¶
Code the approximant and the residual¶
ns = 200
vhat = BasisSpline(ns,0,smax,k=3)
def vhat1(s): return price*s - κ + δ * vhat(h(0))
def vhat0(s): return δ * vhat(h(s))
def resid(c,s=vhat.nodes):
vhat.c = c
return vhat(s) - np.maximum(vhat0(s), vhat1(s))
Solve collocation equation¶
cc = NLP(resid).broyden(vhat.c)
Compute critical biomass¶
scrit = NLP(lambda s: vhat0(s)-vhat1(s)).broyden(0.0)
ANALYSIS¶
Compute refined state grid¶
ss = np.linspace(0,smax,1000)
Plot Conditional Value Functions¶
fig1, ax = plt.subplots()
ax.plot(ss, vhat0(ss), label='Grow')
ax.plot(ss, vhat1(ss), label='Clear-Cut')
ax.set(title='Conditional Value Functions',
xlabel='Biomass',
ylabel='Value of Stand')
ax.legend()
vcrit = vhat(scrit)
ymin = ax.get_ylim()[0]
ax.vlines(scrit, ymin,vcrit,'grey',linestyles='--')
ax.annotate('$s^*$', [scrit,ymin])
ax.plot(scrit,vcrit,'wo')
print(f'Optimal Biomass Harvest Level = {scrit:.4f}')
Optimal Biomass Harvest Level = 0.3067
Plot Value Function Residual¶
fig2, ax = plt.subplots()
ax.set(title='Bellman Equation Residual', xlabel='Biomass', ylabel='Percent Residual')
ax.plot(ss, 100*resid(cc,ss) / vhat(ss))
ax.hlines(0,0,smax,linestyles='--')
<matplotlib.collections.LineCollection at 0x1b6720651c0>
Compute ergodic mean annual harvest¶
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