{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Public Renewable Resource Model\n",
"\n",
"**Randall Romero Aguilar, PhD**\n",
"\n",
"This demo is based on the original Matlab demo accompanying the Computational Economics and Finance 2001 textbook by Mario Miranda and Paul Fackler.\n",
"\n",
"Original (Matlab) CompEcon file: **demdp01.m**\n",
"\n",
"Running this file requires the Python version of CompEcon. This can be installed with pip by running\n",
"\n",
" !pip install compecon --upgrade\n",
"\n",
"Last updated: 2022-Oct-23\n",
"
"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"## About\n",
"\n",
"Welfare maximizing social planner must decide how much of a renewable resource to harvest.\n",
"\n",
"- States\n",
" - s quantity of stock available\n",
"- Actions\n",
" - q quantity of stock harvested\n",
"- Parameters\n",
" - $\\alpha$ growth function parameter\n",
" - $\\beta$ growth function parameter\n",
" - $\\gamma$ relative risk aversion\n",
" - $\\kappa$ unit cost of harvest\n",
" - $\\delta$ discount factor"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from compecon import BasisChebyshev, DPmodel, DPoptions, qnwnorm"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Model parameters"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"α, β, γ, κ, δ = 4.0, 1.0, 0.5, 0.2, 0.9"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Steady-state\n",
"\n",
"The steady-state values for this model are"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Steady-State\n",
"\tStock = 7.3827\n",
"\tHarvest = 4.4938\n"
]
}
],
"source": [
"sstar = (α**2 - 1/δ**2)/(2*β) # steady-state stock\n",
"qstar = sstar - (δ*α-1)/(δ*β) # steady-state action\n",
"\n",
"print('Steady-State')\n",
"print(f\"\\t{'Stock':12s} = {sstar:8.4f}\")\n",
"print(f\"\\t{'Harvest':12s} = {qstar:8.4f}\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Numeric results"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### State space\n",
"The state variable is s=\"Stock\", which we restrict to $s\\in[6, 9]$. \n",
"\n",
"Here, we represent it with a Chebyshev basis, with $n=8$ nodes."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"n, smin, smax = 8, 6, 9\n",
"basis = BasisChebyshev(n, smin, smax, labels=['Stock'])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Action space\n",
"The choice variable q=\"Harvest\" must be nonnegative."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"def bounds(s, i=None, j=None):\n",
" return np.zeros_like(s), s[:]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Reward function\n",
"The reward function is the utility of harvesting $q$ units."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"def reward(s, q, i=None, j=None):\n",
" u = (q**(1-γ))/(1-γ)-κ*q\n",
" ux= q**(-γ)-κ\n",
" uxx = -γ*q**(-γ-1)\n",
" \n",
" return u, ux, uxx"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### State transition function\n",
"Next period, the stock will be equal that is $s' = \\alpha (s-q) - 0.5\\beta(s-q)^2$"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"def transition(s, q, i=None, j=None, in_=None, e=None):\n",
" sq = s-q\n",
" g = α*sq - 0.5*β*sq**2\n",
" gx = -α + β*sq\n",
" gxx = -β*np.ones_like(s)\n",
" return g, gx, gxx"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Model structure\n",
"The value of wealth $s$ satisfies the Bellman equation\n",
"\\begin{equation*}\n",
"V(s) = \\max_q\\left\\{\\frac{q^{1-\\gamma}}{1-\\gamma} -\\kappa q + \\delta V(s') \\right\\}\n",
"\\end{equation*}\n",
"\n",
"To solve and simulate this model,use the CompEcon class `DPmodel`"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"model = DPmodel(basis, reward, transition, bounds,\n",
" x=['Harvest'],\n",
" discount=δ)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Solving the model\n",
"\n",
"Solving the growth model by collocation"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Solving infinite-horizon model collocation equation by Newton's method\n",
"iter change time \n",
"------------------------------\n",
" 0 1.4e+01 0.0049\n",
" 1 1.9e+01 0.0132\n",
" 2 3.7e-01 0.0202\n",
" 3 1.1e-03 0.0257\n",
" 4 1.5e-08 0.0307\n",
" 5 1.7e-14 0.0337\n",
"Elapsed Time = 0.03 Seconds\n"
]
}
],
"source": [
"S = model.solve()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"`DPmodel.solve` returns a pandas `DataFrame` with the following data:"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" Stock | \n",
" value | \n",
" resid | \n",
" Harvest | \n",
"
\n",
" \n",
" | Stock | \n",
" | \n",
" | \n",
" | \n",
" | \n",
"
\n",
" \n",
" \n",
" \n",
" | 6.000000 | \n",
" 6.000000 | \n",
" 32.985613 | \n",
" 6.795958e-09 | \n",
" 3.349254 | \n",
"
\n",
" \n",
" | 6.037975 | \n",
" 6.037975 | \n",
" 32.998722 | \n",
" -1.561702e-09 | \n",
" 3.379828 | \n",
"
\n",
" \n",
" | 6.075949 | \n",
" 6.075949 | \n",
" 33.011736 | \n",
" -5.604697e-09 | \n",
" 3.410456 | \n",
"
\n",
" \n",
" | 6.113924 | \n",
" 6.113924 | \n",
" 33.024658 | \n",
" -6.686079e-09 | \n",
" 3.441138 | \n",
"
\n",
" \n",
" | 6.151899 | \n",
" 6.151899 | \n",
" 33.037489 | \n",
" -5.868451e-09 | \n",
" 3.471872 | \n",
"
\n",
" \n",
"
\n",
"