diff --git a/source/rst/linear_models.rst b/source/rst/linear_models.rst
index d15c07e..a412514 100644
--- a/source/rst/linear_models.rst
+++ b/source/rst/linear_models.rst
@@ -799,7 +799,7 @@ Let's now try with 500,000 observations, showing only the histogram (without rot
     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)
@@ -1006,7 +1006,7 @@ This picture shows cross-sectional distributions for :math:`y` at times
         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/source/rst/perm_income_cons.rst b/source/rst/perm_income_cons.rst
index a21c2af..7dee0f3 100644
--- a/source/rst/perm_income_cons.rst
+++ b/source/rst/perm_income_cons.rst
@@ -382,7 +382,7 @@ First, we create the objects for the optimal linear regulator
     μ_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/source/rst/samuelson.rst b/source/rst/samuelson.rst
index 2d2acc3..6082f60 100644
--- a/source/rst/samuelson.rst
+++ b/source/rst/samuelson.rst
@@ -1348,7 +1348,7 @@ methods and attributes) to add more functions to use
             # 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