diff --git a/.github/workflows/collab.yml b/.github/workflows/collab.yml
index 7f047f794..051fce1ca 100644
--- a/.github/workflows/collab.yml
+++ b/.github/workflows/collab.yml
@@ -10,6 +10,12 @@ jobs:
       - uses: actions/checkout@v4
         with:
           ref: ${{ github.event.pull_request.head.sha }}
+      # Install build software
+      - name: Install Build Software & LaTeX (kalman_2)
+        shell: bash -l {0}
+        run: |
+          pip install jupyter-book==1.0.3 quantecon-book-theme==0.8.2 sphinx-tojupyter==0.3.0 sphinxext-rediraffe==0.2.7 sphinxcontrib-youtube==1.3.0 sphinx-togglebutton==0.3.2 arviz sphinx-proof sphinx-exercise sphinx-reredirects
+          apt-get install dvipng texlive texlive-latex-extra texlive-fonts-recommended cm-super    
       - name: Check nvidia drivers
         shell: bash -l {0}
         run: |
@@ -28,11 +34,6 @@ jobs:
           branch: main
           name: build-cache
           path: _build
-      # Install build software
-      - name: Install Build Software
-        shell: bash -l {0}
-        run: |
-          pip install jupyter-book==1.0.3 quantecon-book-theme==0.8.2 sphinx-tojupyter==0.3.0 sphinxext-rediraffe==0.2.7 sphinxcontrib-youtube==1.3.0 sphinx-togglebutton==0.3.2 arviz sphinx-proof sphinx-exercise sphinx-reredirects
       # Build of HTML (Execution Testing)
       - name: Build HTML
         shell: bash -l {0}
diff --git a/lectures/back_prop.md b/lectures/back_prop.md
index 07bb957cf..6ce9a6dab 100644
--- a/lectures/back_prop.md
+++ b/lectures/back_prop.md
@@ -4,9 +4,9 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.11.5
+    jupytext_version: 1.16.7
 kernelspec:
-  display_name: Python 3
+  display_name: Python 3 (ipykernel)
   language: python
   name: python3
 ---
@@ -22,6 +22,12 @@ kernelspec:
 !pip install --upgrade jax
 ```
 
+```{code-cell} ipython3
+import jax
+## to check that gpu is activated in environment
+print(f"JAX backend: {jax.devices()[0].platform}")
+```
+
 In addition to what's included in base Anaconda, we need to install the following packages
 
 ```{code-cell} ipython3
@@ -604,15 +610,4 @@ Image(fig.to_image(format="png"))
 # notebook locally
 ```
 
-```{code-cell} ipython3
-## to check that gpu is activated in environment
-
-from jax.lib import xla_bridge
-print(xla_bridge.get_backend().platform)
-```
 
-```{note}
-**Cloud Environment:** This lecture site is built in a server environment that doesn't have access to a `gpu`
-If you run this lecture locally this lets you know where your code is being executed, either
-via the `cpu` or the `gpu`
-```
diff --git a/lectures/cass_fiscal.md b/lectures/cass_fiscal.md
index 55b6dcf8f..fbc5b96e0 100644
--- a/lectures/cass_fiscal.md
+++ b/lectures/cass_fiscal.md
@@ -4,7 +4,7 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.16.6
+    jupytext_version: 1.16.7
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
@@ -750,7 +750,7 @@ def plot_results(solution, k_ss, c_ss, shocks, shock_param,
     R_bar_path = compute_R_bar_path(shocks, k_path, model, S)
 
     axes[2].plot(R_bar_path[:T], linestyle=linestyle, label=label)
-    axes[2].set_title('$\overline{R}$')
+    axes[2].set_title(r'$\overline{R}$')
     axes[2].axhline(1 / model.β, linestyle='--', color='black')
     
     η_path = compute_η_path(k_path, model, S=T)
@@ -1041,7 +1041,7 @@ Indeed, {eq}`eq:euler_house` or {eq}`eq:diff_second` indicates that a foreseen i
 crease in $\tau_{ct}$ (i.e., a decrease in $(1+\tau_{ct})$
 $(1+\tau_{ct+1})$) operates like an increase in $\tau_{kt}$.
 
-The following figure portrays the response to a foreseen increase in the consumption tax $\tau_c$. 
+The following figure portrays the response to a foreseen increase in the consumption tax $\tau_c$.
 
 ```{code-cell} ipython3
 shocks = {
@@ -1101,7 +1101,6 @@ The  figure shows that:
 - Transition dynamics push $k_t$ (capital stock) toward a new, lower steady-state level. In the new steady state:
     - Consumption is lower due to reduced output from the lower capital stock.
     - Smoother consumption paths occur when $\gamma = 2$ than when $\gamma = 0.2$.
-    
 
 +++
 
@@ -1111,8 +1110,6 @@ foreseen one-time change in a policy variable (a "pulse").
 
 **Experiment 4: Foreseen one-time increase in $g$ from 0.2 to 0.4 in period 10, after which $g$ returns to 0.2 forever**
 
-
-
 ```{code-cell} ipython3
 g_path = np.repeat(0.2, S + 1)
 g_path[10] = 0.4
@@ -1136,6 +1133,7 @@ The figure indicates how:
     - Before $t = 10$, capital accumulates as interest rate changes induce households to prepare for the anticipated increase in government spending.
     - At $t = 10$, the capital stock sharply decreases as the government consumes part of it.
     - $\bar{R}$ jumps above its steady-state value due to the capital reduction and then gradually declines toward its steady-state level.
+
 +++
 
 ### Method 2: Residual Minimization 
diff --git a/lectures/kalman.md b/lectures/kalman.md
index e5f5a193f..cc468acc0 100644
--- a/lectures/kalman.md
+++ b/lectures/kalman.md
@@ -3,8 +3,10 @@ jupytext:
   text_representation:
     extension: .md
     format_name: myst
+    format_version: 0.13
+    jupytext_version: 1.16.7
 kernelspec:
-  display_name: Python 3
+  display_name: Python 3 (ipykernel)
   language: python
   name: python3
 ---
@@ -29,10 +31,9 @@ kernelspec:
 
 In addition to what's in Anaconda, this lecture will need the following libraries:
 
-```{code-cell} ipython
----
-tags: [hide-output]
----
+```{code-cell} ipython3
+:tags: [hide-output]
+
 !pip install quantecon
 ```
 
@@ -54,9 +55,8 @@ Required knowledge: Familiarity with matrix manipulations, multivariate normal d
 
 We'll need the following imports:
 
-```{code-cell} ipython
+```{code-cell} ipython3
 import matplotlib.pyplot as plt
-plt.rcParams["figure.figsize"] = (11, 5)  #set default figure size
 from scipy import linalg
 import numpy as np
 import matplotlib.cm as cm
@@ -122,10 +122,9 @@ $2 \times 2$ covariance matrix.  In our simulations, we will suppose that
 
 This density $p(x)$ is shown below as a contour map, with the center of the red ellipse being equal to $\hat x$.
 
-```{code-cell} python3
----
-tags: [output_scroll]
----
+```{code-cell} ipython3
+:tags: [output_scroll]
+
 # Set up the Gaussian prior density p
 Σ = [[0.4, 0.3], [0.3, 0.45]]
 Σ = np.matrix(Σ)
@@ -186,7 +185,7 @@ def bivariate_normal(x, y, σ_x=1.0, σ_y=1.0, μ_x=0.0, μ_y=0.0, σ_xy=0.0):
 
 def gen_gaussian_plot_vals(μ, C):
     "Z values for plotting the bivariate Gaussian N(μ, C)"
-    m_x, m_y = float(μ[0]), float(μ[1])
+    m_x, m_y = float(μ[0,0]), float(μ[1,0])
     s_x, s_y = np.sqrt(C[0, 0]), np.sqrt(C[1, 1])
     s_xy = C[0, 1]
     return bivariate_normal(X, Y, s_x, s_y, m_x, m_y, s_xy)
@@ -213,7 +212,7 @@ The good news is that the missile has been located by our sensors, which report
 The next figure shows the original prior $p(x)$ and the new reported
 location $y$
 
-```{code-cell} python3
+```{code-cell} ipython3
 fig, ax = plt.subplots(figsize=(10, 8))
 ax.grid()
 
@@ -221,7 +220,7 @@ Z = gen_gaussian_plot_vals(x_hat, Σ)
 ax.contourf(X, Y, Z, 6, alpha=0.6, cmap=cm.jet)
 cs = ax.contour(X, Y, Z, 6, colors="black")
 ax.clabel(cs, inline=1, fontsize=10)
-ax.text(float(y[0]), float(y[1]), "$y$", fontsize=20, color="black")
+ax.text(float(y[0].item()), float(y[1].item()), "$y$", fontsize=20, color="black")
 
 plt.show()
 ```
@@ -284,7 +283,7 @@ This new density $p(x \,|\, y) = N(\hat x^F, \Sigma^F)$ is shown in the next fig
 
 The original density is left in as contour lines for comparison
 
-```{code-cell} python3
+```{code-cell} ipython3
 fig, ax = plt.subplots(figsize=(10, 8))
 ax.grid()
 
@@ -298,7 +297,7 @@ new_Z = gen_gaussian_plot_vals(x_hat_F, Σ_F)
 cs2 = ax.contour(X, Y, new_Z, 6, colors="black")
 ax.clabel(cs2, inline=1, fontsize=10)
 ax.contourf(X, Y, new_Z, 6, alpha=0.6, cmap=cm.jet)
-ax.text(float(y[0]), float(y[1]), "$y$", fontsize=20, color="black")
+ax.text(float(y[0].item()), float(y[1].item()), "$y$", fontsize=20, color="black")
 
 plt.show()
 ```
@@ -391,7 +390,7 @@ A
 Q = 0.3 * \Sigma
 $$
 
-```{code-cell} python3
+```{code-cell} ipython3
 fig, ax = plt.subplots(figsize=(10, 8))
 ax.grid()
 
@@ -415,7 +414,7 @@ new_Z = gen_gaussian_plot_vals(new_x_hat, new_Σ)
 cs3 = ax.contour(X, Y, new_Z, 6, colors="black")
 ax.clabel(cs3, inline=1, fontsize=10)
 ax.contourf(X, Y, new_Z, 6, alpha=0.6, cmap=cm.jet)
-ax.text(float(y[0]), float(y[1]), "$y$", fontsize=20, color="black")
+ax.text(float(y[0].item()), float(y[1].item()), "$y$", fontsize=20, color="black")
 
 plt.show()
 ```
@@ -577,7 +576,7 @@ Your figure should -- modulo randomness -- look something like this
 :class: dropdown
 ```
 
-```{code-cell} python3
+```{code-cell} ipython3
 # Parameters
 θ = 10  # Constant value of state x_t
 A, C, G, H = 1, 0, 1, 1
@@ -598,7 +597,7 @@ xgrid = np.linspace(θ - 5, θ + 2, 200)
 
 for i in range(N):
     # Record the current predicted mean and variance
-    m, v = [float(z) for z in (kalman.x_hat, kalman.Sigma)]
+    m, v = [float(z) for z in (kalman.x_hat.item(), kalman.Sigma.item())]
     # Plot, update filter
     ax.plot(xgrid, norm.pdf(xgrid, loc=m, scale=np.sqrt(v)), label=f'$t={i}$')
     kalman.update(y[i])
@@ -641,7 +640,7 @@ Your figure should show error erratically declining something like this
 :class: dropdown
 ```
 
-```{code-cell} python3
+```{code-cell} ipython3
 ϵ = 0.1
 θ = 10  # Constant value of state x_t
 A, C, G, H = 1, 0, 1, 1
@@ -657,7 +656,7 @@ y = y.flatten()
 
 for t in range(T):
     # Record the current predicted mean and variance and plot their densities
-    m, v = [float(temp) for temp in (kalman.x_hat, kalman.Sigma)]
+    m, v = [float(temp) for temp in (kalman.x_hat.item(), kalman.Sigma.item())]
 
     f = lambda x: norm.pdf(x, loc=m, scale=np.sqrt(v))
     integral, error = quad(f, θ - ϵ, θ + ϵ)
@@ -745,7 +744,7 @@ Observe how, after an initial learning period, the Kalman filter performs quite
 :class: dropdown
 ```
 
-```{code-cell} python3
+```{code-cell} ipython3
 # Define A, C, G, H
 G = np.identity(2)
 H = np.sqrt(0.5) * np.identity(2)
diff --git a/lectures/kalman_2.md b/lectures/kalman_2.md
index 0689b4538..a1eb32490 100644
--- a/lectures/kalman_2.md
+++ b/lectures/kalman_2.md
@@ -4,7 +4,7 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.14.4
+    jupytext_version: 1.16.7
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
@@ -237,7 +237,7 @@ for t in range(1, T):
     x_hat, Σ = kalman.x_hat, kalman.Sigma
     Σ_t[:, :, t-1] = Σ
     x_hat_t[:, t-1] = x_hat.reshape(-1)
-    y_hat_t[t-1] = worker.G @ x_hat
+    [y_hat_t[t-1]] = worker.G @ x_hat
 
 x_hat_t = np.concatenate((x[:, 1][:, np.newaxis], 
                     x_hat_t), axis=1)
@@ -253,7 +253,6 @@ We also plot $E [u_0 | y^{t-1}]$, which is  the firm inference about  a worker's
 We can  watch as the  firm's inference  $E [u_0 | y^{t-1}]$ of the worker's work ethic converges toward the hidden   $u_0$, which is not directly observed by the firm.
 
 ```{code-cell} ipython3
-
 fig, ax = plt.subplots(1, 2)
 
 ax[0].plot(y_hat_t, label=r'$E[y_t| y^{t-1}]$')
@@ -273,6 +272,7 @@ ax[1].legend()
 fig.tight_layout()
 plt.show()
 ```
+
 ## Some Computational Experiments
 
 Let's look at  $\Sigma_0$ and $\Sigma_T$ in order to see how much the firm learns about the hidden state during the horizon we have set.
@@ -290,7 +290,6 @@ Evidently,  entries in the conditional covariance matrix become smaller over tim
 It is enlightening to  portray how  conditional covariance matrices $\Sigma_t$ evolve by plotting confidence ellipsoides around $E [x_t |y^{t-1}] $ at various $t$'s.
 
 ```{code-cell} ipython3
-
 # Create a grid of points for contour plotting
 h_range = np.linspace(x_hat_t[0, :].min()-0.5*Σ_t[0, 0, 1], 
                       x_hat_t[0, :].max()+0.5*Σ_t[0, 0, 1], 100)
@@ -338,7 +337,6 @@ For example, let's say $h_0 = 0$ and $u_0 = 4$.
 Here is one way to do this.
 
 ```{code-cell} ipython3
-
 # For example, we might want h_0 = 0 and u_0 = 4
 mu_0 = np.array([0.0, 4.0])
 
@@ -361,7 +359,6 @@ print('u_0 =', u_0)
 Another way to accomplish the same goal is to use the following code.
 
 ```{code-cell} ipython3
-
 # If we want to set the initial 
 # h_0 = hhat_0 = 0 and u_0 = uhhat_0 = 4.0:
 worker = create_worker(hhat_0=0.0, uhat_0=4.0)
@@ -398,8 +395,8 @@ for t in range(1, T):
     kalman.update(y[t])
     x_hat, Σ = kalman.x_hat, kalman.Sigma
     Σ_t.append(Σ)
-    y_hat_t[t-1] = worker.G @ x_hat
-    u_hat_t[t-1] = x_hat[1]
+    [y_hat_t[t-1]] = worker.G @ x_hat
+    [u_hat_t[t-1]] = x_hat[1]
 
 
 # Generate plots for y_hat_t and u_hat_t
@@ -426,10 +423,9 @@ plt.show()
 More generally, we can change some or all of the parameters defining a worker in our `create_worker`
 namedtuple.
 
-Here is an example. 
+Here is an example.
 
 ```{code-cell} ipython3
-
 # We can set these parameters when creating a worker -- just like classes!
 hard_working_worker =  create_worker(α=.4, β=.8, 
                         hhat_0=7.0, uhat_0=100, σ_h=2.5, σ_u=3.2)
@@ -475,8 +471,8 @@ def simulate_workers(worker, T, ax, mu_0=None, Sigma_0=None,
         kalman.update(y[i])
         x_hat, Σ = kalman.x_hat, kalman.Sigma
         Σ_t.append(Σ)
-        y_hat_t[i] = worker.G @ x_hat
-        u_hat_t[i] = x_hat[1]
+        [y_hat_t[i]] = worker.G @ x_hat
+        [u_hat_t[i]] = x_hat[1]
 
     if diff == True:
         title = ('Difference between inferred and true work ethic over time' 
@@ -503,7 +499,6 @@ def simulate_workers(worker, T, ax, mu_0=None, Sigma_0=None,
 ```
 
 ```{code-cell} ipython3
-
 num_workers = 3
 T = 50
 fig, ax = plt.subplots(figsize=(7, 7))
@@ -516,7 +511,6 @@ plt.show()
 ```
 
 ```{code-cell} ipython3
-
 # We can also generate plots of u_t:
 
 T = 50
@@ -539,7 +533,6 @@ plt.show()
 ```
 
 ```{code-cell} ipython3
-
 # We can also use exact u_0=1 and h_0=2 for all workers
 
 T = 50
@@ -568,7 +561,6 @@ plt.show()
 ```
 
 ```{code-cell} ipython3
-
 # We can generate a plot for only one of the workers:
 
 T = 50
diff --git a/lectures/kesten_processes.md b/lectures/kesten_processes.md
index 95a6afd30..1fe6e921b 100644
--- a/lectures/kesten_processes.md
+++ b/lectures/kesten_processes.md
@@ -3,8 +3,10 @@ jupytext:
   text_representation:
     extension: .md
     format_name: myst
+    format_version: 0.13
+    jupytext_version: 1.16.7
 kernelspec:
-  display_name: Python 3
+  display_name: Python 3 (ipykernel)
   language: python
   name: python3
 ---
@@ -28,10 +30,9 @@ kernelspec:
 
 In addition to what's in Anaconda, this lecture will need the following libraries:
 
-```{code-cell} ipython
----
-tags: [hide-output]
----
+```{code-cell} ipython3
+:tags: [hide-output]
+
 !pip install quantecon
 !pip install --upgrade yfinance
 ```
@@ -53,9 +54,8 @@ We will discuss these issues as we go along.
 
 Let's start with some imports:
 
-```{code-cell} ipython
+```{code-cell} ipython3
 import matplotlib.pyplot as plt
-plt.rcParams["figure.figsize"] = (11, 5)  #set default figure size
 import numpy as np
 import quantecon as qe
 ```
@@ -63,7 +63,7 @@ import quantecon as qe
 The following two lines are only added to avoid a `FutureWarning` caused by
 compatibility issues between pandas and matplotlib.
 
-```{code-cell} ipython
+```{code-cell} ipython3
 from pandas.plotting import register_matplotlib_converters
 register_matplotlib_converters()
 ```
@@ -105,7 +105,8 @@ For example, consider the following plot of daily returns on the Nasdaq
 Composite Index for the period 1st January 2006 to 1st November 2019.
 
 (ndcode)=
-```{code-cell} python3
+
+```{code-cell} ipython3
 import yfinance as yf
 
 s = yf.download('^IXIC', '2006-1-1', '2019-11-1', auto_adjust=False)['Adj Close']
diff --git a/lectures/linear_models.md b/lectures/linear_models.md
index ac0c0616c..c887aa5db 100644
--- a/lectures/linear_models.md
+++ b/lectures/linear_models.md
@@ -3,8 +3,10 @@ jupytext:
   text_representation:
     extension: .md
     format_name: myst
+    format_version: 0.13
+    jupytext_version: 1.16.7
 kernelspec:
-  display_name: Python 3
+  display_name: Python 3 (ipykernel)
   language: python
   name: python3
 ---
@@ -33,10 +35,9 @@ kernelspec:
 
 In addition to what's in Anaconda, this lecture will need the following libraries:
 
-```{code-cell} ipython
----
-tags: [hide-output]
----
+```{code-cell} ipython3
+:tags: [hide-output]
+
 !pip install quantecon
 ```
 
@@ -65,9 +66,8 @@ Its many applications include:
 
 Let's start with some imports:
 
-```{code-cell} ipython
+```{code-cell} ipython3
 import matplotlib.pyplot as plt
-plt.rcParams["figure.figsize"] = (11, 5)  #set default figure size
 import numpy as np
 from quantecon import LinearStateSpace
 from scipy.stats import norm
@@ -186,7 +186,8 @@ You can confirm that under these definitions, {eq}`st_space_rep` and {eq}`st_ex_
 The next figure shows the dynamics of this process when $\phi_0 = 1.1, \phi_1=0.8, \phi_2 = -0.8, y_0 = y_{-1} = 1$.
 
 (lss_sode_fig)=
-```{code-cell} python3
+
+```{code-cell} ipython3
 def plot_lss(A,
          C,
          G,
@@ -205,7 +206,7 @@ def plot_lss(A,
     plt.show()
 ```
 
-```{code-cell} python3
+```{code-cell} ipython3
 ϕ_0, ϕ_1, ϕ_2 = 1.1, 0.8, -0.8
 
 A = [[1,     0,     0  ],
@@ -269,7 +270,8 @@ y_{-3} = 1
 $$
 
 (lss_uap_fig)=
-```{code-cell} python3
+
+```{code-cell} ipython3
 ϕ_1, ϕ_2, ϕ_3, ϕ_4 = 0.5, -0.2, 0, 0.5
 σ = 0.2
 
@@ -651,7 +653,7 @@ The system in question is the univariate autoregressive model {eq}`eq_ar_rep`.
 
 The values of $y_T$ are represented by black dots in the left-hand figure
 
-```{code-cell} python3
+```{code-cell} ipython3
 def cross_section_plot(A,
                    C,
                    G,
@@ -695,7 +697,7 @@ def cross_section_plot(A,
     plt.show()
 ```
 
-```{code-cell} python3
+```{code-cell} ipython3
 ϕ_1, ϕ_2, ϕ_3, ϕ_4 = 0.5, -0.2, 0, 0.5
 σ = 0.1
 
@@ -716,14 +718,14 @@ that shows relative frequencies from our sample of 20 $y_T$'s.
 
 Here is another figure, this time with 100 observations
 
-```{code-cell} python3
+```{code-cell} ipython3
 t = 100
 cross_section_plot(A_2, C_2, G_2, T=t)
 ```
 
 Let's now try with 500,000 observations, showing only the histogram (without rotation)
 
-```{code-cell} python3
+```{code-cell} ipython3
 T = 100
 ymin=-0.8
 ymax=1.25
@@ -733,7 +735,7 @@ 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, Sigma_yx = ar.stationary_distributions()
-f_y = norm(loc=float(mu_y), scale=float(np.sqrt(Sigma_y)))
+f_y = norm(loc=float(mu_y.item()), scale=float(np.sqrt(Sigma_y.item())))
 y = y.flatten()
 ygrid = np.linspace(ymin, ymax, 150)
 
@@ -776,7 +778,8 @@ The parameters are the same as for the preceding figures,
 and the sample size is relatively small ($I=20$).
 
 (lss_em_fig)=
-```{code-cell} python3
+
+```{code-cell} ipython3
 I = 20
 T = 50
 ymin = -0.5
@@ -800,7 +803,7 @@ m = ar.moment_sequence()
 population_means = []
 for t in range(T):
     μ_x, μ_y, Σ_x, Σ_y = next(m)
-    population_means.append(float(μ_y))
+    population_means.append(float(μ_y.item()))
 
 ax.plot(population_means, color='g', lw=2, alpha=0.8, label=r'$G\mu_t$')
 ax.set_ylim(ymin, ymax)
@@ -904,7 +907,7 @@ Let's look at some more time series from the same model that we analyzed above.
 This picture shows cross-sectional distributions for $y$ at times
 $T, T', T''$
 
-```{code-cell} python3
+```{code-cell} ipython3
 def cross_plot(A,
             C,
             G,
@@ -944,7 +947,7 @@ def cross_plot(A,
     plt.show()
 ```
 
-```{code-cell} python3
+```{code-cell} ipython3
 cross_plot(A_2, C_2, G_2)
 ```
 
@@ -988,7 +991,8 @@ where $\mu_{\infty}$ and $\Sigma_{\infty}$ are fixed points of {eq}`lss_mut_line
 Let's see what happens to the preceding figure if we start $x_0$ at the stationary distribution.
 
 (lss_s_fig)=
-```{code-cell} python3
+
+```{code-cell} ipython3
 cross_plot(A_2, C_2, G_2, steady_state='True')
 ```
 
diff --git a/lectures/markov_perf.md b/lectures/markov_perf.md
index 9e9d51c8f..f4e0a4f24 100644
--- a/lectures/markov_perf.md
+++ b/lectures/markov_perf.md
@@ -3,8 +3,10 @@ jupytext:
   text_representation:
     extension: .md
     format_name: myst
+    format_version: 0.13
+    jupytext_version: 1.16.7
 kernelspec:
-  display_name: Python 3
+  display_name: Python 3 (ipykernel)
   language: python
   name: python3
 ---
@@ -29,10 +31,9 @@ kernelspec:
 
 In addition to what's in Anaconda, this lecture will need the following libraries:
 
-```{code-cell} ipython
----
-tags: [hide-output]
----
+```{code-cell} ipython3
+:tags: [hide-output]
+
 !pip install quantecon
 ```
 
@@ -57,9 +58,8 @@ Other references include chapter 7 of {cite}`Ljungqvist2012`.
 
 Let's start with some standard imports:
 
-```{code-cell} ipython
+```{code-cell} ipython3
 import matplotlib.pyplot as plt
-plt.rcParams["figure.figsize"] = (11, 5)  #set default figure size
 import numpy as np
 import quantecon as qe
 ```
@@ -431,8 +431,10 @@ Consider the previously presented duopoly model with parameter values of:
 
 From these, we compute the infinite horizon MPE using the preceding code
 
-```{code-cell} python3
+```{code-cell} ipython3
 :load: _static/lecture_specific/markov_perf/duopoly_mpe.py
+
+
 ```
 
 Running the code produces the following output.
@@ -443,7 +445,7 @@ In particular, let's take F2 as computed above, plug it into {eq}`eq_mpe_p1p` an
 
 We hope that the resulting policy will agree with F1 as computed above
 
-```{code-cell} python3
+```{code-cell} ipython3
 Λ1 = A - B2 @ F2
 lq1 = qe.LQ(Q1, R1, Λ1, B1, beta=β)
 P1_ih, F1_ih, d = lq1.stationary_values()
@@ -454,7 +456,7 @@ This is close enough for rock and roll, as they say in the trade.
 
 Indeed, np.allclose agrees with our assessment
 
-```{code-cell} python3
+```{code-cell} ipython3
 np.allclose(F1, F1_ih)
 ```
 
@@ -470,7 +472,7 @@ The following program
 * computes the evolution of $x_t$ using {eq}`eq_mpe_cle`.
 * extracts and plots industry output $q_t = q_{1t} + q_{2t}$ and price $p_t = a_0 - a_1 q_t$.
 
-```{code-cell} python3
+```{code-cell} ipython3
 AF = A - B1 @ F1 - B2 @ F2
 n = 20
 x = np.empty((3, n))
@@ -527,7 +529,7 @@ The optimal policy in the monopolist case can be computed using [QuantEcon.py](h
 
 First, let's compute the duopoly MPE under the stated parameters
 
-```{code-cell} python3
+```{code-cell} ipython3
 # == Parameters == #
 a0 = 10.0
 a1 = 2.0
@@ -558,7 +560,7 @@ F1, F2, P1, P2 = qe.nnash(A, B1, B2, R1, R2, Q1,
 Now we evaluate the time path of industry output and prices given
 initial condition $q_{10} = q_{20} = 1$.
 
-```{code-cell} python3
+```{code-cell} ipython3
 AF = A - B1 @ F1 - B2 @ F2
 n = 20
 x = np.empty((3, n))
@@ -600,7 +602,7 @@ in the law of motion $x_{t+1} = A x_t + B u_t$.
 We solve for the optimal policy $u_t = - Fx_t$ and track the
 resulting dynamics of $\{q_t\}$, starting at $q_0 = 2.0$.
 
-```{code-cell} python3
+```{code-cell} ipython3
 R = a1
 Q = γ
 A = B = 1
@@ -613,13 +615,13 @@ x0 = qm[0] - q_bar
 x = x0
 for i in range(1, n):
     x = A * x - B * F * x
-    qm[i] = float(x) + q_bar
+    qm[i] = float(x.item()) + q_bar
 pm = a0 - a1 * qm
 ```
 
 Let's have a look at the different time paths
 
-```{code-cell} python3
+```{code-cell} ipython3
 fig, axes = plt.subplots(2, 1, figsize=(9, 9))
 
 ax = axes[0]
@@ -712,7 +714,7 @@ The exercise is to calculate these matrices and compute the following figures.
 
 The first figure shows the dynamics of inventories for each firm when the parameters are
 
-```{code-cell} python3
+```{code-cell} ipython3
 δ = 0.02
 D = np.array([[-1, 0.5], [0.5, -1]])
 b = np.array([25, 25])
@@ -743,7 +745,7 @@ In this exercise, reproduce the figure when $\delta = 0.02$.
 
 We treat the case $\delta = 0.02$
 
-```{code-cell} python3
+```{code-cell} ipython3
 δ = 0.02
 D = np.array([[-1, 0.5], [0.5, -1]])
 b = np.array([25, 25])
@@ -772,7 +774,7 @@ $$
 
 we set up the matrices as follows:
 
-```{code-cell} python3
+```{code-cell} ipython3
 # ==  Create matrices needed to compute the Nash feedback equilibrium == #
 
 A = np.array([[δ_1,      0,    -δ_1 * b[0]],
@@ -812,7 +814,7 @@ M2 = np.copy(M1)
 
 We can now compute the equilibrium using `qe.nnash`
 
-```{code-cell} python3
+```{code-cell} ipython3
 F1, F2, P1, P2 = qe.nnash(A, B1, B2, R1,
                           R2, Q1, Q2, S1,
                           S2, W1, W2, M1, M2)
@@ -827,7 +829,7 @@ print(F2)
 Now let's look at the dynamics of inventories, and reproduce the graph
 corresponding to $\delta = 0.02$
 
-```{code-cell} python3
+```{code-cell} ipython3
 AF = A - B1 @ F1 - B2 @ F2
 n = 25
 x = np.empty((3, n))
@@ -845,4 +847,4 @@ plt.show()
 ```
 
 ```{solution-end}
-```
+```
\ No newline at end of file
diff --git a/lectures/pandas_panel.md b/lectures/pandas_panel.md
index 42aa72c73..664545b81 100644
--- a/lectures/pandas_panel.md
+++ b/lectures/pandas_panel.md
@@ -3,8 +3,10 @@ jupytext:
   text_representation:
     extension: .md
     format_name: myst
+    format_version: 0.13
+    jupytext_version: 1.16.7
 kernelspec:
-  display_name: Python 3
+  display_name: Python 3 (ipykernel)
   language: python
   name: python3
 ---
@@ -62,11 +64,11 @@ countries and assign it to `realwage`.
 
 The dataset can be accessed with the following link:
 
-```{code-cell} python3
+```{code-cell} ipython3
 url1 = 'https://raw.githubusercontent.com/QuantEcon/lecture-python/master/source/_static/lecture_specific/pandas_panel/realwage.csv'
 ```
 
-```{code-cell} python3
+```{code-cell} ipython3
 import pandas as pd
 
 # Display 6 columns for viewing purposes
@@ -80,7 +82,7 @@ realwage = pd.read_csv(url1)
 
 Let's have a look at what we've got to work with
 
-```{code-cell} python3
+```{code-cell} ipython3
 realwage.head()  # Show first 5 rows
 ```
 
@@ -92,7 +94,7 @@ We will use `pivot_table` to create a wide format panel, with a `MultiIndex` to
 
 By passing a list in columns, we can create a `MultiIndex` in our column axis
 
-```{code-cell} python3
+```{code-cell} ipython3
 realwage = realwage.pivot_table(values='value',
                                 index='Time',
                                 columns=['Country', 'Series', 'Pay period'])
@@ -101,7 +103,7 @@ realwage.head()
 
 To more easily filter our time series data, later on, we will convert the index into a `DateTimeIndex`
 
-```{code-cell} python3
+```{code-cell} ipython3
 realwage.index = pd.to_datetime(realwage.index)
 type(realwage.index)
 ```
@@ -113,18 +115,18 @@ Series > Pay period).
 A `MultiIndex` is the simplest and most flexible way to manage panel
 data in pandas
 
-```{code-cell} python3
+```{code-cell} ipython3
 type(realwage.columns)
 ```
 
-```{code-cell} python3
+```{code-cell} ipython3
 realwage.columns.names
 ```
 
 Like before, we can select the country (the top level of our
 `MultiIndex`)
 
-```{code-cell} python3
+```{code-cell} ipython3
 realwage['United States'].head()
 ```
 
@@ -135,15 +137,15 @@ throughout this lecture to reshape our dataframe into a format we need.
 the row index (`.unstack()` works in the opposite direction - try it
 out)
 
-```{code-cell} python3
-realwage.stack().head()
+```{code-cell} ipython3
+realwage.stack(future_stack=True).head()
 ```
 
 We can also pass in an argument to select the level we would like to
 stack
 
-```{code-cell} python3
-realwage.stack(level='Country').head()
+```{code-cell} ipython3
+realwage.stack(level='Country', future_stack=True).head()
 ```
 
 Using a `DatetimeIndex` makes it easy to select a particular time
@@ -152,8 +154,8 @@ period.
 Selecting one year and stacking the two lower levels of the
 `MultiIndex` creates a cross-section of our panel data
 
-```{code-cell} python3
-realwage.loc['2015'].stack(level=(1, 2)).transpose().head()
+```{code-cell} ipython3
+realwage.loc['2015'].stack(level=(1, 2), future_stack=True).transpose().head()
 ```
 
 For the rest of lecture, we will work with a dataframe of the hourly
@@ -164,7 +166,7 @@ To create our filtered dataframe (`realwage_f`), we can use the `xs`
 method to select values at lower levels in the multiindex, while keeping
 the higher levels (countries in this case)
 
-```{code-cell} python3
+```{code-cell} ipython3
 realwage_f = realwage.xs(('Hourly', 'In 2015 constant prices at 2015 USD exchange rates'),
                          level=('Pay period', 'Series'), axis=1)
 realwage_f.head()
@@ -182,11 +184,11 @@ function.
 
 The dataset can be accessed with the following link:
 
-```{code-cell} python3
+```{code-cell} ipython3
 url2 = 'https://raw.githubusercontent.com/QuantEcon/lecture-python/master/source/_static/lecture_specific/pandas_panel/countries.csv'
 ```
 
-```{code-cell} python3
+```{code-cell} ipython3
 worlddata = pd.read_csv(url2, sep=';')
 worlddata.head()
 ```
@@ -194,7 +196,7 @@ worlddata.head()
 First, we'll select just the country and continent variables from
 `worlddata` and rename the column to 'Country'
 
-```{code-cell} python3
+```{code-cell} ipython3
 worlddata = worlddata[['Country (en)', 'Continent']]
 worlddata = worlddata.rename(columns={'Country (en)': 'Country'})
 worlddata.head()
@@ -209,7 +211,7 @@ Our dataframes will be merged using country names, requiring us to use
 the transpose of `realwage_f` so that rows correspond to country names
 in both dataframes
 
-```{code-cell} python3
+```{code-cell} ipython3
 realwage_f.transpose().head()
 ```
 
@@ -243,7 +245,7 @@ the index, so we set `left_index=True`.
 Our 'right' dataframe (`worlddata`) contains countries in the
 'Country' column, so we set `right_on='Country'`
 
-```{code-cell} python3
+```{code-cell} ipython3
 merged = pd.merge(realwage_f.transpose(), worlddata,
                   how='left', left_index=True, right_on='Country')
 merged.head()
@@ -255,7 +257,7 @@ have `NaN` in the Continent column.
 To check whether this has occurred, we can use `.isnull()` on the
 continent column and filter the merged dataframe
 
-```{code-cell} python3
+```{code-cell} ipython3
 merged[merged['Continent'].isnull()]
 ```
 
@@ -269,7 +271,7 @@ continent from the dictionary.
 
 Notice how countries not in our dictionary are mapped with `NaN`
 
-```{code-cell} python3
+```{code-cell} ipython3
 missing_continents = {'Korea': 'Asia',
                       'Russian Federation': 'Europe',
                       'Slovak Republic': 'Europe'}
@@ -282,7 +284,7 @@ We don't want to overwrite the entire series with this mapping.
 `.fillna()` only fills in `NaN` values in `merged['Continent']`
 with the mapping, while leaving other values in the column unchanged
 
-```{code-cell} python3
+```{code-cell} ipython3
 merged['Continent'] = merged['Continent'].fillna(merged['Country'].map(missing_continents))
 
 # Check for whether continents were correctly mapped
@@ -294,13 +296,12 @@ We will also combine the Americas into a single continent - this will make our v
 
 To do this, we will use `.replace()` and loop through a list of the continent values we want to replace
 
-```{code-cell} python3
+```{code-cell} ipython3
 replace = ['Central America', 'North America', 'South America']
 
 for country in replace:
-    merged['Continent'].replace(to_replace=country,
-                                value='America',
-                                inplace=True)
+    merged.Continent = merged.Continent.replace(to_replace=country,
+                                value='America')
 ```
 
 Now that we have all the data we want in a single `DataFrame`, we will
@@ -311,7 +312,7 @@ can efficiently filter our dataframe later on.
 
 By default, levels will be sorted top-down
 
-```{code-cell} python3
+```{code-cell} ipython3
 merged = merged.set_index(['Continent', 'Country']).sort_index()
 merged.head()
 ```
@@ -319,14 +320,14 @@ merged.head()
 While merging, we lost our `DatetimeIndex`, as we merged columns that
 were not in datetime format
 
-```{code-cell} python3
+```{code-cell} ipython3
 merged.columns
 ```
 
 Now that we have set the merged columns as the index, we can recreate a
 `DatetimeIndex` using `.to_datetime()`
 
-```{code-cell} python3
+```{code-cell} ipython3
 merged.columns = pd.to_datetime(merged.columns)
 merged.columns = merged.columns.rename('Time')
 merged.columns
@@ -335,7 +336,7 @@ merged.columns
 The `DatetimeIndex` tends to work more smoothly in the row axis, so we
 will go ahead and transpose `merged`
 
-```{code-cell} python3
+```{code-cell} ipython3
 merged = merged.transpose()
 merged.head()
 ```
@@ -353,20 +354,20 @@ For example, we can calculate the average real minimum wage for each
 country over the period 2006 to 2016 (the default is to aggregate over
 rows)
 
-```{code-cell} python3
+```{code-cell} ipython3
 merged.mean().head(10)
 ```
 
 Using this series, we can plot the average real minimum wage over the
 past decade for each country in our data set
 
-```{code-cell} ipython
+```{code-cell} ipython3
 import matplotlib.pyplot as plt
 import seaborn as sns
 sns.set_theme()
 ```
 
-```{code-cell} ipython
+```{code-cell} ipython3
 merged.mean().sort_values(ascending=False).plot(kind='bar',
                                                 title="Average real minimum wage 2006 - 2016")
 
@@ -381,13 +382,13 @@ plt.show()
 Passing in `axis=1` to `.mean()` will aggregate over columns (giving
 the average minimum wage for all countries over time)
 
-```{code-cell} python3
+```{code-cell} ipython3
 merged.mean(axis=1).head()
 ```
 
 We can plot this time series as a line graph
 
-```{code-cell} python3
+```{code-cell} ipython3
 merged.mean(axis=1).plot()
 plt.title('Average real minimum wage 2006 - 2016')
 plt.ylabel('2015 USD')
@@ -398,14 +399,14 @@ plt.show()
 We can also specify a level of the `MultiIndex` (in the column axis)
 to aggregate over
 
-```{code-cell} python3
-merged.groupby(level='Continent', axis=1).mean().head()
+```{code-cell} ipython3
+merged.T.groupby(level='Continent').mean().T.head()
 ```
 
 We can plot the average minimum wages in each continent as a time series
 
-```{code-cell} python3
-merged.groupby(level='Continent', axis=1).mean().plot()
+```{code-cell} ipython3
+merged.T.groupby(level='Continent').mean().T.plot()
 plt.title('Average real minimum wage')
 plt.ylabel('2015 USD')
 plt.xlabel('Year')
@@ -414,9 +415,9 @@ plt.show()
 
 We will drop Australia as a continent for plotting purposes
 
-```{code-cell} python3
+```{code-cell} ipython3
 merged = merged.drop('Australia', level='Continent', axis=1)
-merged.groupby(level='Continent', axis=1).mean().plot()
+merged.T.groupby(level='Continent').mean().T.plot()
 plt.title('Average real minimum wage')
 plt.ylabel('2015 USD')
 plt.xlabel('Year')
@@ -426,8 +427,8 @@ plt.show()
 `.describe()` is useful for quickly retrieving a number of common
 summary statistics
 
-```{code-cell} python3
-merged.stack().describe()
+```{code-cell} ipython3
+merged.stack(future_stack=True).describe()
 ```
 
 This is a simplified way to use `groupby`.
@@ -444,8 +445,8 @@ a new `DataFrameGroupBy` object with data split into groups.
 Let's split `merged` by continent again, this time using the
 `groupby` function, and name the resulting object `grouped`
 
-```{code-cell} python3
-grouped = merged.groupby(level='Continent', axis=1)
+```{code-cell} ipython3
+grouped = merged.T.groupby(level='Continent')
 grouped
 ```
 
@@ -457,7 +458,7 @@ each continent using `.size()`.
 
 In this case, our new data structure is a `Series`
 
-```{code-cell} python3
+```{code-cell} ipython3
 grouped.size()
 ```
 
@@ -468,11 +469,11 @@ minimum wages in 2016 for each continent.
 `grouped.groups.keys()` will return the keys from the `groupby`
 object
 
-```{code-cell} python3
+```{code-cell} ipython3
 continents = grouped.groups.keys()
 
 for continent in continents:
-    sns.kdeplot(grouped.get_group(continent).loc['2015'].unstack(), label=continent, fill=True)
+    sns.kdeplot(grouped.get_group(continent).T.loc['2015'].unstack(), label=continent, fill=True)
 
 plt.title('Real minimum wages in 2015')
 plt.xlabel('US dollars')
@@ -500,7 +501,7 @@ in Europe by age and sex from [Eurostat](https://ec.europa.eu/eurostat/data/data
 
 The dataset can be accessed with the following link:
 
-```{code-cell} python3
+```{code-cell} ipython3
 url3 = 'https://raw.githubusercontent.com/QuantEcon/lecture-python/master/source/_static/lecture_specific/pandas_panel/employ.csv'
 ```
 
@@ -519,7 +520,7 @@ Write a program that quickly returns all values in the `MultiIndex`.
 :class: dropdown
 ```
 
-```{code-cell} python3
+```{code-cell} ipython3
 employ = pd.read_csv(url3)
 employ = employ.pivot_table(values='Value',
                             index=['DATE'],
@@ -531,13 +532,13 @@ employ.head()
 This is a large dataset so it is useful to explore the levels and
 variables available
 
-```{code-cell} python3
+```{code-cell} ipython3
 employ.columns.names
 ```
 
 Variables within levels can be quickly retrieved with a loop
 
-```{code-cell} python3
+```{code-cell} ipython3
 for name in employ.columns.names:
     print(name, employ.columns.get_level_values(name).unique())
 ```
@@ -571,7 +572,7 @@ by age group and sex.
 To easily filter by country, swap `GEO` to the top level and sort the
 `MultiIndex`
 
-```{code-cell} python3
+```{code-cell} ipython3
 employ.columns = employ.columns.swaplevel(0,-1)
 employ = employ.sort_index(axis=1)
 ```
@@ -581,7 +582,7 @@ We need to get rid of a few items in `GEO` which are not countries.
 A fast way to get rid of the EU areas is to use a list comprehension to
 find the level values in `GEO` that begin with 'Euro'
 
-```{code-cell} python3
+```{code-cell} ipython3
 geo_list = employ.columns.get_level_values('GEO').unique().tolist()
 countries = [x for x in geo_list if not x.startswith('Euro')]
 employ = employ[countries]
@@ -591,7 +592,7 @@ employ.columns.get_level_values('GEO').unique()
 Select only percentage employed in the active population from the
 dataframe
 
-```{code-cell} python3
+```{code-cell} ipython3
 employ_f = employ.xs(('Percentage of total population', 'Active population'),
                      level=('UNIT', 'INDIC_EM'),
                      axis=1)
@@ -600,11 +601,11 @@ employ_f.head()
 
 Drop the 'Total' value before creating the grouped boxplot
 
-```{code-cell} python3
+```{code-cell} ipython3
 employ_f = employ_f.drop('Total', level='SEX', axis=1)
 ```
 
-```{code-cell} python3
+```{code-cell} ipython3
 box = employ_f.loc['2015'].unstack().reset_index()
 sns.boxplot(x="AGE", y=0, hue="SEX", data=box, palette=("husl"), showfliers=False)
 plt.xlabel('')
@@ -616,4 +617,4 @@ plt.show()
 ```
 
 ```{solution-end}
-```
+ ```
diff --git a/lectures/perm_income_cons.md b/lectures/perm_income_cons.md
index 2f46724e9..1732f6ba0 100644
--- a/lectures/perm_income_cons.md
+++ b/lectures/perm_income_cons.md
@@ -3,8 +3,10 @@ jupytext:
   text_representation:
     extension: .md
     format_name: myst
+    format_version: 0.13
+    jupytext_version: 1.16.7
 kernelspec:
-  display_name: Python 3
+  display_name: Python 3 (ipykernel)
   language: python
   name: python3
 ---
@@ -29,10 +31,9 @@ kernelspec:
 
 In addition to what's in Anaconda, this lecture will need the following libraries:
 
-```{code-cell} ipython
----
-tags: [hide-output]
----
+```{code-cell} ipython3
+:tags: [hide-output]
+
 !pip install quantecon
 ```
 
@@ -74,9 +75,8 @@ The model will prove useful for illustrating concepts such as
 
 Let's start with some imports:
 
-```{code-cell} ipython
+```{code-cell} ipython3
 import matplotlib.pyplot as plt
-plt.rcParams["figure.figsize"] = (11, 5)  #set default figure size
 import quantecon as qe
 import numpy as np
 import scipy.linalg as la
@@ -329,7 +329,7 @@ In what follows we set it equal to unity.
 
 First, we create the objects for the optimal linear regulator
 
-```{code-cell} python3
+```{code-cell} ipython3
 # Set parameters
 α, β, ρ1, ρ2, σ = 10.0, 0.95, 0.9, 0.0, 1.0
 
@@ -364,7 +364,7 @@ sxbewley = sxo
 
 The next step is to create the matrices for the LQ system
 
-```{code-cell} python3
+```{code-cell} ipython3
 A12 = np.zeros((3,1))
 ALQ_l = np.hstack([A, A12])
 ALQ_r = np.array([[0, -R, 0, R]])
@@ -383,7 +383,7 @@ CLQ = np.array([0., σ, 0., 0.]).reshape(4,1)
 
 Let's print these out and have a look at them
 
-```{code-cell} python3
+```{code-cell} ipython3
 print(f"A = \n {ALQ}")
 print(f"B = \n {BLQ}")
 print(f"R = \n {RLQ}")
@@ -392,14 +392,14 @@ print(f"Q = \n {QLQ}")
 
 Now create the appropriate instance of an LQ model
 
-```{code-cell} python3
+```{code-cell} ipython3
 lqpi = qe.LQ(QLQ, RLQ, ALQ, BLQ, C=CLQ, beta=β_LQ)
 ```
 
 We'll save the implied optimal policy function soon compare them with what we get by
 employing an alternative solution method
 
-```{code-cell} python3
+```{code-cell} ipython3
 P, F, d = lqpi.stationary_values()  # Compute value function and decision rule
 ABF = ALQ - BLQ @ F  #  Form closed loop system
 ```
@@ -428,7 +428,7 @@ $$
 
 Now we'll apply the formulas in this system
 
-```{code-cell} python3
+```{code-cell} ipython3
 # Use the above formulas to create the optimal policies for b_{t+1} and c_t
 b_pol = G @ la.inv(np.eye(3, 3) - β * A) @ (A - np.eye(3, 3))
 c_pol = (1 - β) * G @ la.inv(np.eye(3, 3) - β * A)
@@ -453,13 +453,13 @@ G_LSS = np.hstack([G_LSS1, G_LSS2])
 
 `A_LSS` calculated as we have here should equal `ABF` calculated above using the LQ model
 
-```{code-cell} python3
+```{code-cell} ipython3
 ABF - A_LSS
 ```
 
 Now compare pertinent elements of `c_pol` and `F`
 
-```{code-cell} python3
+```{code-cell} ipython3
 print(c_pol, "\n", -F)
 ```
 
@@ -501,7 +501,7 @@ A second graph  plots a collection of simulations against the population distrib
 
 Comparing sample paths with population distributions at each date $t$ is a useful exercise---see {ref}`our discussion <lln_mr>` of the laws of large numbers
 
-```{code-cell} python3
+```{code-cell} ipython3
 lss = qe.LinearStateSpace(A_LSS, C_LSS, G_LSS, mu_0=μ_0, Sigma_0=Σ_0)
 ```
 
@@ -514,7 +514,7 @@ In the code below, we use the [LinearStateSpace](https://github.com/QuantEcon/Qu
 - simulate a group of 25 consumers and plot sample paths on the same
   graph as the population distribution.
 
-```{code-cell} python3
+```{code-cell} ipython3
 def income_consumption_debt_series(A, C, G, μ_0, Σ_0, T=150, npaths=25):
     """
     This function takes initial conditions (μ_0, Σ_0) and uses the
@@ -545,8 +545,8 @@ def income_consumption_debt_series(A, C, G, μ_0, Σ_0, T=150, npaths=25):
     debt_var = np.empty(T)
     for t in range(T):
         μ_x, μ_y, Σ_x, Σ_y = next(moment_generator)
-        cons_mean[t], cons_var[t] = μ_y[1], Σ_y[1, 1]
-        debt_mean[t], debt_var[t] = μ_x[3], Σ_x[3, 3]
+        cons_mean[t], cons_var[t] = μ_y[1,0], Σ_y[1, 1]
+        debt_mean[t], debt_var[t] = μ_x[3,0], Σ_x[3, 3]
 
     return bsim, csim, ysim, cons_mean, cons_var, debt_mean, debt_var
 
@@ -622,7 +622,7 @@ def consumption_debt_fanchart(csim, cons_mean, cons_var,
 
 Now let's create figures with initial conditions of zero for $y_0$ and $b_0$
 
-```{code-cell} python3
+```{code-cell} ipython3
 out = income_consumption_debt_series(A_LSS, C_LSS, G_LSS, μ_0, Σ_0)
 bsim0, csim0, ysim0 = out[:3]
 cons_mean0, cons_var0, debt_mean0, debt_var0 = out[3:]
@@ -632,7 +632,7 @@ consumption_income_debt_figure(bsim0, csim0, ysim0)
 plt.show()
 ```
 
-```{code-cell} python3
+```{code-cell} ipython3
 consumption_debt_fanchart(csim0, cons_mean0, cons_var0,
                           bsim0, debt_mean0, debt_var0)
 
@@ -698,7 +698,7 @@ behavior early in the sample.
 
 By altering initial conditions, we shall remove this transient in our second example to be presented below
 
-```{code-cell} python3
+```{code-cell} ipython3
 def cointegration_figure(bsim, csim):
     """
     Plots the cointegration
@@ -713,7 +713,7 @@ def cointegration_figure(bsim, csim):
     return fig
 ```
 
-```{code-cell} python3
+```{code-cell} ipython3
 cointegration_figure(bsim0, csim0)
 plt.show()
 ```
@@ -756,7 +756,7 @@ There is no need for foreigners to lend to our group.
 
 Let's have a look at the corresponding figures
 
-```{code-cell} python3
+```{code-cell} ipython3
 out = income_consumption_debt_series(A_LSS, C_LSS, G_LSS, mxbewley, sxbewley)
 bsimb, csimb, ysimb = out[:3]
 cons_meanb, cons_varb, debt_meanb, debt_varb = out[3:]
@@ -766,7 +766,7 @@ consumption_income_debt_figure(bsimb, csimb, ysimb)
 plt.show()
 ```
 
-```{code-cell} python3
+```{code-cell} ipython3
 consumption_debt_fanchart(csimb, cons_meanb, cons_varb,
                           bsimb, debt_meanb, debt_varb)
 
@@ -785,7 +785,7 @@ But now there is some initial dispersion because there is *ex-ante* heterogeneit
 
 Let's have a look at the cointegration figure
 
-```{code-cell} python3
+```{code-cell} ipython3
 cointegration_figure(bsimb, csimb)
 plt.show()
 ```
diff --git a/lectures/status.md b/lectures/status.md
index fd319bf08..ca1d9dd2e 100644
--- a/lectures/status.md
+++ b/lectures/status.md
@@ -33,7 +33,15 @@ and the following package versions
 !conda list
 ```
 
-This lecture series also has access to the following GPU
+You can check the backend used by JAX using:
+
+```{code-cell} ipython3
+import jax
+# Check if JAX is using GPU
+print(f"JAX backend: {jax.devices()[0].platform}")
+```
+
+and this lecture series also has access to the following GPU
 
 ```{code-cell} ipython
 !nvidia-smi