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: demdp07.m
Running this file requires the Python version of CompEcon. This can be installed with pip by running
!pip install compecon --upgrade
Last updated: 2022-Oct-23
Welfare maximizing social planner must decide how much society should consume and invest. Unlike the deterministic model, this model allows arbitrary constant relative risk aversion, capital depreciation, and stochastic production shocks. It lacks a known closed-form solution.
import numpy as np
import matplotlib.pyplot as plt
from compecon import BasisChebyshev, DPmodel, DPoptions, qnwlogn
import seaborn as sns
import pandas as pd
α, β, γ, σ, δ = 0.2, 0.5, 0.9, 0.1, 0.9
The deterministic steady-state values for this model are
kstar = ((1 - δ*γ)/(δ*β))**(1/(β-1)) # steady-state capital investment
sstar = γ*kstar + kstar**β # steady-state wealth
The state variable is s="Wealth", which we restrict to $s\in[5, 10]$.
Here, we represent it with a Chebyshev basis, with $n=10$ nodes.
n, smin, smax = 10, 5, 10
basis = BasisChebyshev(n, smin, smax, labels=['Wealth'])
m = 5
e, w = qnwlogn(m, -σ**2/2,σ**2)
The choice variable k="Investment" must be nonnegative.
def bounds(s, i=None, j=None):
return np.zeros_like(s), 0.99*s
The reward function is the utility of consumption=$s-k$.
def reward(s, k, i=None, j=None):
sk = s - k
u = sk**(1-α) / (1-α)
ux= -sk **(-α)
uxx = -α * sk**(-1-α)
return u, ux, uxx
Next period, wealth will be equal to production from available initial capital $k$, that is $s' = k^\beta$
def transition(s, k, i=None, j=None, in_=None, e=None):
g = γ * k + e * k**β
gx = γ + β*e * k**(β-1)
gxx = β*(β-1)*e * k**(β-2)
return g, gx, gxx
The value of wealth $s$ satisfies the Bellman equation \begin{equation*} V(s) = \max_k\left\{\log(s-k) + \delta V(k^\beta) \right\} \end{equation*}
To solve and simulate this model,use the CompEcon class DPmodel
growth_model = DPmodel(basis, reward, transition, bounds,e=e,w=w,
x=['Investment'],
discount=δ )
Solving the growth model by collocation.
S = growth_model.solve()
Solving infinite-horizon model collocation equation by Newton's method iter change time ------------------------------ 0 8.2e+00 0.0099 1 9.6e+00 0.0932 2 2.2e+00 0.1444 3 3.8e-02 0.1826 4 7.0e-06 0.2109 5 7.1e-15 0.2169 Elapsed Time = 0.22 Seconds
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).
To analyze the dynamics of the model, it also helps to compute the optimal change of wealth.
S['shadow price'] = growth_model.Value(S['Wealth'],1)
S['D.Wealth'] = transition(S['Wealth'], S['Investment'],e=1)[0] - S['Wealth']
S.head()
Wealth | value | resid | Investment | shadow price | D.Wealth | |
---|---|---|---|---|---|---|
Wealth | ||||||
5.000000 | 5.000000 | 17.793592 | 4.324104e-09 | 3.997670 | 0.999535 | 0.597321 |
5.050505 | 5.050505 | 17.843995 | -1.837417e-09 | 4.032514 | 0.996440 | 0.586869 |
5.101010 | 5.101010 | 17.894243 | -4.078490e-09 | 4.067300 | 0.993391 | 0.576314 |
5.151515 | 5.151515 | 17.944339 | -3.954334e-09 | 4.102028 | 0.990386 | 0.565657 |
5.202020 | 5.202020 | 17.994283 | -2.583509e-09 | 4.136700 | 0.987425 | 0.554898 |
The DPmodel.lqapprox solves the linear-quadratic approximation, in this case arround the steady-state. It returns a LQmodel which works similar to the DPmodel object. We also compute the shadow price and the approximation error to compare these results to the collocation results.
growth_lq = growth_model.lqapprox(sstar, kstar)
L = growth_lq.solution(basis.nodes)
L['shadow price'] = L['value_Wealth']
L['D.Wealth'] = L['Wealth_next']- L['Wealth']
L.head()
Wealth | Investment | value | value_Wealth | Wealth_next | shadow price | D.Wealth | |
---|---|---|---|---|---|---|---|
Wealth | |||||||
5.030779 | 5.030779 | 3.461913 | 17.923286 | 0.911808 | 5.030781 | 0.911808 | 1.397911e-06 |
5.272484 | 5.272484 | 3.679447 | 18.143387 | 0.909432 | 5.272485 | 0.909432 | 1.256542e-06 |
5.732233 | 5.732233 | 4.093221 | 18.560459 | 0.904912 | 5.732234 | 0.904912 | 9.876423e-07 |
6.365024 | 6.365024 | 4.662732 | 19.131111 | 0.898692 | 6.365024 | 0.898692 | 6.175332e-07 |
7.108914 | 7.108914 | 5.332233 | 19.796920 | 0.891380 | 7.108914 | 0.891380 | 1.824439e-07 |
fig1, ax = plt.subplots()
ax.set(title='Optimal Investment Policy', xlabel='Wealth', ylabel='Investment')
ax.plot(S.Investment, label='Chebychev Collocation')
ax.plot(L.Investment, label='L-Q Approximation')
ax.plot(sstar, kstar,'wo')
ax.annotate(f"$s^*$ = {sstar:.2f}\n$k^*$ = {kstar:.2f}",[sstar, kstar],va='top')
ax.legend(loc= 'upper left');
fig2, ax = plt.subplots()
ax.set(title='Value Function', xlabel='Wealth', ylabel='Value')
ax.plot(S.value, label='Chebychev Collocation')
ax.plot(L.value, label='L-Q Approximation')
ax.legend(loc= 'upper left');
fig3, ax = plt.subplots()
ax.set(title='Shadow Price Function', xlabel='Wealth', ylabel='Shadow Price')
ax.plot(S['shadow price'], label='Chebychev Collocation')
ax.plot(L['shadow price'], label='L-Q Approximation')
ax.legend(loc= 'upper right');
fig4, ax = plt.subplots()
ax.set(title='Bellman Equation Residual', xlabel='Wealth', ylabel='Residual')
ax.axhline(0, color='w')
ax.plot(S[['resid']])
ax.scatter(basis.nodes[0], np.zeros_like(basis.nodes[0]), color='C1')
ax.ticklabel_format(style='sci', axis='y', scilimits=(-1,1))
Notice how the steady-state is stable in the Chebyshev collocation solution, but unstable in the linear-quadratic approximation. In particular, simulated paths of wealth in the L-Q approximation will converge to zero, unless the initial states is within a small neighborhood of the steady-state.
fig5, ax = plt.subplots()
ax.set(title='Wealth dynamics', xlabel='Wealth', ylabel='Wealth change')
ax.plot(S['D.Wealth'], label='Chebychev Collocation')
ax.plot(L['D.Wealth'], label='L-Q Approximation')
ax.axhline(0, linestyle=':')
ax.plot(sstar, 0, 'wo')
ax.annotate(f"$s^*$ = {sstar:.2f}\n$\Delta s^*$ = 0.00",[sstar, 0], va='bottom')
plt.legend(loc= 'lower left');
We simulate 21 periods of the model starting from $s=s_{\min}$
T = 21
nrep = 50_000
data = growth_model.simulate(T, np.tile(smin,nrep))
subdata = data.query('_rep < 100')
fig6, ax = plt.subplots()
ax.set(title='Simulated and Expected Wealth',
xlabel='Period',
ylabel='Wealth',
xlim=[0, T + 0.5])
ax.plot(data[['time','Wealth']].groupby('time').mean())
ax.plot(subdata.pivot('time','_rep','Wealth'), lw=1, color='C0', alpha=0.2)
ax.annotate(f'steady-state wealth\n = {sstar:.2f}', [T, sstar-0.03], color='C0', va='top', ha='right')
Text(21, 7.386897506925211, 'steady-state wealth\n = 7.42')
fig7, ax = plt.subplots()
ax.set(title='Simulated and Expected Investment',
xlabel='Period',
ylabel='Investment',
xlim=[0, T + 0.5],
xticks=range(0,T+1,5))
ax.plot(data[['time','Investment']].groupby('time').mean(), color='C1')
ax.plot(subdata.pivot('time','_rep','Investment'),lw=1, color='C1', alpha=0.2)
ax.annotate(f'steady-state investment\n = {sstar:.2f}', [T, kstar-0.03], color='C1', va='top', ha='right')
Text(21, 5.579418282548479, 'steady-state investment\n = 7.42')
data.query("time == @T")
time | _rep | Wealth | Investment | |
---|---|---|---|---|
1050000 | 21 | 0 | 7.326797 | 5.551655 |
1050001 | 21 | 1 | 6.874412 | 5.256752 |
1050002 | 21 | 2 | 7.258543 | 5.507353 |
1050003 | 21 | 3 | 7.396451 | 5.596796 |
1050004 | 21 | 4 | 7.534192 | 5.685864 |
... | ... | ... | ... | ... |
1099995 | 21 | 49995 | 7.495991 | 5.661188 |
1099996 | 21 | 49996 | 7.283651 | 5.523658 |
1099997 | 21 | 49997 | 7.849644 | 5.888873 |
1099998 | 21 | 49998 | 8.163080 | 6.089300 |
1099999 | 21 | 49999 | 7.691887 | 5.787515 |
50000 rows × 4 columns
subdata = data.query("time == @T")[['Wealth','Investment']]
stats =pd.DataFrame({'Deterministic Steady-State': [sstar, kstar],
'Ergodic Means': subdata.mean(),
'Ergodic Standard Deviations': subdata.std()}).T
stats
Wealth | Investment | |
---|---|---|
Deterministic Steady-State | 7.416898 | 5.609418 |
Ergodic Means | 7.411443 | 5.605683 |
Ergodic Standard Deviations | 0.342352 | 0.221515 |
fig8, ax = plt.subplots()
ax.set(title='Ergodic Wealth and Investment Distribution',
xlabel='Wealth',
ylabel='Probability',
xlim=[4.5, 9.5])
sns.kdeplot(subdata['Wealth'], ax=ax, label='Wealth')
sns.kdeplot(subdata['Investment'], ax=ax, label='Investment')
ax.legend();