from __future__ import division
import numpy as np
import quantecon as qe
Exercises using Hamilton's Markov chain:
P = [[0.971, 0.029, 0 ],
[0.145, 0.778, 0.077],
[0 , 0.508, 0.492]]
mc_H = qe.MarkovChain(P)
mc_H.P
array([[ 0.971, 0.029, 0. ], [ 0.145, 0.778, 0.077], [ 0. , 0.508, 0.492]])
states = {'NG': 0, 'MR': 1, 'SR': 2}
Let us use the following function from the previous exercise set:
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
psi = [0, 0, 1]
psi_10 = cross_sectional_dist(mc_H, init=psi, T=10)
print (psi_10)
[ 0.57 0.41 0.02]
psi_10[states['NG']]
0.56999999999999995
profits = [1000, 0, -1000]
t = 5
exp_profits_vec = np.linalg.matrix_power(mc_H.P, t).dot(profits)
print (exp_profits_vec)
[ 885.34763268 388.68267371 217.74304877]
exp_profits_vec[states['SR']]
217.74304876607994
t = 1000
exp_profits_vec = np.linalg.matrix_power(mc_H.P, t).dot(profits)
print (exp_profits_vec)
[ 788.16 788.16 788.16]
t = 5
psi=[0.2,0.2,0.6]
prob_t_ini=np.dot(psi,np.linalg.matrix_power(mc_H.P, t))
exp_profits = prob_t_ini.dot(profits)
print (exp_profits)
385.451890538
def path_prob(mc, init, X):
"""
calucualtion of the probablility of a path x
"""
prob=init[X[0]]
for t in range(len(X)-1):
prob=prob*mc[X[t],X[t+1]]
return prob
Path=[0,1,0]
path_010=path_prob(mc_H.P, init=psi, X=Path)
print(path_010)
0.000841
Y_0=[1,1,1]
Y_1=[1,2,1]
Y_2=[1,1,2]
Y_3=[1,2,2]
Y_4=[2,1,1]
Y_5=[2,1,2]
Y_6=[2,2,1]
Y_7=[2,2,2]
Y=[Y_0,Y_1,Y_2,Y_3,Y_4,Y_5,Y_6,Y_7]
path_recession= np.empty(len(Y))
for i in range(len(Y)):
path_recession[i]=path_prob(mc_H.P, init=psi, X=Y[i])
print(path_recession)
I=[1,1,1,1,1,1,1,1]
Pr_recession=np.dot(path_recession,I)
print(Pr_recession)
[ 0.1210568 0.0078232 0.0119812 0.0075768 0.2371344 0.0234696 0.1499616 0.1452384] 0.704242
Feasible_Path_0=[0,0,0]
Feasible_Path_1=[0,0,1]
Feasible_Path_2=[0,0,2]
Feasible_Path_3=[0,1,0]
Feasible_Path_4=[0,1,1]
Feasible_Path_5=[0,1,2]
Feasible_Path_6=[0,2,0]
Feasible_Path_7=[0,2,1]
Feasible_Path_8=[0,2,2]
Feasible_Path_9=[1,0,0]
Feasible_Path_10=[1,0,1]
Feasible_Path_11=[1,0,2]
Feasible_Path_12=[1,1,0]
Feasible_Path_13=[1,1,1]
Feasible_Path_14=[1,1,2]
Feasible_Path_15=[1,2,0]
Feasible_Path_16=[1,2,1]
Feasible_Path_17=[1,2,2]
Feasible_Path_18=[2,0,0]
Feasible_Path_19=[2,0,1]
Feasible_Path_20=[2,0,2]
Feasible_Path_21=[2,1,0]
Feasible_Path_22=[2,1,1]
Feasible_Path_23=[2,1,2]
Feasible_Path_24=[2,2,0]
Feasible_Path_25=[2,2,1]
Feasible_Path_26=[2,2,2]
Feasible_Path=[Feasible_Path_0,Feasible_Path_1,Feasible_Path_2,Feasible_Path_3
,Feasible_Path_4,Feasible_Path_5,Feasible_Path_6,Feasible_Path_7
,Feasible_Path_8,Feasible_Path_9,Feasible_Path_10,Feasible_Path_11
,Feasible_Path_12,Feasible_Path_13,Feasible_Path_14,Feasible_Path_15
,Feasible_Path_16,Feasible_Path_17,Feasible_Path_18,Feasible_Path_19
,Feasible_Path_20,Feasible_Path_21,Feasible_Path_22,Feasible_Path_23
,Feasible_Path_24,Feasible_Path_25,Feasible_Path_26,]
rho=1/(1+0.05)
T=2
profits_path=np.empty(len(Feasible_Path))
for j in range(len(Feasible_Path)):
for i in range(T):
profits_path[j]=profits_path[j]+rho**i*profits[Feasible_Path[j][i]]
Feasible_path_prob=np.empty(len(Feasible_Path))
for j in range(len(Feasible_Path)):
Feasible_path_prob[j]=path_prob(mc_H.P, init=psi, X=Feasible_Path[j])
print(profits_path)
print(Feasible_path_prob)
NPV=np.dot(Feasible_path_prob,profits_path)
print(NPV)
[ 1952.38095238 1952.38095238 1952.38095238 1000. 1000. 1000. 47.61904762 47.61904762 47.61904762 952.38095238 952.38095238 952.38095238 0. 0. 0. -952.38095238 -952.38095238 -952.38095238 -47.61904762 -47.61904762 -47.61904762 -1000. -1000. -1000. -1952.38095238 -1952.38095238 -1952.38095238] [ 0.1885682 0.0056318 0. 0.000841 0.0045124 0.0004466 0. 0. 0. 0.028159 0.000841 0. 0.022562 0.1210568 0.0119812 0. 0.0078232 0.0075768 0. 0. 0. 0.044196 0.2371344 0.0234696 0. 0.1499616 0.1452384] -483.238095238
T=2
NPV=0
prob_t_ini=np.dot(psi,np.linalg.matrix_power(mc_H.P, t))
exp_profits = prob_t_ini.dot(profits)
for t in range(T):
NPV=NPV+rho**t*exp_profits
print(NPV)
-170.638095238
vecs = mc_H.stationary_distributions
print(vecs)
exp_profits_st=vecs.dot(profits)
print(exp_profits_st)
[[ 0.8128 0.16256 0.02464]] [ 788.16]
from numpy import ones, identity,transpose
from numpy.linalg import solve
I=identity(3)
Q,b=ones((3,3)),ones((3,1))
A=transpose(I-mc_H.P+Q)
print (solve(A,b))
psi = solve(A,b)
psi_10 = cross_sectional_dist(mc_H, init=psi, T=20)
print (psi_10)
[[ 0.8128 ] [ 0.16256] [ 0.02464]] [ 0.83 0.15 0.02]