# Mine management model

**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: **demddp03.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

A mine operator must decide how much ore to extract from a mine that will be shut down and abandoned after $T$ years of operation. The price of extracted ore is $p$ dollars per ton, and the total cost of extracting $x$ tons of ore in any year, given that the mine contains $s$ tons at the beginning of the year, is $c(s, x)$ dollars. The mine currently contains $\bar{s}$ tons of ore. Assuming the amount of ore extracted in any year must be an integer number of tons, what extraction schedule maximizes profits?

## Initial tasks



```python
import numpy as np
import pandas as pd
from compecon import DDPmodel, getindex
```

## Model Parameters

Assuming a market price $p=1$, initial stock of ore $\bar{s}=100$, and annual discount factor $\delta = 0.9$


```python
price = 1
sbar  = 100
delta = 0.9
```

### State Space

This is a finite horizon, deterministic model with time $t$ measured in years. The state
variable $s \in \{0, 1, 2, \dots, \bar{s}\}$ is the amount of ore remaining in the mine at the beginning of the year, measured in tons.


```python
S = np.arange(sbar + 1)      # vector of states
n = S.size                   # number of states
```

### Action Space

The action variable $x \in \{0, 1, 2, \dots, s\}$ is the amount of ore extracted over the year, measured in tons.


```python
X = np.arange(sbar + 1)      # vector of actions
m = X.size                   # number of actions
```

### Reward Function

The reward function is $f(s, x) = px − c(s, x)$. The cost of extraction is $c(s, x) = \frac{x^2}{1+s}$.

Here, the reward is set to negative infinity if the extraction level exceeds the available stock in order to preclude the choice of an infeasible
action:


```python
f = np.full((m, n), -np.inf)
for c, s in enumerate(S):
    for r, x in enumerate(X):
        if x <= s:
            f[r, c] = price * x - (x ** 2) / (1 + s)
```

### State Transition Function

The state transition function is $g(s, x) = s − x$

Here, the routine `getindex` is used to find the index of the following period’s state.


```python
g = np.empty_like(f)
for r, x in enumerate(X):
    snext = S - x
    g[r] = getindex(snext, S)
```

## Model Structure

The value of the mine, given that it contains $s$ tons of ore at the beginning of year $t$, satisfies the Bellman equation

\begin{equation} V_t(s) = max_{x\in\{0,1,\dots,s\}} \left\{px−c(s, x) + \delta V_{t+1}(s−x)\right\} \end{equation}

subject to the terminal condition $V_{T+1}(s) = 0$

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


```python
model = DDPmodel(f, g, delta)
model.solve()
```




    A deterministic discrete state, discrete action, dynamic model.
    There are 101 possible actions over 101 possible states



By default, the `solve()` method uses Newton's algorithm. This and other default settings can be changed when solving the model. For example,
```python
model.solve(algorithm='funcit', print=True)
```
solves the model by function iteration, printing a summary of each iteration to screen.

In either case, `solve()` updates the model itself, storing the $n$ vector of values `.value`, the $n$ vector of optimal actions `.policy`, and the $n\times n$ controlled state `.transition` probability.



```python
solution = pd.DataFrame({
    'Stock': S,
    'Extraction': X[model.policy], 
    'Value': model.value}).set_index('Stock')
```

## Analysis

### Simulate Model
The path followed by the stock level over time is computed by the `simulate()` method. Here, the simulation assumes an initial stock level of 100 and 15 years in duration.


```python
sinit = S.max()
nyrs = 15
t = np.arange(nyrs + 1)
spath, xpath = model.simulate(sinit, nyrs)

simul = pd.DataFrame({
    'Year': t,
    'Stock': S[spath],
    'Extraction': X[xpath]}).set_index('Year')
simul
```




<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>Stock</th>
      <th>Extraction</th>
    </tr>
    <tr>
      <th>Year</th>
      <th></th>
      <th></th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>100</td>
      <td>24</td>
    </tr>
    <tr>
      <th>1</th>
      <td>76</td>
      <td>18</td>
    </tr>
    <tr>
      <th>2</th>
      <td>58</td>
      <td>14</td>
    </tr>
    <tr>
      <th>3</th>
      <td>44</td>
      <td>11</td>
    </tr>
    <tr>
      <th>4</th>
      <td>33</td>
      <td>8</td>
    </tr>
    <tr>
      <th>5</th>
      <td>25</td>
      <td>6</td>
    </tr>
    <tr>
      <th>6</th>
      <td>19</td>
      <td>5</td>
    </tr>
    <tr>
      <th>7</th>
      <td>14</td>
      <td>3</td>
    </tr>
    <tr>
      <th>8</th>
      <td>11</td>
      <td>3</td>
    </tr>
    <tr>
      <th>9</th>
      <td>8</td>
      <td>2</td>
    </tr>
    <tr>
      <th>10</th>
      <td>6</td>
      <td>2</td>
    </tr>
    <tr>
      <th>11</th>
      <td>4</td>
      <td>1</td>
    </tr>
    <tr>
      <th>12</th>
      <td>3</td>
      <td>1</td>
    </tr>
    <tr>
      <th>13</th>
      <td>2</td>
      <td>1</td>
    </tr>
    <tr>
      <th>14</th>
      <td>1</td>
      <td>1</td>
    </tr>
    <tr>
      <th>15</th>
      <td>0</td>
      <td>0</td>
    </tr>
  </tbody>
</table>
</div>



### Plot Optimal Policy


```python
ax = solution['Extraction'].plot(title='Optimal Extraction Policy')
ax.set(ylabel='Extraction');
```

### Plot Value Function
The value of the firm is very nearly proportional to the stock level.


```python
ax = solution['Value'].plot(title='Optimal Value Function')
ax.set(ylabel='Value');
```

### Plot State and Extraction Path
As seen in this figure, the content of the mine is optimally exhausted in 15 years.


```python
ax = simul.plot(subplots=True, layout=[1,2], legend=None)
for i, lab in enumerate(simul.columns):
    ax[0, i].set(ylabel=lab);
```

eac0b012054a4d0031b6c3517b7ed5120fef2e44
ipynb
Jupyter Notebook
notebooks/ddp/01 Mine management model.ipynb
randall-romero/CompEcon-python
c7a75f57f8472c972fddcace8ff7b86fee049d29
MIT
2016-12-14T13:21:27.000Z
2020-08-23T21:04:34.000Z
notebooks/ddp/01 Mine management model.ipynb
randall-romero/CompEcon
c7a75f57f8472c972fddcace8ff7b86fee049d29
MIT
2017-09-10T04:48:54.000Z
2018-03-31T01:36:46.000Z
_build/html/_sources/notebooks/ddp/01 Mine management model.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