Skip to content
Merged
Show file tree
Hide file tree
Changes from 6 commits
Commits
Show all changes
19 commits
Select commit Hold shift + click to select a range
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions quantecon/_compute_fp.py
Original file line number Diff line number Diff line change
Expand Up @@ -262,11 +262,11 @@ def _compute_fixed_point_ig(T, v, max_iter, verbose, print_skip, is_approx_fp,
_, rho = _get_mixed_actions(tableaux_curr, bases_curr)

if Y.ndim <= 2:
x_new = rho.dot(Y[:m])
x_new = rho @ Y[:m]
else:
shape_Y = Y.shape
Y_2d = Y.reshape(shape_Y[0], np.prod(shape_Y[1:]))
x_new = rho.dot(Y_2d[:m]).reshape(shape_Y[1:])
x_new = (rho @ Y_2d[:m]).reshape(shape_Y[1:])

if verbose == 2:
error = np.max(np.abs(y_new - x_new))
Expand Down
146 changes: 73 additions & 73 deletions quantecon/_dle.py
Original file line number Diff line number Diff line change
Expand Up @@ -76,11 +76,11 @@ def __init__(self, information, technology, preferences):
uc = np.hstack((np.eye(self.nc), np.zeros((self.nc, self.ng))))
ug = np.hstack((np.zeros((self.ng, self.nc)), np.eye(self.ng)))
phiin = np.linalg.inv(np.hstack((self.phic, self.phig)))
phiinc = uc.dot(phiin)
b11 = - self.thetah.dot(phiinc).dot(self.phii)
a1 = self.thetah.dot(phiinc).dot(self.gamma)
a12 = np.vstack((self.thetah.dot(phiinc).dot(
self.ud), np.zeros((self.nk, self.nz))))
phiinc = uc @ phiin
b11 = - self.thetah @ phiinc @ self.phii
a1 = self.thetah @ phiinc @ self.gamma
a12 = np.vstack((self.thetah @ phiinc @
self.ud, np.zeros((self.nk, self.nz))))

# === Creation of the A Matrix for the state transition of the LQ problem === #

Expand All @@ -100,11 +100,11 @@ def __init__(self, information, technology, preferences):

# === Define R,W and Q for the payoff function of the LQ problem === #

self.H = np.hstack((self.llambda, self.pih.dot(uc).dot(phiin).dot(self.gamma), self.pih.dot(
uc).dot(phiin).dot(self.ud) - self.ub, -self.pih.dot(uc).dot(phiin).dot(self.phii)))
self.G = ug.dot(phiin).dot(
self.H = np.hstack((self.llambda, self.pih @ uc @ phiin @ self.gamma, self.pih @
uc @ phiin @ self.ud - self.ub, -self.pih @ uc @ phiin @ self.phii))
self.G = (ug @ phiin @
np.hstack((np.zeros((self.nd, self.nh)), self.gamma, self.ud, -self.phii)))
self.S = (self.G.T.dot(self.G) + self.H.T.dot(self.H)) / 2
self.S = (self.G.T @ self.G + self.H.T @ self.H) / 2

self.nx = self.nh + self.nk + self.nz
self.n = self.ni + self.nh + self.nk + self.nz
Expand All @@ -122,7 +122,7 @@ def __init__(self, information, technology, preferences):

# === Construct output matrices for our economy using the solution to the LQ problem === #

self.A0 = self.A - self.B.dot(self.F)
self.A0 = self.A - self.B @ self.F

self.Sh = self.A0[0:self.nh, 0:self.nx]
self.Sk = self.A0[self.nh:self.nh + self.nk, 0:self.nx]
Expand All @@ -131,12 +131,12 @@ def __init__(self, information, technology, preferences):
self.Si = -self.F
self.Sd = np.hstack((np.zeros((self.nd, self.nh + self.nk)), self.ud))
self.Sb = np.hstack((np.zeros((self.nb, self.nh + self.nk)), self.ub))
self.Sc = uc.dot(phiin).dot(-self.phii.dot(self.Si) +
self.gamma.dot(self.Sk1) + self.Sd)
self.Sg = ug.dot(phiin).dot(-self.phii.dot(self.Si) +
self.gamma.dot(self.Sk1) + self.Sd)
self.Ss = self.llambda.dot(np.hstack((np.eye(self.nh), np.zeros(
(self.nh, self.nk + self.nz))))) + self.pih.dot(self.Sc)
self.Sc = uc @ phiin @ (-self.phii @ self.Si +
self.gamma @ self.Sk1 + self.Sd)
self.Sg = ug @ phiin @ (-self.phii @ self.Si +
self.gamma @ self.Sk1 + self.Sd)
self.Ss = self.llambda @ np.hstack((np.eye(self.nh), np.zeros(
(self.nh, self.nk + self.nz)))) + self.pih @ self.Sc

# === Calculate eigenvalues of A0 === #
self.A110 = self.A0[0:self.nh + self.nk, 0:self.nh + self.nk]
Expand All @@ -146,14 +146,14 @@ def __init__(self, information, technology, preferences):
# === Construct matrices for Lagrange Multipliers === #

self.Mk = -2 * self.beta.item() * (np.hstack((np.zeros((self.nk, self.nh)), np.eye(
self.nk), np.zeros((self.nk, self.nz))))).dot(self.P).dot(self.A0)
self.nk), np.zeros((self.nk, self.nz))))) @ self.P @ self.A0
self.Mh = -2 * self.beta.item() * (np.hstack((np.eye(self.nh), np.zeros(
(self.nh, self.nk)), np.zeros((self.nh, self.nz))))).dot(self.P).dot(self.A0)
(self.nh, self.nk)), np.zeros((self.nh, self.nz))))) @ self.P @ self.A0
self.Ms = -(self.Sb - self.Ss)
self.Md = -(np.linalg.inv(np.vstack((self.phic.T, self.phig.T))).dot(
np.vstack((self.thetah.T.dot(self.Mh) + self.pih.T.dot(self.Ms), -self.Sg))))
self.Mc = -(self.thetah.T.dot(self.Mh) + self.pih.T.dot(self.Ms))
self.Mi = -(self.thetak.T.dot(self.Mk))
self.Md = -(np.linalg.inv(np.vstack((self.phic.T, self.phig.T))) @
np.vstack((self.thetah.T @ self.Mh + self.pih.T @ self.Ms, -self.Sg)))
self.Mc = -(self.thetah.T @ self.Mh + self.pih.T @ self.Ms)
self.Mi = -(self.thetak.T @ self.Mk)

def compute_steadystate(self, nnc=2):
"""
Expand All @@ -168,13 +168,13 @@ def compute_steadystate(self, nnc=2):
zx = np.eye(self.A0.shape[0])-self.A0
self.zz = nullspace(zx)
self.zz /= self.zz[nnc]
self.css = self.Sc.dot(self.zz)
self.sss = self.Ss.dot(self.zz)
self.iss = self.Si.dot(self.zz)
self.dss = self.Sd.dot(self.zz)
self.bss = self.Sb.dot(self.zz)
self.kss = self.Sk.dot(self.zz)
self.hss = self.Sh.dot(self.zz)
self.css = self.Sc @ self.zz
self.sss = self.Ss @ self.zz
self.iss = self.Si @ self.zz
self.dss = self.Sd @ self.zz
self.bss = self.Sb @ self.zz
self.kss = self.Sk @ self.zz
self.hss = self.Sh @ self.zz

def compute_sequence(self, x0, ts_length=None, Pay=None):
"""
Expand All @@ -195,14 +195,14 @@ def compute_sequence(self, x0, ts_length=None, Pay=None):
lq = LQ(self.Q, self.R, self.A, self.B,
self.C, N=self.W, beta=self.beta)
xp, up, wp = lq.compute_sequence(x0, ts_length)
self.h = self.Sh.dot(xp)
self.k = self.Sk.dot(xp)
self.i = self.Si.dot(xp)
self.b = self.Sb.dot(xp)
self.d = self.Sd.dot(xp)
self.c = self.Sc.dot(xp)
self.g = self.Sg.dot(xp)
self.s = self.Ss.dot(xp)
self.h = self.Sh @ xp
self.k = self.Sk @ xp
self.i = self.Si @ xp
self.b = self.Sb @ xp
self.d = self.Sd @ xp
self.c = self.Sc @ xp
self.g = self.Sg @ xp
self.s = self.Ss @ xp

# === Value of J-period risk-free bonds === #
# === See p.145: Equation (7.11.2) === #
Expand All @@ -212,12 +212,12 @@ def compute_sequence(self, x0, ts_length=None, Pay=None):
self.R2_Price = np.empty((ts_length + 1, 1))
self.R5_Price = np.empty((ts_length + 1, 1))
for i in range(ts_length + 1):
self.R1_Price[i, 0] = self.beta * e1.dot(self.Mc).dot(np.linalg.matrix_power(
self.A0, 1)).dot(xp[:, i]) / e1.dot(self.Mc).dot(xp[:, i])
self.R2_Price[i, 0] = self.beta**2 * e1.dot(self.Mc).dot(
np.linalg.matrix_power(self.A0, 2)).dot(xp[:, i]) / e1.dot(self.Mc).dot(xp[:, i])
self.R5_Price[i, 0] = self.beta**5 * e1.dot(self.Mc).dot(
np.linalg.matrix_power(self.A0, 5)).dot(xp[:, i]) / e1.dot(self.Mc).dot(xp[:, i])
self.R1_Price[i, 0] = self.beta * e1 @ self.Mc @ np.linalg.matrix_power(
self.A0, 1) @ xp[:, i] / (e1 @ self.Mc @ xp[:, i])
self.R2_Price[i, 0] = self.beta**2 * (e1 @ self.Mc @
np.linalg.matrix_power(self.A0, 2) @ xp[:, i]) / (e1 @ self.Mc @ xp[:, i])
self.R5_Price[i, 0] = self.beta**5 * (e1 @ self.Mc @
np.linalg.matrix_power(self.A0, 5) @ xp[:, i]) / (e1 @ self.Mc @ xp[:, i])

# === Gross rates of return on 1-period risk-free bonds === #
self.R1_Gross = 1 / self.R1_Price
Expand All @@ -231,20 +231,20 @@ def compute_sequence(self, x0, ts_length=None, Pay=None):
# === Value of asset whose payout vector is Pay*xt === #
# See p.145: Equation (7.11.1)
if isinstance(Pay, np.ndarray) == True:
self.Za = Pay.T.dot(self.Mc)
self.Za = Pay.T @ self.Mc
self.Q = solve_discrete_lyapunov(
self.A0.T * self.beta**0.5, self.Za)
self.q = self.beta / (1 - self.beta) * \
np.trace(self.C.T.dot(self.Q).dot(self.C))
np.trace(self.C.T @ self.Q @ self.C)
self.Pay_Price = np.empty((ts_length + 1, 1))
self.Pay_Gross = np.empty((ts_length + 1, 1))
self.Pay_Gross[0, 0] = np.nan
for i in range(ts_length + 1):
self.Pay_Price[i, 0] = (xp[:, i].T.dot(self.Q).dot(
xp[:, i]) + self.q) / e1.dot(self.Mc).dot(xp[:, i])
self.Pay_Price[i, 0] = (xp[:, i].T @ self.Q @
xp[:, i] + self.q) / (e1 @ self.Mc @ xp[:, i])
for i in range(ts_length):
self.Pay_Gross[i + 1, 0] = self.Pay_Price[i + 1,
0] / (self.Pay_Price[i, 0] - Pay.dot(xp[:, i]))
0] / (self.Pay_Price[i, 0] - Pay @ xp[:, i])
return

def irf(self, ts_length=100, shock=None):
Expand Down Expand Up @@ -276,22 +276,22 @@ def irf(self, ts_length=100, shock=None):
self.b_irf = np.empty((ts_length, self.nb))

for i in range(ts_length):
self.c_irf[i, :] = self.Sc.dot(
np.linalg.matrix_power(self.A0, i)).dot(self.C).dot(shock).T
self.s_irf[i, :] = self.Ss.dot(
np.linalg.matrix_power(self.A0, i)).dot(self.C).dot(shock).T
self.i_irf[i, :] = self.Si.dot(
np.linalg.matrix_power(self.A0, i)).dot(self.C).dot(shock).T
self.k_irf[i, :] = self.Sk.dot(
np.linalg.matrix_power(self.A0, i)).dot(self.C).dot(shock).T
self.h_irf[i, :] = self.Sh.dot(
np.linalg.matrix_power(self.A0, i)).dot(self.C).dot(shock).T
self.g_irf[i, :] = self.Sg.dot(
np.linalg.matrix_power(self.A0, i)).dot(self.C).dot(shock).T
self.d_irf[i, :] = self.Sd.dot(
np.linalg.matrix_power(self.A0, i)).dot(self.C).dot(shock).T
self.b_irf[i, :] = self.Sb.dot(
np.linalg.matrix_power(self.A0, i)).dot(self.C).dot(shock).T
self.c_irf[i, :] = (self.Sc @
np.linalg.matrix_power(self.A0, i) @ self.C @ shock).T
self.s_irf[i, :] = (self.Ss @
np.linalg.matrix_power(self.A0, i) @ self.C @ shock).T
self.i_irf[i, :] = (self.Si @
np.linalg.matrix_power(self.A0, i) @ self.C @ shock).T
self.k_irf[i, :] = (self.Sk @
np.linalg.matrix_power(self.A0, i) @ self.C @ shock).T
self.h_irf[i, :] = (self.Sh @
np.linalg.matrix_power(self.A0, i) @ self.C @ shock).T
self.g_irf[i, :] = (self.Sg @
np.linalg.matrix_power(self.A0, i) @ self.C @ shock).T
self.d_irf[i, :] = (self.Sd @
np.linalg.matrix_power(self.A0, i) @ self.C @ shock).T
self.b_irf[i, :] = (self.Sb @
np.linalg.matrix_power(self.A0, i) @ self.C @ shock).T

return

Expand All @@ -305,13 +305,13 @@ def canonical(self):
Ac2 = np.hstack((np.zeros((self.nz, self.nh)), self.a22))
Ac = np.vstack((Ac1, Ac2))
Bc = np.vstack((self.thetah, np.zeros((self.nz, self.nc))))
Rc1 = np.hstack((self.llambda.T.dot(self.llambda), -
self.llambda.T.dot(self.ub)))
Rc2 = np.hstack((-self.ub.T.dot(self.llambda), self.ub.T.dot(self.ub)))
Rc1 = np.hstack((self.llambda.T @ self.llambda, -
self.llambda.T @ self.ub))
Rc2 = np.hstack((-self.ub.T @ self.llambda, self.ub.T @ self.ub))
Rc = np.vstack((Rc1, Rc2))
Qc = self.pih.T.dot(self.pih)
Qc = self.pih.T @ self.pih
Nc = np.hstack(
(self.pih.T.dot(self.llambda), -self.pih.T.dot(self.ub)))
(self.pih.T @ self.llambda, -self.pih.T @ self.ub))

lq_aux = LQ(Qc, Rc, Ac, Bc, N=Nc, beta=self.beta)

Expand All @@ -320,9 +320,9 @@ def canonical(self):
self.F_b = F1[:, 0:self.nh]
self.F_f = F1[:, self.nh:]

self.pihat = np.linalg.cholesky(self.pih.T.dot(
self.pih) + self.beta.dot(self.thetah.T).dot(P1[0:self.nh, 0:self.nh]).dot(self.thetah)).T
self.llambdahat = self.pihat.dot(self.F_b)
self.ubhat = - self.pihat.dot(self.F_f)
self.pihat = np.linalg.cholesky((self.pih.T @
self.pih) + (self.beta @ self.thetah.T @ P1[0:self.nh, 0:self.nh] @ self.thetah)).T
self.llambdahat = self.pihat @ self.F_b
self.ubhat = - self.pihat @ self.F_f

return
44 changes: 22 additions & 22 deletions quantecon/_kalman.py
Original file line number Diff line number Diff line change
Expand Up @@ -162,9 +162,9 @@ def whitener_lss(self):
A, C, G, H = self.ss.A, self.ss.C, self.ss.G, self.ss.H

Atil = np.vstack([np.hstack([A, np.zeros((n, n)), np.zeros((n, l))]),
np.hstack([np.dot(K, G),
A-np.dot(K, G),
np.dot(K, H)]),
np.hstack([K @ G,
A - (K @ G),
K @ H]),
np.zeros((l, 2*n + l))])

Ctil = np.vstack([np.hstack([C, np.zeros((n, l))]),
Expand Down Expand Up @@ -200,16 +200,16 @@ def prior_to_filtered(self, y):
"""
# === simplify notation === #
G, H = self.ss.G, self.ss.H
R = np.dot(H, H.T)
R = H @ H.T

# === and then update === #
y = np.atleast_2d(y)
y.shape = self.ss.k, 1
E = np.dot(self.Sigma, G.T)
F = np.dot(np.dot(G, self.Sigma), G.T) + R
M = np.dot(E, inv(F))
self.x_hat = self.x_hat + np.dot(M, (y - np.dot(G, self.x_hat)))
self.Sigma = self.Sigma - np.dot(M, np.dot(G, self.Sigma))
E = self.Sigma @ G.T
F = (G @ self.Sigma @ G.T) + R
M = E @ inv(F)
self.x_hat = self.x_hat + (M @ (y - (G @ self.x_hat)))
self.Sigma = self.Sigma - (M @ G @ self.Sigma)

def filtered_to_forecast(self):
"""
Expand All @@ -220,11 +220,11 @@ def filtered_to_forecast(self):
"""
# === simplify notation === #
A, C = self.ss.A, self.ss.C
Q = np.dot(C, C.T)
Q = C @ C.T

# === and then update === #
self.x_hat = np.dot(A, self.x_hat)
self.Sigma = np.dot(A, np.dot(self.Sigma, A.T)) + Q
self.x_hat = A @ self.x_hat
self.Sigma = (A @ self.Sigma @ A.T) + Q

def update(self, y):
"""
Expand Down Expand Up @@ -265,13 +265,13 @@ def stationary_values(self, method='doubling'):

# === simplify notation === #
A, C, G, H = self.ss.A, self.ss.C, self.ss.G, self.ss.H
Q, R = np.dot(C, C.T), np.dot(H, H.T)
Q, R = C @ C.T, H @ H.T

# === solve Riccati equation, obtain Kalman gain === #
Sigma_infinity = solve_discrete_riccati(A.T, G.T, Q, R, method=method)
temp1 = np.dot(np.dot(A, Sigma_infinity), G.T)
temp2 = inv(np.dot(G, np.dot(Sigma_infinity, G.T)) + R)
K_infinity = np.dot(temp1, temp2)
temp1 = (A @ Sigma_infinity @ G.T)
temp2 = inv((G @ Sigma_infinity @ G.T) + R)
K_infinity = temp1 @ temp2

# == record as attributes and return == #
self._Sigma_infinity, self._K_infinity = Sigma_infinity, K_infinity
Expand Down Expand Up @@ -302,21 +302,21 @@ def stationary_coefficients(self, j, coeff_type='ma'):
P_mat = A
P = np.identity(self.ss.n) # Create a copy
elif coeff_type == 'var':
coeffs.append(np.dot(G, K_infinity))
P_mat = A - np.dot(K_infinity, G)
coeffs.append(G @ K_infinity)
P_mat = A - (K_infinity @ G)
P = np.copy(P_mat) # Create a copy
else:
raise ValueError("Unknown coefficient type")
while i <= j:
coeffs.append(np.dot(np.dot(G, P), K_infinity))
P = np.dot(P, P_mat)
coeffs.append((G @ P) @ K_infinity)
P = P @ P_mat
i += 1
return coeffs

def stationary_innovation_covar(self):
# == simplify notation == #
H, G = self.ss.H, self.ss.G
R = np.dot(H, H.T)
R = H @ H.T
Sigma_infinity = self.Sigma_infinity

return np.dot(G, np.dot(Sigma_infinity, G.T)) + R
return (G @ Sigma_infinity @ G.T) + R
Loading
Loading