diff --git a/lectures/dynamic_programming_squared/amss.md b/lectures/dynamic_programming_squared/amss.md
index d5dbb1ef..5ae14997 100644
--- a/lectures/dynamic_programming_squared/amss.md
+++ b/lectures/dynamic_programming_squared/amss.md
@@ -44,11 +44,6 @@ In this lecture, we
 We begin with an introduction to the model.
 
 
-
-```{code-cell} julia
-using LinearAlgebra, Statistics
-```
-
 ## Competitive Equilibrium with Distorting Taxes
 
 Many but not all features of the economy are identical to those of {doc}`the Lucas-Stokey economy <../dynamic_programming_squared/opt_tax_recur>`.
@@ -390,13 +385,15 @@ on optimal taxation with state-contingent debt  sequential allocation implementa
 ---
 tags: [remove-cell]
 ---
-using Test, Random
+using Test
 ```
 
 ```{code-cell} julia
 ---
 tags: [output_scroll]
 ---
+
+using LinearAlgebra, Statistics, Random
 using QuantEcon, NLsolve, NLopt
 
 import QuantEcon.simulate
@@ -417,7 +414,6 @@ mutable struct Model{TF <: AbstractFloat,
     n_less_than_one::Bool
 end
 
-
 struct SequentialAllocation{TP <: Model,
                             TI <: Integer,
                             TV <: AbstractVector}
@@ -430,7 +426,6 @@ struct SequentialAllocation{TP <: Model,
     zFB::TV
 end
 
-
 function SequentialAllocation(model::Model)
     beta, Pi, G, Theta = model.beta, model.Pi, model.G, model.Theta
     mc = MarkovChain(Pi)
@@ -962,6 +957,7 @@ The recursive formulation is implemented as follows
 
 ```{code-cell} julia
 
+# Interpolations.jl doesn't support irregular grids for splines
 using DataInterpolations
 
 mutable struct BellmanEquation_Recursive{TP <: Model, TI <: Integer, TR <: Real}
@@ -976,11 +972,10 @@ mutable struct BellmanEquation_Recursive{TP <: Model, TI <: Integer, TR <: Real}
     zFB::Vector{Vector{TR}}
 end
 
-
 struct RecursiveAllocation{TP <: Model,
-                        TI <: Integer,
-                        TVg <: AbstractVector,
-                        TT <: Tuple}
+                           TI <: Integer,
+                           TVg <: AbstractVector,
+                           TT <: Tuple}
     model::TP
     mc::MarkovChain
     S::TI
@@ -991,7 +986,6 @@ struct RecursiveAllocation{TP <: Model,
     policies::TT
 end
 
-
 function RecursiveAllocation(model::Model, mugrid::AbstractArray)
     G = model.G
     S = size(model.Pi, 1)             # number of states
@@ -1002,8 +996,8 @@ function RecursiveAllocation(model::Model, mugrid::AbstractArray)
     return RecursiveAllocation(model, mc, S, T, mugrid, xgrid, Vf, policies)
 end
 
-
-function solve_time1_bellman(model::Model{TR}, mugrid::AbstractArray) where {TR <: Real}
+function solve_time1_bellman(model::Model{TR},
+                             mugrid::AbstractArray) where {TR <: Real}
     Pi = model.Pi
     S = size(model.Pi, 1)
 
@@ -1012,8 +1006,8 @@ function solve_time1_bellman(model::Model{TR}, mugrid::AbstractArray) where {TR
     PP = SequentialAllocation(model)
 
     function incomplete_allocation(PP::SequentialAllocation,
-                                mu_::AbstractFloat,
-                                s_::Integer)
+                                   mu_::AbstractFloat,
+                                   s_::Integer)
         c, n, x, V = time1_value(PP, mu_)
         return c, n, dot(Pi[s_, :], x), dot(Pi[s_, :], V)
     end
@@ -1030,20 +1024,24 @@ function solve_time1_bellman(model::Model{TR}, mugrid::AbstractArray) where {TR
         x = Array{TR}(undef, length(mugrid))
         V = Array{TR}(undef, length(mugrid))
         for (i_mu, mu) in enumerate(mugrid)
-            c[i_mu, :], n[i_mu, :], x[i_mu], V[i_mu] =
-                incomplete_allocation(PP, mu, s_)
+            c[i_mu, :], n[i_mu, :], x[i_mu], V[i_mu] = incomplete_allocation(PP,
+                                                                             mu,
+                                                                             s_)
         end
         xprimes = repeat(x, 1, S)
         xgrid[s_, :] = x
-        for sprime = 1:S
-            splc = CubicSpline(c[:, sprime][end:-1:1], x[end:-1:1];extrapolate=true)
-            spln = CubicSpline(n[:, sprime][end:-1:1], x[end:-1:1];extrapolate=true)
-            splx = CubicSpline(xprimes[:, sprime][end:-1:1], x[end:-1:1];extrapolate=true)
+        for sprime in 1:S
+            splc = CubicSpline(c[:, sprime][end:-1:1], x[end:-1:1];
+                               extrapolate = true)
+            spln = CubicSpline(n[:, sprime][end:-1:1], x[end:-1:1];
+                               extrapolate = true)
+            splx = CubicSpline(xprimes[:, sprime][end:-1:1], x[end:-1:1];
+                               extrapolate = true)
             cf[s_, sprime] = y -> splc(y)
             nf[s_, sprime] = y -> spln(y)
             xprimef[s_, sprime] = y -> splx(y)
         end
-        splV = CubicSpline(V[end:-1:1], x[end:-1:1];extrapolate=true)
+        splV = CubicSpline(V[end:-1:1], x[end:-1:1]; extrapolate = true)
         Vf[s_] = y -> splV(y)
     end
 
@@ -1061,9 +1059,11 @@ function solve_time1_bellman(model::Model{TR}, mugrid::AbstractArray) where {TR
         Vfnew, policies = fit_policy_function(T, PF, xgrid)
 
         diff = 0.0
-        for s=1:S
-            diff = max(diff, maximum(abs, (Vf[s].(xgrid) - Vfnew[s].(xgrid)) ./
-                                            Vf[s].(xgrid)))
+        for s in 1:S
+            diff = max(diff,
+                       maximum(abs,
+                               (Vf[s].(xgrid) - Vfnew[s].(xgrid)) ./
+                               Vf[s].(xgrid)))
         end
 
         println("diff = $diff")
@@ -1073,10 +1073,11 @@ function solve_time1_bellman(model::Model{TR}, mugrid::AbstractArray) where {TR
     return Vf, policies, T, xgrid
 end
 
-
 function fit_policy_function(T::BellmanEquation_Recursive,
                              PF::Function,
-                             xgrid::AbstractVector{TF}) where {TF <: AbstractFloat}
+                             xgrid::AbstractVector{TF}) where {
+                                                               TF <:
+                                                               AbstractFloat}
     S = T.S
     # preallocation
     PFvec = Array{TF}(undef, 4S + 1, length(xgrid))
@@ -1090,9 +1091,9 @@ function fit_policy_function(T::BellmanEquation_Recursive,
         for (i_x, x) in enumerate(xgrid)
             PFvec[:, i_x] = PF(i_x, x, s_)
         end
-        splV = CubicSpline(PFvec[1,:], xgrid)
+        splV = CubicSpline(PFvec[1, :], xgrid)
         Vf[s_] = y -> splV(y)
-        for sprime=1:S
+        for sprime in 1:S
             splc = CubicSpline(PFvec[1 + sprime, :], xgrid)
             spln = CubicSpline(PFvec[1 + S + sprime, :], xgrid)
             splxprime = CubicSpline(PFvec[1 + 2S + sprime, :], xgrid)
@@ -1107,18 +1108,16 @@ function fit_policy_function(T::BellmanEquation_Recursive,
     return Vf, policies
 end
 
-
 function Tau(pab::RecursiveAllocation,
-            c::AbstractArray,
-        n::AbstractArray)
+             c::AbstractArray,
+             n::AbstractArray)
     model = pab.model
     Uc, Un = model.Uc(c, n), model.Un(c, n)
-    return 1. .+ Un ./ (model.Theta .* Uc)
+    return 1.0 .+ Un ./ (model.Theta .* Uc)
 end
 
 Tau(pab::RecursiveAllocation, c::Real, n::Real) = Tau(pab, [c], [n])
 
-
 function time0_allocation(pab::RecursiveAllocation, B_::Real, s0::Integer)
     T, Vf = pab.T, pab.Vf
     xbar = T.xbar
@@ -1128,10 +1127,12 @@ function time0_allocation(pab::RecursiveAllocation, B_::Real, s0::Integer)
     return c0, n0, xprime0, T0
 end
 
-
 function simulate(pab::RecursiveAllocation,
                   B_::TF, s_0::Integer, T::Integer,
-                  sHist::Vector=simulate(pab.mc, T, init=s_0)) where {TF <: AbstractFloat}
+                  sHist::Vector = simulate(pab.mc, T, init = s_0)) where {
+                                                                          TF <:
+                                                                          AbstractFloat
+                                                                          }
     model, mc, Vf, S = pab.model, pab.mc, pab.Vf, pab.S
     Pi, Uc = model.Pi, model.Uc
     cf, nf, xprimef, TTf = pab.policies
@@ -1145,22 +1146,24 @@ function simulate(pab::RecursiveAllocation,
     muHist = Array{TF}(undef, T)
 
     #time0
-    cHist[1], nHist[1], xHist[1], THist[1]  = time0_allocation(pab, B_, s_0)
+    cHist[1], nHist[1], xHist[1], THist[1] = time0_allocation(pab, B_, s_0)
     TauHist[1] = Tau(pab, cHist[1], nHist[1])[s_0]
     Bhist[1] = B_
     muHist[1] = Vf[s_0](xHist[1])
 
     #time 1 onward
     for t in 2:T
-        s_, x, s = sHist[t-1], xHist[t-1], sHist[t]
+        s_, x, s = sHist[t - 1], xHist[t - 1], sHist[t]
         c = Array{TF}(undef, S)
         n = Array{TF}(undef, S)
         xprime = Array{TF}(undef, S)
         TT = Array{TF}(undef, S)
-        for sprime=1:S
-            c[sprime], n[sprime], xprime[sprime], TT[sprime] =
-                cf[s_, sprime](x), nf[s_, sprime](x),
-                xprimef[s_, sprime](x), TTf[s_, sprime](x)
+        for sprime in 1:S
+            c[sprime], n[sprime], xprime[sprime], TT[sprime] = cf[s_, sprime](x),
+                                                               nf[s_, sprime](x),
+                                                               xprimef[s_,
+                                                                       sprime](x),
+                                                               TTf[s_, sprime](x)
         end
 
         Tau_val = Tau(pab, c, n)[s]
@@ -1169,17 +1172,15 @@ function simulate(pab::RecursiveAllocation,
 
         muHist[t] = Vf[s](xprime[s])
 
-        cHist[t], nHist[t], Bhist[t], TauHist[t] = c[s], n[s], x/Eu_c, Tau_val
+        cHist[t], nHist[t], Bhist[t], TauHist[t] = c[s], n[s], x / Eu_c, Tau_val
         xHist[t], THist[t] = xprime[s], TT[s]
     end
     return cHist, nHist, Bhist, xHist, TauHist, THist, muHist, sHist
 end
 
-
 function BellmanEquation_Recursive(model::Model{TF},
                                    xgrid::AbstractVector{TF},
                                    policies0::Array) where {TF <: AbstractFloat}
-
     S = size(model.Pi, 1) # number of states
     xbar = [minimum(xgrid), maximum(xgrid)]
     time_0 = false
@@ -1190,18 +1191,18 @@ function BellmanEquation_Recursive(model::Model{TF},
             cs = Array{TF}(undef, S)
             ns = Array{TF}(undef, S)
             xprimes = Array{TF}(undef, S)
-            for j = 1:S
-                cs[j], ns[j], xprimes[j] = cf[s, j](x), nf[s, j](x), xprimef[s, j](x)
+            for j in 1:S
+                cs[j], ns[j], xprimes[j] = cf[s, j](x), nf[s, j](x),
+                                           xprimef[s, j](x)
             end
             z0[i_x, s] = vcat(cs, ns, xprimes, zeros(S))
         end
     end
     cFB, nFB, IFB, xFB, zFB = find_first_best(model, S, 2)
-    return BellmanEquation_Recursive(model, S, xbar, time_0, z0, cFB, nFB, xFB, zFB)
+    return BellmanEquation_Recursive(model, S, xbar, time_0, z0, cFB, nFB, xFB,
+                                     zFB)
 end
 
-
-
 function get_policies_time1(T::BellmanEquation_Recursive,
                             i_x::Integer,
                             x::Real,
@@ -1210,28 +1211,28 @@ function get_policies_time1(T::BellmanEquation_Recursive,
                             xbar::AbstractVector)
     model, S = T.model, T.S
     beta, Theta, G, Pi = model.beta, model.Theta, model.G, model.Pi
-    U,Uc,Un = model.U, model.Uc, model.Un
+    U, Uc, Un = model.U, model.Uc, model.Un
 
-    S_possible = sum(Pi[s_, :].>0)
-    sprimei_possible = findall(Pi[s_, :].>0)
+    S_possible = sum(Pi[s_, :] .> 0)
+    sprimei_possible = findall(Pi[s_, :] .> 0)
 
     function objf(z, grad)
-        c, xprime = z[1:S_possible], z[S_possible+1:2S_possible]
+        c, xprime = z[1:S_possible], z[(S_possible + 1):(2S_possible)]
         n = (c .+ G[sprimei_possible]) ./ Theta[sprimei_possible]
         Vprime = [Vf[sprimei_possible[si]](xprime[si]) for si in 1:S_possible]
         return -dot(Pi[s_, sprimei_possible], U.(c, n) + beta * Vprime)
     end
 
     function cons(out, z, grad)
-        c, xprime, TT =
-            z[1:S_possible], z[S_possible + 1:2S_possible], z[2S_possible + 1:3S_possible]
+        c, xprime, TT = z[1:S_possible], z[(S_possible + 1):(2S_possible)],
+                        z[(2S_possible + 1):(3S_possible)]
         n = (c .+ G[sprimei_possible]) ./ Theta[sprimei_possible]
         u_c = Uc.(c, n)
         Eu_c = dot(Pi[s_, sprimei_possible], u_c)
-        out .= x * u_c/Eu_c - u_c .* (c - TT) - Un(c, n) .* n - beta * xprime
+        out .= x * u_c / Eu_c - u_c .* (c - TT) - Un(c, n) .* n - beta * xprime
     end
     function cons_no_trans(out, z, grad)
-        c, xprime = z[1:S_possible], z[S_possible + 1:2S_possible]
+        c, xprime = z[1:S_possible], z[(S_possible + 1):(2S_possible)]
         n = (c .+ G[sprimei_possible]) ./ Theta[sprimei_possible]
         u_c = Uc.(c, n)
         Eu_c = dot(Pi[s_, sprimei_possible], u_c)
@@ -1239,14 +1240,15 @@ function get_policies_time1(T::BellmanEquation_Recursive,
     end
 
     if model.transfers == true
-        lb = vcat(zeros(S_possible), ones(S_possible)*xbar[1], zeros(S_possible))
+        lb = vcat(zeros(S_possible), ones(S_possible) * xbar[1],
+                  zeros(S_possible))
         if model.n_less_than_one == true
             ub = vcat(ones(S_possible) - G[sprimei_possible],
-                    ones(S_possible) * xbar[2], ones(S_possible))
+                      ones(S_possible) * xbar[2], ones(S_possible))
         else
             ub = vcat(100 * ones(S_possible),
-                    ones(S_possible) * xbar[2],
-                    100 * ones(S_possible))
+                      ones(S_possible) * xbar[2],
+                      100 * ones(S_possible))
         end
         init = vcat(T.z0[i_x, s_][sprimei_possible],
                     T.z0[i_x, s_][2S .+ sprimei_possible],
@@ -1254,9 +1256,10 @@ function get_policies_time1(T::BellmanEquation_Recursive,
         opt = Opt(:LN_COBYLA, 3S_possible)
         equality_constraint!(opt, cons, zeros(S_possible))
     else
-        lb = vcat(zeros(S_possible), ones(S_possible)*xbar[1])
+        lb = vcat(zeros(S_possible), ones(S_possible) * xbar[1])
         if model.n_less_than_one == true
-            ub = vcat(ones(S_possible)-G[sprimei_possible], ones(S_possible)*xbar[2])
+            ub = vcat(ones(S_possible) - G[sprimei_possible],
+                      ones(S_possible) * xbar[2])
         else
             ub = vcat(ones(S_possible), ones(S_possible) * xbar[2])
         end
@@ -1279,15 +1282,16 @@ function get_policies_time1(T::BellmanEquation_Recursive,
     (minf, minx, ret) = NLopt.optimize(opt, init)
 
     if ret != :SUCCESS && ret != :ROUNDOFF_LIMITED && ret != :MAXEVAL_REACHED &&
-        ret != :FTOL_REACHED && ret != :MAXTIME_REACHED
+       ret != :FTOL_REACHED && ret != :MAXTIME_REACHED
         error("optimization failed: ret = $ret")
     end
 
     T.z0[i_x, s_][sprimei_possible] = minx[1:S_possible]
-    T.z0[i_x, s_][S .+ sprimei_possible] = minx[1:S_possible] .+ G[sprimei_possible]
-    T.z0[i_x, s_][2S .+ sprimei_possible] = minx[S_possible .+ 1:2S_possible]
+    T.z0[i_x, s_][S .+ sprimei_possible] = minx[1:S_possible] .+
+                                           G[sprimei_possible]
+    T.z0[i_x, s_][2S .+ sprimei_possible] = minx[(S_possible .+ 1):(2S_possible)]
     if model.transfers == true
-        T.z0[i_x, s_][3S .+ sprimei_possible] = minx[2S_possible + 1:3S_possible]
+        T.z0[i_x, s_][3S .+ sprimei_possible] = minx[(2S_possible + 1):(3S_possible)]
     else
         T.z0[i_x, s_][3S .+ sprimei_possible] = zeros(S)
     end
@@ -1295,7 +1299,6 @@ function get_policies_time1(T::BellmanEquation_Recursive,
     return vcat(-minf, T.z0[i_x, s_])
 end
 
-
 function get_policies_time0(T::BellmanEquation_Recursive,
                             B_::Real,
                             s0::Integer,
@@ -1311,7 +1314,7 @@ function get_policies_time0(T::BellmanEquation_Recursive,
         return -(U(c, n) + beta * Vf[s0](xprime))
     end
 
-    function cons(z,grad)
+    function cons(z, grad)
         c, xprime, TT = z[1], z[2], z[3]
         n = (c + G[s0]) / Theta[s0]
         return -Uc(c, n) * (c - B_ - TT) - Un(c, n) * n - beta * xprime
@@ -1326,13 +1329,13 @@ function get_policies_time0(T::BellmanEquation_Recursive,
             ub = [100.0, xbar[2], 100.0]
         end
         init = vcat(T.zFB[s0][1], T.zFB[s0][3], T.zFB[s0][4])
-        init = [0.95124922, -1.15926816,  0.0]
+        init = [0.95124922, -1.15926816, 0.0]
         opt = Opt(:LN_COBYLA, 3)
         equality_constraint!(opt, cons)
     else
         lb = [0.0, xbar[1]]
         if model.n_less_than_one == true
-            ub = [1-G[s0], xbar[2]]
+            ub = [1 - G[s0], xbar[2]]
         else
             ub = [100, xbar[2]]
         end
@@ -1344,7 +1347,6 @@ function get_policies_time0(T::BellmanEquation_Recursive,
     init[init .> ub] = ub[init .> ub]
     init[init .< lb] = lb[init .< lb]
 
-
     min_objective!(opt, objf)
     lower_bounds!(opt, lb)
     upper_bounds!(opt, ub)
@@ -1354,14 +1356,14 @@ function get_policies_time0(T::BellmanEquation_Recursive,
     (minf, minx, ret) = NLopt.optimize(opt, init)
 
     if ret != :SUCCESS && ret != :ROUNDOFF_LIMITED && ret != :MAXEVAL_REACHED &&
-        ret != :FTOL_REACHED
-            error("optimization failed: ret = $ret")
+       ret != :FTOL_REACHED
+        error("optimization failed: ret = $ret")
     end
 
     if model.transfers == true
-        return -minf, minx[1], minx[1]+G[s0], minx[2], minx[3]
+        return -minf, minx[1], minx[1] + G[s0], minx[2], minx[3]
     else
-        return -minf, minx[1], minx[1]+G[s0], minx[2], 0
+        return -minf, minx[1], minx[1] + G[s0], minx[2], 0
     end
 end
 ```
@@ -1429,31 +1431,30 @@ We assume the same utility parameters as in the {doc}`Lucas-Stokey economy <../d
 This utility function is implemented in the following constructor
 
 ```{code-cell} julia
-function crra_utility(;
-    beta = 0.9,
-    sigma = 2.0,
-    gamma = 2.0,
-    Pi = 0.5 * ones(2, 2),
-    G = [0.1, 0.2],
-    Theta = ones(Float64, 2),
-    transfers = false
-    )
+function CRRAModel(;
+                   beta = 0.9,
+                   sigma = 2.0,
+                   gamma = 2.0,
+                   Pi = 0.5 * ones(2, 2),
+                   G = [0.1, 0.2],
+                   Theta = ones(Float64, 2),
+                   transfers = false)
     function U(c, n)
         if sigma == 1.0
             U = log(c)
         else
-            U = (c.^(1.0 - sigma) - 1.0) / (1.0 - sigma)
+            U = (c .^ (1.0 - sigma) - 1.0) / (1.0 - sigma)
         end
-        return U - n.^(1 + gamma) / (1 + gamma)
+        return U - n .^ (1 + gamma) / (1 + gamma)
     end
     # Derivatives of utility function
-    Uc(c,n) =  c.^(-sigma)
-    Ucc(c,n) = -sigma * c.^(-sigma - 1.0)
-    Un(c,n) = -n.^gamma
-    Unn(c,n) = -gamma * n.^(gamma - 1.0)
+    Uc(c, n) = c .^ (-sigma)
+    Ucc(c, n) = -sigma * c .^ (-sigma - 1.0)
+    Un(c, n) = -n .^ gamma
+    Unn(c, n) = -gamma * n .^ (gamma - 1.0)
     n_less_than_one = false
     return Model(beta, Pi, G, Theta, transfers,
-                U, Uc, Ucc, Un, Unn, n_less_than_one)
+                 U, Uc, Ucc, Un, Unn, n_less_than_one)
 end
 ```
 
@@ -1468,10 +1469,10 @@ Paths with circles are histories in which there is peace, while those with
 triangle denote war.
 
 ```{code-cell} julia
-time_example = crra_utility(G=[0.1, 0.1, 0.1, 0.2, 0.1, 0.1],
-                            Theta = ones(6)) # Theta can in principle be random
+time_example = CRRAModel(;G = [0.1, 0.1, 0.1, 0.2, 0.1, 0.1],
+                         Theta = ones(6)) # Theta can in principle be random
 
-time_example.Pi = [ 0.0 1.0 0.0 0.0 0.0 0.0;
+time_example.Pi = [0.0 1.0 0.0 0.0 0.0 0.0;
                    0.0 0.0 1.0 0.0 0.0 0.0;
                    0.0 0.0 0.0 0.5 0.5 0.0;
                    0.0 0.0 0.0 0.0 0.0 1.0;
@@ -1489,10 +1490,10 @@ time_bellman = RecursiveAllocation(time_example, mugrid)
 sHist_h = [1, 2, 3, 4, 6, 6, 6]
 sHist_l = [1, 2, 3, 5, 6, 6, 6]
 
-sim_seq_h = simulate(time_sequential, 1., 1, 7, sHist_h)
-sim_bel_h = simulate(time_bellman, 1., 1, 7, sHist_h)
-sim_seq_l = simulate(time_sequential, 1., 1, 7, sHist_l)
-sim_bel_l = simulate(time_bellman, 1., 1, 7, sHist_l)
+sim_seq_h = simulate(time_sequential, 1.0, 1, 7, sHist_h)
+sim_bel_h = simulate(time_bellman, 1.0, 1, 7, sHist_h)
+sim_seq_l = simulate(time_sequential, 1.0, 1, 7, sHist_l)
+sim_bel_l = simulate(time_bellman, 1.0, 1, 7, sHist_l)
 
 using Plots
 
@@ -1510,14 +1511,18 @@ sim_seq_h_plot = hcat(sim_seq_h[1:3]..., sim_seq_h[4],
 sim_bel_h_plot = hcat(sim_bel_h[1:3]..., sim_bel_h[5],
                       time_example.G[sHist_h],
                       time_example.Theta[sHist_h] .* sim_bel_h[2])
-p = plot(size = (920, 750), layout =(3, 2),
- xaxis=(0:6), grid=false, titlefont=Plots.font("sans-serif", 10))
+p = plot(size = (920, 750), layout = (3, 2),
+         xaxis = (0:6), grid = false, titlefont = Plots.font("sans-serif", 10))
 plot!(p, title = titles)
-for i=1:6
-    plot!(p[i], 0:6, sim_seq_l_plot[:, i], marker=:circle, color=:black, lab="")
-    plot!(p[i], 0:6, sim_bel_l_plot[:, i], marker=:circle, color=:red, lab="")
-    plot!(p[i], 0:6, sim_seq_h_plot[:, i], marker=:utriangle, color=:black, lab="")
-    plot!(p[i], 0:6, sim_bel_h_plot[:, i], marker=:utriangle, color=:red, lab="")
+for i in 1:6
+    plot!(p[i], 0:6, sim_seq_l_plot[:, i], marker = :circle, color = :black,
+          lab = "")
+    plot!(p[i], 0:6, sim_bel_l_plot[:, i], marker = :circle, color = :red,
+          lab = "")
+    plot!(p[i], 0:6, sim_seq_h_plot[:, i], marker = :utriangle, color = :black,
+          lab = "")
+    plot!(p[i], 0:6, sim_bel_h_plot[:, i], marker = :utriangle, color = :red,
+          lab = "")
 end
 p
 ```
@@ -1574,21 +1579,21 @@ $$
 In accordance, we will re-define our utility function
 
 ```{code-cell} julia
-function log_utility(;beta = 0.9,
-                    psi = 0.69,
-                    Pi = 0.5 * ones(2, 2),
-                    G = [0.1, 0.2],
-                    Theta = ones(2),
-                    transfers = false)
+function log_utility(; beta = 0.9,
+                     psi = 0.69,
+                     Pi = 0.5 * ones(2, 2),
+                     G = [0.1, 0.2],
+                     Theta = ones(2),
+                     transfers = false)
     # Derivatives of utility function
-    U(c,n) = log(c) + psi * log(1 - n)
-    Uc(c,n) = 1 ./ c
-    Ucc(c,n) = -c.^(-2.0)
-    Un(c,n) = -psi ./ (1.0 .- n)
-    Unn(c,n) = -psi ./ (1.0 .- n).^2.0
+    U(c, n) = log(c) + psi * log(1 - n)
+    Uc(c, n) = 1 ./ c
+    Ucc(c, n) = -c .^ (-2.0)
+    Un(c, n) = -psi ./ (1.0 .- n)
+    Unn(c, n) = -psi ./ (1.0 .- n) .^ 2.0
     n_less_than_one = true
     return Model(beta, Pi, G, Theta, transfers,
-                U, Uc, Ucc, Un, Unn, n_less_than_one)
+                 U, Uc, Ucc, Un, Unn, n_less_than_one)
 end
 ```
 
@@ -1614,22 +1619,26 @@ sim_seq = simulate(log_sequential, 0.5, 1, T, sHist)
 sim_bel = simulate(log_bellman, 0.5, 1, T, sHist)
 
 sim_seq_plot = hcat(sim_seq[1:3]...,
-            sim_seq[4], log_example.G[sHist], log_example.Theta[sHist] .* sim_seq[2])
+                    sim_seq[4], log_example.G[sHist],
+                    log_example.Theta[sHist] .* sim_seq[2])
 sim_bel_plot = hcat(sim_bel[1:3]...,
-            sim_bel[5], log_example.G[sHist], log_example.Theta[sHist] .* sim_bel[2])
+                    sim_bel[5], log_example.G[sHist],
+                    log_example.Theta[sHist] .* sim_bel[2])
 
 #plot policies
 p = plot(size = (920, 750), layout = grid(3, 2),
-         xaxis=(0:T), grid=false, titlefont=Plots.font("sans-serif", 10))
+         xaxis = (0:T), grid = false, titlefont = Plots.font("sans-serif", 10))
 labels = fill(("", ""), 6)
 labels[3] = ("Complete Market", "Incomplete Market")
 plot!(p, title = titles)
-for i = vcat(collect(1:4), 6)
-    plot!(p[i], sim_seq_plot[:, i], marker=:circle, color=:black, lab=labels[i][1])
-    plot!(p[i], sim_bel_plot[:, i], marker=:utriangle, color=:blue, lab=labels[i][2],
-          legend=:bottomright)
+for i in vcat(collect(1:4), 6)
+    plot!(p[i], sim_seq_plot[:, i], marker = :circle, color = :black,
+          lab = labels[i][1])
+    plot!(p[i], sim_bel_plot[:, i], marker = :utriangle, color = :blue,
+          lab = labels[i][2],
+          legend = :bottomright)
 end
-plot!(p[5], sim_seq_plot[:, 5], marker=:circle, color=:blue, lab="")
+plot!(p[5], sim_seq_plot[:, 5], marker = :circle, color = :blue, lab = "")
 ```
 
 When the government experiences a prolonged period of peace, it is able to reduce
@@ -1653,38 +1662,29 @@ Random.seed!(42)
 ```{code-cell} julia
 T_long = 200
 sim_seq_long = simulate(log_sequential, 0.5, 1, T_long)
-sHist_long = sim_seq_long[end-2]
+sHist_long = sim_seq_long[end - 2]
 sim_bel_long = simulate(log_bellman, 0.5, 1, T_long, sHist_long)
 sim_seq_long_plot = hcat(sim_seq_long[1:4]...,
-             log_example.G[sHist_long], log_example.Theta[sHist_long] .* sim_seq_long[2])
+                         log_example.G[sHist_long],
+                         log_example.Theta[sHist_long] .* sim_seq_long[2])
 sim_bel_long_plot = hcat(sim_bel_long[1:3]..., sim_bel_long[5],
-             log_example.G[sHist_long], log_example.Theta[sHist_long] .* sim_bel_long[2])
+                         log_example.G[sHist_long],
+                         log_example.Theta[sHist_long] .* sim_bel_long[2])
 
-p = plot(size = (920, 750), layout = (3, 2), xaxis=(0:50:T_long), grid=false,
-         titlefont=Plots.font("sans-serif", 10))
+p = plot(size = (920, 750), layout = (3, 2), xaxis = (0:50:T_long),
+         grid = false,
+         titlefont = Plots.font("sans-serif", 10))
 plot!(p, title = titles)
-for i = 1:6
-    plot!(p[i], sim_seq_long_plot[:, i], color=:black, linestyle=:solid, lab=labels[i][1])
-    plot!(p[i], sim_bel_long_plot[:, i], color=:blue, linestyle=:dot, lab=labels[i][2],
-          legend=:bottomright)
+for i in 1:6
+    plot!(p[i], sim_seq_long_plot[:, i], color = :black, linestyle = :solid,
+          lab = labels[i][1])
+    plot!(p[i], sim_bel_long_plot[:, i], color = :blue, linestyle = :dot,
+          lab = labels[i][2],
+          legend = :bottomright)
 end
 p
 ```
 
-```{code-cell} julia
----
-tags: [remove-cell]
----
-@testset begin
-  # @test sim_seq_long_plot[50, 3] ≈ 0.3951985593686047
-  # @test sim_bel_long_plot[50, 3] ≈ 0.05684753244006188 atol = 1e-2
-  # @test sim_seq_long_plot[100, 4] ≈ 0.340233842670859
-  # @test sim_bel_long_plot[100, 4] ≈ 0.2093423366870517 atol = 1e-3
-  # @test sim_seq_long_plot[200, 2] ≈ 0.5839693539786998
-  # @test sim_bel_long_plot[200, 2] ≈ 0.6324036099550768 atol = 1e-3
-end
-```
-
 [^fn_a]: In an allocation that solves the Ramsey problem and that levies distorting
 taxes on labor, why would the government ever want to hand revenues back
 to the private sector? It would not in an economy with state-contingent debt, since
diff --git a/lectures/dynamic_programming_squared/dyn_stack.md b/lectures/dynamic_programming_squared/dyn_stack.md
index 36012503..4527a9d8 100644
--- a/lectures/dynamic_programming_squared/dyn_stack.md
+++ b/lectures/dynamic_programming_squared/dyn_stack.md
@@ -914,7 +914,8 @@ dynamic program as outlined above
 
 
 ```{code-cell} julia
-using LaTeXStrings, QuantEcon, Plots, LinearAlgebra, Statistics, Parameters, Random
+using LaTeXStrings, QuantEcon, Plots, LinearAlgebra, Statistics,
+      Random
 
 ```
 
@@ -929,10 +930,12 @@ We define named tuples and default values for the model and solver settings, and
 instantiate one copy of each
 
 ```{code-cell} julia
-model(;a0 = 10, a1 = 2, beta = 0.96, gamma = 120., n = 300) = (; a0, a1, beta, gamma, n)
+function model(; a0 = 10, a1 = 2, beta = 0.96, gamma = 120.0, n = 300)
+    return (; a0, a1, beta, gamma, n)
+end
 
 # things like tolerances, etc.
-settings(;tol0 = 1e-8,tol1 = 1e-16,tol2 = 1e-2) = (;tol0, tol1, tol2)
+settings(; tol0 = 1e-8, tol1 = 1e-16, tol2 = 1e-2) = (; tol0, tol1, tol2)
 
 defaultModel = model();
 defaultSettings = settings();
@@ -941,21 +944,26 @@ defaultSettings = settings();
 Now we can compute the actual policy using the LQ routine from QuantEcon.jl
 
 ```{code-cell} julia
-(;a0, a1, beta, gamma, n) = defaultModel
-(;tol0, tol1, tol2) = defaultSettings
+(; a0, a1, beta, gamma, n) = defaultModel
+(; tol0, tol1, tol2) = defaultSettings
 
-betas = [beta^x for x = 0:n-1]
+betas = [beta^x for x in 0:(n - 1)]
 Alhs = I + zeros(4, 4);
-Alhs[4, :] = [beta * a0 / (2 * gamma), -beta * a1 / (2 * gamma), -beta * a1 / gamma, beta] # Euler equation coefficients
+Alhs[4, :] = [
+    beta * a0 / (2 * gamma),
+    -beta * a1 / (2 * gamma),
+    -beta * a1 / gamma,
+    beta,
+] # Euler equation coefficients
 Arhs = I + zeros(4, 4);
-Arhs[3, 4] = 1.;
+Arhs[3, 4] = 1.0;
 Alhsinv = inv(Alhs);
 
 A = Alhsinv * Arhs;
-B = Alhsinv * [0, 1, 0, 0,];
+B = Alhsinv * [0, 1, 0, 0];
 R = [0 -a0/2 0 0; -a0/2 a1 a1/2 0; 0 a1/2 0 0; 0 0 0 0];
 Q = gamma;
-lq = QuantEcon.LQ(Q, R, A, B, bet=beta);
+lq = QuantEcon.LQ(Q, R, A, B, bet = beta);
 P, F, d = stationary_values(lq)
 
 P22 = P[4:end, 4:end];
@@ -974,7 +982,7 @@ yt, ut = compute_sequence(lq, y0, n);
 pi_matrix = R + F' * Q * F;
 
 for t in 1:n
-    pi_leader[t] = -(yt[:, t]' * pi_matrix * yt[:, t]);
+    pi_leader[t] = -(yt[:, t]' * pi_matrix * yt[:, t])
 end
 
 println("Computed policy for Stackelberg leader: $F")
@@ -984,7 +992,7 @@ println("Computed policy for Stackelberg leader: $F")
 ---
 tags: [remove-cell]
 ---
-#test F ≈ [-1.580044538772657, 0.29461312747031404, 0.6748093760774972, 6.539705936147515]'
+@test F ≈ [-1.580044538772657, 0.29461312747031404, 0.6748093760774972, 6.539705936147515]'
 ```
 
 ### Implied Time Series for Price and Quantities
@@ -995,12 +1003,12 @@ The following code plots the price and quantities
 q_leader = yt[2, 1:end];
 q_follower = yt[3, 1:end];
 q = q_leader + q_follower;
-p = a0 .- a1*q;
+p = a0 .- a1 * q;
 
-plot(1:n+1, [q_leader, q_follower, p],
-    title = "Output and Prices, Stackelberg Duopoly",
-    labels = ["leader output" "follower output" "price"],
-    xlabel = L"t")
+plot(1:(n + 1), [q_leader, q_follower, p],
+     title = "Output and Prices, Stackelberg Duopoly",
+     labels = ["leader output" "follower output" "price"],
+     xlabel = L"t")
 ```
 
 ### Value of Stackelberg Leader
@@ -1022,8 +1030,8 @@ println("v_leader_direct is $v_leader_direct")
 ---
 tags: [remove-cell]
 ---
-#test v_leader_forward ≈ 150.0316212532547
-#test v_leader_direct ≈ 150.03237147548967
+@test v_leader_forward ≈ 150.0316212532547
+@test v_leader_direct ≈ 150.03237147548967
 ```
 
 ```{code-cell} julia
@@ -1036,14 +1044,14 @@ all(P - P_next .< tol0)
 ---
 tags: [remove-cell]
 ---
-#test all(P - P_next .< tol0)
+@test all(P - P_next .< tol0)
 ```
 
 ```{code-cell} julia
 # manually checks whether two different ways of computing the
 # value function give approximately the same answer
 v_expanded = -((y0' * R * y0 + ut[:, 1]' * Q * ut[:, 1] +
-           beta * (y0' * (A - B * F)' * P * (A - B * F) * y0)));
+                beta * (y0' * (A - B * F)' * P * (A - B * F) * y0)));
 (v_leader_direct - v_expanded < tol0)[1, 1]
 ```
 
@@ -1051,7 +1059,7 @@ v_expanded = -((y0' * R * y0 + ut[:, 1]' * Q * ut[:, 1] +
 ---
 tags: [remove-cell]
 ---
-#test (v_leader_direct - v_expanded < tol0)[1, 1]
+@test (v_leader_direct - v_expanded < tol0)[1, 1]
 ```
 
 ## Exhibiting Time Inconsistency of Stackelberg Plan
@@ -1079,11 +1087,13 @@ for t in 1:n
     vt_reset_leader[t] = -yt_reset[:, t]' * P * yt_reset[:, t]
 end
 
-p1 = plot(1:n+1, [(-F * yt)', (-F * yt_reset)'], labels = ["Stackelberg Leader" L"Continuation Leader at $t$"],
-        title = "Leader Control Variable", xlabel = L"t");
-p2 = plot(1:n+1, [yt[4, :], yt_reset[4, :]], title = "Follower Control Variable", xlabel = L"t", legend = false);
+p1 = plot(1:(n + 1), [(-F * yt)', (-F * yt_reset)'],
+          labels = ["Stackelberg Leader" L"Continuation Leader at $t$"],
+          title = "Leader Control Variable", xlabel = L"t");
+p2 = plot(1:(n + 1), [yt[4, :], yt_reset[4, :]],
+          title = "Follower Control Variable", xlabel = L"t", legend = false);
 p3 = plot(1:n, [vt_leader, vt_reset_leader], legend = false,
-            xlabel = L"t", title = "Leader Value Function");
+          xlabel = L"t", title = "Leader Value Function");
 plot(p1, p2, p3, layout = (3, 1), size = (800, 600))
 ```
 
@@ -1102,7 +1112,7 @@ R_tilde = [0 0 0 0 -a0/2; 0 0 0 0 a1/2; 0 0 0 0 0; 0 0 0 0 0; -a0/2 a1/2 0 0 a1]
 Q_tilde = Q;
 B_tilde = [0, 0, 0, 0, 1];
 
-lq_tilde = QuantEcon.LQ(Q_tilde, R_tilde, A_tilde, B_tilde, bet=beta);
+lq_tilde = QuantEcon.LQ(Q_tilde, R_tilde, A_tilde, B_tilde, bet = beta);
 P_tilde, F_tilde, d_tilde = stationary_values(lq_tilde);
 y0_tilde = vcat(y0, y0[3]);
 yt_tilde = compute_sequence(lq_tilde, y0_tilde, n)[1];
@@ -1111,7 +1121,7 @@ yt_tilde = compute_sequence(lq_tilde, y0_tilde, n)[1];
 ```{code-cell} julia
 # checks that the recursive formulation of the follower's problem gives
 # the same solution as the original Stackelberg problem
-plot(1:n+1, [yt_tilde[5, :], yt_tilde[3, :]], labels = [L"\tilde{q}" L"q"])
+plot(1:(n + 1), [yt_tilde[5, :], yt_tilde[3, :]], labels = [L"\tilde{q}" L"q"])
 ```
 
 Note: Variables with `_tilde` are obtained from solving the follower's
@@ -1126,7 +1136,7 @@ max(abs(yt_tilde[5] - yt_tilde[3]))
 ---
 tags: [remove-cell]
 ---
-#test max(abs(yt_tilde[5] - yt_tilde[3])) ≈ 0. atol = 1e-15
+@test max(abs(yt_tilde[5] - yt_tilde[3])) ≈ 0. atol = 1e-15
 ```
 
 ```{code-cell} julia
@@ -1138,7 +1148,7 @@ yt[:, 1][end] - (yt_tilde[:, 2] - yt_tilde[:, 1])[end] < tol0
 ---
 tags: [remove-cell]
 ---
-#test yt[:, 1][end] - (yt_tilde[:, 2] - yt_tilde[:, 1])[end] < tol0
+@test yt[:, 1][end] - (yt_tilde[:, 2] - yt_tilde[:, 1])[end] < tol0
 ```
 
 ### Explanation of Alignment
@@ -1173,7 +1183,7 @@ all((P - ((R + F' * Q * F) + beta * (A - B * F)' * P * (A - B * F)) .< tol0))
 ---
 tags: [remove-cell]
 ---
-#test all((P  - ((R + F' * Q * F) + beta * (A - B * F)' * P * (A - B * F)) .< tol0))
+@test all((P  - ((R + F' * Q * F) + beta * (A - B * F)' * P * (A - B * F)) .< tol0))
 ```
 
 ```{code-cell} julia
@@ -1184,7 +1194,7 @@ P_guess = zeros(5, 5);
 for i in 1:1000
     P_guess = ((R_tilde + F_tilde_star' * Q_tilde * F_tilde_star) +
                beta * (A_tilde - B_tilde * F_tilde_star)' * P_guess
-               * (A_tilde - B_tilde * F_tilde_star));
+               * (A_tilde - B_tilde * F_tilde_star))
 end
 ```
 
@@ -1215,7 +1225,7 @@ tags: [remove-cell]
 ```{code-cell} julia
 # c policy using policy iteration algorithm
 F_iter = (beta * inv(Q + beta * B_tilde' * P_guess * B_tilde)
-      * B_tilde' * P_guess * A_tilde);
+          * B_tilde' * P_guess * A_tilde);
 P_iter = zeros(5, 5);
 dist_vec = zeros(5, 5);
 
@@ -1223,27 +1233,30 @@ for i in 1:100
     # compute P_iter
     dist_vec = similar(P_iter)
     for j in 1:1000
-        P_iter = (R_tilde + F_iter' * Q * F_iter) + beta *
-                  (A_tilde - B_tilde * F_iter)' * P_iter *
-                  (A_tilde - B_tilde * F_iter);
+        P_iter = (R_tilde + F_iter' * Q * F_iter) +
+                 beta *
+                 (A_tilde - B_tilde * F_iter)' * P_iter *
+                 (A_tilde - B_tilde * F_iter)
 
         # update F_iter
         F_iter = beta * inv(Q + beta * B_tilde' * P_iter * B_tilde) *
-                    B_tilde' * P_iter * A_tilde;
+                 B_tilde' * P_iter * A_tilde
 
         dist_vec = P_iter - ((R_tilde + F_iter' * Q * F_iter) +
-                beta * (A_tilde - B_tilde * F_iter)' * P_iter *
-                (A_tilde - B_tilde * F_iter));
+                    beta * (A_tilde - B_tilde * F_iter)' * P_iter *
+                    (A_tilde - B_tilde * F_iter))
     end
 end
 
 if maximum(abs.(dist_vec)) < 1e-8
-    dist_vec2 = F_iter - (beta * inv(Q + beta * B_tilde' * P_iter * B_tilde) * B_tilde' * P_iter * A_tilde)
-        if maximum(abs.(dist_vec2)) < 1e-8
-            @show F_iter
-        else
-            println("The policy didn't converge: try increasing the number of outer loop iterations")
-        end
+    dist_vec2 = F_iter -
+                (beta * inv(Q + beta * B_tilde' * P_iter * B_tilde) * B_tilde' *
+                 P_iter * A_tilde)
+    if maximum(abs.(dist_vec2)) < 1e-8
+        @show F_iter
+    else
+        println("The policy didn't converge: try increasing the number of outer loop iterations")
+    end
 else
     println("The policy didn't converge: try increasing the number of inner loop iterations")
 end
@@ -1253,22 +1266,23 @@ end
 yt_tilde_star = zeros(n, 5);
 yt_tilde_star[1, :] = y0_tilde;
 
-for t in 1:n-1
-    yt_tilde_star[t+1, :] = (A_tilde - B_tilde * F_tilde_star) * yt_tilde_star[t, :];
+for t in 1:(n - 1)
+    yt_tilde_star[t + 1, :] = (A_tilde - B_tilde * F_tilde_star) *
+                              yt_tilde_star[t, :]
 end
 
 plot([yt_tilde_star[:, 5], yt_tilde[3, :]], labels = [L"\tilde{q}" L"q"])
 ```
 
 ```{code-cell} julia
-maximum(abs.(yt_tilde_star[:, 5] - yt_tilde[3, 1:end-1]))
+maximum(abs.(yt_tilde_star[:, 5] - yt_tilde[3, 1:(end - 1)]))
 ```
 
 ```{code-cell} julia
 ---
 tags: [remove-cell]
 ---
-#test maximum(abs.(yt_tilde_star[:, 5] - yt_tilde[3, 1:end-1])) < 1e-15
+@test maximum(abs.(yt_tilde_star[:, 5] - yt_tilde[3, 1:end-1])) < 1e-15
 ```
 
 ## Markov Perfect Equilibrium
@@ -1312,20 +1326,20 @@ B2 = [0, 1, 0];
 R1 = [0 0 -a0/2; 0 0 a1/2; -a0/2 a1/2 a1];
 R2 = [0 -a0/2 0; -a0/2 a1 a1/2; 0 a1/2 0];
 Q1 = Q2 = gamma;
-S1 = S2 = W1 = W2 = M1 = M2 = 0.;
+S1 = S2 = W1 = W2 = M1 = M2 = 0.0;
 
 # solve using nnash from QE
 F1, F2, P1, P2 = nnash(A, B1, B2, R1, R2, Q1, Q2,
-                        S1, S2, W1, W2, M1, M2,
-                        beta = beta,
-                        tol = tol1);
+                       S1, S2, W1, W2, M1, M2,
+                       beta = beta,
+                       tol = tol1);
 
 # simulate forward
 AF = A - B1 * F1 - B2 * F2;
 z = zeros(3, n);
 z[:, 1] .= 1;
-for t in 1:n-1
-    z[:, t+1] = AF * z[:, t]
+for t in 1:(n - 1)
+    z[:, t + 1] = AF * z[:, t]
 end
 
 println("Policy for F1 is $F1")
@@ -1336,8 +1350,8 @@ println("Policy for F2 is $F2")
 ---
 tags: [remove-cell]
 ---
-#test round(F1[1], digits = 4) == -0.227
-#test round(F2[2], digits = 4) == 0.0945
+@test round(F1[1], digits = 4) == -0.227
+@test round(F2[2], digits = 4) == 0.0945
 ```
 
 ```{code-cell} julia
@@ -1345,7 +1359,8 @@ q1 = z[2, :];
 q2 = z[3, :];
 q = q1 + q2; # total output, MPE
 p = a0 .- a1 * q; # total price, MPE
-plot([q, p], labels = ["total ouput" "total price"], title = "Output and prices, duopoly MPE", xlabel = L"t")
+plot([q, p], labels = ["total ouput" "total price"],
+     title = "Output and prices, duopoly MPE", xlabel = L"t")
 ```
 
 ```{code-cell} julia
@@ -1364,8 +1379,8 @@ tags: [remove-cell]
 # compute values
 u1 = -F1 * z;
 u2 = -F2 * z;
-pi_1 = (p .* q1)' - gamma * u1.^2;
-pi_2 = (p .* q2)' - gamma * u2.^2;
+pi_1 = (p .* q1)' - gamma * u1 .^ 2;
+pi_2 = (p .* q2)' - gamma * u2 .^ 2;
 
 v1_forward = pi_1 * betas;
 v2_forward = pi_2 * betas;
@@ -1381,10 +1396,10 @@ println("Firm 2: Direct is $v2_direct, Forward is $(v2_forward[1])");
 ---
 tags: [remove-cell]
 ---
-#test round(v1_direct, digits = 4) == 133.3296
-#test round(v2_direct, digits = 4) == 133.3296
-#test round(v1_forward[1], digits = 4) == 133.3303
-#test round(v2_forward[1], digits = 4) == 133.3303
+@test round(v1_direct, digits = 3) == 133.329
+@test round(v2_direct, digits = 3) == 133.329
+@test round(v1_forward[1], digits = 3) == 133.330
+@test round(v2_forward[1], digits = 3) == 133.330
 ```
 
 ```{code-cell} julia
@@ -1404,14 +1419,15 @@ vt_MPE = zeros(n);
 vt_follower = zeros(n);
 
 for t in 1:n
-    vt_MPE[t] = -z[:, t]' * P1 * z[:, t];
-    vt_follower[t] = -yt_tilde[:, t]' * P_tilde * yt_tilde[:, t];
+    vt_MPE[t] = -z[:, t]' * P1 * z[:, t]
+    vt_follower[t] = -yt_tilde[:, t]' * P_tilde * yt_tilde[:, t]
 end
 
-plot([vt_MPE, vt_leader, vt_follower], labels = ["MPE" "Stackelberg leader" "Stackelberg follower"], 
-        title = "MPE vs Stackelberg Values",
-        xlabel = L"t", 
-        legend = :outertopright)
+plot([vt_MPE, vt_leader, vt_follower],
+     labels = ["MPE" "Stackelberg leader" "Stackelberg follower"],
+     title = "MPE vs Stackelberg Values",
+     xlabel = L"t",
+     legend = :outertopright)
 ```
 
 ```{code-cell} julia
@@ -1423,7 +1439,7 @@ println("vt_MPE(y0) = $(vt_MPE[1])");
 
 ```{code-cell} julia
 # total difference in value b/t Stackelberg and MPE
-vt_leader[1] + vt_follower[1] - 2*vt_MPE[1]
+vt_leader[1] + vt_follower[1] - 2 * vt_MPE[1]
 ```
 
 ```{code-cell} julia
diff --git a/lectures/dynamic_programming_squared/lqramsey.md b/lectures/dynamic_programming_squared/lqramsey.md
index 3fefa1b7..84a22866 100644
--- a/lectures/dynamic_programming_squared/lqramsey.md
+++ b/lectures/dynamic_programming_squared/lqramsey.md
@@ -582,8 +582,7 @@ using Test
 ```
 
 ```{code-cell} julia
-using LaTeXStrings, QuantEcon, Plots, LinearAlgebra, Parameters
-
+using LaTeXStrings, QuantEcon, Plots, LinearAlgebra
 
 abstract type AbstractStochProcess end
 
@@ -592,7 +591,6 @@ struct ContStochProcess{TF <: AbstractFloat} <: AbstractStochProcess
     C::Matrix{TF}
 end
 
-
 struct DiscreteStochProcess{TF <: AbstractFloat} <: AbstractStochProcess
     P::Matrix{TF}
     x_vals::Matrix{TF}
@@ -621,9 +619,8 @@ function compute_exog_sequences(econ, x)
     return g, d, b, s, Sm
 end
 
-
 function compute_allocation(econ, Sm, nu, x, b)
-    (;Sg, Sd, Sb, Ss) = econ
+    (; Sg, Sd, Sb, Ss) = econ
     # solve for the allocation given nu and x
     Sc = 0.5 .* (Sb + Sd - Sg - nu .* Sm)
     Sl = 0.5 .* (Sb - Sd + Sg - nu .* Sm)
@@ -636,12 +633,11 @@ function compute_allocation(econ, Sm, nu, x, b)
     return Sc, Sl, c, l, p, tau, rvn
 end
 
-
 function compute_nu(a0, b0)
     disc = a0^2 - 4a0 * b0
 
     if disc >= 0
-        nu = 0.5 *(a0 - sqrt(disc)) / a0
+        nu = 0.5 * (a0 - sqrt(disc)) / a0
     else
         println("There is no Ramsey equilibrium for these parameters.")
         error("Government spending (economy.g) too low")
@@ -656,21 +652,21 @@ function compute_nu(a0, b0)
     return nu
 end
 
-
-function compute_Pi(B, R, rvn, g,xi)
-    pi = B[2:end] - R[1:end-1] .* B[1:end-1] - rvn[1:end-1] + g[1:end-1]
-    Pi = cumsum(pi .*xi)
+function compute_Pi(B, R, rvn, g, xi)
+    pi = B[2:end] - R[1:(end - 1)] .* B[1:(end - 1)] - rvn[1:(end - 1)] +
+         g[1:(end - 1)]
+    Pi = cumsum(pi .* xi)
     return pi, Pi
 end
 
-
-function compute_paths(econ::Economy{<:AbstractFloat, <:DiscreteStochProcess}, T)
+function compute_paths(econ::Economy{<:AbstractFloat, <:DiscreteStochProcess},
+                       T)
     # simplify notation
-    (;beta, Sg, Sd, Sb, Ss) = econ
-    (;P, x_vals) = econ.proc
+    (; beta, Sg, Sd, Sb, Ss) = econ
+    (; P, x_vals) = econ.proc
 
     mc = MarkovChain(P)
-    state = simulate(mc, T, init=1)
+    state = simulate(mc, T, init = 1)
     x = x_vals[:, state]
 
     # Compute exogenous sequence
@@ -678,9 +674,9 @@ function compute_paths(econ::Economy{<:AbstractFloat, <:DiscreteStochProcess}, T
 
     # compute a0, b0
     ns = size(P, 1)
-    F = I - beta.*P
-    a0 = (F \ ((Sm * x_vals)'.^2))[1] ./ 2
-    H = ((Sb - Sd + Sg) * x_vals) .* ((Sg - Ss)*x_vals)
+    F = I - beta .* P
+    a0 = (F \ ((Sm * x_vals)' .^ 2))[1] ./ 2
+    H = ((Sb - Sd + Sg) * x_vals) .* ((Sg - Ss) * x_vals)
     b0 = (F \ H')[1] ./ 2
 
     # compute lagrange multiplier
@@ -690,24 +686,24 @@ function compute_paths(econ::Economy{<:AbstractFloat, <:DiscreteStochProcess}, T
     Sc, Sl, c, l, p, tau, rvn = compute_allocation(econ, Sm, nu, x, b)
 
     # compute remaining variables
-    H = ((Sb - Sc) * x_vals) .* ((Sl - Sg) * x_vals) - (Sl * x_vals).^2
+    H = ((Sb - Sc) * x_vals) .* ((Sl - Sg) * x_vals) - (Sl * x_vals) .^ 2
     temp = dropdims(F * H', dims = 2)
     B = temp[state] ./ p
     H = dropdims(P[state, :] * ((Sb - Sc) * x_vals)', dims = 2)
     R = p ./ (beta .* H)
-    temp = dropdims(P[state, :] *((Sb - Sc) * x_vals)', dims = 2)
-   xi = p[2:end] ./ temp[1:end-1]
+    temp = dropdims(P[state, :] * ((Sb - Sc) * x_vals)', dims = 2)
+    xi = p[2:end] ./ temp[1:(end - 1)]
 
     # compute pi
     pi, Pi = compute_Pi(B, R, rvn, g, xi)
 
-    return (;g, d, b, s, c, l, p, tau, rvn, B, R, pi, Pi, xi)
+    return (; g, d, b, s, c, l, p, tau, rvn, B, R, pi, Pi, xi)
 end
 
 function compute_paths(econ::Economy{<:AbstractFloat, <:ContStochProcess}, T)
     # simplify notation
-    (;beta, Sg, Sd, Sb, Ss) = econ
-    (;A, C) = econ.proc
+    (; beta, Sg, Sd, Sb, Ss) = econ
+    (; A, C) = econ.proc
 
     # generate an initial condition x0 satisfying x0 = A x0
     nx, nx = size(A)
@@ -722,7 +718,7 @@ function compute_paths(econ::Economy{<:AbstractFloat, <:ContStochProcess}, T)
     w = randn(nw, T)
     x[:, 1] = x0
     for t in 2:T
-        x[:, t] = A *x[:, t-1] + C * w[:, t]
+        x[:, t] = A * x[:, t - 1] + C * w[:, t]
     end
 
     # compute exogenous sequence
@@ -731,7 +727,7 @@ function compute_paths(econ::Economy{<:AbstractFloat, <:ContStochProcess}, T)
     # compute a0 and b0
     H = Sm'Sm
     a0 = 0.5 * var_quadratic_sum(A, C, H, beta, x0)
-    H = (Sb - Sd + Sg)'*(Sg + Ss)
+    H = (Sb - Sd + Sg)' * (Sg + Ss)
     b0 = 0.5 * var_quadratic_sum(A, C, H, beta, x0)
 
     # compute lagrange multiplier
@@ -741,51 +737,50 @@ function compute_paths(econ::Economy{<:AbstractFloat, <:ContStochProcess}, T)
     Sc, Sl, c, l, p, tau, rvn = compute_allocation(econ, Sm, nu, x, b)
 
     # compute remaining variables
-    H = Sl'Sl - (Sb - Sc)' *(Sl - Sg)
+    H = Sl'Sl - (Sb - Sc)' * (Sl - Sg)
     L = zeros(T)
     for t in eachindex(L)
         L[t] = var_quadratic_sum(A, C, H, beta, x[:, t])
     end
     B = L ./ p
-    Rinv = dropdims(beta .* (Sb- Sc)*A*x, dims = 1) ./ p
+    Rinv = dropdims(beta .* (Sb - Sc) * A * x, dims = 1) ./ p
     R = 1 ./ Rinv
     AF1 = (Sb - Sc) * x[:, 2:end]
-    AF2 = (Sb - Sc) * A * x[:, 1:end-1]
-   xi =  AF1 ./ AF2
-   xi = dropdims(xi, dims = 1)
+    AF2 = (Sb - Sc) * A * x[:, 1:(end - 1)]
+    xi = AF1 ./ AF2
+    xi = dropdims(xi, dims = 1)
 
     # compute pi
     pi, Pi = compute_Pi(B, R, rvn, g, xi)
 
-    return(;g, d, b, s, c, l, p, tau, rvn, B, R, pi, Pi, xi)
+    return (; g, d, b, s, c, l, p, tau, rvn, B, R, pi, Pi, xi)
 end
 
 function gen_fig_1(path)
     T = length(path.c)
 
-    plt_1 = plot(path.rvn, lw=2, label = L"\tau_t l_t")
-    plot!(plt_1, path.g, lw=2, label= L"g_t")
-    plot!(plt_1, path.c, lw=2, label= L"c_t")
-    plot!(xlabel="Time", grid=true)
+    plt_1 = plot(path.rvn, lw = 2, label = L"\tau_t l_t")
+    plot!(plt_1, path.g, lw = 2, label = L"g_t")
+    plot!(plt_1, path.c, lw = 2, label = L"c_t")
+    plot!(xlabel = "Time", grid = true)
 
-    plt_2 = plot(path.rvn, lw=2, label=L"\tau_t l_t")
-    plot!(plt_2, path.g, lw=2, label=L"g_t")
-    plot!(plt_2, path.B[2:end], lw=2, label=L"B_{t+1}")
-    plot!(xlabel="Time", grid=true)
+    plt_2 = plot(path.rvn, lw = 2, label = L"\tau_t l_t")
+    plot!(plt_2, path.g, lw = 2, label = L"g_t")
+    plot!(plt_2, path.B[2:end], lw = 2, label = L"B_{t+1}")
+    plot!(xlabel = "Time", grid = true)
 
-    plt_3 = plot(path.R, lw=2, label=L"R_{t-1}")
-    plot!(plt_3, xlabel="Time", grid=true)
+    plt_3 = plot(path.R, lw = 2, label = L"R_{t-1}")
+    plot!(plt_3, xlabel = "Time", grid = true)
 
-    plt_4 = plot(path.rvn, lw=2, label=L"\tau_t l_t")
-    plot!(plt_4, path.g, lw=2, label=L"g_t")
-    plot!(plt_4, path.pi, lw=2, label=L"\pi_t")
-    plot!(plt_4, xlabel="Time", grid=true)
+    plt_4 = plot(path.rvn, lw = 2, label = L"\tau_t l_t")
+    plot!(plt_4, path.g, lw = 2, label = L"g_t")
+    plot!(plt_4, path.pi, lw = 2, label = L"\pi_t")
+    plot!(plt_4, xlabel = "Time", grid = true)
 
-    plot(plt_1, plt_2, plt_3, plt_4, layout=(2,2), size = (800,600))
+    plot(plt_1, plt_2, plt_3, plt_4, layout = (2, 2), size = (800, 600))
 end
 
 function gen_fig_2(path)
-
     T = length(path.c)
 
     paths = [path.xi, path.Pi]
@@ -795,9 +790,10 @@ function gen_fig_2(path)
     plots = [plt_1, plt_2]
 
     for (plot, path, label) in zip(plots, paths, labels)
-        plot!(plot, 2:T, path, lw=2, label=label, xlabel="Time", grid=true)
+        plot!(plot, 2:T, path, lw = 2, label = label, xlabel = "Time",
+              grid = true)
     end
-    plot(plt_1, plt_2, layout=(2,1), size = (600,500))
+    plot(plt_1, plt_2, layout = (2, 1), size = (600, 500))
 end
 ```
 
@@ -842,9 +838,9 @@ Random.seed!(42)
 
 # parameters
 beta = 1 / 1.05
-rho, mg = .7, .35
+rho, mg = 0.7, 0.35
 A = [rho mg*(1 - rho); 0.0 1.0]
-C = [sqrt(1 - rho^2) * mg / 10 0.0; 0 0]
+C = [sqrt(1 - rho^2) * mg/10 0.0; 0 0]
 Sg = [1.0 0.0]
 Sd = [0.0 0.0]
 Sb = [0 2.135]
@@ -858,22 +854,6 @@ path = compute_paths(econ, T)
 gen_fig_1(path)
 ```
 
-```{code-cell} julia
----
-tags: [remove-cell]
----
-@testset begin
-  #test path.p[3] ≈ 1.5395294981420302 atol = 1e-3 # randomness check.
-  #test path.g[31] ≈ 0.366487070014081 atol = 1e-3 # stuff we plot
-  #test path.c[36] ≈ 0.6291101011610297 atol = 1e-3
-  #test path.B[9] ≈ 0.07442403655989423 atol = 1e-3
-  #test path.rvn[27] ≈ 0.35013269833342753 atol = 1e-3
-  #test path.pi[31] ≈ -0.05846215388377568 atol = 1e-3
-  #test path.R[43] ≈ 1.0437715852385672 atol = 1e-3
-  #test path.xi[43] ≈  1.001895202392805 atol = 1e-3
-  #test path.Pi[43] ≈ -0.4185282208457552 atol = 1e-3 # plot tests
-end
-```
 
 The legends on the figures indicate the variables being tracked.
 
@@ -907,16 +887,16 @@ Random.seed!(42);
 # Parameters
 beta = 1 / 1.05
 P = [0.8 0.2 0.0
-    0.0 0.5 0.5
-    0.0 0.0 1.0]
+     0.0 0.5 0.5
+     0.0 0.0 1.0]
 
 # Possible states of the world
 # Each column is a state of the world. The rows are [g d b s 1]
 x_vals = [0.5 0.5 0.25;
-        0.0 0.0  0.0;
-        2.2 2.2  2.2;
-        0.0 0.0  0.0;
-        1.0 1.0  1.0]
+          0.0 0.0 0.0;
+          2.2 2.2 2.2;
+          0.0 0.0 0.0;
+          1.0 1.0 1.0]
 Sg = [1.0 0.0 0.0 0.0 0.0]
 Sd = [0.0 1.0 0.0 0.0 0.0]
 Sb = [0.0 0.0 1.0 0.0 0.0]
@@ -930,23 +910,6 @@ path = compute_paths(econ, T)
 gen_fig_1(path)
 ```
 
-```{code-cell} julia
----
-tags: [remove-cell]
----
-@testset begin
-  #test path.p[3] ≈ 1.5852129146694405
-  #test path.B[13] ≈ 0.003279632025474284
-  #@test path.g ≈ [0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25,
-  #                 0.25, 0.25]
-  #test path.rvn[7] ≈ 0.3188722725349599
-  #test path.c[2] ≈ 0.6147870853305598
-  #@test path.R ≈ [1.05, 1.05, 1.05, 1.05, 1.05, 1.0930974212983846, 1.05, 1.05, 1.05, 1.05,
-  #                1.05, 1.05, 1.05, 1.05, 1.05]
-  #@test path.xi ≈ [1.0, 1.0, 1.0, 1.0, 1.0, 0.9589548368586813, 1.0, 1.0, 1.0, 1.0, 1.0,
-  #                1.0, 1.0, 1.0]
-end
-```
 
 The call `gen_fig_2(path)` generates
 
@@ -1011,21 +974,6 @@ T = 50
 path = compute_paths(econ, T)
 ```
 
-```{code-cell} julia
----
-tags: [remove-cell]
----
-@testset begin
-  #test path.p[3]  ≈ 1.524261187305079
-  #test path.B[13]   ≈ -0.053219518947408805
-  #test path.g[7] ≈ 0.36908804521710115
-  #test path.rvn[7]  ≈ 0.35146870025913474
-  #test path.c[2]  ≈ 0.6259521929536346
-  #test path.R[5][1]  ≈ 1.0501742289013196
-  #test path.xi[10] ≈ 1.002202281639002
-end
-```
-
 ```{code-cell} julia
 gen_fig_1(path)
 ```
diff --git a/lectures/dynamic_programming_squared/opt_tax_recur.md b/lectures/dynamic_programming_squared/opt_tax_recur.md
index a8e0efb2..b58ab291 100644
--- a/lectures/dynamic_programming_squared/opt_tax_recur.md
+++ b/lectures/dynamic_programming_squared/opt_tax_recur.md
@@ -1490,30 +1490,29 @@ This utility function is implemented in the type CRRAutility
 
 ```{code-cell} julia
 function crra_utility(;
-    beta = 0.9,
-    sigma = 2.0,
-    gamma = 2.0,
-    Pi = 0.5 * ones(2, 2),
-    G = [0.1, 0.2],
-    Theta = ones(2),
-    transfers = false
-    )
+                      beta = 0.9,
+                      sigma = 2.0,
+                      gamma = 2.0,
+                      Pi = 0.5 * ones(2, 2),
+                      G = [0.1, 0.2],
+                      Theta = ones(2),
+                      transfers = false)
     function U(c, n)
         if sigma == 1.0
             U = log(c)
         else
-            U = (c.^(1.0 .- sigma) .- 1.0) / (1.0 - sigma)
+            U = (c .^ (1.0 .- sigma) .- 1.0) / (1.0 - sigma)
         end
-        return U .- n.^(1 + gamma) / (1 + gamma)
+        return U .- n .^ (1 + gamma) / (1 + gamma)
     end
     # Derivatives of utility function
-    Uc(c,n) =  c.^(-sigma)
-    Ucc(c,n) = -sigma * c.^(-sigma - 1.0)
-    Un(c,n) = -n.^gamma
-    Unn(c,n) = -gamma * n.^(gamma - 1.0)
+    Uc(c, n) = c .^ (-sigma)
+    Ucc(c, n) = -sigma * c .^ (-sigma - 1.0)
+    Un(c, n) = -n .^ gamma
+    Unn(c, n) = -gamma * n .^ (gamma - 1.0)
     n_less_than_one = false
     return Model(beta, Pi, G, Theta, transfers,
-                U, Uc, Ucc, Un, Unn, n_less_than_one)
+                 U, Uc, Ucc, Un, Unn, n_less_than_one)
 end
 ```
 
@@ -1528,15 +1527,15 @@ We can now plot the Ramsey tax  under both realizations of time $t = 3$ governme
 using Random
 Random.seed!(42) # For reproducible results.
 
-M_time_example = crra_utility(G=[0.1, 0.1, 0.1, 0.2, 0.1, 0.1],
-                              Theta=ones(6))            # Theta can in principle be random
+M_time_example = crra_utility(G = [0.1, 0.1, 0.1, 0.2, 0.1, 0.1],
+                              Theta = ones(6))            # Theta can in principle be random
 
 M_time_example.Pi = [0.0 1.0 0.0 0.0 0.0 0.0;
-                    0.0 0.0 1.0 0.0 0.0 0.0;
-                    0.0 0.0 0.0 0.5 0.5 0.0;
-                    0.0 0.0 0.0 0.0 0.0 1.0;
-                    0.0 0.0 0.0 0.0 0.0 1.0;
-                    0.0 0.0 0.0 0.0 0.0 1.0]
+                     0.0 0.0 1.0 0.0 0.0 0.0;
+                     0.0 0.0 0.0 0.5 0.5 0.0;
+                     0.0 0.0 0.0 0.0 0.0 1.0;
+                     0.0 0.0 0.0 0.0 0.0 1.0;
+                     0.0 0.0 0.0 0.0 0.0 1.0]
 
 PP_seq_time = SequentialAllocation(M_time_example)  # Solve sequential problem
 
@@ -1556,38 +1555,23 @@ titles = hcat("Consumption",
               "Output")
 
 sim_seq_l_plot = [sim_seq_l[1:4]..., M_time_example.G[sHist_l],
-                  M_time_example.Theta[sHist_l].*sim_seq_l[2]]
+    M_time_example.Theta[sHist_l] .* sim_seq_l[2]]
 sim_seq_h_plot = [sim_seq_h[1:4]..., M_time_example.G[sHist_h],
-                  M_time_example.Theta[sHist_h].*sim_seq_h[2]]
-
+    M_time_example.Theta[sHist_h] .* sim_seq_h[2]]
 
 #plots = plot(layout=(3,2), size=(800,600))
 plots = [plot(), plot(), plot(), plot(), plot(), plot()]
-for i = 1:6
-    plot!(plots[i], sim_seq_l_plot[i], color=:black, lw=2,
-          marker=:circle, markersize=2, label="")
-    plot!(plots[i], sim_seq_h_plot[i], color=:red, lw=2,
-          marker=:circle, markersize=2, label="")
-    plot!(plots[i], title=titles[i], grid=true)
+for i in 1:6
+    plot!(plots[i], sim_seq_l_plot[i], color = :black, lw = 2,
+          marker = :circle, markersize = 2, label = "")
+    plot!(plots[i], sim_seq_h_plot[i], color = :red, lw = 2,
+          marker = :circle, markersize = 2, label = "")
+    plot!(plots[i], title = titles[i], grid = true)
 end
-plot(plots[1], plots[2], plots[3], plots[4], plots[5], plots[6], layout=(3,2), size=(800,600))
+plot(plots[1], plots[2], plots[3], plots[4], plots[5], plots[6],
+     layout = (3, 2), size = (800, 600))
 ```
 
-```{code-cell} julia
----
-tags: [remove-cell]
----
-@testset begin
-  @test M_time_example.G[sHist_l] ≈ [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1]
-  #test M_time_example.Theta[sHist_l] .* sim_seq_l[2] ≈ [1.026385289423105, 0.9945696863679917,
-  #                                                  0.9945696863679917, 0.9945696863679917,
-  #                                                  0.9945696863679917, 0.9945696863679917,
-  #                                                  0.9945696863679917]
-  @test M_time_example.G[sHist_h] ≈ [0.1, 0.1, 0.1, 0.2, 0.1, 0.1, 0.1]
-  #test sim_seq_l[end] ≈ [1.0361020796451619, 1.111111111111111, 1.052459380877434,
-  #                       1.111111111111111, 1.111111111111111, 1.111111111111111]
-end
-```
 
 **Tax smoothing**
 
@@ -1628,11 +1612,11 @@ The following plot illustrates how the government lowers the interest rate at
 time 0 by raising consumption
 
 ```{code-cell} julia
-plot(sim_seq_l[end], color=:black, lw=2,
-        marker=:circle, markersize=2, label="")
-plot!(sim_seq_h[end], color=:red, lw=2,
-        marker=:circle, markersize=2, label="")
-plot!(title="Gross Interest Rate", grid=true)
+plot(sim_seq_l[end], color = :black, lw = 2,
+     marker = :circle, markersize = 2, label = "")
+plot!(sim_seq_h[end], color = :red, lw = 2,
+      marker = :circle, markersize = 2, label = "")
+plot!(title = "Gross Interest Rate", grid = true)
 ```
 
 ### Government Saving
@@ -1701,35 +1685,25 @@ Random.seed!(42); # For reproducible results.
 ```
 
 ```{code-cell} julia
-M2 = crra_utility(G=[0.15], Pi=ones(1, 1), Theta=[1.0])
+M2 = crra_utility(G = [0.15], Pi = ones(1, 1), Theta = [1.0])
 
 PP_seq_time0 = SequentialAllocation(M2) # solve sequential problem
 
 B_vec = range(-1.5, 1.0, length = 100)
 taxpolicy = Matrix(hcat([simulate(PP_seq_time0, B_, 1, 2)[4] for B_ in B_vec]...)')
-interest_rate = Matrix(hcat([simulate(PP_seq_time0, B_, 1, 3)[end] for B_ in B_vec]...)')
+interest_rate = Matrix(hcat([simulate(PP_seq_time0, B_, 1, 3)[end]
+                             for B_ in B_vec]...)')
 
 titles = ["Tax Rate" "Gross Interest Rate"]
 labels = [[L"Time , $t = 0$" L"Time , $t \geq 0$"], ""]
-plots = plot(layout=(2,1), size =(700,600))
+plots = plot(layout = (2, 1), size = (700, 600))
 for (i, series) in enumerate((taxpolicy, interest_rate))
-    plot!(plots[i], B_vec, series, linewidth=2, label=labels[i])
-    plot!(plots[i], title=titles[i], grid=true, legend=:topleft)
+    plot!(plots[i], B_vec, series, linewidth = 2, label = labels[i])
+    plot!(plots[i], title = titles[i], grid = true, legend = :topleft)
 end
 plot(plots)
 ```
 
-```{code-cell} julia
----
-tags: [remove-cell]
----
-@testset begin
-  #test B_vec[3] ≈ -1.4494949494949494
-  #test taxpolicy[2, 2] ≈ 0.0020700125847712414
-  #test interest_rate[3, 1] ≈ 1.113064964490116
-end
-```
-
 The figure indicates  that if the government enters with  positive debt, it sets
 a tax rate at $t=0$ that is less than all later tax rates.
 
@@ -1781,9 +1755,10 @@ B1_vec = hcat([simulate(PP_seq_time0, B_, 1, 2)[3][2] for B_ in B_vec]...)'
 # Compute the optimal policy if the government could reset
 tau1_reset = Matrix(hcat([simulate(PP_seq_time0, B1, 1, 1)[4] for B1 in B1_vec]...)')
 
-plot(B_vec, taxpolicy[:, 2], linewidth=2, label=L"\tau_1")
-plot!(B_vec, tau1_reset, linewidth=2, label=L"\tau_1^R")
-plot!(title="Tax Rate", xlabel="Initial Government Debt", legend=:topleft, grid=true)
+plot(B_vec, taxpolicy[:, 2], linewidth = 2, label = L"\tau_1")
+plot!(B_vec, tau1_reset, linewidth = 2, label = L"\tau_1^R")
+plot!(title = "Tax Rate", xlabel = "Initial Government Debt", legend = :topleft,
+      grid = true)
 ```
 
 The tax rates in the figure are equal  for only two values of initial government debt.
@@ -1819,21 +1794,21 @@ $$
 We will write a new constructor LogUtility to represent this utility function
 
 ```{code-cell} julia
-function log_utility(;beta = 0.9,
-                      psi = 0.69,
-                      Pi = 0.5 * ones(2, 2),
-                      G = [0.1, 0.2],
-                      Theta = ones(2),
-                      transfers = false)
+function log_utility(; beta = 0.9,
+                     psi = 0.69,
+                     Pi = 0.5 * ones(2, 2),
+                     G = [0.1, 0.2],
+                     Theta = ones(2),
+                     transfers = false)
     # Derivatives of utility function
-    U(c,n) = log(c) + psi * log(1 - n)
-    Uc(c,n) = 1 ./ c
-    Ucc(c,n) = -c.^(-2.0)
-    Un(c,n) = -psi ./ (1.0 .- n)
-    Unn(c,n) = -psi ./ (1.0 .- n).^2.0
+    U(c, n) = log(c) + psi * log(1 - n)
+    Uc(c, n) = 1 ./ c
+    Ucc(c, n) = -c .^ (-2.0)
+    Un(c, n) = -psi ./ (1.0 .- n)
+    Unn(c, n) = -psi ./ (1.0 .- n) .^ 2.0
     n_less_than_one = true
     return Model(beta, Pi, G, Theta, transfers,
-                U, Uc, Ucc, Un, Unn, n_less_than_one)
+                 U, Uc, Ucc, Un, Unn, n_less_than_one)
 end
 ```
 
@@ -1853,7 +1828,7 @@ Random.seed!(42); # For reproducible results.
 
 ```{code-cell} julia
 M1 = log_utility()
-mu_grid = range(-0.6,  0.0, length = 200)
+mu_grid = range(-0.6, 0.0, length = 200)
 PP_seq = SequentialAllocation(M1)         # Solve sequential problem
 PP_bel = RecursiveAllocation(M1, mu_grid)  # Solve recursive problem
 
@@ -1865,8 +1840,8 @@ sim_seq = simulate(PP_seq, 0.5, 1, T, sHist)
 sim_bel = simulate(PP_bel, 0.5, 1, T, sHist)
 
 # Plot policies
-sim_seq_plot = [sim_seq[1:4]..., M1.G[sHist], M1.Theta[sHist].*sim_seq[2]]
-sim_bel_plot = [sim_bel[1:4]..., M1.G[sHist], M1.Theta[sHist].*sim_bel[2]]
+sim_seq_plot = [sim_seq[1:4]..., M1.G[sHist], M1.Theta[sHist] .* sim_seq[2]]
+sim_bel_plot = [sim_bel[1:4]..., M1.G[sHist], M1.Theta[sHist] .* sim_bel[2]]
 
 titles = hcat("Consumption",
               "Labor Supply",
@@ -1874,36 +1849,26 @@ titles = hcat("Consumption",
               "Tax Rate",
               "Government Spending",
               "Output")
-labels = [["Sequential", "Recursive"], ["",""], ["",""], ["",""], ["",""], ["",""]]
-plots=plot(layout=(3,2), size=(850,780))
-
-for i = 1:6
-    plot!(plots[i], sim_seq_plot[i], color=:black, lw=2, marker=:circle,
-          markersize=2, label=labels[i][1])
-    plot!(plots[i], sim_bel_plot[i], color=:blue, lw=2, marker=:xcross,
-          markersize=2, label=labels[i][2])
-    plot!(plots[i], title=titles[i], grid=true, legend=:topright)
+labels = [
+    ["Sequential", "Recursive"],
+    ["", ""],
+    ["", ""],
+    ["", ""],
+    ["", ""],
+    ["", ""],
+]
+plots = plot(layout = (3, 2), size = (850, 780))
+
+for i in 1:6
+    plot!(plots[i], sim_seq_plot[i], color = :black, lw = 2, marker = :circle,
+          markersize = 2, label = labels[i][1])
+    plot!(plots[i], sim_bel_plot[i], color = :blue, lw = 2, marker = :xcross,
+          markersize = 2, label = labels[i][2])
+    plot!(plots[i], title = titles[i], grid = true, legend = :topright)
 end
 plot(plots)
 ```
 
-```{code-cell} julia
----
-tags: [remove-cell]
----
-@testset begin
-  #test sim_seq_plot[1][14] ≈ 0.38396935397869975
-  #test sim_seq_plot[2][14] ≈ 0.5839693539786998
-  #test sim_seq_plot[3][14] ≈ 0.3951985593686047
-  #test sim_seq_plot[4][14] ≈ 0.3631746680706347
-  #test sim_seq_plot[5][14] ≈ 0.2
-  #test sim_seq_plot[6][14] ≈ 0.5839693539786998
-  #test sim_bel_plot[3][5] ≈ 0.5230509296608254 atol = 1e-3
-  #test sim_bel_plot[5][7] ≈ 0.1
-  #test sim_bel_plot[2][3] ≈ 0.5402933557593538 atol = 1e-3
-end
-```
-
 As should be expected, the recursive and sequential solutions produce almost
 identical allocations.