Just to replicate Figure 3, p.1025, given the estimated parameter values.
In [1]:
%matplotlib inline
In [2]:
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
In [3]:
def operator_iteration(T, w_0, tol, max_iter):
w = w_0
for i in range(max_iter):
new_w = T(w)
if np.abs(new_w - w).max() < tol:
w = new_w
break
w = new_w
num_iter = i + 1
return w, num_iter
In [4]:
class Zurcher(object):
def __init__(self, num_states, beta,
c=lambda x, theta_1: 0.001 * theta_1[0] * x):
self.num_states = num_states # 0, 5000, ..., 5000*(num_states-1)
self.num_actions = 2 # 0: keep, 1: replace
self.beta = beta
self.c = c
def set_params(self, RC, theta_1, theta_3):
self.RC, self.theta_1, self.theta_3 = RC, theta_1, theta_3
# Flow utility
self.u = np.empty((self.num_states, self.num_actions))
self.u[:, 0] = -self.c(np.arange(self.num_states)+1, self.theta_1)
self.u[:, 1] = -self.RC - self.c(0, self.theta_1)
# Transition probabilities
n = self.num_states
supp_size = len(theta_3)
p = np.zeros((n, n))
for i in range(n-supp_size+1):
p[i, i:i+supp_size] = theta_3
for i in range(n-supp_size+1, n):
p[i, i:] = theta_3[:n-i]/theta_3[:n-i].sum()
self.p = p
def continuation_values(self, ev):
"""
u + beta * ev
"""
v0 = self.u[:, 0] + self.beta * ev
v1 = self.u[0, 1] + self.beta * ev[0]
return v0, v1
def T(self, ev):
v0, v1 = self.continuation_values(ev)
return self.p.dot(np.logaddexp(v0, v1))
def replacement_prob(self, ev):
v0, v1 = self.continuation_values(ev)
return 1/(1 + np.exp(v0-v1))
def compute_EV(self, w_0=None, tol=1e-5, max_iter=10**5):
if w_0 is None:
w_0 = np.zeros(self.num_states)
EV, num_iter = operator_iteration(self.T, w_0=w_0, tol=tol, max_iter=max_iter)
P_replace = self.replacement_prob(EV)
return EV, num_iter, P_replace
Consider Model 11.
In [5]:
num_states = 90
Estimated parameter values for Group 4, Table IX, p.1021
In [6]:
# Transition probabilities
supp_size = 3
theta_3 = np.empty(supp_size)
theta_3[:-1] = .3919, .5953
theta_3[-1] = 1 - theta_3[:-1].sum()
In [7]:
betas = np.empty(2)
RCs = np.empty(2)
theta_1s = np.empty((2, 1))
In [8]:
# Estimated parameter values for beta = 0
betas[0] = 0
# Replacement cost
RCs[0] = 7.6358
# Cost function parameter
theta_1s[0, 0] = 71.5133
In [9]:
# Estimated parameter values for beta = .9999
betas[1] = .9999
# Replacement cost
RCs[1] = 10.0750
# Cost function parameter
theta_1s[1, 0] = 2.2930
Create Zurcher instances, one for each of $\beta = 0$ and $\beta = .9999$:
In [10]:
zurchers = [Zurcher(num_states=num_states, beta=betas[i]) for i in range(2)]
In [11]:
for i, zurcher in enumerate(zurchers):
zurcher.set_params(RC=RCs[i], theta_1=theta_1s[i], theta_3=theta_3)
Solve for $\mathit{EV}$:
In [12]:
EVs = np.empty((2, num_states))
num_iters = np.empty(2, dtype=int)
P_replaces = np.empty((2, num_states))
In [13]:
for i, zurcher in enumerate(zurchers):
EVs[i], num_iters[i], P_replaces[i] = zurcher.compute_EV()
In [14]:
num_iters
Out[14]:
In [15]:
EVs
Out[15]:
In [16]:
P_replaces
Out[16]:
In [17]:
# V for beta = .9999
V, _ = zurchers[1].continuation_values(EVs[1])
In [18]:
V
Out[18]:
Plot $c(x) - \beta \mathit{EV}(x)$ and $- \beta \mathit{EV}(x)$ for $\beta = .9999$ (TODO: find out how to "normalize" the values):
In [19]:
fig, ax = plt.subplots(figsize=(8,5))
ax.plot(-zurchers[1].beta * EVs[1], label=r'$c - \beta \mathit{EV}$')
ax.plot(-V, label=r'$- \beta \mathit{EV}$')
#for zurcher in zurchers:
# ax.plot(-zurcher.u[:, 0], label=r'$c$, $\beta = {0}$'.format(zurcher.beta))
ax.set_title('Estimated Value Functions')
ax.set_xlabel('State')
ax.set_ylabel(r'Non-normalized value $\times\ (-1)$')
plt.legend(loc=2)
plt.show()
Plot $P(1|x)$ for $\beta = 0$ and $\beta = .9999$:
In [20]:
fig, ax = plt.subplots(figsize=(8,5))
for i in range(2):
ax.plot(P_replaces[i], label=r'$\beta = {0}$'.format(betas[i]))
ax.set_ylim(0, 0.16)
ax.set_yticks(np.linspace(0, 0.16, 17, endpoint=True))
ax.set_title('Estimated Hazard Functions')
ax.set_xlabel('State')
ax.set_ylabel('Probability of engine replacement')
plt.legend(loc=2)
plt.show()
Timing the code (for $\beta = .9999$):
In [21]:
%timeit zurchers[1].compute_EV()
In [21]: