In [1]:
from __future__ import division
import numpy as np
import quantecon as qe
Exercises using Hamilton's Markov chain:
In [2]:
P = [[0.971, 0.029, 0 ],
[0.145, 0.778, 0.077],
[0 , 0.508, 0.492]]
mc_H = qe.MarkovChain(P)
In [3]:
mc_H.P
Out[3]:
In [4]:
states = {'NG': 0, 'MR': 1, 'SR': 2}
Let us use the following function from the previous exercise set:
In [5]:
def cross_sectional_dist(mc, init, T, num_reps=100):
"""
Return a distribution of visits at time T across num_reps sample paths
of mc with an initial state init.
"""
x = np.empty(num_reps, dtype=int)
for i in range(num_reps):
x[i] = mc.simulate(init=init, sample_size=T+1)[-1]
bins = np.arange(mc.n+1)
hist, bin_edges = np.histogram(x, bins=bins)
dist = hist/len(x)
return dist
In [6]:
psi = [0, 0, 1]
psi_10 = cross_sectional_dist(mc_H, init=psi, T=10)
print psi_10
In [7]:
psi_10[states['NG']]
Out[7]:
In [8]:
profits = [1000, 0, -1000]
In [9]:
t = 5
exp_profits_vec = np.linalg.matrix_power(mc_H.P, t).dot(profits)
print exp_profits_vec
The expected profit when the economy starts in NG:
In [10]:
exp_profits_vec[states['NG']]
Out[10]:
The change in profits when the economy starts in SR:
In [11]:
exp_profits_vec[states['SR']] - exp_profits_vec[states['NG']]
Out[11]:
In [12]:
t = 1000
exp_profits_vec_4213 = np.linalg.matrix_power(mc_H.P, t).dot(profits)
print exp_profits_vec_4213
In [13]:
psi = np.array([0.2, 0.2, 0.6])
t = 5
exp_profits = psi.dot(np.linalg.matrix_power(mc_H.P, t)).dot(profits)
print exp_profits
In [14]:
def path_prob(P, psi, X):
prob = psi[X[0]]
for t in range(len(X)-1):
prob *= P[X[t], X[t+1]]
return prob
In [15]:
psi = [0.2, 0.2, 0.6]
X = [states['NG'], states['MR'], states['NG']]
path_prob(mc_H.P, psi, X)
Out[15]:
A recursion version of the path_prob function:
In [16]:
def path_prob_recursion(P, psi, X):
if len(X) == 1:
return psi[X[0]]
return path_prob_recursion(P, psi, X[:-1]) * P[X[-2], X[-1]]
Note that, given a list X = [x^0, ..., x^{T-1}, x^T],
the path [x^0, ..., x^{T-1}] can be obtained by X[:-1]
(or equivalently, X[0:-1]).
In [17]:
path_prob_recursion(mc_H.P, psi, X)
Out[17]:
In [18]:
recession_states = [states[s] for s in states.keys() if s not in ['NG']]
print recession_states
We want to obtain the Cartesian product
recession_states x recession_states x recession_states.
One way is to use
cartesian
from quantecon:
In [19]:
from quantecon.cartesian import cartesian
paths = cartesian([recession_states for t in range(3)])
print paths
(Refer to this issue.)
Another is to use product
from itertools:
In [20]:
from itertools import product
paths = product(*[recession_states for t in range(3)])
for path in paths:
print path
In [21]:
paths = cartesian([recession_states for t in range(3)]).astype(int)
# paths = product(*[recession_states for t in range(3)])
prob = 0
for path in paths:
prob += path_prob(mc_H.P, psi, path)
print prob
Suppose that the path was [2, 0, 1]:
In [22]:
path = [2, 0, 1]
You can determine whether the path contains state NG by the not in operator:
In [23]:
states['NG'] not in path
Out[23]:
Count the number of True's out of 10,000 observations:
In [24]:
sample_size = 3
num_obs = 10000
is_recession = np.empty(num_obs, dtype=bool)
for i in range(num_obs):
path = mc_H.simulate(init=psi, sample_size=sample_size)
is_recession[i] = (states['NG'] not in path)
frequency = is_recession.sum()/num_obs
print frequency
You may use the mean method instead:
In [25]:
is_recession.mean()
Out[25]:
Recall:
In [26]:
profits
Out[26]:
In [27]:
psi
Out[27]:
In [28]:
def discounted_profit(h, r, X):
profit = 0
for t in range(len(X)):
profit += (1/(1+r)**t) * h[X[t]]
return profit
In [29]:
T = 2
r = 0.05
paths = cartesian([states.values() for t in range(T+1)]).astype(int)
NPV = 0
for path in paths:
NPV += discounted_profit(profits, r, path) * path_prob(mc_H.P, psi, path)
print NPV
In [30]:
num_obs = 10000
discounted_profit_values = np.empty(num_obs)
for i in range(num_obs):
path = mc_H.simulate(init=psi, sample_size=sample_size)
discounted_profit_values[i] = discounted_profit(profits, r, path)
avg_profit = discounted_profit_values.sum()/num_obs
print avg_profit
You may use the mean method instead:
In [31]:
discounted_profit_values.mean()
Out[31]:
In [32]:
def NPV_matrix_power(P, psi, h, r, T):
psi = np.asarray(psi) # To use the method dot from NumPy
NPV = 0
values = np.empty(T+1)
for t in range(T+1):
values[t] = \
(1/(1+r)**t) * psi.dot(np.linalg.matrix_power(P, t)).dot(h)
return values.sum()
In [33]:
NPV_matrix_power(mc_H.P, psi, profits, r, T)
Out[33]:
Another implementation starategy is to use recursion:
In [34]:
def NPV_recursion(P, psi, h, r, T):
if T == 0:
return np.dot(psi, h)
NPV = np.dot(psi, h) + \
(1/(1+r)) * NPV_recursion(P, np.dot(psi, P), h, r, T-1)
return NPV
In [35]:
NPV_recursion(mc_H.P, psi, profits, r, T)
Out[35]:
Plot NPVs against T:
In [36]:
%matplotlib inline
import matplotlib.pyplot as plt
In [37]:
ts = 7
NPV_values = np.empty(ts)
for t in range(ts):
NPV_values[t] = NPV_recursion(mc_H.P, psi, profits, r, t)
In [38]:
fig, ax = plt.subplots()
ax.plot(NPV_values)
ax.set_xlabel(r'$T$')
ax.set_ylabel('NPV')
ax.axhline(color='k', linewidth=0.5)
ax.grid(True)
plt.show()
In [39]:
mc_H.num_recurrent_classes
Out[39]:
mc_H has only one recurrent class, and hence has only one stationary distribution.
In [40]:
psi_star = mc_H.stationary_distributions[0]
In [41]:
print psi_star
In [42]:
extected_profit_stationary = psi_star.dot(profits)
print extected_profit_stationary
Recall from Exercise 4.2.13:
In [43]:
print exp_profits_vec_4213
In [44]:
mc_H.is_aperiodic
Out[44]:
mc_H is aperiodic, so that the marginal distribution converges to a stationary distribution,
which is unique.
In [45]:
psi_20 = cross_sectional_dist(mc_H, init=psi_star, T=20, num_reps=1000)
print psi_20
In [46]:
def sample_return_time(mc, state):
t = 0
x = state
while True:
t += 1
x = mc.simulate(init=x, sample_size=2)[-1]
if x == state:
break
return t
Or, alternatively:
In [47]:
def sample_return_time2(mc, state):
t = 1
x = mc.simulate(init=state, sample_size=2)[-1]
while (x != state):
t += 1
x = mc.simulate(init=x, sample_size=2)[-1]
return t
In [48]:
def simulated_mean_return_time(mc, state, num_reps):
return_times = np.empty(num_reps)
for i in range(num_reps):
return_times[i] = sample_return_time(mc, state)
return return_times.mean()
In [49]:
mc_H.is_irreducible
Out[49]:
mc_H is irreducible, so that every state is recurrent, and therefore for any initial state, the chain returns to it in finite time with probability one.
In [50]:
nums = [10**i for i in range(1, 5)]
mean_return_time_vec = np.empty(mc_H.n)
for num_reps in nums:
for state in range(mc_H.n):
mean_return_time_vec[state] = \
simulated_mean_return_time(mc_H, state, num_reps)
print 'num_reps = {0}:'.format(num_reps)
print ' mean return times = {0}'.format(mean_return_time_vec)
Compare:
In [51]:
psi_star ** (-1)
Out[51]:
In [52]:
h = np.zeros(mc_H.n)
h[states['NG']] = 1
print h
In [53]:
def h_mean(mc, init, h, sample_size):
path = mc.simulate(init=init, sample_size=sample_size)
h_sum = 0
for t in range(sample_size):
h_sum += h[path[t]]
return h_sum/sample_size
In [54]:
sample_size = 50000
init = 0
h_mean(mc_H, init=init, h=h, sample_size=sample_size)
Out[54]:
In [55]:
psi_star.dot(h)
Out[55]:
Recall:
In [56]:
profits
Out[56]:
In [57]:
h_mean(mc_H, init=init, h=profits, sample_size=sample_size)
Out[57]:
Recall from Exercise 4.3.9:
In [58]:
extected_profit_stationary
Out[58]:
In [58]: