CompEcon Toolbox:
DemDdp01
Mine management model
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.

Last updated: 2020-Sep-09

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

In [ ]:
if 'google.colab' in str(get_ipython()):
    print("This notebook is running on Google Colab. Installing the compecon package.")
    !pip install compecon
In [ ]:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

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$

In [ ]:
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.

In [ ]:
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.

In [ ]:
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:

In [ ]:
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.

In [ ]:
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.

In [ ]:
model = DDPmodel(f, g, delta)
model.solve()

By default, the solve() method uses Newton's algorithm. This and other default settings can be changed when solving the model. For example,

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.

In [ ]:
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.

In [ ]:
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

Plot Optimal Policy

In [ ]:
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.

In [ ]:
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.

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