<jupyter_start><jupyter_text>Example: Water Management Miranda and Fackler, Applied Computational Economics and Finance, 2002,Section 7.6.5<jupyter_code>%matplotlib inline
from __future__ import division
import numpy as np
import itertools
import matplotlib.pyplot as plt
from quantecon.markov import DiscreteDP
M = 30
n = M + 1  # Number of states
m = n  # Number of actions

a1, b1 = 14, 0.8
a2, b2 = 10, 0.4
F = lambda s, x: a1 * x**b1
G = lambda s, x: a2 * (s - x)**b2

beta = 0.9
# Reward array
R = np.empty((n, m))
for s, x in itertools.product(range(n), range(m)):
    R[s, x] = F(s, x) + G(s, x) if x <= s else -np.inf
# Probability transition array
probs = np.array([0.1, 0.2, 0.4, 0.2, 0.1])
size_supp = len(probs)
Q = np.zeros((n, m, n))
for s, x in itertools.product(range(n), range(m)):
    if x <= s:
        for j in range(size_supp):
            Q[s, x, np.minimum(s-x+j, n-1)] += probs[j]
# Create a DiscreteDP
ddp = DiscreteDP(R, Q, beta)
# Solve the dynamic optimization problem (by policy iteration)
res = ddp.solve()
# Number of iterations
res.num_iter
# Optimal policy
res.sigma
# Optimal value function
res.v
# Simulate the controlled Markov chain for num_rep times
# and compute the average
init = 0
ts_length = 50+1
num_rep = 10**4
ave_path = np.zeros(ts_length)
for i in range(num_rep):
    path = res.mc.simulate(ts_length=ts_length, init=init)
    ave_path = (i/(i+1)) * ave_path + (1/(i+1)) * path
ave_path
# Stationary distribution of the Markov chain
stationary_dist = res.mc.stationary_distributions[0]
stationary_dist
# Plot sigma, v, ave_path, stationary_dist
hspace = 0.3
fig, ax = plt.subplots(2, 2, figsize=(12, 8+hspace))
fig.subplots_adjust(hspace=hspace)

ax[0, 0].plot(res.sigma, '*')
ax[0, 0].set_xlim(-1, 31)
ax[0, 0].set_ylim(-0.5, 5.5)
ax[0, 0].set_xlabel('Water Level')
ax[0, 0].set_ylabel('Irrigation')
ax[0, 0].set_title('Optimal Irrigation Policy')

ax[0, 1].plot(res.v)
ax[0, 1].set_xlabel('Water Level')
ax[0, 1].set_ylabel('Value')
ax[0, 1].set_title('Optimal Value Function')

ax[1, 0].plot(ave_path)
ax[1, 0].set_xlabel('Year')
ax[1, 0].set_ylabel('Water Level')
ax[1, 0].set_title('Average Optimal State Path')

ax[1, 1].bar(range(n), stationary_dist, align='center')
ax[1, 1].set_xlim(-1, n)
ax[1, 1].set_xlabel('Water Level')
ax[1, 1].set_ylabel('Probability')
ax[1, 1].set_title('Stationary Distribution')

plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16174636174636176, 'fraction_numerical': 0.04282744282744283, 'mean_word_length': 4.714964370546318, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 5, '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': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '436180', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 1034, 'n_tokens_neox': 971, 'n_words': 263, 'path': 'ddp_ex_MF_7_6_5_py.ipynb', 'repo_name': 'tyleransom/QuantEcon.notebooks', 'size': '2405'}