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.
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
WARNING This demo is not running. Problem with dpmodel.
TODO: Fix error in dpmodel.
Last updated: 2022-Oct-23
Profit maximizing entrepeneur must decide how much to produce, subject to production adjustment costs.
i market price (discrete)
s lagged production (continuous)
x current production
import numpy as np
import matplotlib.pyplot as plt
from compecon import BasisSpline, DPmodel, DPoptions, qnwlogn, BasisChebyshev
import seaborn as sns
import pandas as pd
α, β0, β1, pbar = 0.01, 0.8, 0.03, 1.0
σ, δ = 0.2, 0.9
μ = np.log(pbar) - σ**2 / 2
m = 3 #number of market price shocks
p, w = qnwlogn(m, μ, σ**2)
q = np.repeat(w,3).reshape(3,3).T
The state variable is s="lagged production", which we restrict to $s\in[0, 20]$.
Here, we represent it with a cubic spline basis, with $n=50$ nodes.
n, smin, smax = 5, 0.0, 20.0
basis = BasisChebyshev(n, smin, smax, labels=['lagged production'])
The discrete state is given by the price p
prices = ['p_low', 'p_mid', 'p_high']
The choice variable x="current production" must be nonnegative.
def bounds(s, i, j=None):
return np.zeros_like(s), np.inf + np.zeros_like(s)
The reward function is
def reward(s, q, i, j=None):
u = p[i]*q - (β0*q + 0.5*β1*q**2) - 0.5*α*((q-s)**2)
ux = p[i] - β0 - β1*q - α*(q-s)
uxx = (-β1-α)*np.ones_like(s)
return u, ux, uxx
Next period, reservoir level wealth will be equal to current level minus irrigation plus random rainfall:
def transition(s, q, i, j=None, in_=None, e=None):
g = q
gx = np.ones_like(q)
gxx = np.zeros_like(q)
return g, gx, gxx
The value of wealth $s$ satisfies the Bellman equation \begin{equation*} 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\} \end{equation*}
To solve and simulate this model,use the CompEcon class DPmodel
firm = DPmodel(basis, reward, transition, bounds,q=q,
i=prices, x=['Production'],discount=δ )
firm
A CONTINUOUS STATE, CONTINUOUS ACTION DYNAMIC MODEL. * Continuous states: 0 : lagged production --> 5 nodes in [0.00, 20.00] * Continuous actions 0 : Production * Discrete states 0 : p_low 1 : p_mid 2 : p_high
Solving the growth model by collocation.
S = firm.solve()
S.head()
Solving infinite-horizon model collocation equation by Newton's method iter change time ------------------------------
--------------------------------------------------------------------------- AttributeError Traceback (most recent call last) c:\Users\randa\OneDrive\Documents\Python\CompEcon\notebooks\dp\12 Production Management Model.ipynb Cell 23 in <cell line: 1>() ----> <a href='vscode-notebook-cell:/c%3A/Users/randa/OneDrive/Documents/Python/CompEcon/notebooks/dp/12%20Production%20Management%20Model.ipynb#X32sZmlsZQ%3D%3D?line=0'>1</a> S = firm.solve() <a href='vscode-notebook-cell:/c%3A/Users/randa/OneDrive/Documents/Python/CompEcon/notebooks/dp/12%20Production%20Management%20Model.ipynb#X32sZmlsZQ%3D%3D?line=1'>2</a> S.head() File c:\ProgramData\Anaconda3\lib\site-packages\compecon\dpmodel.py:458, in DPmodel.solve(self, v, x, nr, **kwargs) 456 self.__solve_by_function_iteration() 457 elif self.options.algorithm == 'newton': --> 458 self.__solve_by_Newton_method() 459 else: 460 raise ValueError('Unknown solution algorithm') File c:\ProgramData\Anaconda3\lib\site-packages\compecon\dpmodel.py:829, in DPmodel.__solve_by_Newton_method(self) 827 cold = self.Value.c.copy().flatten() 828 # print('\ncold', cold) --> 829 self.Value_j[:], vc = self.vmax(s, x, self.Value, True) 830 self.make_discrete_choice() 831 step = - SOLVE(Phik - vc, Phik @ cold - self.Value.y.flatten()) File c:\ProgramData\Anaconda3\lib\site-packages\compecon\dpmodel.py:896, in DPmodel.vmax(self, s, x, Value, dVc) 894 snext = self.transition(s[:, is_], x[i, j, :, is_], i , j, in_, ee[:, is_]) #fixme need to know number of output arguments!!! 895 prob = w[k] * q[j, i, in_,] --> 896 vc[is_, i, :, in_] += prob * Value.Phi(snext).toarray().reshape((is_.sum(), ms), order='F') #fixme I can't find the proper way to index this 898 vc = vc.reshape((ns*ni,ms*ni),order='F') 899 else: AttributeError: 'numpy.ndarray' object has no attribute 'toarray'
firm.Policy_j(firm.Policy.nodes,dropdim=True).shape
DPmodel.solve
returns a pandas DataFrame
with the following data:
We are also interested in the shadow price of wealth (the first derivative of the value function).
S['shadow price'] = water_model.Value(S['Reservoir'],1)
S.head()
fig1 = demo.figure('Optimal Irrigation Policy', 'Reservoir Level', 'Irrigation')
plt.plot(S['Irrigation'])
demo.annotate(sstar, xstar,f'$s^*$ = {sstar:.2f}\n$x^*$ = {xstar:.2f}','bo', (10, -6),ms=10,fs=11)
fig2 = demo.figure('Value Function', 'Reservoir Level', 'Value')
plt.plot(S['value'])
fig3 = demo.figure('Shadow Price Function', 'Reservoir Level', 'Shadow Price')
plt.plot(S['shadow price'])
fig4 = demo.figure('Bellman Equation Residual', 'Reservoir Level', 'Residual')
plt.hlines(0,smin,smax,'k',linestyles='--')
plt.plot(S[['resid']])
plt.ticklabel_format(style='sci', axis='y', scilimits=(-1,1))
We simulate 21 periods of the model starting from $s=s_{\min}$
T = 31
nrep = 100_000
data = water_model.simulate(T, np.tile(smin,nrep))
subdata = data[data['_rep'].isin(range(3))]
opts = dict(spec='r*', offset=(0, -15), fs=11, ha='right')
fig6 = demo.figure('Simulated and Expected Reservoir Level','Year', 'Reservoir Level',[0, T + 0.5])
plt.plot(data[['time','Reservoir']].groupby('time').mean())
plt.plot(subdata.pivot('time','_rep','Reservoir'),lw=1)
demo.annotate(T, sstar, f'steady-state reservoir\n = {sstar:.2f}', **opts)
fig7 = demo.figure('Simulated and Expected Irrigation','Year', 'Irrigation',[0, T + 0.5])
plt.plot(data[['time','Irrigation']].groupby('time').mean())
plt.plot(subdata.pivot('time','_rep','Irrigation'),lw=1)
demo.annotate(T, xstar, f'steady-state irrigation\n = {xstar:.2f}', **opts)
subdata = data[data['time']==T][['Reservoir','Irrigation']]
stats =pd.DataFrame({'Deterministic Steady-State': [sstar, xstar],
'Ergodic Means': subdata.mean(),
'Ergodic Standard Deviations': subdata.std()}).T
stats
fig8 = demo.figure('Ergodic Reservoir and Irrigation Distribution','Wealth','Probability')
sns.kdeplot(subdata['Reservoir'])
sns.kdeplot(subdata['Irrigation'])
#demo.savefig([fig1,fig2,fig3,fig4,fig5,fig6,fig7,fig8])