From 81346f8729728d6289c030303ea77cae56fdf1b5 Mon Sep 17 00:00:00 2001
From: Matt McKay <mamckay@gmail.com>
Date: Mon, 13 Aug 2018 10:24:29 +1000
Subject: [PATCH 01/15] initial commit of dle code from Sebastian Graves

---
 quantecon/dle.py | 237 +++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 237 insertions(+)
 create mode 100644 quantecon/dle.py

diff --git a/quantecon/dle.py b/quantecon/dle.py
new file mode 100644
index 000000000..97245f977
--- /dev/null
+++ b/quantecon/dle.py
@@ -0,0 +1,237 @@
+"""
+Filename: DynLinEcon.py
+
+Author: Sebastian Graves
+
+Provides a class called DLE to convert and solve dynamic linear economics (as set out in Hansen & Sargent (2013)) as LQ problems.
+"""
+
+import numpy as np
+from quantecon import LQ
+from quantecon import solve_discrete_lyapunov
+from sympy import Matrix
+from pylab import array
+
+
+class DLE(object):
+    """
+    This class is for analysing economies that are characterised by matrices defining information and shocks {a22, c2, ub, ud}, 
+    production technology {phic, phig, phii, gamma, deltak, thetak},
+    household technology {deltah, thetah, lambda, pih},
+    and the scalar {beta}.
+
+    Section 5.5 of HS2013 describes how to map these matrices into those of a LQ problem. 
+    
+    This class solves the associated LQ problem for such economies, and provides methods to: 
+    find the non-stochastic steady state,
+    simulate quantities and prices,
+    find impulse response functions,
+    and create canonical preference representations.
+    """
+
+    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
+        junk,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.dot(phiin)
+        phiing = ug.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))))
+
+        # === 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.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(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.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.dot(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.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)
+
+        # ===  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*np.asscalar(self.beta)*(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.Mh = -2*np.asscalar(self.beta)*(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.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))
+        #Note to check: I had to put a minus in front of everything to get same as the Matlab code
+        
+    def compute_steadystate(self, nnc=2): #nnc is the position of the constant in the state vector
+        """ Find non-stochastic steady-state of the economy """
+        zx = Matrix(np.eye(self.A0.shape[0])-self.A0)
+        self.zz = zx.nullspace()
+        self.zz = array(self.zz)
+        self.zz = self.zz.T
+        self.zz = zz = self.zz/self.zz[nnc] 
+        self.css = self.Sc.dot(self.zz).astype(float)
+        self.sss = self.Ss.dot(self.zz).astype(float)
+        self.iss = self.Si.dot(self.zz).astype(float)
+        self.dss = self.Sd.dot(self.zz).astype(float)
+        self.bss = self.Sb.dot(self.zz).astype(float)
+        self.kss = self.Sk.dot(self.zz).astype(float)
+        self.hss = self.Sh.dot(self.zz).astype(float)
+
+    def compute_sequence(self, x0, ts_length=None, Pay=None):
+        """ Simulate quantities and prices for our economy """
+        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)
+
+        # === 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.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])
+
+        # === 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.dot(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))
+            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])
+            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]))
+        return
+
+    def irf(self, ts_length=100, shock=None):
+        """ Create Impulse Response Functions """
+        
+        if type(shock) != np.ndarray:
+            shock = np.vstack((np.ones((1,1)),np.zeros((self.nw-1,1)))) #Default is to select first element of w
+
+        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.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
+
+        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))))
+        Cc = np.vstack((np.zeros((self.nh,self.nw)),self.c2))
+        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)))
+        Rc = np.vstack((Rc1,Rc2))
+        Qc = self.pih.T.dot(self.pih)
+        Nc = np.hstack((self.pih.T.dot(self.llambda),-self.pih.T.dot(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.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)
+        
+        return 
\ No newline at end of file

From b37a902b18d3c99c170de676b5cd856c7a57e026 Mon Sep 17 00:00:00 2001
From: Matt McKay <mamckay@gmail.com>
Date: Mon, 13 Aug 2018 10:41:07 +1000
Subject: [PATCH 02/15] adjustments for package

---
 quantecon/__init__.py |  1 +
 quantecon/dle.py      | 15 ++++++---------
 2 files changed, 7 insertions(+), 9 deletions(-)

diff --git a/quantecon/__init__.py b/quantecon/__init__.py
index 1428a7f35..728a879dd 100644
--- a/quantecon/__init__.py
+++ b/quantecon/__init__.py
@@ -17,6 +17,7 @@
 #-Objects-#
 from .compute_fp import compute_fixed_point
 from .discrete_rv import DiscreteRV
+from .dle import DLE
 from .ecdf import ECDF
 from .estspec import smooth, periodogram, ar_periodogram
 # from .game_theory import <objects-here> 							#Place Holder if we wish to promote any general objects to the qe namespace.
diff --git a/quantecon/dle.py b/quantecon/dle.py
index 97245f977..a0aa4fb1b 100644
--- a/quantecon/dle.py
+++ b/quantecon/dle.py
@@ -1,21 +1,18 @@
 """
-Filename: DynLinEcon.py
-
-Author: Sebastian Graves
-
-Provides a class called DLE to convert and solve dynamic linear economics (as set out in Hansen & Sargent (2013)) as LQ problems.
+Provides a class called DLE to convert and solve dynamic linear economics 
+(as set out in Hansen & Sargent (2013)) as LQ problems.
 """
 
 import numpy as np
-from quantecon import LQ
-from quantecon import solve_discrete_lyapunov
 from sympy import Matrix
 from pylab import array
-
+from .lqcontrol import LQ
+from .matrix_eqn import solve_discrete_lyapunov
 
 class DLE(object):
     """
-    This class is for analysing economies that are characterised by matrices defining information and shocks {a22, c2, ub, ud}, 
+    This class is for analysing economies that are characterised by matrices 
+    defining information and shocks {a22, c2, ub, ud}, 
     production technology {phic, phig, phii, gamma, deltak, thetak},
     household technology {deltah, thetah, lambda, pih},
     and the scalar {beta}.

From a314e9b0e242ef1e9bf66e97d01595de93333f6a Mon Sep 17 00:00:00 2001
From: Matt McKay <mamckay@gmail.com>
Date: Thu, 23 Aug 2018 12:36:26 +1000
Subject: [PATCH 03/15] update code with docstrings

---
 quantecon/dle.py | 326 ++++++++++++++++++++++++++++++-----------------
 1 file changed, 210 insertions(+), 116 deletions(-)

diff --git a/quantecon/dle.py b/quantecon/dle.py
index a0aa4fb1b..b78ad1fd6 100644
--- a/quantecon/dle.py
+++ b/quantecon/dle.py
@@ -9,21 +9,49 @@
 from .lqcontrol import LQ
 from .matrix_eqn import solve_discrete_lyapunov
 
+
 class DLE(object):
-    """
-    This class is for analysing economies that are characterised by matrices 
-    defining information and shocks {a22, c2, ub, ud}, 
-    production technology {phic, phig, phii, gamma, deltak, thetak},
-    household technology {deltah, thetah, lambda, pih},
-    and the scalar {beta}.
+    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. 
-    
-    This class solves the associated LQ problem for such economies, and provides methods to: 
-    find the non-stochastic steady state,
-    simulate quantities and prices,
-    find impulse response functions,
-    and create canonical preference representations.
+
+    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 matrices \Delta_h, \Theta_h, \Lambda, \Pi, and the scalar \beta
+
     """
 
     def __init__(self, information, technology, preferences):
@@ -32,54 +60,60 @@ def __init__(self, information, technology, 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
-        junk,self.ng = self.phig.shape
-        self.nk,self.ni = self.thetak.shape
+        self.nb, self.nh = self.llambda.shape
+        self.nd, self.nc = self.phic.shape
+        self.nz, self.nw = self.c2.shape
+        junk, 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)))
+        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)
         phiing = ug.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))))
+        a12 = np.vstack((self.thetah.dot(phiinc).dot(
+            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))))
+        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))))
+        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))
+        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.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(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.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(
+            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.nx = self.nh+self.nk+self.nz
-        self.n = self.ni+self.nh+self.nk+self.nz
+        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]
+        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)
+        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()
 
@@ -87,38 +121,52 @@ def __init__(self, information, technology, preferences):
 
         self.A0 = self.A - self.B.dot(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.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.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.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)
 
         # ===  Calculate eigenvalues of A0 === #
-        self.A110 = self.A0[0:self.nh+self.nk,0:self.nh+self.nk]
+        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*np.asscalar(self.beta)*(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.Mh = -2*np.asscalar(self.beta)*(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.Mk = -2 * np.asscalar(self.beta) * (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.Mh = -2 * np.asscalar(self.beta) * (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.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.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))
-        #Note to check: I had to put a minus in front of everything to get same as the Matlab code
-        
-    def compute_steadystate(self, nnc=2): #nnc is the position of the constant in the state vector
-        """ Find non-stochastic steady-state of the economy """
-        zx = Matrix(np.eye(self.A0.shape[0])-self.A0)
+
+    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 = Matrix(np.eye(self.A0.shape[0]) - self.A0)
         self.zz = zx.nullspace()
         self.zz = array(self.zz)
         self.zz = self.zz.T
-        self.zz = zz = self.zz/self.zz[nnc] 
+        self.zz = zz = self.zz / self.zz[nnc]
         self.css = self.Sc.dot(self.zz).astype(float)
         self.sss = self.Ss.dot(self.zz).astype(float)
         self.iss = self.Si.dot(self.zz).astype(float)
@@ -128,8 +176,23 @@ def compute_steadystate(self, nnc=2): #nnc is the position of the constant in th
         self.hss = self.Sh.dot(self.zz).astype(float)
 
     def compute_sequence(self, x0, ts_length=None, Pay=None):
-        """ Simulate quantities and prices for our economy """
-        lq = LQ(self.Q, self.R, self.A, self.B, self.C, N=self.W, beta=self.beta)
+        """
+        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.dot(xp)
         self.k = self.Sk.dot(xp)
@@ -142,64 +205,92 @@ def compute_sequence(self, x0, ts_length=None, Pay=None):
 
         # === 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.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])
+        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.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])
 
         # === Gross rates of return on 1-period risk-free bonds === #
-        self.R1_Gross = 1/self.R1_Price
+        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
+        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:
+        if isinstance(Pay, np.ndarray) == True:
             self.Za = Pay.T.dot(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))
-            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.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))
+            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])
             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]))
+                self.Pay_Gross[i + 1, 0] = self.Pay_Price[i + 1,
+                                                          0] / (self.Pay_Price[i, 0] - Pay.dot(xp[:, i]))
         return
 
     def irf(self, ts_length=100, shock=None):
-        """ Create Impulse Response Functions """
-        
-        if type(shock) != np.ndarray:
-            shock = np.vstack((np.ones((1,1)),np.zeros((self.nw-1,1)))) #Default is to select first element of w
+        """
+        Create Impulse Response Functions
+
+        Parameters
+        ----------
+
+        ts_length : scalar(int)
+            Number of periods to calculate IRF
 
-        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))
+        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.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.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
 
         return
 
@@ -209,26 +300,29 @@ def canonical(self):
         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))))
-        Cc = np.vstack((np.zeros((self.nh,self.nw)),self.c2))
-        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)))
-        Rc = np.vstack((Rc1,Rc2))
+        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))))
+        Cc = np.vstack((np.zeros((self.nh, self.nw)), self.c2))
+        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)))
+        Rc = np.vstack((Rc1, Rc2))
         Qc = self.pih.T.dot(self.pih)
-        Nc = np.hstack((self.pih.T.dot(self.llambda),-self.pih.T.dot(self.ub)))
-        
-        lq_aux = LQ(Qc, Rc, Ac, Bc, N = Nc, beta=self.beta)
+        Nc = np.hstack(
+            (self.pih.T.dot(self.llambda), -self.pih.T.dot(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.dot(self.pih) + self.beta.dot(self.thetah.T).dot(P1[0:self.nh,0:self.nh]).dot(self.thetah)).T
+
+        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)
-        
-        return 
\ No newline at end of file
+
+        return

From 49674074d7f87718c7042f9c46d432b2ae574d92 Mon Sep 17 00:00:00 2001
From: Matt McKay <mamckay@gmail.com>
Date: Thu, 23 Aug 2018 13:04:18 +1000
Subject: [PATCH 04/15] remove pylab imports

---
 quantecon/dle.py | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/quantecon/dle.py b/quantecon/dle.py
index b78ad1fd6..be4c3b20b 100644
--- a/quantecon/dle.py
+++ b/quantecon/dle.py
@@ -4,8 +4,6 @@
 """
 
 import numpy as np
-from sympy import Matrix
-from pylab import array
 from .lqcontrol import LQ
 from .matrix_eqn import solve_discrete_lyapunov
 
@@ -164,7 +162,7 @@ def compute_steadystate(self, nnc=2):
         """
         zx = Matrix(np.eye(self.A0.shape[0]) - self.A0)
         self.zz = zx.nullspace()
-        self.zz = array(self.zz)
+        self.zz = np.array(self.zz)
         self.zz = self.zz.T
         self.zz = zz = self.zz / self.zz[nnc]
         self.css = self.Sc.dot(self.zz).astype(float)

From 3fc1f95852ab66f743e3336dae0b28054b62b791 Mon Sep 17 00:00:00 2001
From: Matt McKay <mamckay@gmail.com>
Date: Thu, 23 Aug 2018 13:13:12 +1000
Subject: [PATCH 05/15] move from pylab Matrix to np.matrix

---
 quantecon/dle.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/quantecon/dle.py b/quantecon/dle.py
index be4c3b20b..f3fb22451 100644
--- a/quantecon/dle.py
+++ b/quantecon/dle.py
@@ -160,7 +160,7 @@ def compute_steadystate(self, nnc=2):
             nnc is the location of the constant in the state vector x_t
 
         """
-        zx = Matrix(np.eye(self.A0.shape[0]) - self.A0)
+        zx =np.matrix(np.eye(self.A0.shape[0]) - self.A0)
         self.zz = zx.nullspace()
         self.zz = np.array(self.zz)
         self.zz = self.zz.T

From eff4806723bb433fed379dced847dc8bb95851ff Mon Sep 17 00:00:00 2001
From: Matt McKay <mamckay@gmail.com>
Date: Thu, 23 Aug 2018 14:05:54 +1000
Subject: [PATCH 06/15] update to use nullspace from quantecon

---
 quantecon/dle.py | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/quantecon/dle.py b/quantecon/dle.py
index f3fb22451..8635625d7 100644
--- a/quantecon/dle.py
+++ b/quantecon/dle.py
@@ -6,6 +6,7 @@
 import numpy as np
 from .lqcontrol import LQ
 from .matrix_eqn import solve_discrete_lyapunov
+from .rank_nullspace import nullspace
 
 
 class DLE(object):
@@ -160,8 +161,9 @@ def compute_steadystate(self, nnc=2):
             nnc is the location of the constant in the state vector x_t
 
         """
-        zx =np.matrix(np.eye(self.A0.shape[0]) - self.A0)
-        self.zz = zx.nullspace()
+        # zx =Matrix(np.eye(self.A0.shape[0]) - self.A0)  #remove once confirmed can use nullspace from QE
+        # self.zz = zx.nullspace()
+        self.zz = nullspace(A0)
         self.zz = np.array(self.zz)
         self.zz = self.zz.T
         self.zz = zz = self.zz / self.zz[nnc]

From b62fe6f2652424acd2eea971c56db995fbb25eeb Mon Sep 17 00:00:00 2001
From: Matt McKay <mamckay@gmail.com>
Date: Thu, 23 Aug 2018 14:43:26 +1000
Subject: [PATCH 07/15] correct algebra for use of quantecon nullspace method

---
 quantecon/dle.py | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/quantecon/dle.py b/quantecon/dle.py
index 8635625d7..81e3218f2 100644
--- a/quantecon/dle.py
+++ b/quantecon/dle.py
@@ -161,9 +161,8 @@ def compute_steadystate(self, nnc=2):
             nnc is the location of the constant in the state vector x_t
 
         """
-        # zx =Matrix(np.eye(self.A0.shape[0]) - self.A0)  #remove once confirmed can use nullspace from QE
-        # self.zz = zx.nullspace()
-        self.zz = nullspace(A0)
+        zx = np.eye(self.A0.shape[0]) - self.A0
+        self.zz = nullspace(zx)
         self.zz = np.array(self.zz)
         self.zz = self.zz.T
         self.zz = zz = self.zz / self.zz[nnc]

From a65a3039704d4da26a68f1f6b55cbc6c7def3ba2 Mon Sep 17 00:00:00 2001
From: Matt McKay <mamckay@gmail.com>
Date: Fri, 7 Sep 2018 16:04:36 +1000
Subject: [PATCH 08/15] add test for transformation of problem to LQ type
 problem using known solutions verified by matlab programs

---
 quantecon/tests/test_dle.py | 93 +++++++++++++++++++++++++++++++++++++
 1 file changed, 93 insertions(+)
 create mode 100644 quantecon/tests/test_dle.py

diff --git a/quantecon/tests/test_dle.py b/quantecon/tests/test_dle.py
new file mode 100644
index 000000000..e6ddf12eb
--- /dev/null
+++ b/quantecon/tests/test_dle.py
@@ -0,0 +1,93 @@
+"""
+Tests for dle.py file
+"""
+
+import sys
+import os
+import unittest
+import numpy as np
+from scipy.linalg import LinAlgError
+from numpy.testing import assert_allclose
+from numpy import dot
+from quantecon.dle import DLE
+
+
+class TestDLE(unittest.TestCase):
+
+    def setUp(self):
+        """
+        Given LQ control is tested we will test the transformation
+        to alter the problem into a form suitable to solve using LQ
+        """
+        # Initial Values
+        gam = 0
+        gamma = np.array([[gam], [0]])
+        phic = np.array([[1], [0]])
+        phig = np.array([[0], [1]])
+        phi1 = 1e-4
+        phii = np.array([[0], [-phi1]])
+        deltak = np.array([[.95]])
+        thetak = np.array([[1]])
+        beta = np.array([[1 / 1.05]])
+        ud = np.array([[5, 1, 0], [0, 0, 0]])
+        a22 = np.array([[1, 0, 0], [0, 0.8, 0], [0, 0, 0.5]])
+        c2 = np.array([[0, 1, 0], [0, 0, 1]]).T
+        llambda = np.array([[0]])
+        pih = np.array([[1]])
+        deltah = np.array([[.9]])
+        thetah = np.array([[1]]) - deltah
+        ub = np.array([[30, 0, 0]])
+        x0 = np.array([[5, 150,1,0,0]]).T
+
+        information = (a22, c2, ub, ud)
+        technology = (phic, phig, phii, gamma, deltak, thetak)
+        preferences = (beta, llambda, pih, deltah, thetah)        
+
+        self.dle = DLE(information, technology, preferences)
+
+    def tearDown(self):
+        del self.dle
+
+    def test_transformation_Q(self):
+        Q_solution = np.array([[5.e-09]])
+        assert_allclose(Q_solution, self.dle.Q)
+
+    def test_transformation_R(self):
+        R_solution = np.array([[0.,   0.,   0.,   0.,   0.],
+                               [0.,   0.,   0.,   0.,   0.],
+                               [0.,   0., 312.5, -12.5,   0.],
+                               [0.,   0., -12.5,   0.5,   0.],
+                               [0.,   0.,   0.,   0.,   0.]])
+        assert_allclose(R_solution, self.dle.R)
+
+    def test_transformation_A(self):
+        A_solution = np.array([[0.9, 0., 0.5, 0.1, 0.],
+                               [0., 0.95, 0., 0., 0.],
+                               [0., 0., 1., 0., 0.],
+                               [0., 0., 0., 0.8, 0.],
+                               [0., 0., 0., 0., 0.5]])
+        assert_allclose(A_solution, self.dle.A)
+
+    def test_transformation_B(self):
+        B_solution = np.array([[-0.],
+                               [1.],
+                               [0.],
+                               [0.],
+                               [0.]])
+        assert_allclose(B_solution, self.dle.B)
+
+    def test_transformation_C(self):
+        C_solution = np.array([[0., 0.],
+                               [0., 0.],
+                               [0., 0.],
+                               [1., 0.],
+                               [0., 1.]])
+        assert_allclose(C_solution, self.dle.C)
+
+    def test_transformation_W(self):
+        W_solution = np.array([[0., 0., 0., 0., 0.]])
+        assert_allclose(W_solution, self.dle.W)
+
+if __name__ == '__main__':
+    suite = unittest.TestLoader().loadTestsFromTestCase(TestDLE)
+    unittest.TextTestRunner(verbosity=2, stream=sys.stderr).run(suite)

From 9cccc38b505009539637ce2b1878b04f198e76fa Mon Sep 17 00:00:00 2001
From: Matt McKay <mamckay@gmail.com>
Date: Fri, 7 Sep 2018 16:22:38 +1000
Subject: [PATCH 09/15] reduce line lengths in docstrings

---
 quantecon/dle.py | 20 +++++++++++++-------
 1 file changed, 13 insertions(+), 7 deletions(-)

diff --git a/quantecon/dle.py b/quantecon/dle.py
index 81e3218f2..fe2217c28 100644
--- a/quantecon/dle.py
+++ b/quantecon/dle.py
@@ -31,25 +31,31 @@ class DLE(object):
 
     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. 
+    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
+    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
+    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.
+    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
+        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 matrices \Delta_h, \Theta_h, \Lambda, \Pi, and the scalar \beta
+        Preferences is a tuple containing the matrices \Delta_h, \Theta_h, \Lambda, 
+        \Pi, and the scalar \beta
 
     """
 

From e0ac1a917b386ac5feef7690d7ae68b9b690d53b Mon Sep 17 00:00:00 2001
From: Matt McKay <mamckay@gmail.com>
Date: Fri, 7 Sep 2018 16:30:31 +1000
Subject: [PATCH 10/15] revert use of quantecon nullspace as representations
 are not exactly the same ... need to investigate

---
 quantecon/dle.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/quantecon/dle.py b/quantecon/dle.py
index fe2217c28..22f1e9aa7 100644
--- a/quantecon/dle.py
+++ b/quantecon/dle.py
@@ -7,7 +7,7 @@
 from .lqcontrol import LQ
 from .matrix_eqn import solve_discrete_lyapunov
 from .rank_nullspace import nullspace
-
+from sympy import Matrix
 
 class DLE(object):
     r"""
@@ -167,8 +167,8 @@ def compute_steadystate(self, nnc=2):
             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)
+        zx = Matrix(np.eye(self.A0.shape[0])-self.A0)
+        self.zz = zx.nullspace()
         self.zz = np.array(self.zz)
         self.zz = self.zz.T
         self.zz = zz = self.zz / self.zz[nnc]

From 4f30abd8053ae4a22a5b505cc1401af60c226d6a Mon Sep 17 00:00:00 2001
From: Matt McKay <mamckay@gmail.com>
Date: Tue, 11 Sep 2018 10:45:22 +1000
Subject: [PATCH 11/15] migrate to using quantecon nullspace method

---
 quantecon/dle.py            | 23 ++++++++++-------------
 quantecon/tests/test_dle.py | 13 +++++++++++++
 2 files changed, 23 insertions(+), 13 deletions(-)

diff --git a/quantecon/dle.py b/quantecon/dle.py
index 22f1e9aa7..f2f7c6b7c 100644
--- a/quantecon/dle.py
+++ b/quantecon/dle.py
@@ -7,7 +7,6 @@
 from .lqcontrol import LQ
 from .matrix_eqn import solve_discrete_lyapunov
 from .rank_nullspace import nullspace
-from sympy import Matrix
 
 class DLE(object):
     r"""
@@ -167,18 +166,16 @@ def compute_steadystate(self, nnc=2):
             nnc is the location of the constant in the state vector x_t
 
         """
-        zx = Matrix(np.eye(self.A0.shape[0])-self.A0)
-        self.zz = zx.nullspace()
-        self.zz = np.array(self.zz)
-        self.zz = self.zz.T
-        self.zz = zz = self.zz / self.zz[nnc]
-        self.css = self.Sc.dot(self.zz).astype(float)
-        self.sss = self.Ss.dot(self.zz).astype(float)
-        self.iss = self.Si.dot(self.zz).astype(float)
-        self.dss = self.Sd.dot(self.zz).astype(float)
-        self.bss = self.Sb.dot(self.zz).astype(float)
-        self.kss = self.Sk.dot(self.zz).astype(float)
-        self.hss = self.Sh.dot(self.zz).astype(float)
+        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)
 
     def compute_sequence(self, x0, ts_length=None, Pay=None):
         """
diff --git a/quantecon/tests/test_dle.py b/quantecon/tests/test_dle.py
index e6ddf12eb..33ec8c6ba 100644
--- a/quantecon/tests/test_dle.py
+++ b/quantecon/tests/test_dle.py
@@ -88,6 +88,19 @@ def test_transformation_W(self):
         W_solution = np.array([[0., 0., 0., 0., 0.]])
         assert_allclose(W_solution, self.dle.W)
 
+    def test_compute_steadystate(self):
+        solutions = {
+            'css' : np.array([[5.]]),
+            'sss' : np.array([[5.]]),
+            'iss' : np.array([[0.]]),
+            'dss' : np.array([[5.], [0.]]),
+            'bss' : np.array([[30.]]),
+            'kss' : np.array([[0.]]),
+            'hss' : np.array([[5.]]),
+        }
+        for item in solutions.keys():
+            assert_allclose(self.dle.__dict__[item], solutions[item])
+
 if __name__ == '__main__':
     suite = unittest.TestLoader().loadTestsFromTestCase(TestDLE)
     unittest.TextTestRunner(verbosity=2, stream=sys.stderr).run(suite)

From 584fb22a1d254bb5dd2ce5509c39160c9bdf8f84 Mon Sep 17 00:00:00 2001
From: Matt McKay <mamckay@gmail.com>
Date: Tue, 11 Sep 2018 10:59:13 +1000
Subject: [PATCH 12/15] fix test for steadystate

---
 quantecon/tests/test_dle.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/quantecon/tests/test_dle.py b/quantecon/tests/test_dle.py
index 33ec8c6ba..8efc2ba05 100644
--- a/quantecon/tests/test_dle.py
+++ b/quantecon/tests/test_dle.py
@@ -98,6 +98,7 @@ def test_compute_steadystate(self):
             'kss' : np.array([[0.]]),
             'hss' : np.array([[5.]]),
         }
+        self.dle.compute_steadystate()
         for item in solutions.keys():
             assert_allclose(self.dle.__dict__[item], solutions[item])
 

From 89c98a55fda7548412736f27b927283632ccb4cf Mon Sep 17 00:00:00 2001
From: Matt McKay <mamckay@gmail.com>
Date: Tue, 11 Sep 2018 11:28:33 +1000
Subject: [PATCH 13/15] add atol for comparing small values

---
 quantecon/tests/test_dle.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/quantecon/tests/test_dle.py b/quantecon/tests/test_dle.py
index 8efc2ba05..8451e7c75 100644
--- a/quantecon/tests/test_dle.py
+++ b/quantecon/tests/test_dle.py
@@ -11,6 +11,7 @@
 from numpy import dot
 from quantecon.dle import DLE
 
+ATOL = 1e-10
 
 class TestDLE(unittest.TestCase):
 
@@ -100,7 +101,7 @@ def test_compute_steadystate(self):
         }
         self.dle.compute_steadystate()
         for item in solutions.keys():
-            assert_allclose(self.dle.__dict__[item], solutions[item])
+            assert_allclose(self.dle.__dict__[item], solutions[item], atol=ATOL)
 
 if __name__ == '__main__':
     suite = unittest.TestLoader().loadTestsFromTestCase(TestDLE)

From 97749df49156b652d12104b030e7ec7cc327d123 Mon Sep 17 00:00:00 2001
From: Matt McKay <mamckay@gmail.com>
Date: Tue, 11 Sep 2018 17:32:38 +1000
Subject: [PATCH 14/15] add test for canonical

---
 quantecon/tests/test_dle.py | 11 +++++++++++
 1 file changed, 11 insertions(+)

diff --git a/quantecon/tests/test_dle.py b/quantecon/tests/test_dle.py
index 8451e7c75..526bd1120 100644
--- a/quantecon/tests/test_dle.py
+++ b/quantecon/tests/test_dle.py
@@ -103,6 +103,17 @@ def test_compute_steadystate(self):
         for item in solutions.keys():
             assert_allclose(self.dle.__dict__[item], solutions[item], atol=ATOL)
 
+    def test_canonical(self):
+        solutions = {
+            'pihat': np.array([[1.]]),
+            'llambdahat': np.array([[-1.48690584e-19]]),
+            'ubhat': np.array([[30., -0., -0.]])
+        }
+        self.dle.canonical()
+        for item in solutions.keys():
+            assert_allclose(self.dle.__dict__[item], solutions[item], atol=ATOL)
+
+
 if __name__ == '__main__':
     suite = unittest.TestLoader().loadTestsFromTestCase(TestDLE)
     unittest.TextTestRunner(verbosity=2, stream=sys.stderr).run(suite)

From 54be06aa0a365cb65968565b00e17fe04f2a8812 Mon Sep 17 00:00:00 2001
From: Matt McKay <mamckay@gmail.com>
Date: Tue, 11 Sep 2018 17:40:58 +1000
Subject: [PATCH 15/15] add dle to documentation

---
 docs/source/tools.rst     | 1 +
 docs/source/tools/dle.rst | 7 +++++++
 2 files changed, 8 insertions(+)
 create mode 100644 docs/source/tools/dle.rst

diff --git a/docs/source/tools.rst b/docs/source/tools.rst
index 7b8456710..ad6f60592 100644
--- a/docs/source/tools.rst
+++ b/docs/source/tools.rst
@@ -9,6 +9,7 @@ Tools
    tools/compute_fp
    tools/discrete_rv
    tools/distributions
+   tools/dle
    tools/ecdf
    tools/estspec
    tools/filter
diff --git a/docs/source/tools/dle.rst b/docs/source/tools/dle.rst
new file mode 100644
index 000000000..4e1ecc6e0
--- /dev/null
+++ b/docs/source/tools/dle.rst
@@ -0,0 +1,7 @@
+dle
+===
+
+.. automodule:: quantecon.dle
+    :members:
+    :undoc-members:
+    :show-inheritance: