<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 DiscreteDP Example
Step1:  Setup
Step2:  Continuous-state benchmark
Let us compute the value function of the continuous-state version
as described in equations (2.22) and (2.23) in Section 2.3.
Step3:  The optimal stopping boundary $\gamma$ for the contiuous-state version, given by (2.23)
Step4:  The value function for the continuous-state version, given by (2.24)
Step5:  Solving the problem with DiscreteDP
Construct a DiscreteDP instance for the disrete-state version
Step6:  Let us solve the decision problem by
(0) value iteration,
(1) value iteration with span-based termination
(equivalent to modified policy iteration with step $k = 0$),
(2) policy iteration,
(3) modified policy iteration.
Following Rust (1996), we set
Step7:  The numbers of iterations
Step8:  Policy iteration gives the optimal policy
Step9:  Takes action 1 ("replace") if and only if $s \geq \bar{\gamma}$, where $\bar{\gamma}$ is equal to
Step10:  Check that the other methods gave the correct answer
Step11:  The deviations of the returned value function from the continuous-state benchmark
Step12:  In the following we try to reproduce Table 14.1 in Rust (1996), p.660,
although the precise definitions and procedures there are not very clear.
The maximum absolute differences of $v$ from that by policy iteration
Step13:  Compute $\lVert v - T(v)\rVert$
Step14:  Next we compute $\overline{b} - \underline{b}$
for the three methods other than policy iteration, where
$I$ is the number of iterations required to fulfill the termination condition, and
$$
\begin{aligned}
\underline{b} &= \frac{\beta}{1-\beta} \min\left[T(v^{I-1}) - v^{I-1}\right], \\
\overline{b}  &= \frac{\beta}{1-\beta} \max\left[T(v^{I-1}) - v^{I-1}\right].
\end{aligned}
$$
Step15:  For policy iteration, while it does not seem really relevant,
we compute $\overline{b} - \underline{b}$ with the returned value of $v$
in place of $v^{I-1}$
Step16:  Last, time each algorithm
Step17:  Notes
It appears that our value iteration with span-based termination is different in some details
from the corresponding algorithm (successive approximation with error bounds) in Rust.
In returing the value function, our algorithm returns
$T(v^{I-1}) + (\overline{b} + \underline{b})/2$,
while Rust's seems to return $v^{I-1} + (\overline{b} + \underline{b})/2$.
In fact
Step18:  $\lVert v - v_{\mathrm{pi}}\rVert$
Step19:  $\lVert v - T(v)\rVert$
Step20:  Compare the Table in Rust.
Convergence of trajectories
Let us plot the convergence of $v^i$ for the four algorithms;
see also Figure 14.2 in Rust.
Value iteration
Step21:  Value iteration with span-based termination
Step22:  Policy iteration
Step23:  Modified policy iteration
Step24:  Increasing the discount factor
Let us consider the case with a discount factor closer to $1$, $\beta = 0.9999$.
Step25:  The numbers of iterations
Step26:  Policy iteration gives the optimal policy
Step27:  Takes action 1 ("replace") if and only if $s \geq \bar{\gamma}$, where $\bar{\gamma}$ is equal to
Step28:  Check that the other methods gave the correct answer
Step29:  $\lVert v - v_{\mathrm{pi}}\rVert$
Step30:  $\lVert v - T(v)\rVert$
Step31:  $\overline{b} - \underline{b}$
Step32:  For policy iteration
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
import numpy as np
import itertools
import scipy.optimize
import matplotlib.pyplot as plt
import pandas as pd
from quantecon.markov import DiscreteDP
# matplotlib settings
plt.rcParams['axes.autolimit_mode'] = 'round_numbers'
plt.rcParams['axes.xmargin'] = 0
plt.rcParams['axes.ymargin'] = 0
plt.rcParams['patch.force_edgecolor'] = True
from cycler import cycler
plt.rcParams['axes.prop_cycle'] = cycler(color='bgrcmyk')
Explanation: DiscreteDP Example: Automobile Replacement
Daisuke Oyama
Faculty of Economics, University of Tokyo
We study the finite-state version of the automobile replacement problem as considered in
Rust (1996, Section 4.2.2).
J. Rust, "Numerical Dynamic Programming in Economics",
  <i>Handbook of Computational Economics</i>, Volume 1, 619-729, 1996.
End of explanation
lambd = 0.5        # Exponential distribution parameter
c = 200            # (Constant) marginal cost of maintainance
net_price = 10**5  # Replacement cost
n = 100  # Number of states; s = 0, ..., n-1: level of utilization of the asset
m = 2    # Number of actions; 0: keep, 1: replace
# Reward array
R = np.empty((n, m))
R[:, 0] = -c * np.arange(n)   # Costs for maintainance
R[:, 1] = -net_price - c * 0  # Costs for replacement
# Transition probability array
# For each state s, s' distributes over
# s, s+1, ..., min{s+supp_size-1, n-1} if a = 0
# 0, 1, ..., supp_size-1 if a = 1
# according to the (discretized and truncated) exponential distribution
# with parameter lambd
supp_size = 12
probs = np.empty(supp_size)
probs[0] = 1 - np.exp(-lambd * 0.5)
for j in range(1, supp_size-1):
    probs[j] = np.exp(-lambd * (j - 0.5)) - np.exp(-lambd * (j + 0.5))
probs[supp_size-1] = 1 - np.sum(probs[:-1])
Q = np.zeros((n, m, n))
# a = 0
for i in range(n-supp_size):
    Q[i, 0, i:i+supp_size] = probs
for k in range(supp_size):
    Q[n-supp_size+k, 0, n-supp_size+k:] = probs[:supp_size-k]/probs[:supp_size-k].sum()
# a = 1
for i in range(n):
    Q[i, 1, :supp_size] = probs
# Discount factor
beta = 0.95
Explanation: Setup
End of explanation
def f(x, s):
    return (c/(1-beta)) * \
        ((x-s) - (beta/(lambd*(1-beta))) * (1 - np.exp(-lambd*(1-beta)*(x-s))))
Explanation: Continuous-state benchmark
Let us compute the value function of the continuous-state version
as described in equations (2.22) and (2.23) in Section 2.3.
End of explanation
gamma = scipy.optimize.brentq(lambda x: f(x, 0) - net_price, 0, 100)
print(gamma)
Explanation: The optimal stopping boundary $\gamma$ for the contiuous-state version, given by (2.23):
End of explanation
def value_func_cont_time(s):
    return -c*gamma/(1-beta) + (s < gamma) * f(gamma, s)
v_cont = value_func_cont_time(np.arange(n))
Explanation: The value function for the continuous-state version, given by (2.24):
End of explanation
ddp = DiscreteDP(R, Q, beta)
Explanation: Solving the problem with DiscreteDP
Construct a DiscreteDP instance for the disrete-state version:
End of explanation
v_init = np.zeros(ddp.num_states)
epsilon = 1164
methods = ['vi', 'mpi', 'pi', 'mpi']
labels = ['Value iteration', 'Value iteration with span-based termination',
          'Policy iteration', 'Modified policy iteration']
results = {}
for i in range(4):
    k = 20 if labels[i] == 'Modified policy iteration' else 0
    results[labels[i]] = \
        ddp.solve(method=methods[i], v_init=v_init, epsilon=epsilon, k=k)
columns = [
    'Iterations', 'Time (second)', r'$\lVert v - v_{\mathrm{pi}} \rVert$',
    r'$\overline{b} - \underline{b}$', r'$\lVert v - T(v)\rVert$'
]
df = pd.DataFrame(index=labels, columns=columns)
Explanation: Let us solve the decision problem by
(0) value iteration,
(1) value iteration with span-based termination
(equivalent to modified policy iteration with step $k = 0$),
(2) policy iteration,
(3) modified policy iteration.
Following Rust (1996), we set:
$\varepsilon = 1164$ (for value iteration and modified policy iteration),
$v^0 \equiv 0$,
the number of iteration for iterative policy evaluation $k = 20$.
End of explanation
for label in labels:
    print(results[label].num_iter, '\t' + '(' + label + ')')
    df[columns[0]].loc[label] = results[label].num_iter
Explanation: The numbers of iterations:
End of explanation
print(results['Policy iteration'].sigma)
Explanation: Policy iteration gives the optimal policy:
End of explanation
(1-results['Policy iteration'].sigma).sum()
Explanation: Takes action 1 ("replace") if and only if $s \geq \bar{\gamma}$, where $\bar{\gamma}$ is equal to:
End of explanation
for result in results.values():
    if result != results['Policy iteration']:
        print(np.array_equal(result.sigma, results['Policy iteration'].sigma))
Explanation: Check that the other methods gave the correct answer:
End of explanation
diffs_cont = {}
for label in labels:
    diffs_cont[label] = np.abs(results[label].v - v_cont).max()
    print(diffs_cont[label], '\t' + '(' + label + ')')
label = 'Policy iteration'
fig, ax = plt.subplots(figsize=(8,5))
ax.plot(-v_cont, label='Continuous-state')
ax.plot(-results[label].v, label=label)
ax.set_title('Comparison of discrete vs. continuous value functions')
ax.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
ax.set_xlabel('State')
ax.set_ylabel(r'Value $\times\ (-1)$')
plt.legend(loc=4)
plt.show()
Explanation: The deviations of the returned value function from the continuous-state benchmark:
End of explanation
for label in labels:
    diff_pi = \
        np.abs(results[label].v - results['Policy iteration'].v).max()
    print(diff_pi, '\t' + '(' + label + ')')
    df[columns[2]].loc[label] = diff_pi
Explanation: In the following we try to reproduce Table 14.1 in Rust (1996), p.660,
although the precise definitions and procedures there are not very clear.
The maximum absolute differences of $v$ from that by policy iteration:
End of explanation
for label in labels:
    v = results[label].v
    diff_max = \
        np.abs(v - ddp.bellman_operator(v)).max()
    print(diff_max, '\t' + '(' + label + ')')
    df[columns[4]].loc[label] = diff_max
Explanation: Compute $\lVert v - T(v)\rVert$:
End of explanation
for i in range(4):
    if labels[i] != 'Policy iteration':
        k = 20 if labels[i] == 'Modified policy iteration' else 0
        res = ddp.solve(method=methods[i], v_init=v_init, k=k,
                        max_iter=results[labels[i]].num_iter-1)
        diff = ddp.bellman_operator(res.v) - res.v
        diff_span = (diff.max() - diff.min()) * ddp.beta / (1 - ddp.beta)
        print(diff_span, '\t' + '(' + labels[i] + ')')
        df[columns[3]].loc[labels[i]] = diff_span
Explanation: Next we compute $\overline{b} - \underline{b}$
for the three methods other than policy iteration, where
$I$ is the number of iterations required to fulfill the termination condition, and
$$
\begin{aligned}
\underline{b} &= \frac{\beta}{1-\beta} \min\left[T(v^{I-1}) - v^{I-1}\right], \\
\overline{b}  &= \frac{\beta}{1-\beta} \max\left[T(v^{I-1}) - v^{I-1}\right].
\end{aligned}
$$
End of explanation
label = 'Policy iteration'
v = results[label].v
diff = ddp.bellman_operator(v) - v
diff_span = (diff.max() - diff.min()) * ddp.beta / (1 - ddp.beta)
print(diff_span, '\t' + '(' + label + ')')
df[columns[3]].loc[label] = diff_span
Explanation: For policy iteration, while it does not seem really relevant,
we compute $\overline{b} - \underline{b}$ with the returned value of $v$
in place of $v^{I-1}$:
End of explanation
for i in range(4):
    k = 20 if labels[i] == 'Modified policy iteration' else 0
    print(labels[i])
    t = %timeit -o ddp.solve(method=methods[i], v_init=v_init, epsilon=epsilon, k=k)
    df[columns[1]].loc[labels[i]] = t.best
df
Explanation: Last, time each algorithm:
End of explanation
i = 1
k = 0
res = ddp.solve(method=methods[i], v_init=v_init, k=k,
                max_iter=results[labels[i]].num_iter-1)
diff = ddp.bellman_operator(res.v) - res.v
v = res.v + (diff.max() + diff.min()) * ddp.beta / (1 - ddp.beta) / 2
Explanation: Notes
It appears that our value iteration with span-based termination is different in some details
from the corresponding algorithm (successive approximation with error bounds) in Rust.
In returing the value function, our algorithm returns
$T(v^{I-1}) + (\overline{b} + \underline{b})/2$,
while Rust's seems to return $v^{I-1} + (\overline{b} + \underline{b})/2$.
In fact:
End of explanation
np.abs(v - results['Policy iteration'].v).max()
Explanation: $\lVert v - v_{\mathrm{pi}}\rVert$:
End of explanation
np.abs(v - ddp.bellman_operator(v)).max()
Explanation: $\lVert v - T(v)\rVert$:
End of explanation
label = 'Value iteration'
iters = [2, 20, 40, 80]
v = np.zeros(ddp.num_states)
fig, ax = plt.subplots(figsize=(8,5))
for i in range(iters[-1]):
    v = ddp.bellman_operator(v)
    if i+1 in iters:
        ax.plot(-v, label='Iteration {0}'.format(i+1))
ax.plot(-results['Policy iteration'].v, label='Fixed Point')
ax.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
ax.set_ylim(0, 2.4e5)
ax.set_yticks([0.4e5 * i for i in range(7)])
ax.set_title(label)
ax.set_xlabel('State')
ax.set_ylabel(r'Value $\times\ (-1)$')
plt.legend(loc=(0.7, 0.2))
plt.show()
Explanation: Compare the Table in Rust.
Convergence of trajectories
Let us plot the convergence of $v^i$ for the four algorithms;
see also Figure 14.2 in Rust.
Value iteration
End of explanation
label = 'Value iteration with span-based termination'
iters = [1, 10, 15, 20]
v = np.zeros(ddp.num_states)
fig, ax = plt.subplots(figsize=(8,5))
for i in range(iters[-1]):
    u = ddp.bellman_operator(v)
    if i+1 in iters:
        diff = u - v
        w = u + ((diff.max() + diff.min()) / 2) * ddp.beta / (1 - ddp.beta)
        ax.plot(-w, label='Iteration {0}'.format(i+1))
    v = u
ax.plot(-results['Policy iteration'].v, label='Fixed Point')
ax.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
ax.set_ylim(1.0e5, 2.4e5)
ax.set_yticks([1.0e5+0.2e5 * i for i in range(8)])
ax.set_title(label)
ax.set_xlabel('State')
ax.set_ylabel(r'Value $\times\ (-1)$')
plt.legend(loc=(0.7, 0.2))
plt.show()
Explanation: Value iteration with span-based termination
End of explanation
label = 'Policy iteration'
iters = [1, 2, 3]
v_init = np.zeros(ddp.num_states)
fig, ax = plt.subplots(figsize=(8,5))
sigma = ddp.compute_greedy(v_init)
for i in range(iters[-1]):
    # Policy evaluation
    v_sigma = ddp.evaluate_policy(sigma)
    if i+1 in iters:
        ax.plot(-v_sigma, label='Iteration {0}'.format(i+1))
    # Policy improvement
    new_sigma = ddp.compute_greedy(v_sigma)
    sigma = new_sigma
ax.plot(-results['Policy iteration'].v, label='Fixed Point')
ax.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
ax.set_ylim(1e5, 4.2e5)
ax.set_yticks([1e5 + 0.4e5 * i for i in range(9)])
ax.set_title(label)
ax.set_xlabel('State')
ax.set_ylabel(r'Value $\times\ (-1)$')
plt.legend(loc=4)
plt.show()
Explanation: Policy iteration
End of explanation
label = 'Modified policy iteration'
iters = [1, 2, 3, 4]
v = np.zeros(ddp.num_states)
k = 20 #- 1
fig, ax = plt.subplots(figsize=(8,5))
for i in range(iters[-1]):
    # Policy improvement
    sigma = ddp.compute_greedy(v)
    u = ddp.bellman_operator(v)
    if i == results[label].num_iter-1:
        diff = u - v
        break
    # Partial policy evaluation with k=20 iterations
    for j in range(k):
        u = ddp.T_sigma(sigma)(u)
    v = u
    if i+1 in iters:
        ax.plot(-v, label='Iteration {0}'.format(i+1))
ax.plot(-results['Policy iteration'].v, label='Fixed Point')
ax.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
ax.set_ylim(0, 2.8e5)
ax.set_yticks([0.4e5 * i for i in range(8)])
ax.set_title(label)
ax.set_xlabel('State')
ax.set_ylabel(r'Value $\times\ (-1)$')
plt.legend(loc=4)
plt.show()
Explanation: Modified policy iteration
End of explanation
ddp.beta = 0.9999
v_init = np.zeros(ddp.num_states)
epsilon = 1164
ddp.max_iter = 10**5 * 2
results_9999 = {}
for i in range(4):
    k = 20 if labels[i] == 'Modified policy iteration' else 0
    results_9999[labels[i]] = \
        ddp.solve(method=methods[i], v_init=v_init, epsilon=epsilon, k=k)
df_9999 = pd.DataFrame(index=labels, columns=columns)
Explanation: Increasing the discount factor
Let us consider the case with a discount factor closer to $1$, $\beta = 0.9999$.
End of explanation
for label in labels:
    print(results_9999[label].num_iter, '\t' + '(' + label + ')')
    df_9999[columns[0]].loc[label] = results_9999[label].num_iter
Explanation: The numbers of iterations:
End of explanation
print(results_9999['Policy iteration'].sigma)
Explanation: Policy iteration gives the optimal policy:
End of explanation
(1-results_9999['Policy iteration'].sigma).sum()
Explanation: Takes action 1 ("replace") if and only if $s \geq \bar{\gamma}$, where $\bar{\gamma}$ is equal to:
End of explanation
for result in results_9999.values():
    if result != results_9999['Policy iteration']:
        print(np.array_equal(result.sigma, results_9999['Policy iteration'].sigma))
Explanation: Check that the other methods gave the correct answer:
End of explanation
for label in labels:
    diff_pi = \
        np.abs(results_9999[label].v - results_9999['Policy iteration'].v).max()
    print(diff_pi, '\t' + '(' + label + ')')
    df_9999[columns[2]].loc[label] = diff_pi
Explanation: $\lVert v - v_{\mathrm{pi}}\rVert$:
End of explanation
for label in labels:
    v = results_9999[label].v
    diff_max = \
        np.abs(v - ddp.bellman_operator(v)).max()
    print(diff_max, '\t' + '(' + label + ')')
    df_9999[columns[4]].loc[label] = diff_max
Explanation: $\lVert v - T(v)\rVert$:
End of explanation
for i in range(4):
    if labels[i] != 'Policy iteration':
        k = 20 if labels[i] == 'Modified policy iteration' else 0
        res = ddp.solve(method=methods[i], v_init=v_init, k=k,
                        max_iter=results_9999[labels[i]].num_iter-1)
        diff = ddp.bellman_operator(res.v) - res.v
        diff_span = (diff.max() - diff.min()) * ddp.beta / (1 - ddp.beta)
        print(diff_span, '\t' + '(' + labels[i] + ')')
        df_9999[columns[3]].loc[labels[i]] = diff_span
Explanation: $\overline{b} - \underline{b}$:
End of explanation
label = 'Policy iteration'
v = results_9999[label].v
diff = ddp.bellman_operator(v) - v
diff_span = (diff.max() - diff.min()) * ddp.beta / (1 - ddp.beta)
print(diff_span, '\t' + '(' + label + ')')
df_9999[columns[3]].loc[label] = diff_span
for i in range(4):
    k = 20 if labels[i] == 'Modified policy iteration' else 0
    print(labels[i])
    t = %timeit -o ddp.solve(method=methods[i], v_init=v_init, epsilon=epsilon, k=k)
    df_9999[columns[1]].loc[labels[i]] = t.best
df_9999
df_time = pd.DataFrame(index=labels)
df_time[r'$\beta = 0.95$'] = df[columns[1]]
df_time[r'$\beta = 0.9999$'] = df_9999[columns[1]]
second_max = df_time[r'$\beta = 0.9999$'][1:].max()
for xlim in [None, (0, second_max*1.2)]:
    ax = df_time.loc[reversed(labels)][df_time.columns[::-1]].plot(
        kind='barh', legend='reverse', xlim=xlim, figsize=(8,5)
    )
    ax.set_xlabel('Time (second)')
import platform
print(platform.platform())
import sys
print(sys.version)
print(np.__version__)
Explanation: For policy iteration:
End of explanation
<END_TASK>