{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Production Management 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",
"**WARNING** This demo is not running. Problem with dpmodel.\n",
"\n",
"TODO: Fix error in dpmodel.\n",
"\n",
"Last updated: 2022-Oct-23\n",
"
"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"## About\n",
"\n",
"Profit maximizing entrepeneur must decide how much to produce, subject to production adjustment costs.\n",
"\n",
"- States\n",
" - i market price (discrete)\n",
" - s lagged production (continuous)\n",
"- Actions\n",
" - x current production\n",
"- Parameters\n",
" - $\\alpha$ -- marginal adjustment cost\n",
" - $\\beta$ -- marginal production cost parameters\n",
" - pbar -- long-run average market price \n",
" - $\\mu$ -- mean log price\n",
" - $\\sigma$ -- market price shock standard deviation\n",
" - $\\delta$ -- discount factor\n",
" \n",
" "
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from compecon import BasisSpline, DPmodel, DPoptions, qnwlogn, BasisChebyshev\n",
"import seaborn as sns\n",
"import pandas as pd"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Model parameters\n"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"α, β0, β1, pbar = 0.01, 0.8, 0.03, 1.0 \n",
"σ, δ = 0.2, 0.9\n",
"μ = np.log(pbar) - σ**2 / 2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Continuous state shock distribution"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"m = 3 #number of market price shocks\n",
"p, w = qnwlogn(m, μ, σ**2) \n",
"q = np.repeat(w,3).reshape(3,3).T"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### State space\n",
"The state variable is s=\"lagged production\", which we restrict to $s\\in[0, 20]$. \n",
"\n",
"Here, we represent it with a cubic spline basis, with $n=50$ nodes."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"n, smin, smax = 5, 0.0, 20.0\n",
"basis = BasisChebyshev(n, smin, smax, labels=['lagged production'])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The discrete state is given by the price *p*"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"prices = ['p_low', 'p_mid', 'p_high']"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Action space\n",
"The choice variable x=\"current production\" must be nonnegative."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"def bounds(s, i, j=None):\n",
" return np.zeros_like(s), np.inf + np.zeros_like(s)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Reward function\n",
"The reward function is "
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"def reward(s, q, i, j=None):\n",
" u = p[i]*q - (β0*q + 0.5*β1*q**2) - 0.5*α*((q-s)**2)\n",
" ux = p[i] - β0 - β1*q - α*(q-s)\n",
" uxx = (-β1-α)*np.ones_like(s) \n",
" return u, ux, uxx"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### State transition function\n",
"Next period, reservoir level wealth will be equal to current level minus irrigation plus random rainfall:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"def transition(s, q, i, j=None, in_=None, e=None):\n",
" g = q\n",
" gx = np.ones_like(q)\n",
" gxx = np.zeros_like(q)\n",
" return g, gx, gxx"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Model structure\n",
"# TODO: CORREGIR ESTA ECUACION\n",
"\n",
"The value of wealth $s$ satisfies the Bellman equation \n",
"\\begin{equation*}\n",
"V(s) = \\max_x\\left\\{\\left(\\frac{a_0}{1+a_1}\\right)x^{1+a1} + \\left(\\frac{b_0}{1+b_1}\\right)(s-x)^{1+b1}+ \\delta V(s-x+e) \\right\\}\n",
"\\end{equation*}\n",
"\n",
"To solve and simulate this model,use the CompEcon class `DPmodel`"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"firm = DPmodel(basis, reward, transition, bounds,q=q, \n",
" i=prices, x=['Production'],discount=δ )"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"A CONTINUOUS STATE, CONTINUOUS ACTION DYNAMIC MODEL.\n",
"\n",
"\t* Continuous states:\n",
"\t\t0 : lagged production --> 5 nodes in [0.00, 20.00]\n",
"\n",
"\t* Continuous actions\n",
"\t\t0 : Production\n",
"\n",
"\t* Discrete states\n",
"\t\t0 : p_low\n",
"\t\t1 : p_mid\n",
"\t\t2 : p_high"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"firm"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Solving the model"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Solving the growth model by collocation. "
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Solving infinite-horizon model collocation equation by Newton's method\n",
"iter change time \n",
"------------------------------\n"
]
},
{
"ename": "AttributeError",
"evalue": "'numpy.ndarray' object has no attribute 'toarray'",
"output_type": "error",
"traceback": [
"\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[1;31mAttributeError\u001b[0m Traceback (most recent call last)",
"\u001b[1;32mc:\\Users\\randa\\OneDrive\\Documents\\Python\\CompEcon\\notebooks\\dp\\12 Production Management Model.ipynb Cell 23\u001b[0m in \u001b[0;36m\u001b[1;34m()\u001b[0m\n\u001b[1;32m----> 1\u001b[0m S \u001b[39m=\u001b[39m firm\u001b[39m.\u001b[39;49msolve()\n\u001b[0;32m 2\u001b[0m S\u001b[39m.\u001b[39mhead()\n",
"File \u001b[1;32mc:\\ProgramData\\Anaconda3\\lib\\site-packages\\compecon\\dpmodel.py:458\u001b[0m, in \u001b[0;36mDPmodel.solve\u001b[1;34m(self, v, x, nr, **kwargs)\u001b[0m\n\u001b[0;32m 456\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m__solve_by_function_iteration()\n\u001b[0;32m 457\u001b[0m \u001b[39melif\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39moptions\u001b[39m.\u001b[39malgorithm \u001b[39m==\u001b[39m \u001b[39m'\u001b[39m\u001b[39mnewton\u001b[39m\u001b[39m'\u001b[39m:\n\u001b[1;32m--> 458\u001b[0m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49m__solve_by_Newton_method()\n\u001b[0;32m 459\u001b[0m \u001b[39melse\u001b[39;00m:\n\u001b[0;32m 460\u001b[0m \u001b[39mraise\u001b[39;00m \u001b[39mValueError\u001b[39;00m(\u001b[39m'\u001b[39m\u001b[39mUnknown solution algorithm\u001b[39m\u001b[39m'\u001b[39m)\n",
"File \u001b[1;32mc:\\ProgramData\\Anaconda3\\lib\\site-packages\\compecon\\dpmodel.py:829\u001b[0m, in \u001b[0;36mDPmodel.__solve_by_Newton_method\u001b[1;34m(self)\u001b[0m\n\u001b[0;32m 827\u001b[0m cold \u001b[39m=\u001b[39m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mValue\u001b[39m.\u001b[39mc\u001b[39m.\u001b[39mcopy()\u001b[39m.\u001b[39mflatten()\n\u001b[0;32m 828\u001b[0m \u001b[39m# print('\\ncold', cold)\u001b[39;00m\n\u001b[1;32m--> 829\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mValue_j[:], vc \u001b[39m=\u001b[39m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mvmax(s, x, \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mValue, \u001b[39mTrue\u001b[39;49;00m)\n\u001b[0;32m 830\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mmake_discrete_choice()\n\u001b[0;32m 831\u001b[0m step \u001b[39m=\u001b[39m \u001b[39m-\u001b[39m SOLVE(Phik \u001b[39m-\u001b[39m vc, Phik \u001b[39m@\u001b[39m cold \u001b[39m-\u001b[39m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mValue\u001b[39m.\u001b[39my\u001b[39m.\u001b[39mflatten())\n",
"File \u001b[1;32mc:\\ProgramData\\Anaconda3\\lib\\site-packages\\compecon\\dpmodel.py:896\u001b[0m, in \u001b[0;36mDPmodel.vmax\u001b[1;34m(self, s, x, Value, dVc)\u001b[0m\n\u001b[0;32m 894\u001b[0m snext \u001b[39m=\u001b[39m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mtransition(s[:, is_], x[i, j, :, is_], i , j, in_, ee[:, is_]) \u001b[39m#fixme need to know number of output arguments!!!\u001b[39;00m\n\u001b[0;32m 895\u001b[0m prob \u001b[39m=\u001b[39m w[k] \u001b[39m*\u001b[39m q[j, i, in_,]\n\u001b[1;32m--> 896\u001b[0m vc[is_, i, :, in_] \u001b[39m+\u001b[39m\u001b[39m=\u001b[39m prob \u001b[39m*\u001b[39m Value\u001b[39m.\u001b[39;49mPhi(snext)\u001b[39m.\u001b[39;49mtoarray()\u001b[39m.\u001b[39mreshape((is_\u001b[39m.\u001b[39msum(), ms), order\u001b[39m=\u001b[39m\u001b[39m'\u001b[39m\u001b[39mF\u001b[39m\u001b[39m'\u001b[39m) \u001b[39m#fixme I can't find the proper way to index this\u001b[39;00m\n\u001b[0;32m 898\u001b[0m vc \u001b[39m=\u001b[39m vc\u001b[39m.\u001b[39mreshape((ns\u001b[39m*\u001b[39mni,ms\u001b[39m*\u001b[39mni),order\u001b[39m=\u001b[39m\u001b[39m'\u001b[39m\u001b[39mF\u001b[39m\u001b[39m'\u001b[39m)\n\u001b[0;32m 899\u001b[0m \u001b[39melse\u001b[39;00m:\n",
"\u001b[1;31mAttributeError\u001b[0m: 'numpy.ndarray' object has no attribute 'toarray'"
]
}
],
"source": [
"S = firm.solve()\n",
"S.head()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"firm.Policy_j(firm.Policy.nodes,dropdim=True).shape"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"`DPmodel.solve` returns a pandas `DataFrame` with the following data:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We are also interested in the shadow price of wealth (the first derivative of the value function)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"S['shadow price'] = water_model.Value(S['Reservoir'],1)\n",
"S.head()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Plotting the results"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Optimal Policy"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fig1 = demo.figure('Optimal Irrigation Policy', 'Reservoir Level', 'Irrigation')\n",
"plt.plot(S['Irrigation'])\n",
"demo.annotate(sstar, xstar,f'$s^*$ = {sstar:.2f}\\n$x^*$ = {xstar:.2f}','bo', (10, -6),ms=10,fs=11)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Value Function"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fig2 = demo.figure('Value Function', 'Reservoir Level', 'Value')\n",
"plt.plot(S['value'])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Shadow Price Function"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fig3 = demo.figure('Shadow Price Function', 'Reservoir Level', 'Shadow Price')\n",
"plt.plot(S['shadow price'])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Chebychev Collocation Residual"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fig4 = demo.figure('Bellman Equation Residual', 'Reservoir Level', 'Residual')\n",
"plt.hlines(0,smin,smax,'k',linestyles='--')\n",
"plt.plot(S[['resid']])\n",
"plt.ticklabel_format(style='sci', axis='y', scilimits=(-1,1))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Simulating the model\n",
"\n",
"We simulate 21 periods of the model starting from $s=s_{\\min}$"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"T = 31\n",
"nrep = 100_000\n",
"data = water_model.simulate(T, np.tile(smin,nrep))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Simulated State and Policy Paths"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"subdata = data[data['_rep'].isin(range(3))]\n",
"opts = dict(spec='r*', offset=(0, -15), fs=11, ha='right')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fig6 = demo.figure('Simulated and Expected Reservoir Level','Year', 'Reservoir Level',[0, T + 0.5])\n",
"plt.plot(data[['time','Reservoir']].groupby('time').mean())\n",
"plt.plot(subdata.pivot('time','_rep','Reservoir'),lw=1)\n",
"demo.annotate(T, sstar, f'steady-state reservoir\\n = {sstar:.2f}', **opts)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"fig7 = demo.figure('Simulated and Expected Irrigation','Year', 'Irrigation',[0, T + 0.5])\n",
"plt.plot(data[['time','Irrigation']].groupby('time').mean())\n",
"plt.plot(subdata.pivot('time','_rep','Irrigation'),lw=1)\n",
"demo.annotate(T, xstar, f'steady-state irrigation\\n = {xstar:.2f}', **opts)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Ergodic Wealth Distribution"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"subdata = data[data['time']==T][['Reservoir','Irrigation']]\n",
"stats =pd.DataFrame({'Deterministic Steady-State': [sstar, xstar],\n",
" 'Ergodic Means': subdata.mean(),\n",
" 'Ergodic Standard Deviations': subdata.std()}).T\n",
"stats"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"fig8 = demo.figure('Ergodic Reservoir and Irrigation Distribution','Wealth','Probability')\n",
"sns.kdeplot(subdata['Reservoir'])\n",
"sns.kdeplot(subdata['Irrigation'])"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"#demo.savefig([fig1,fig2,fig3,fig4,fig5,fig6,fig7,fig8])"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3.9.7 ('base')",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.9.7"
},
"vscode": {
"interpreter": {
"hash": "ad2bdc8ecc057115af97d19610ffacc2b4e99fae6737bb82f5d7fb13d2f2c186"
}
}
},
"nbformat": 4,
"nbformat_minor": 4
}
|