# Timber Harvesting Model - Cubic Spline Approximation¶

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

In [1]:
%matplotlib inline
from warnings import simplefilter
simplefilter('ignore')
from compecon import DPmodel, BasisSpline
from demos.setup import np, plt, demo


### 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'}\}$.

### 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 [4]:
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 [5]:
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 [6]:
model = DPmodel(basis, reward, transition,
discount=delta,
j=['keep', 'clear-cut'])

S = model.solve()

Solving infinite-horizon model collocation equation by Newton's method
iter change       time
------------------------------
0       3.1e-01    0.0100
1       1.1e-01    0.0180
2       1.1e-02    0.0260
3       4.6e-03    0.0325
4       1.9e-16    0.0405
Elapsed Time =    0.04 Seconds


The solve method retuns a pandas DataFrame, which can easily be used with the ggplot package.

## Analysis¶

### Plot the Value Function¶

In [7]:
demo.qplot('biomass', 'value_j','j',
data=S,
main='Action-Contingent Value functions',
ylab='Value of Stand'
)

Out[7]:
<ggplot: (152050177505)>

### Plot Residuals¶

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

demo.qplot('biomass', 'resid2',
data=S,
geom='line',
main='Bellman Equation Residual',
ylab='Percent Residual'
)

Out[8]:
<ggplot: (152050202323)>

### Simulation¶

The path followed by the biomass is computed by the simulate() method. Here we simulate 31 periods starting with a biomass level $s_0 = 0$.

In [9]:
H = model.simulate(31, 0.0)
demo.qplot('time','biomass',geom='line', data=H)

Out[9]:
<ggplot: (-9223371884804430041)>