DEMDP01a

Timber Harvesting -- using 2 nodes

In [1]:
import numpy as np
from compecon import NLP, demo
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]:
snodes = np.array([0.2, 0.4])
In [4]:
def h(s): return s + gamma*(smax - s)

SOLUTION

Code the approximant and the residual

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

def resid(c,s=snodes): return vhat(c,s) - np.maximum(vhat0(c,s), vhat1(c,s))

Solve collocation equation

In [6]:
cc = NLP(resid).broyden(np.zeros(2))

Compute critical biomass

In [7]:
scrit = NLP(lambda s: vhat0(cc,s)-vhat1(cc,s)).broyden(0.0)[0]
In [8]:
scrit
Out[8]:
0.3619999999996372

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(cc,ss),label='Grow')
plt.plot(ss,vhat1(cc,ss),label='Clear-Cut')
plt.legend()

vcrit = vhat(cc,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.3620

Plot Value Function Residual

In [11]:
fig2 = demo.figure('Bellman Equation Residual', 'Biomass', 'Percent Residual')
plt.plot(ss, 100*resid(cc,ss) / vhat(cc,ss))
plt.hlines(0,0,smax,linestyles='--')
plt.plot(snodes,resid(cc),'ro')
Out[11]:
[<matplotlib.lines.Line2D at 0x203e554d648>]

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.0311 after 13 iterations
In [13]:
#demo.savefig([fig1,fig2])