#!/usr/bin/env python
# coding: utf-8

# # Example: Water Management

# Miranda and Fackler, <i>Applied Computational Economics and Finance</i>, 2002,
# Section 7.6.5

# In[1]:


get_ipython().run_line_magic('matplotlib', 'inline')


# In[2]:


from __future__ import division
import numpy as np
import itertools
import matplotlib.pyplot as plt
from quantecon.markov import DiscreteDP


# In[3]:


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


# In[4]:


# 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


# In[5]:


# 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]


# In[6]:


# Create a DiscreteDP
ddp = DiscreteDP(R, Q, beta)


# In[7]:


# Solve the dynamic optimization problem (by policy iteration)
res = ddp.solve()


# In[8]:


# Number of iterations
res.num_iter


# In[9]:


# Optimal policy
res.sigma


# In[10]:


# Optimal value function
res.v


# In[11]:


# 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


# In[12]:


ave_path


# In[13]:


# Stationary distribution of the Markov chain
stationary_dist = res.mc.stationary_distributions[0]


# In[14]:


stationary_dist


# In[15]:


# 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()


# In[ ]: