diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index a4cc2c646..64de906d0 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -31,7 +31,7 @@ jobs:
       - name: Build HTML
         shell: bash -l {0}
         run: |
-          jb build lectures --path-output ./
+          jb build -W --keep-going lectures --path-output ./
       - name: Save Build as Artifact
         uses: actions/upload-artifact@v1
         with:
diff --git a/environment.yml b/environment.yml
index b2a7a09ed..38cdc94bd 100644
--- a/environment.yml
+++ b/environment.yml
@@ -8,7 +8,7 @@ dependencies:
   - pip:
     - jupyter-book
     - sphinx-multitoc-numbering
-    - git+https://github.com/quantecon/quantecon-book-theme
+    - quantecon-book-theme
     - sphinx-tojupyter
     - sphinxext-rediraffe
     - sphinx-exercise
diff --git a/lectures/linear_models.md b/lectures/linear_models.md
index aeaeb9aa9..cc48009d2 100644
--- a/lectures/linear_models.md
+++ b/lectures/linear_models.md
@@ -734,7 +734,7 @@ sample_size = 500_000
 ar = LinearStateSpace(A_2, C_2, G_2, mu_0=np.ones(4))
 fig, ax = plt.subplots()
 x, y = ar.simulate(sample_size)
-mu_x, mu_y, Sigma_x, Sigma_y = ar.stationary_distributions()
+mu_x, mu_y, Sigma_x, Sigma_y, Sigma_yx = ar.stationary_distributions()
 f_y = norm(loc=float(mu_y), scale=float(np.sqrt(Sigma_y)))
 y = y.flatten()
 ygrid = np.linspace(ymin, ymax, 150)
@@ -920,7 +920,7 @@ def cross_plot(A,
     ar = LinearStateSpace(A, C, G, mu_0=np.ones(4))
 
     if steady_state == 'True':
-        μ_x, μ_y, Σ_x, Σ_y = ar.stationary_distributions()
+        μ_x, μ_y, Σ_x, Σ_y, Σ_yx = ar.stationary_distributions()
         ar_state = LinearStateSpace(A, C, G, mu_0=μ_x, Sigma_0=Σ_x)
 
     ymin, ymax = -0.6, 0.6
diff --git a/lectures/perm_income_cons.md b/lectures/perm_income_cons.md
index 2fefa8fb9..2c5def65d 100644
--- a/lectures/perm_income_cons.md
+++ b/lectures/perm_income_cons.md
@@ -345,7 +345,7 @@ G = np.array([[0., 1., 0.]])
 μ_z0 = np.array([[1.0], [0.0], [0.0]])
 Σ_z0 = np.zeros((3, 3))
 Lz = qe.LinearStateSpace(A, C, G, mu_0=μ_z0, Sigma_0=Σ_z0)
-μ_z, μ_y, Σ_z, Σ_y = Lz.stationary_distributions()
+μ_z, μ_y, Σ_z, Σ_y, Σ_yx = Lz.stationary_distributions()
 
 # Mean vector of state for the savings problem
 mxo = np.vstack([μ_z, 0.0])
diff --git a/lectures/samuelson.md b/lectures/samuelson.md
index 244658e52..c1918f4db 100644
--- a/lectures/samuelson.md
+++ b/lectures/samuelson.md
@@ -1316,7 +1316,7 @@ class SamuelsonLSS(LinearStateSpace):
         # values for simulation
         if stationary == True:
             try:
-                self.μ_x, self.μ_y, self.σ_x, self.σ_y = \
+                self.μ_x, self.μ_y, self.σ_x, self.σ_y, self.σ_yx = \
                     self.stationary_distributions()
                 self.μ_0 = self.μ_y
                 self.Σ_0 = self.σ_y