17227928141

Committed 26 Aug 2025 04:28AM UTC coverage: 92.626%. Remained the same

Build # 17227928141

Build Type

Pull #787

github

Committed by

web-flow

Commit Message

Merge 2d53230c0 into 4a889ad1b

Pull Request Pull Request #787: Migrate np.dot and .dot() method calls to @ operator in library code only

Run Details

113 of 158 new or added lines in 14 files covered. (71.52%)

1 existing line in 1 file now uncovered.

7512 of 8110 relevant lines covered (92.63%)

2.78 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

60.99

/quantecon/_dle.py

"""
Provides a class called DLE to convert and solve dynamic linear economies
(as set out in Hansen & Sargent (2013)) as LQ problems.
"""

import numpy as np
from ._lqcontrol import LQ
from ._matrix_eqn import solve_discrete_lyapunov
from ._rank_nullspace import nullspace

class DLE(object):
    r"""
    This class is for analyzing dynamic linear economies, as set out in Hansen & Sargent (2013).
    The planner's problem is to choose \{c_t, s_t, i_t, h_t, k_t, g_t\}_{t=0}^\infty to maximize

        \max -(1/2) \mathbb{E} \sum_{t=0}^{\infty} \beta^t [(s_t - b_t).(s_t-b_t) + g_t.g_t]

    subject to the linear constraints

        \Phi_c c_t + \Phi_g g_t + \Phi_i i_t = \Gamma k_{t-1} + d_t
        k_t = \Delta_k k_{t-1} + \Theta_k i_t
        h_t = \Delta_h h_{t-1} + \Theta_h c_t
        s_t = \Lambda h_{t-1} + \Pi c_t

    and

        z_{t+1} = A_{22} z_t + C_2 w_{t+1}
        b_t = U_b z_t
        d_t = U_d z_t

    where h_{-1}, k_{-1}, and z_0 are given as initial conditions.

    Section 5.5 of HS2013 describes how to map these matrices into those of
    a LQ problem.

    HS2013 sort the matrices defining the problem into three groups:

    Information: A_{22}, C_2, U_b , and U_d characterize the motion of information
    sets and of taste and technology shocks

    Technology: \Phi_c, \Phi_g, \Phi_i, \Gamma, \Delta_k, and \Theta_k determine the
    technology for producing consumption goods

    Preferences: \Delta_h, \Theta_h, \Lambda, and \Pi determine the technology for
    producing consumption services from consumer goods. A scalar discount factor \beta
    determines the preference ordering over consumption services.

    Parameters
    ----------
    Information : tuple
        Information is a tuple containing the matrices A_{22}, C_2, U_b, and U_d
    Technology : tuple
        Technology is a tuple containing the matrices \Phi_c, \Phi_g, \Phi_i, \Gamma,
        \Delta_k, and \Theta_k
    Preferences : tuple
        Preferences is a tuple containing the scalar \beta and the
        matrices \Lambda, \Pi, \Delta_h, and \Theta_h

    """

    def __init__(self, information, technology, preferences):

        # === Unpack the tuples which define information, technology and preferences === #
        self.a22, self.c2, self.ub, self.ud = information
        self.phic, self.phig, self.phii, self.gamma, self.deltak, self.thetak = technology
        self.beta, self.llambda, self.pih, self.deltah, self.thetah = preferences

        # === Computation of the dimension of the structural parameter matrices === #
        self.nb, self.nh = self.llambda.shape
        self.nd, self.nc = self.phic.shape
        self.nz, self.nw = self.c2.shape
        _, self.ng = self.phig.shape
        self.nk, self.ni = self.thetak.shape

        # === Creation of various useful matrices === #
        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 @ 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 === #

        a11 = np.vstack((np.hstack((self.deltah, a1)), np.hstack(
            (np.zeros((self.nk, self.nh)), self.deltak))))
        self.A = np.vstack((np.hstack((a11, a12)), np.hstack(
            (np.zeros((self.nz, self.nk + self.nh)), self.a22))))

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

        b1 = np.vstack((b11, self.thetak))
        self.B = np.vstack((b1, np.zeros((self.nz, self.ni))))

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

        self.C = np.vstack((np.zeros((self.nk + self.nh, self.nw)), self.c2))

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

        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 @ 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

        self.R = self.S[0:self.nx, 0:self.nx]
        self.W = self.S[self.nx:self.n, 0:self.nx]
        self.Q = self.S[self.nx:self.n, self.nx:self.n]

        # === Use quantecon's LQ code to solve our LQ problem === #

        lq = LQ(self.Q, self.R, self.A, self.B,
                self.C, N=self.W, beta=self.beta)

        self.P, self.F, self.d = lq.stationary_values()

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

        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]
        self.Sk1 = np.hstack((np.zeros((self.nk, self.nh)), np.eye(
            self.nk), np.zeros((self.nk, self.nz))))
        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 @ 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]
        self.endo = np.linalg.eigvals(self.A110)
        self.exo = np.linalg.eigvals(self.a22)

        # === 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))))) @ 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))))) @ self.P @ self.A0
        self.Ms = -(self.Sb - self.Ss)
        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):
        """
        Computes the non-stochastic steady-state of the economy.

        Parameters
        ----------
        nnc : array_like(float)
            nnc is the location of the constant in the state vector x_t

        """
        zx = np.eye(self.A0.shape[0])-self.A0
        self.zz = nullspace(zx)
        self.zz /= self.zz[nnc]
        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):
        """
        Simulate quantities and prices for the economy

        Parameters
        ----------
        x0 : array_like(float)
            The initial state

        ts_length : scalar(int)
            Length of the simulation

        Pay : array_like(float)
            Vector to price an asset whose payout is Pay*xt

        """
        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 @ 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) === #
        e1 = np.zeros((1, self.nc))
        e1[0, 0] = 1
        self.R1_Price = np.empty((ts_length + 1, 1))
        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 @ 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

        # === Net rates of return on J-period risk-free bonds === #
        # === See p.148: log of gross rate of return, divided by j === #
        self.R1_Net = np.log(1 / self.R1_Price) / 1
        self.R2_Net = np.log(1 / self.R2_Price) / 2
        self.R5_Net = np.log(1 / self.R5_Price) / 5

        # === 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 @ 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 @ 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 @ 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 @ xp[:, i])
        return

    def irf(self, ts_length=100, shock=None):
        """
        Create Impulse Response Functions

        Parameters
        ----------

        ts_length : scalar(int)
            Number of periods to calculate IRF

        Shock : array_like(float)
            Vector of shocks to calculate IRF to. Default is first element of w

        """

        if type(shock) != np.ndarray:
            # Default is to select first element of w
            shock = np.vstack((np.ones((1, 1)), np.zeros((self.nw - 1, 1))))

        self.c_irf = np.empty((ts_length, self.nc))
        self.s_irf = np.empty((ts_length, self.nb))
        self.i_irf = np.empty((ts_length, self.ni))
        self.k_irf = np.empty((ts_length, self.nk))
        self.h_irf = np.empty((ts_length, self.nh))
        self.g_irf = np.empty((ts_length, self.ng))
        self.d_irf = np.empty((ts_length, self.nd))
        self.b_irf = np.empty((ts_length, self.nb))

        for i in range(ts_length):
            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

    def canonical(self):
        """
        Compute canonical preference representation
        Uses auxiliary problem of 9.4.2, with the preference shock process reintroduced
        Calculates pihat, llambdahat and ubhat for the equivalent canonical household technology
        """
        Ac1 = np.hstack((self.deltah, np.zeros((self.nh, self.nz))))
        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 @ 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 @ self.pih
        Nc = np.hstack(
            (self.pih.T @ self.llambda, -self.pih.T @ self.ub))

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

        P1, F1, d1 = lq_aux.stationary_values()

        self.F_b = F1[:, 0:self.nh]
        self.F_f = F1[:, self.nh:]

        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

1	"""
2	Provides a class called DLE to convert and solve dynamic linear economies
3	(as set out in Hansen & Sargent (2013)) as LQ problems.
4	"""
5
6	import numpy as np	3✔
7	from ._lqcontrol import LQ	3✔
8	from ._matrix_eqn import solve_discrete_lyapunov	3✔
9	from ._rank_nullspace import nullspace	3✔
10
11	class DLE(object):	3✔
12	r"""
13	This class is for analyzing dynamic linear economies, as set out in Hansen & Sargent (2013).
14	The planner's problem is to choose \{c_t, s_t, i_t, h_t, k_t, g_t\}_{t=0}^\infty to maximize
15
16	\max -(1/2) \mathbb{E} \sum_{t=0}^{\infty} \beta^t [(s_t - b_t).(s_t-b_t) + g_t.g_t]
17
18	subject to the linear constraints
19
20	\Phi_c c_t + \Phi_g g_t + \Phi_i i_t = \Gamma k_{t-1} + d_t
21	k_t = \Delta_k k_{t-1} + \Theta_k i_t
22	h_t = \Delta_h h_{t-1} + \Theta_h c_t
23	s_t = \Lambda h_{t-1} + \Pi c_t
24
25	and
26
27	z_{t+1} = A_{22} z_t + C_2 w_{t+1}
28	b_t = U_b z_t
29	d_t = U_d z_t
30
31	where h_{-1}, k_{-1}, and z_0 are given as initial conditions.
32
33	Section 5.5 of HS2013 describes how to map these matrices into those of
34	a LQ problem.
35
36	HS2013 sort the matrices defining the problem into three groups:
37
38	Information: A_{22}, C_2, U_b , and U_d characterize the motion of information
39	sets and of taste and technology shocks
40
41	Technology: \Phi_c, \Phi_g, \Phi_i, \Gamma, \Delta_k, and \Theta_k determine the
42	technology for producing consumption goods
43
44	Preferences: \Delta_h, \Theta_h, \Lambda, and \Pi determine the technology for
45	producing consumption services from consumer goods. A scalar discount factor \beta
46	determines the preference ordering over consumption services.
47
48	Parameters
49	----------
50	Information : tuple
51	Information is a tuple containing the matrices A_{22}, C_2, U_b, and U_d
52	Technology : tuple
53	Technology is a tuple containing the matrices \Phi_c, \Phi_g, \Phi_i, \Gamma,
54	\Delta_k, and \Theta_k
55	Preferences : tuple
56	Preferences is a tuple containing the scalar \beta and the
57	matrices \Lambda, \Pi, \Delta_h, and \Theta_h
58
59	"""
60
61	def __init__(self, information, technology, preferences):	3✔
62
63	# === Unpack the tuples which define information, technology and preferences === #
64	self.a22, self.c2, self.ub, self.ud = information	3✔
65	self.phic, self.phig, self.phii, self.gamma, self.deltak, self.thetak = technology	3✔
66	self.beta, self.llambda, self.pih, self.deltah, self.thetah = preferences	3✔
67
68	# === Computation of the dimension of the structural parameter matrices === #
69	self.nb, self.nh = self.llambda.shape	3✔
70	self.nd, self.nc = self.phic.shape	3✔
71	self.nz, self.nw = self.c2.shape	3✔
72	_, self.ng = self.phig.shape	3✔
73	self.nk, self.ni = self.thetak.shape	3✔
74
75	# === Creation of various useful matrices === #
76	uc = np.hstack((np.eye(self.nc), np.zeros((self.nc, self.ng))))	3✔
77	ug = np.hstack((np.zeros((self.ng, self.nc)), np.eye(self.ng)))	3✔
78	phiin = np.linalg.inv(np.hstack((self.phic, self.phig)))	3✔
79	phiinc = uc @ phiin	3✔
80	b11 = - self.thetah @ phiinc @ self.phii	3✔
81	a1 = self.thetah @ phiinc @ self.gamma	3✔
82	a12 = np.vstack((self.thetah @ phiinc @	3✔
83	self.ud, np.zeros((self.nk, self.nz))))
84
85	# === Creation of the A Matrix for the state transition of the LQ problem === #
86
87	a11 = np.vstack((np.hstack((self.deltah, a1)), np.hstack(	3✔
88	(np.zeros((self.nk, self.nh)), self.deltak))))
89	self.A = np.vstack((np.hstack((a11, a12)), np.hstack(	3✔
90	(np.zeros((self.nz, self.nk + self.nh)), self.a22))))
91
92	# === Creation of the B Matrix for the state transition of the LQ problem === #
93
94	b1 = np.vstack((b11, self.thetak))	3✔
95	self.B = np.vstack((b1, np.zeros((self.nz, self.ni))))	3✔
96
97	# === Creation of the C Matrix for the state transition of the LQ problem === #
98
99	self.C = np.vstack((np.zeros((self.nk + self.nh, self.nw)), self.c2))	3✔
100
101	# === Define R,W and Q for the payoff function of the LQ problem === #
102
103	self.H = np.hstack((self.llambda, self.pih @ uc @ phiin @ self.gamma, self.pih @	3✔
104	uc @ phiin @ self.ud - self.ub, -self.pih @ uc @ phiin @ self.phii))
105	self.G = (ug @ phiin @	3✔
106	np.hstack((np.zeros((self.nd, self.nh)), self.gamma, self.ud, -self.phii)))
107	self.S = (self.G.T @ self.G + self.H.T @ self.H) / 2	3✔
108
109	self.nx = self.nh + self.nk + self.nz	3✔
110	self.n = self.ni + self.nh + self.nk + self.nz	3✔
111
112	self.R = self.S[0:self.nx, 0:self.nx]	3✔
113	self.W = self.S[self.nx:self.n, 0:self.nx]	3✔
114	self.Q = self.S[self.nx:self.n, self.nx:self.n]	3✔
115
116	# === Use quantecon's LQ code to solve our LQ problem === #
117
118	lq = LQ(self.Q, self.R, self.A, self.B,	3✔
119	self.C, N=self.W, beta=self.beta)
120
121	self.P, self.F, self.d = lq.stationary_values()	3✔
122
123	# === Construct output matrices for our economy using the solution to the LQ problem === #
124
125	self.A0 = self.A - self.B @ self.F	3✔
126
127	self.Sh = self.A0[0:self.nh, 0:self.nx]	3✔
128	self.Sk = self.A0[self.nh:self.nh + self.nk, 0:self.nx]	3✔
129	self.Sk1 = np.hstack((np.zeros((self.nk, self.nh)), np.eye(	3✔
130	self.nk), np.zeros((self.nk, self.nz))))
131	self.Si = -self.F	3✔
132	self.Sd = np.hstack((np.zeros((self.nd, self.nh + self.nk)), self.ud))	3✔
133	self.Sb = np.hstack((np.zeros((self.nb, self.nh + self.nk)), self.ub))	3✔
134	self.Sc = uc @ phiin @ (-self.phii @ self.Si +	3✔
135	self.gamma @ self.Sk1 + self.Sd)
136	self.Sg = ug @ phiin @ (-self.phii @ self.Si +	3✔
137	self.gamma @ self.Sk1 + self.Sd)
138	self.Ss = self.llambda @ np.hstack((np.eye(self.nh), np.zeros(	3✔
139	(self.nh, self.nk + self.nz)))) + self.pih @ self.Sc
140
141	# === Calculate eigenvalues of A0 === #
142	self.A110 = self.A0[0:self.nh + self.nk, 0:self.nh + self.nk]	3✔
143	self.endo = np.linalg.eigvals(self.A110)	3✔
144	self.exo = np.linalg.eigvals(self.a22)	3✔
145
146	# === Construct matrices for Lagrange Multipliers === #
147
148	self.Mk = -2 * self.beta.item() * (np.hstack((np.zeros((self.nk, self.nh)), np.eye(	3✔
149	self.nk), np.zeros((self.nk, self.nz))))) @ self.P @ self.A0
150	self.Mh = -2 * self.beta.item() * (np.hstack((np.eye(self.nh), np.zeros(	3✔
151	(self.nh, self.nk)), np.zeros((self.nh, self.nz))))) @ self.P @ self.A0
152	self.Ms = -(self.Sb - self.Ss)	3✔
153	self.Md = -(np.linalg.inv(np.vstack((self.phic.T, self.phig.T))) @	3✔
154	np.vstack((self.thetah.T @ self.Mh + self.pih.T @ self.Ms, -self.Sg)))
155	self.Mc = -(self.thetah.T @ self.Mh + self.pih.T @ self.Ms)	3✔
156	self.Mi = -(self.thetak.T @ self.Mk)	3✔
157
158	def compute_steadystate(self, nnc=2):	3✔
159	"""
160	Computes the non-stochastic steady-state of the economy.
161
162	Parameters
163	----------
164	nnc : array_like(float)
165	nnc is the location of the constant in the state vector x_t
166
167	"""
168	zx = np.eye(self.A0.shape[0])-self.A0	3✔
169	self.zz = nullspace(zx)	3✔
170	self.zz /= self.zz[nnc]	3✔
171	self.css = self.Sc @ self.zz	3✔
172	self.sss = self.Ss @ self.zz	3✔
173	self.iss = self.Si @ self.zz	3✔
174	self.dss = self.Sd @ self.zz	3✔
175	self.bss = self.Sb @ self.zz	3✔
176	self.kss = self.Sk @ self.zz	3✔
177	self.hss = self.Sh @ self.zz	3✔
178
179	def compute_sequence(self, x0, ts_length=None, Pay=None):	3✔
180	"""
181	Simulate quantities and prices for the economy
182
183	Parameters
184	----------
185	x0 : array_like(float)
186	The initial state
187
188	ts_length : scalar(int)
189	Length of the simulation
190
191	Pay : array_like(float)
192	Vector to price an asset whose payout is Pay*xt
193
194	"""
195	lq = LQ(self.Q, self.R, self.A, self.B,	×
196	self.C, N=self.W, beta=self.beta)
197	xp, up, wp = lq.compute_sequence(x0, ts_length)	×
NEW 198	self.h = self.Sh @ xp	×
NEW 199	self.k = self.Sk @ xp	×
NEW 200	self.i = self.Si @ xp	×
NEW 201	self.b = self.Sb @ xp	×
NEW 202	self.d = self.Sd @ xp	×
NEW 203	self.c = self.Sc @ xp	×
NEW 204	self.g = self.Sg @ xp	×
NEW 205	self.s = self.Ss @ xp	×
206
207	# === Value of J-period risk-free bonds === #
208	# === See p.145: Equation (7.11.2) === #
209	e1 = np.zeros((1, self.nc))	×
210	e1[0, 0] = 1	×
211	self.R1_Price = np.empty((ts_length + 1, 1))	×
212	self.R2_Price = np.empty((ts_length + 1, 1))	×
213	self.R5_Price = np.empty((ts_length + 1, 1))	×
214	for i in range(ts_length + 1):	×
NEW 215	self.R1_Price[i, 0] = self.beta * e1 @ self.Mc @ np.linalg.matrix_power(	×
216	self.A0, 1) @ xp[:, i] / (e1 @ self.Mc @ xp[:, i])
NEW 217	self.R2_Price[i, 0] = self.beta*2 (e1 @ self.Mc @	×
218	np.linalg.matrix_power(self.A0, 2) @ xp[:, i]) / (e1 @ self.Mc @ xp[:, i])
NEW 219	self.R5_Price[i, 0] = self.beta*5 (e1 @ self.Mc @	×
220	np.linalg.matrix_power(self.A0, 5) @ xp[:, i]) / (e1 @ self.Mc @ xp[:, i])
221
222	# === Gross rates of return on 1-period risk-free bonds === #
223	self.R1_Gross = 1 / self.R1_Price	×
224
225	# === Net rates of return on J-period risk-free bonds === #
226	# === See p.148: log of gross rate of return, divided by j === #
227	self.R1_Net = np.log(1 / self.R1_Price) / 1	×
228	self.R2_Net = np.log(1 / self.R2_Price) / 2	×
229	self.R5_Net = np.log(1 / self.R5_Price) / 5	×
230
231	# === Value of asset whose payout vector is Pay*xt === #
232	# See p.145: Equation (7.11.1)
233	if isinstance(Pay, np.ndarray) == True:	×
NEW 234	self.Za = Pay.T @ self.Mc	×
235	self.Q = solve_discrete_lyapunov(	×
236	self.A0.T * self.beta**0.5, self.Za)
237	self.q = self.beta / (1 - self.beta) * \	×
238	np.trace(self.C.T @ self.Q @ self.C)
239	self.Pay_Price = np.empty((ts_length + 1, 1))	×
240	self.Pay_Gross = np.empty((ts_length + 1, 1))	×
241	self.Pay_Gross[0, 0] = np.nan	×
242	for i in range(ts_length + 1):	×
NEW 243	self.Pay_Price[i, 0] = (xp[:, i].T @ self.Q @	×
244	xp[:, i] + self.q) / (e1 @ self.Mc @ xp[:, i])
245	for i in range(ts_length):	×
246	self.Pay_Gross[i + 1, 0] = self.Pay_Price[i + 1,	×
247	0] / (self.Pay_Price[i, 0] - Pay @ xp[:, i])
248	return	×
249
250	def irf(self, ts_length=100, shock=None):	3✔
251	"""
252	Create Impulse Response Functions
253
254	Parameters
255	----------
256
257	ts_length : scalar(int)
258	Number of periods to calculate IRF
259
260	Shock : array_like(float)
261	Vector of shocks to calculate IRF to. Default is first element of w
262
263	"""
264
265	if type(shock) != np.ndarray:	×
266	# Default is to select first element of w
267	shock = np.vstack((np.ones((1, 1)), np.zeros((self.nw - 1, 1))))	×
268
269	self.c_irf = np.empty((ts_length, self.nc))	×
270	self.s_irf = np.empty((ts_length, self.nb))	×
271	self.i_irf = np.empty((ts_length, self.ni))	×
272	self.k_irf = np.empty((ts_length, self.nk))	×
273	self.h_irf = np.empty((ts_length, self.nh))	×
274	self.g_irf = np.empty((ts_length, self.ng))	×
275	self.d_irf = np.empty((ts_length, self.nd))	×
276	self.b_irf = np.empty((ts_length, self.nb))	×
277
278	for i in range(ts_length):	×
NEW 279	self.c_irf[i, :] = (self.Sc @	×
280	np.linalg.matrix_power(self.A0, i) @ self.C @ shock).T
NEW 281	self.s_irf[i, :] = (self.Ss @	×
282	np.linalg.matrix_power(self.A0, i) @ self.C @ shock).T
NEW 283	self.i_irf[i, :] = (self.Si @	×
284	np.linalg.matrix_power(self.A0, i) @ self.C @ shock).T
NEW 285	self.k_irf[i, :] = (self.Sk @	×
286	np.linalg.matrix_power(self.A0, i) @ self.C @ shock).T
NEW 287	self.h_irf[i, :] = (self.Sh @	×
288	np.linalg.matrix_power(self.A0, i) @ self.C @ shock).T
NEW 289	self.g_irf[i, :] = (self.Sg @	×
290	np.linalg.matrix_power(self.A0, i) @ self.C @ shock).T
NEW 291	self.d_irf[i, :] = (self.Sd @	×
292	np.linalg.matrix_power(self.A0, i) @ self.C @ shock).T
NEW 293	self.b_irf[i, :] = (self.Sb @	×
294	np.linalg.matrix_power(self.A0, i) @ self.C @ shock).T
295
296	return	×
297
298	def canonical(self):	3✔
299	"""
300	Compute canonical preference representation
301	Uses auxiliary problem of 9.4.2, with the preference shock process reintroduced
302	Calculates pihat, llambdahat and ubhat for the equivalent canonical household technology
303	"""
304	Ac1 = np.hstack((self.deltah, np.zeros((self.nh, self.nz))))	3✔
305	Ac2 = np.hstack((np.zeros((self.nz, self.nh)), self.a22))	3✔
306	Ac = np.vstack((Ac1, Ac2))	3✔
307	Bc = np.vstack((self.thetah, np.zeros((self.nz, self.nc))))	3✔
308	Rc1 = np.hstack((self.llambda.T @ self.llambda, -	3✔
309	self.llambda.T @ self.ub))
310	Rc2 = np.hstack((-self.ub.T @ self.llambda, self.ub.T @ self.ub))	3✔
311	Rc = np.vstack((Rc1, Rc2))	3✔
312	Qc = self.pih.T @ self.pih	3✔
313	Nc = np.hstack(	3✔
314	(self.pih.T @ self.llambda, -self.pih.T @ self.ub))
315
316	lq_aux = LQ(Qc, Rc, Ac, Bc, N=Nc, beta=self.beta)	3✔
317
318	P1, F1, d1 = lq_aux.stationary_values()	3✔
319
320	self.F_b = F1[:, 0:self.nh]	3✔
321	self.F_f = F1[:, self.nh:]	3✔
322
323	self.pihat = np.linalg.cholesky((self.pih.T @	3✔
324	self.pih) + (self.beta @ self.thetah.T @ P1[0:self.nh, 0:self.nh] @ self.thetah)).T
325	self.llambdahat = self.pihat @ self.F_b	3✔
326	self.ubhat = - self.pihat @ self.F_f	3✔
327
328	return	3✔

QuantEcon / QuantEcon.py / 17227928141

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous