diff --git a/examples/paths_and_hist.py b/examples/paths_and_hist.py
index ff7c828c7..b3f5c4bb1 100644
--- a/examples/paths_and_hist.py
+++ b/examples/paths_and_hist.py
@@ -1,7 +1,7 @@
 
 import numpy as np
 import matplotlib.pyplot as plt
-from quantecon import LSS
+from quantecon import LinearStateSpace
 import random
 
 phi_1, phi_2, phi_3, phi_4 = 0.5, -0.2, 0, 0.5
@@ -15,7 +15,7 @@
 G = [1, 0, 0, 0]
 
 T = 30
-ar = LSS(A, C, G, mu_0=np.ones(4))
+ar = LinearStateSpace(A, C, G, mu_0=np.ones(4))
 
 ymin, ymax = -0.8, 1.25
 
diff --git a/examples/paths_and_stationarity.py b/examples/paths_and_stationarity.py
index fd42b9bd5..84f0c6d8d 100644
--- a/examples/paths_and_stationarity.py
+++ b/examples/paths_and_stationarity.py
@@ -1,7 +1,7 @@
 
 import numpy as np
 import matplotlib.pyplot as plt
-from quantecon import LSS
+from quantecon import LinearStateSpace
 import random
 
 phi_1, phi_2, phi_3, phi_4 = 0.5, -0.2, 0, 0.5
@@ -19,7 +19,7 @@
 T2 = 75
 T4 = 100
 
-ar = LSS(A, C, G, mu_0=np.ones(4))
+ar = LinearStateSpace(A, C, G, mu_0=np.ones(4))
 ymin, ymax = -0.8, 1.25
 
 fig, ax = plt.subplots(figsize=(8, 5))
diff --git a/examples/tsh_hg.py b/examples/tsh_hg.py
index 624668283..182e7b101 100644
--- a/examples/tsh_hg.py
+++ b/examples/tsh_hg.py
@@ -2,7 +2,7 @@
 import numpy as np
 import matplotlib.pyplot as plt
 from scipy.stats import norm
-from quantecon import LSS
+from quantecon import LinearStateSpace
 
 phi_1, phi_2, phi_3, phi_4 = 0.5, -0.2, 0, 0.5
 sigma = 0.1
@@ -15,7 +15,7 @@
 G = [1, 0, 0, 0]
 
 T = 30
-ar = LSS(A, C, G)
+ar = LinearStateSpace(A, C, G)
 
 ymin, ymax = -0.8, 1.25
 
diff --git a/quantecon/__init__.py b/quantecon/__init__.py
index 1da79174b..69d37431a 100644
--- a/quantecon/__init__.py
+++ b/quantecon/__init__.py
@@ -15,7 +15,7 @@
 from .arma import ARMA
 from .lqcontrol import LQ
 from .lqnash import nnash
-from .lss import LSS
+from .lss import LinearStateSpace
 from .matrix_eqn import solve_discrete_lyapunov, solve_discrete_riccati
 from .mc_tools import MarkovChain, mc_compute_stationary, mc_sample_path
 from .quadsums import var_quadratic_sum, m_quadratic_sum
diff --git a/quantecon/lss.py b/quantecon/lss.py
index bc150188b..d72056acb 100644
--- a/quantecon/lss.py
+++ b/quantecon/lss.py
@@ -1,16 +1,8 @@
 """
-Authors: Thomas J. Sargent, John Stachurski
-
 Filename: lss.py
+Reference: http://quant-econ.net/py/linear_models.html
 
-Computes quantities associated with the Gaussian linear state space model
-
-    x_{t+1} = A x_t + C w_{t+1}
-
-        y_t = G x_t
-
-The shocks {w_t} are iid and N(0, I)
-
+Computes quantities associated with the Gaussian linear state space model.
 """
 from textwrap import dedent
 import numpy as np
@@ -18,31 +10,29 @@
 from scipy.linalg import solve
 
 
-class LSS(object):
+class LinearStateSpace(object):
     """
     A class that describes a Gaussian linear state space model of the
     form:
 
       x_{t+1} = A x_t + C w_{t+1}
 
-      y_t = G x_t
+      y_t = G x_t + H v_t
 
-    where {w_t} are iid and N(0, I).  If the initial conditions mu_0
-    and Sigma_0 for x_0 ~ N(mu_0, Sigma_0) are not supplied, both
-    are set to zero. When Sigma_0=0, the draw of x_0 is exactly
-    mu_0.
+    where {w_t} and {v_t} are independent and standard normal with dimensions
+    k and l respectively.  The initial conditions are mu_0 and Sigma_0 for x_0
+    ~ N(mu_0, Sigma_0).  When Sigma_0=0, the draw of x_0 is exactly mu_0.
 
     Parameters
     ----------
     A : array_like or scalar(float)
-        This is part of the state transition equation.  It should be
-        `n x n`
+        Part of the state transition equation.  It should be `n x n`
     C : array_like or scalar(float)
-        This is part of the state transition equation.  It should be
-        `n x m`
+        Part of the state transition equation.  It should be `n x m`
     G : array_like or scalar(float)
-        This describes the relation between y_t and x_t and should
-        be `k x n`
+        Part of the observation equation.  It should be `k x n`
+    H : array_like or scalar(float), optional(default=None)
+        Part of the observation equation.  It should be `k x l`
     mu_0 : array_like or scalar(float), optional(default=None)
         This is the mean of initial draw and is `n x 1`
     Sigma_0 : array_like or scalar(float), optional(default=None)
@@ -51,25 +41,30 @@ class LSS(object):
 
     Attributes
     ----------
-    A, C, G, mu_0, Sigma_0 : see Parameters
-    k, n, m : scalar(int)
-        The matrix dimensions
+    A, C, G, H, mu_0, Sigma_0 : see Parameters
+    n, k, m, l : scalar(int)
+        The dimensions of x_t, y_t, w_t and v_t respectively
 
     """
 
-    def __init__(self, A, C, G, mu_0=None, Sigma_0=None):
+    def __init__(self, A, C, G, H=None, mu_0=None, Sigma_0=None):
         self.A, self.G, self.C = list(map(self.convert, (A, G, C)))
         self.k, self.n = self.G.shape
         self.m = self.C.shape[1]
-        # == Default initial conditions == #
+        if H is None:
+            self.H = None
+            self.l = None
+        else:
+            self.H = self.convert(H)
+            self.l = self.H.shape[1]
         if mu_0 is None:
             self.mu_0 = np.zeros((self.n, 1))
         else:
-            self.mu_0 = np.asarray(mu_0)
+            self.mu_0 = self.convert(mu_0)
         if Sigma_0 is None:
             self.Sigma_0 = np.zeros((self.n, self.n))
         else:
-            self.Sigma_0 = Sigma_0
+            self.Sigma_0 = self.convert(Sigma_0)
 
     def __repr__(self):
         return self.__str__()
@@ -116,7 +111,11 @@ def simulate(self, ts_length=100):
         w = np.random.randn(self.m, ts_length-1)
         for t in range(ts_length-1):
             x[:, t+1] = self.A.dot(x[:, t]) + self.C.dot(w[:, t])
-        y = self.G.dot(x)
+        if self.H is not None:
+            v = np.random.randn(self.l, ts_length)
+            y = self.G.dot(x) + self.H.dot(v)
+        else:
+            y = self.G.dot(x)
 
         return x, y
 
@@ -147,7 +146,11 @@ def replicate(self, T=10, num_reps=100):
         for j in range(num_reps):
             x_T, _ = self.simulate(ts_length=T+1)
             x[:, j] = x_T[:, -1]
-        y = self.G.dot(x)
+        if self.H is not None:
+            v = np.random.randn(self.l, num_reps)
+            y = self.G.dot(x) + self.H.dot(v)
+        else:
+            y = self.G.dot(x)
 
         return x, y
 
@@ -174,12 +177,19 @@ def moment_sequence(self):
 
         """
         # == Simplify names == #
-        A, C, G = self.A, self.C, self.G
+        A, C, G, H = self.A, self.C, self.G, self.H
         # == Initial moments == #
         mu_x, Sigma_x = self.mu_0, self.Sigma_0
+
         while 1:
-            mu_y, Sigma_y = G.dot(mu_x), G.dot(Sigma_x).dot(G.T)
+            mu_y = G.dot(mu_x)
+            if H is None:
+                Sigma_y = G.dot(Sigma_x).dot(G.T)
+            else:
+                Sigma_y = G.dot(Sigma_x).dot(G.T) + H.dot(H.T)
+
             yield mu_x, mu_y, Sigma_x, Sigma_y
+
             # == Update moments of x == #
             mu_x = A.dot(mu_x)
             Sigma_x = A.dot(Sigma_x).dot(A.T) + C.dot(C.T)
@@ -216,7 +226,8 @@ def stationary_distributions(self, max_iter=200, tol=1e-5):
         mu_x, mu_y, Sigma_x, Sigma_y = next(m)
         i = 0
         error = tol + 1
-        # == Loop until convergence or failuer == #
+
+        # == Loop until convergence or failure == #
         while error > tol:
 
             if i > max_iter:
diff --git a/quantecon/tests/test_lss.py b/quantecon/tests/test_lss.py
index 0509b0a17..9c7301e16 100644
--- a/quantecon/tests/test_lss.py
+++ b/quantecon/tests/test_lss.py
@@ -11,7 +11,7 @@
 import unittest
 import numpy as np
 from numpy.testing import assert_allclose
-from quantecon.lss import LSS
+from quantecon.lss import LinearStateSpace
 
 
 class TestLinearStateSpace(unittest.TestCase):
@@ -23,7 +23,7 @@ def setUp(self):
         G = 1.
         mu_0 = .75
 
-        self.ss = LSS(A, C, G, mu_0)
+        self.ss = LinearStateSpace(A, C, G, mu_0=mu_0)
 
     def tearDown(self):
         del self.ss
diff --git a/solutions/lss_solutions.ipynb b/solutions/lss_solutions.ipynb
index 462c722a5..da709b58c 100644
--- a/solutions/lss_solutions.ipynb
+++ b/solutions/lss_solutions.ipynb
@@ -1,7 +1,7 @@
 {
  "metadata": {
   "name": "",
-  "signature": "sha256:74d6017db9fff83f5e920a3f6f30592ea24c605b29c54965059aeefcb540a595"
+  "signature": "sha256:b081069f474bcd3cae9777f77b493b432b34a77e8ae68c38b6293505b6d38bba"
  },
  "nbformat": 3,
  "nbformat_minor": 0,
@@ -50,7 +50,7 @@
      "input": [
       "import numpy as np\n",
       "import matplotlib.pyplot as plt\n",
-      "from quantecon import LSS"
+      "from quantecon import LinearStateSpace"
      ],
      "language": "python",
      "metadata": {},
@@ -77,7 +77,7 @@
       "C = np.zeros((3, 1))\n",
       "G = [0, 1, 0]\n",
       "\n",
-      "ar = LSS(A, C, G, mu_0=np.ones(3))\n",
+      "ar = LinearStateSpace(A, C, G, mu_0=np.ones(3))\n",
       "x, y = ar.simulate(ts_length=50)\n",
       "\n",
       "fig, ax = plt.subplots(figsize=(8, 4.6))\n",
@@ -127,7 +127,7 @@
       "     [0]]\n",
       "G = [1, 0, 0, 0]\n",
       "\n",
-      "ar = LSS(A, C, G, mu_0=np.ones(4))\n",
+      "ar = LinearStateSpace(A, C, G, mu_0=np.ones(4))\n",
       "x, y = ar.simulate(ts_length=200)\n",
       "\n",
       "fig, ax = plt.subplots(figsize=(8, 4.6))\n",
@@ -183,7 +183,7 @@
       "\n",
       "I = 20\n",
       "T = 50\n",
-      "ar = LSS(A, C, G, mu_0=np.ones(4))\n",
+      "ar = LinearStateSpace(A, C, G, mu_0=np.ones(4))\n",
       "ymin, ymax = -0.5, 1.15\n",
       "\n",
       "fig, ax = plt.subplots(figsize=(8, 5))\n",
@@ -255,7 +255,7 @@
       "T2 = 75\n",
       "T4 = 100\n",
       "\n",
-      "ar = LSS(A, C, G, mu_0=np.ones(4))\n",
+      "ar = LinearStateSpace(A, C, G, mu_0=np.ones(4))\n",
       "ymin, ymax = -0.6, 0.6\n",
       "\n",
       "fig, ax = plt.subplots(figsize=(8, 5))\n",