<jupyter_start><jupyter_text>Public Renewable Resource Model**<NAME>, PhD**This demo is based on the original Matlab demo accompanying the  Computational Economics and Finance 2001 textbook by <NAME> and <NAME>.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 --upgradeLast updated: 2021-Oct-01  AboutWelfare maximizing social planner must decide how much of a renewable resource to harvest.- States    -     s       quantity of stock available- Actions    -     q       quantity of stock harvested- Parameters    - $\alpha$     growth function parameter    - $\beta$      growth function parameter    - $\gamma$    relative risk aversion    - $\kappa$     unit cost of harvest    - $\delta$     discount factor<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
from compecon import BasisChebyshev, DPmodel, DPoptions, qnwnorm, demo<jupyter_output><empty_output><jupyter_text>Model parameters<jupyter_code>α, β, γ, κ, δ = 4.0, 1.0, 0.5, 0.2, 0.9<jupyter_output><empty_output><jupyter_text>Steady-stateThe steady-state values for this model are<jupyter_code>sstar = (α**2 - 1/δ**2)/(2*β)     # steady-state stock
qstar = sstar - (δ*α-1)/(δ*β) # steady-state action

print('Steady-State')
print(f"\t{'Stock':12s} = {sstar:8.4f}")
print(f"\t{'Harvest':12s} = {qstar:8.4f}")<jupyter_output>Steady-State
	Stock        =   7.3827
	Harvest      =   4.4938<jupyter_text>Numeric results  State spaceThe state variable is s="Stock", which we restrict to $s\in[6, 9]$. Here, we represent it with a Chebyshev basis, with $n=8$ nodes.<jupyter_code>n, smin, smax = 8, 6, 9
basis = BasisChebyshev(n, smin, smax, labels=['Stock'])<jupyter_output><empty_output><jupyter_text>Action spaceThe choice variable q="Harvest" must be nonnegative.<jupyter_code>def bounds(s, i=None, j=None):
    return np.zeros_like(s), s[:]<jupyter_output><empty_output><jupyter_text>Reward functionThe reward function is the utility of harvesting $q$ units.<jupyter_code>def reward(s, q, i=None, j=None):
    u = (q**(1-γ))/(1-γ)-κ*q
    ux= q**(-γ)-κ
    uxx = -γ*q**(-γ-1)
    
    return u, ux, uxx<jupyter_output><empty_output><jupyter_text>State transition functionNext period, the stock will be equal that is $s' = \alpha (s-q) - 0.5\beta(s-q)^2$<jupyter_code>def transition(s, q, i=None, j=None, in_=None, e=None):
    sq = s-q
    g = α*sq - 0.5*β*sq**2
    gx = -α + β*sq
    gxx = -β*np.ones_like(s)
    return g, gx, gxx<jupyter_output><empty_output><jupyter_text>Model structureThe value of wealth $s$ satisfies the Bellman equation\begin{equation*}V(s) = \max_q\left\{\frac{q^{1-\gamma}}{1-\gamma} -\kappa q + \delta V(s')  \right\}\end{equation*}To solve and simulate this model,use the CompEcon class `DPmodel`<jupyter_code>model = DPmodel(basis, reward, transition, bounds,
                       x=['Harvest'],
                       discount=δ)<jupyter_output><empty_output><jupyter_text>Solving the modelSolving the growth model by collocation<jupyter_code>S = model.solve()<jupyter_output>Solving infinite-horizon model collocation equation by Newton's method
iter change       time    
------------------------------
   0       1.4e+01    0.0040
   1       1.9e+01    0.0120
   2       3.7e-01    0.0180
   3       1.1e-03    0.0219
   4       1.5e-08    0.0249
   5       1.7e-14    0.0259
Elapsed Time =    0.03 Seconds<jupyter_text>`DPmodel.solve` returns a pandas `DataFrame` with the following data:<jupyter_code>S.head()<jupyter_output><empty_output><jupyter_text>We are also interested in the shadow price of wealth (the first derivative of the value function).<jupyter_code>S['shadow price'] = model.Value(S['Stock'],1)<jupyter_output><empty_output><jupyter_text>Plotting the results  Optimal Policy<jupyter_code>fig1 = demo.figure('Optimal Harvest Policy','Stock','Harvest')
plt.plot(S['Harvest'])
demo.annotate(sstar, qstar,f'$s^*$ = {sstar:.2f}\n$q^*$ = {qstar:.2f}','bo', (8, -8),ms=10,fs=11)<jupyter_output><empty_output><jupyter_text>Value Function<jupyter_code>fig2 = demo.figure('Value Function', 'Stock', 'Value')
plt.plot(S.value)<jupyter_output><empty_output><jupyter_text>Shadow Price Function<jupyter_code>fig3 = demo.figure('Shadow Price Function', 'Stock', 'Shadow Price')
plt.plot(S['shadow price'])<jupyter_output><empty_output><jupyter_text>Residual<jupyter_code>fig4 = demo.figure('Bellman Equation Residua', 'Stock', 'Residual')
plt.hlines(0,smin,smax,'k',linestyles='--')
plt.plot(S['resid'])
plt.ticklabel_format(style='sci', axis='y', scilimits=(-1,1))<jupyter_output><empty_output><jupyter_text>Simulating the modelWe simulate 16 periods of the model starting from $s=s_{\min}$<jupyter_code>T = 16
data = model.simulate(T, smin)<jupyter_output><empty_output><jupyter_text>Simulated State and Policy Paths<jupyter_code>opts = dict(spec='r*', offset=(-5, -5), fs=11, ha='right')

fig5 = demo.figure('State and Policy Paths','Period', 'Stock/Harvest',[0, T + 0.5])
plt.plot(data[['Stock', 'Harvest']])
demo.annotate(T, sstar, f'steady-state stock\n = {sstar:.2f}', **opts)
demo.annotate(T, qstar, f'steady-state harvest\n = {qstar:.2f}', **opts)
#demo.savefig([fig1,fig2,fig3,fig4,fig5])<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1527065527065527, 'fraction_numerical': 0.030389363722697058, 'mean_word_length': 4.650214592274678, 'pattern_counts': {'":': 0, '<': 70, '<?xml version=': 0, '>': 70, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '176042', 'license': 'MIT', 'n_tokens_mistral': 2037, 'n_tokens_neox': 1908, 'n_words': 435, 'path': '_build/jupyter_execute/notebooks/dp/08 Public Renewable Resource Model.ipynb', 'repo_name': 'randall-romero/CompEcon-python', 'size': '5294'}