<jupyter_start><jupyter_text>DEMDP01b <NAME> -- using cubic splines<jupyter_code>import numpy as np
from compecon import NLP, demo, BasisSpline
import matplotlib.pyplot as plt
price = 1.0  # price of biomass
kappa = 0.2  # clearcut-replant cost
smax  = 0.5  # stand carrying capacity
gamma = 0.1  # biomass growth parameter
delta = 0.9  # discount factor<jupyter_output><empty_output><jupyter_text>Code the growth function<jupyter_code>def h(s): return np.array(s + gamma*(smax - s))<jupyter_output><empty_output><jupyter_text>SOLUTION Code the approximant and the residual<jupyter_code>ns = 200
vhat = BasisSpline(ns,0,smax,k=3)
def vhat1(s): return price*s - kappa + delta * vhat(h(0))
def vhat0(s): return delta * vhat(h(s))
def resid(c,s=vhat.nodes):
    vhat.c = c
    return vhat(s) - np.maximum(vhat0(s), vhat1(s))<jupyter_output><empty_output><jupyter_text>Solve collocation equation<jupyter_code>cc = NLP(resid).broyden(vhat.c)<jupyter_output><empty_output><jupyter_text>Compute critical biomass<jupyter_code>scrit = NLP(lambda s: vhat0(s)-vhat1(s)).broyden(0.0)[0]<jupyter_output><empty_output><jupyter_text>ANALYSIS  Compute refined state grid<jupyter_code>ss = np.linspace(0,smax,1000)<jupyter_output><empty_output><jupyter_text>Plot Conditional Value Functions<jupyter_code>fig1 =demo.figure('Conditional Value Functions','Biomass','Value of Stand')
plt.plot(ss,vhat0(ss),label='Grow')
plt.plot(ss,vhat1(ss),label='Clear-Cut')
plt.legend()

vcrit = vhat(scrit)
ymin = plt.ylim()[0]
plt.vlines(scrit, ymin,vcrit,'grey',linestyles='--')
demo.annotate(scrit,ymin,'$s^*$',ms=10)
demo.bullet(scrit,vcrit)
print(f'Optimal Biomass Harvest Level = {scrit:.4f}')<jupyter_output>Optimal Biomass Harvest Level = 0.3067<jupyter_text>Plot Value Function Residual<jupyter_code>fig2 = demo.figure('Bellman Equation Residual', 'Biomass', 'Percent Residual')
plt.plot(ss, 100*resid(cc,ss) / vhat(ss))
plt.hlines(0,0,smax,linestyles='--')<jupyter_output><empty_output><jupyter_text>Compute ergodic mean annual harvest<jupyter_code>s = h(0)
for n in range(100):
    if s > scrit: break
    s = h(s)
    
print(f'Ergodic Mean Annual Harvest = {s/n:.4f} after {n+1} iterations') 
#demo.savefig([fig1,fig2])<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1490125673249551, 'fraction_numerical': 0.026032315978456014, 'mean_word_length': 7.443181818181818, 'pattern_counts': {'":': 0, '<': 37, '<?xml version=': 0, '>': 38, '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': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '69552', 'license': 'MIT', 'n_tokens_mistral': 962, 'n_tokens_neox': 904, 'n_words': 156, 'path': 'notebooks/dp/01b Timber Harvesting -- cubic spline.ipynb', 'repo_name': 'snowdj/CompEcon-python', 'size': '2239'}