# Timber Harvesting Model - Cubic Spline Approximation

**Randall Romero Aguilar, PhD**

This demo is based on the original Matlab demo accompanying the  <a href="https://mitpress.mit.edu/books/applied-computational-economics-and-finance">Computational Economics and Finance</a> 2001 textbook by Mario Miranda and Paul Fackler.

Original (Matlab) CompEcon file: **demdp01.m**

Running this file requires the Python version of CompEcon. This can be installed with pip by running

    !pip install compecon --upgrade

<i>Last updated: 2021-Oct-01</i>
<hr>

## About

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


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


```python
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.



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



```python
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.


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



```python
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

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

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


```python
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.0160
       1       1.1e-01    0.0209
       2       1.1e-02    0.0259
       3       4.6e-03    0.0309
       4       1.6e-16    0.0389
    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


```python
S.head()
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>biomass</th>
      <th>value</th>
      <th>resid</th>
      <th>j*</th>
      <th>value[keep]</th>
      <th>value[clear-cut]</th>
    </tr>
    <tr>
      <th>biomass</th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
      <th></th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0.000000</th>
      <td>0.000000</td>
      <td>0.067287</td>
      <td>-1.387779e-17</td>
      <td>keep</td>
      <td>0.067287</td>
      <td>-0.132713</td>
    </tr>
    <tr>
      <th>0.000250</th>
      <td>0.000250</td>
      <td>0.067320</td>
      <td>-7.491474e-09</td>
      <td>keep</td>
      <td>0.067320</td>
      <td>-0.132462</td>
    </tr>
    <tr>
      <th>0.000500</th>
      <td>0.000500</td>
      <td>0.067353</td>
      <td>-5.953236e-09</td>
      <td>keep</td>
      <td>0.067353</td>
      <td>-0.132212</td>
    </tr>
    <tr>
      <th>0.000750</th>
      <td>0.000750</td>
      <td>0.067387</td>
      <td>-1.332984e-09</td>
      <td>keep</td>
      <td>0.067387</td>
      <td>-0.131962</td>
    </tr>
    <tr>
      <th>0.001001</th>
      <td>0.001001</td>
      <td>0.067420</td>
      <td>7.309628e-10</td>
      <td>keep</td>
      <td>0.067420</td>
      <td>-0.131712</td>
    </tr>
  </tbody>
</table>
</div>



## Analysis

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


```python
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


```python
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([])
```

### Plot Residuals


```python
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')
```

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


```python
H = model.simulate(32, 0.0)

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

```


```python
#demo.savefig([fig1, fig2, fig3])
```

d8ac9d40d9d8e76ebfafb6b26507441bcc337b65
ipynb
Jupyter Notebook
_build/jupyter_execute/notebooks/dp/01 Timber Harvesting Model - Cubic Spline Approximation.ipynb
randall-romero/CompEcon-python
c7a75f57f8472c972fddcace8ff7b86fee049d29
MIT
2016-12-14T13:21:27.000Z
2020-08-23T21:04:34.000Z
_build/jupyter_execute/notebooks/dp/01 Timber Harvesting Model - Cubic Spline Approximation.ipynb
randall-romero/CompEcon
c7a75f57f8472c972fddcace8ff7b86fee049d29
MIT
2017-09-10T04:48:54.000Z
2018-03-31T01:36:46.000Z
_build/jupyter_execute/notebooks/dp/01 Timber Harvesting Model - Cubic Spline Approximation.ipynb
randall-romero/CompEcon-python
c7a75f57f8472c972fddcace8ff7b86fee049d29
MIT
2017-02-25T08:10:38.000Z
2020-05-15T09:49:16.000Z
Qwen/Qwen-72B
1. YES
2. YES