diff --git a/HARK/distribution.py b/HARK/distribution.py
index 274765b2d..eb3f7ab38 100644
--- a/HARK/distribution.py
+++ b/HARK/distribution.py
@@ -963,11 +963,11 @@ def _approx_equiprobable(
             )
 
         if np.array_equal(self.Sigma, np.diag(np.diag(self.Sigma))):
-            ind_atoms = np.empty((self.M, N))
+            ind_atoms = np.empty((self.M, N + 2 * tail_N))
 
             for i in range(self.M):
                 if self.Sigma[i, i] == 0.0:
-                    x_atoms = np.repeat(np.exp(self.mu[i]), N)
+                    x_atoms = np.repeat(np.exp(self.mu[i]), N + 2 * tail_N)
                     ind_atoms[i] = x_atoms
                 else:
                     x_atoms = (
@@ -983,7 +983,22 @@ def _approx_equiprobable(
             atoms = np.stack(
                 [ar.flatten() for ar in list(np.meshgrid(*atoms_list))], axis=1
             ).T
-            pmv = np.repeat(1 / (N**self.M), N**self.M)
+
+            interiors = np.empty([self.M, (N + 2 * tail_N) ** (self.M)])
+
+            inners = np.zeros(N + 2 * tail_N)
+
+            if tail_N > 0:
+                inners[:tail_N] = [(tail_N - i) for i in range(tail_N)]
+                inners[-tail_N:] = [(i + 1) for i in range(tail_N)]
+
+            for i in range(self.M):
+                inners_i = [inners for _ in range((N + 2 * tail_N) ** i)]
+
+                interiors[i] = np.repeat(
+                    [*inners_i], (N + 2 * tail_N) ** (self.M - (i + 1))
+                )
+
         else:
             if tail_bound is not None:
                 if type(tail_bound) is float:
@@ -1024,10 +1039,10 @@ def eval(params, z):
                 excl = []
 
                 for j in range(len(z)):
-                    if z[j, 0] != z[j, 1]:
-                        inds.append(j)
-                    else:
+                    if z[j, 0] == z[j, 1]:
                         excl.append(j)
+                    elif params[j] != 0.0:
+                        inds.append(j)
 
                 dim = len(inds)
 
@@ -1111,21 +1126,21 @@ def eval(params, z):
 
                     atoms[i] = xi_atoms
 
-            max_locs = np.argmax(np.abs(interiors), axis=0)
+        max_locs = np.argmax(np.abs(interiors), axis=0)
 
-            max_inds = np.stack([max_locs, np.arange(len(max_locs))], axis=1)
+        max_inds = np.stack([max_locs, np.arange(len(max_locs))], axis=1)
 
-            prob_locs = interiors[max_inds[:, 0], max_inds[:, 1]]
+        prob_locs = interiors[max_inds[:, 0], max_inds[:, 1]]
 
-            def prob_assign(x):
-                if x == 0:
-                    return 1 / (N**self.M)
-                else:
-                    return 0.0
+        def prob_assign(x):
+            if x == 0:
+                return 1 / (N**self.M)
+            else:
+                return 0.0
 
-            prob_vec = np.vectorize(prob_assign)
+        prob_vec = np.vectorize(prob_assign)
 
-            pmv = prob_vec(prob_locs)
+        pmv = prob_vec(prob_locs)
 
         limit = {
             "dist": self,
diff --git a/binder/postBuild b/binder/postBuild
deleted file mode 100755
index 4e4b22864..000000000
--- a/binder/postBuild
+++ /dev/null
@@ -1,10 +0,0 @@
-#!/bin/sh
-jupyter contrib nbextension install --user
-jupyter nbextension enable codefolding/main
-python -m cite2c.install
-result=$(python <<EOF
-from notebook.services.config.manager import ConfigManager
-cm = ConfigManager()
-cm.update('cite2c', {'zotero':{"user_id": "5043554","username": "econ-ark","access_token": "XZpH9NsoAZmDMmjLKiy8xMXX"}})
-EOF
-)
diff --git a/binder/postBuild.bat b/binder/postBuild.bat
deleted file mode 100644
index b95f9ed64..000000000
--- a/binder/postBuild.bat
+++ /dev/null
@@ -1,3 +0,0 @@
-jupyter contrib nbextension install --user
-jupyter nbextension enable codefolding/main
-python -m cite2c.install
diff --git a/binder/requirements.txt b/binder/requirements.txt
index 78add7d76..ed0b97cfd 100644
--- a/binder/requirements.txt
+++ b/binder/requirements.txt
@@ -1,8 +1,2 @@
-cite2c
 git+https://github.com/econ-ark/hark@master
-ipywidgets
-jupyter_contrib_nbextensions
-matplotlib
-numpy
-scipy
-seaborn
+sequence_jacobian # required for KS-HARK-presentation example
diff --git a/examples/ConsNewKeynesianModel/KS-HARK-presentation.ipynb b/examples/ConsNewKeynesianModel/KS-HARK-presentation.ipynb
index c9de4c1a8..8da2bb30d 100644
--- a/examples/ConsNewKeynesianModel/KS-HARK-presentation.ipynb
+++ b/examples/ConsNewKeynesianModel/KS-HARK-presentation.ipynb
@@ -10,31 +10,6 @@
     "By William Du"
    ]
   },
-  {
-   "cell_type": "code",
-   "execution_count": 1,
-   "id": "46e98a1e",
-   "metadata": {
-    "execution": {
-     "iopub.execute_input": "2024-07-11T15:18:12.929098Z",
-     "iopub.status.busy": "2024-07-11T15:18:12.928854Z",
-     "iopub.status.idle": "2024-07-11T15:18:14.164584Z",
-     "shell.execute_reply": "2024-07-11T15:18:14.163965Z"
-    }
-   },
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "Note: you may need to restart the kernel to use updated packages.\n"
-     ]
-    }
-   ],
-   "source": [
-    "%pip install sequence_jacobian -q"
-   ]
-  },
   {
    "cell_type": "code",
    "execution_count": 2,
@@ -60,6 +35,20 @@
     "from HARK.ConsumptionSaving.ConsNewKeynesianModel import NewKeynesianConsumerType"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "28d24bc5",
+   "metadata": {},
+   "source": [
+    "# Calibration and setup\n",
+    "\n",
+    "This notebook uses the HARK toolkit to solve [Krusell and Smith (1998)](https://www.journals.uchicago.edu/doi/abs/10.1086/250034), applying the [Sequence Space Jacobian](https://onlinelibrary.wiley.com/doi/abs/10.3982/ECTA17434) and [SSJ toolkit](https://github.com/shade-econ/sequence-jacobian). \n",
+    "\n",
+    "## Firm setup\n",
+    "\n",
+    "Collect calibration for the production economy in a dictionary."
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": 3,
@@ -75,43 +64,23 @@
    "outputs": [],
    "source": [
     "calibration = {\n",
-    "    \"eis\": 1,\n",
-    "    \"delta\": 0.025,\n",
-    "    \"alpha\": 0.11,\n",
-    "    \"L\": 1.0,\n",
-    "    \"K\": 1.0,\n",
-    "    \"Y\": 1.0,\n",
-    "    \"r\": 0.01,\n",
+    "    \"eis\": 1,  # Elasticity of intertemporal substitution\n",
+    "    \"delta\": 0.025,  # Depreciation rate\n",
+    "    \"alpha\": 0.11,  # Capital share of income\n",
+    "    \"L_ss\": 1.0,  # Steady state labor\n",
+    "    \"Y_ss\": 1.0,  # Steady state output\n",
+    "    \"r_ss\": 0.01,  # Steady state real interest rate\n",
     "}"
    ]
   },
-  {
-   "cell_type": "code",
-   "execution_count": 4,
-   "id": "2a97701f",
-   "metadata": {
-    "execution": {
-     "iopub.execute_input": "2024-07-11T15:18:16.176291Z",
-     "iopub.status.busy": "2024-07-11T15:18:16.176051Z",
-     "iopub.status.idle": "2024-07-11T15:18:16.178499Z",
-     "shell.execute_reply": "2024-07-11T15:18:16.178018Z"
-    }
-   },
-   "outputs": [],
-   "source": [
-    "L_ss = 1.0  # Steady state labor\n",
-    "r_ss = 0.01  # steady state interest rate\n",
-    "Y_ss = 1.0  # steady state output"
-   ]
-  },
   {
    "cell_type": "markdown",
    "id": "1bac99bf",
    "metadata": {},
    "source": [
-    "## \n",
+    "### Steady State Capital \n",
     "\n",
-    "Given these steady state choices, we will need to find $K, Z$ to clear the firm first order conditions.\n"
+    "Find steady state capital (`K_ss`) and productivity  (`Z`) implied by the firm's first order condition, aggregate resource constraint and the above steady state values of labor (`L_ss`), the real interest rate (`r_ss`) and output (`Y_ss`). "
    ]
   },
   {
@@ -132,15 +101,18 @@
     "\n",
     "\n",
     "def your_funcs(X):\n",
-    "    L = calibration[\"L\"]\n",
+    "    L_ss = calibration[\"L_ss\"]\n",
     "    alpha = calibration[\"alpha\"]\n",
     "    delta = calibration[\"delta\"]\n",
+    "    r_ss = calibration[\"r_ss\"]\n",
+    "    Y_ss = calibration[\"Y_ss\"]\n",
     "\n",
     "    K, Z = X\n",
-    "    # all RHS have to be 0\n",
+    "\n",
+    "    # Firms first order condition and aggregate resource constraint\n",
     "    f = [\n",
-    "        alpha * Z * (K / L) ** (alpha - 1) - delta - r_ss,  # r = MPK\n",
-    "        Z * K**alpha * L ** (1 - alpha) - Y_ss,  # Y = Z*F(K,L)\n",
+    "        alpha * Z * (K / L_ss) ** (alpha - 1) - delta - r_ss,  # r_ss = MPK\n",
+    "        Z * K**alpha * L_ss ** (1 - alpha) - Y_ss,  # Y = Z*F(K,L)\n",
     "    ]\n",
     "\n",
     "    return f\n",
@@ -148,7 +120,10 @@
     "\n",
     "sol = root(your_funcs, [1.0, 1.0])  # find roots\n",
     "\n",
-    "K_ss, Z_ss = sol.x"
+    "K_ss, Z_ss = sol.x\n",
+    "\n",
+    "calibration[\"K_ss\"] = K_ss\n",
+    "calibration[\"Z_ss\"] = Z_ss"
    ]
   },
   {
@@ -191,7 +166,7 @@
    "id": "922dce5d",
    "metadata": {},
    "source": [
-    " Let's double check the roots we find produce our chosen steady state values for $ r, Y , L$."
+    "Double check the roots we find produce our chosen steady state values. "
    ]
   },
   {
@@ -211,17 +186,19 @@
     "def firm(\n",
     "    K,\n",
     "    Z,\n",
-    "    L=calibration[\"L\"],\n",
+    "    L_ss=calibration[\"L_ss\"],\n",
     "    alpha=calibration[\"alpha\"],\n",
     "    delta=calibration[\"delta\"],\n",
     "):\n",
-    "    r = alpha * Z * (K / L) ** (alpha - 1) - delta\n",
-    "    w = (1 - alpha) * Z * (K / L) ** alpha\n",
-    "    Y = Z * K**alpha * L ** (1 - alpha)\n",
+    "    r = alpha * Z * (K / L_ss) ** (alpha - 1) - delta\n",
+    "    w = (1 - alpha) * Z * (K / L_ss) ** alpha\n",
+    "    Y = Z * K**alpha * L_ss ** (1 - alpha)\n",
     "    return r, w, Y\n",
     "\n",
     "\n",
-    "r_ss, w_ss, Y_ss = firm(sol.x[0], sol.x[1])"
+    "r_ss, w_ss, Y_ss = firm(sol.x[0], sol.x[1])\n",
+    "\n",
+    "calibration[\"w_ss\"] = w_ss"
    ]
   },
   {
@@ -246,7 +223,30 @@
     }
    ],
    "source": [
-    "print(r_ss, w_ss, Y_ss)"
+    "print(\"--------------------------------------+------------\")\n",
+    "print(f\"Steady state capital                  | {K_ss:.4f}\")\n",
+    "print(f\"Steady state output                   | {Y_ss:.4f}\")\n",
+    "print(f\"Steady state real interest rate       | {r_ss:.4f}\")\n",
+    "print(f\"Steady state wage                     | {w_ss:.4f}\")\n",
+    "print(\"--------------------------------------+------------\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "7ca88d16",
+   "metadata": {},
+   "source": [
+    "## HARK agent\n",
+    "\n",
+    "HARK represents agent problems where aggregate variables can affect an individual income process as instances of `NewKeynesianConsumerType`, a subclass of `AgentType` (see [A Gentle Introduction to HARK](https://docs.econ-ark.org/examples/Gentle-Intro/Gentle-Intro-To-HARK.html)).\n",
+    "\n",
+    "Instances of the `NewKeynesianConsumerType` class contain functions that solve policy functions of optimizing agents and compute Jacobian matrices in response to a policy shocks using the [SSJ toolkit](https://github.com/shade-econ/sequence-jacobian). \n",
+    "\n",
+    "**_NOTE:_**  `NewKeynesianConsumerType` is still a microeconomic agent solving an income fluctuation problem. The only difference from [`IndShockConsumerType`](https://docs.econ-ark.org/examples/HowWeSolveIndShockConsumerType/HowWeSolveIndShockConsumerType.html) is that additional aggregate labour income variables are passed to the agent's income process. This allows the agent's problem to be defined within a general equilibrium framework such as (but not restricted to) HANK. \n",
+    "\n",
+    "**_NOTE:_** Researchers can also create their own agent types by subclassing `AgentType`. See [HARK's documentation](https://docs.econ-ark.org/examples/HowWeSolveIndShockConsumerType/HowWeSolveIndShockConsumerType.html) for more information.\n",
+    "\n",
+    "To create an instance of `NewKeynesianConsumerType`, first specify parameters in a dictionary. To start, we use a discount factor of 0.98."
    ]
   },
   {
@@ -263,40 +263,42 @@
    },
    "outputs": [],
    "source": [
+    "L_ss = calibration[\"L_ss\"]\n",
+    "w_ss = calibration[\"w_ss\"]\n",
+    "r_ss = calibration[\"r_ss\"]\n",
+    "\n",
     "HANK_Dict = {\n",
-    "    # Parameters shared with the perfect foresight model\n",
+    "    # Individual agent 'preferences' (shared with perfect foresight model)\n",
     "    \"CRRA\": calibration[\"eis\"],  # Coefficient of relative risk aversion\n",
-    "    \"Rfree\": (1 + r_ss),  # Interest factor on assets\n",
     "    \"DiscFac\": 0.98,  # Intertemporal discount factor\n",
     "    \"LivPrb\": [0.99375],  # Survival probability\n",
+    "    # Individual lifcycle income process parameters\n",
     "    \"PermGroFac\": [1.00],  # Permanent income growth factor\n",
-    "    # Parameters that specify the income distribution over the lifecycle\n",
-    "    # Standard deviation of log permanent shocks to income\n",
-    "    \"PermShkStd\": [0.06],\n",
+    "    \"PermShkStd\": [0.06],  # Standard deviation of log permanent shocks to income\n",
     "    \"PermShkCount\": 5,  # Number of points in discrete approximation to permanent income shocks\n",
-    "    # Standard deviation of log transitory shocks to income\n",
-    "    \"TranShkStd\": [0.2],\n",
-    "    \"TranShkCount\": 5,\n",
-    "    # HANK params\n",
-    "    \"tax_rate\": [\n",
-    "        0,\n",
-    "    ],  # set to 0.0 because we are going to assume that labor here is actually after tax income\n",
-    "    \"labor\": [L_ss],\n",
-    "    \"wage\": [w_ss],\n",
-    "    # Number of points in discrete approximation to transitory income shocks\n",
+    "    \"TranShkStd\": [0.2],  # Standard deviation of log transitory shocks to income\n",
+    "    \"TranShkCount\": 5,  # Number of points in discrete approximation to transitory income shocks\n",
+    "    # Parameters related to unemployment and retirement\n",
     "    \"UnempPrb\": 0.0,  # Probability of unemployment while working\n",
     "    \"IncUnemp\": 0.0,  # Unemployment benefits replacement rate\n",
     "    \"UnempPrbRet\": 0.0000,  # Probability of \"unemployment\" while retired\n",
     "    \"IncUnempRet\": 0.0,  # \"Unemployment\" benefits when retired\n",
     "    \"T_retire\": 0.0,  # Period of retirement (0 --> no retirement)\n",
-    "    # Parameters for constructing the \"assets above minimum\" grid\n",
+    "    # Aggregates affecting the agent's decision\n",
+    "    \"Rfree\": (1 + r_ss),  # Interest factor for assets faced by agents\n",
+    "    \"wage\": [w_ss],  # Wage rate faced by agents\n",
+    "    \"tax_rate\": [\n",
+    "        0\n",
+    "    ],  # set to 0.0 because we are going to assume that labor here is actually after tax income\n",
+    "    \"labor\": [L_ss],  # Aggregate (mean) labor supply\n",
+    "    # Parameters for constructing \"assets above minimum\" grid\n",
     "    \"aXtraMin\": 0.0001,  # Minimum end-of-period \"assets above minimum\" value\n",
     "    \"aXtraMax\": 2000,  # Maximum end-of-period \"assets above minimum\" value\n",
     "    \"aXtraCount\": 200,  # Number of points in the base grid of \"assets above minimum\"\n",
     "    # Exponential nesting factor when constructing \"assets above minimum\" grid\n",
     "    \"aXtraNestFac\": 3,\n",
     "    \"aXtraExtra\": None,  # Additional values to add to aXtraGrid\n",
-    "    # Transition Matrix simulation parameters\n",
+    "    # Transition matrix simulation parameters\n",
     "    \"mCount\": 200,\n",
     "    \"mMax\": 2000,\n",
     "    \"mMin\": 0.0001,\n",
@@ -309,7 +311,7 @@
    "id": "3be9593e",
    "metadata": {},
    "source": [
-    "# Create HARK agent"
+    "Let's create a temporary instance of a `NewKeynesianConsumerType` agent."
    ]
   },
   {
@@ -326,7 +328,7 @@
    },
    "outputs": [],
    "source": [
-    "Agent = NewKeynesianConsumerType(**HANK_Dict)"
+    "tempAgent = NewKeynesianConsumerType(**HANK_Dict)"
    ]
   },
   {
@@ -334,10 +336,28 @@
    "id": "fa266888",
    "metadata": {},
    "source": [
-    "# Find Steady state discount factor clear asset market\n",
+    "Given a discount factor, the agent will supply a steady state capital stock `A_ss`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ebe1792e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "A_ss = tempAgent.compute_steady_state()[0]\n",
+    "print(f\"Steady state agent asset supply for beta = 0.98 is: {A_ss:.3f}\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "911805ab",
+   "metadata": {},
+   "source": [
+    "### Steady State Assets \n",
     "\n",
-    "We will estimate the discount factor to ensure that asset supply equals the steady state capital we found earlier. \n",
-    "\n"
+    "Since we are interested in computing a steady state equilibrium, we find `discFac` such that `A_ss`  clears the asset market given that `K_ss` is the steady state firm capital demand. "
    ]
   },
   {
@@ -365,7 +385,10 @@
     "def A_ss_func(beta):\n",
     "    HANK_Dict[\"DiscFac\"] = beta\n",
     "\n",
+    "    # For a given value of beta, we need to re-create the Agent\n",
     "    Agent_func = NewKeynesianConsumerType(**HANK_Dict, verbose=False)\n",
+    "\n",
+    "    # And then solve for the steady state supply\n",
     "    A_ss = Agent_func.compute_steady_state()[0]\n",
     "\n",
     "    return A_ss\n",
@@ -380,7 +403,15 @@
     "\n",
     "DiscFac = optimize.brentq(ss_dif, 0.8, 0.9999)\n",
     "\n",
-    "print(\"Time taken to solve for steady state\", time.time() - start)"
+    "print(f\"Time taken to solve for steady state {time.time() - start:.3f} secs.\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "01be6a8a",
+   "metadata": {},
+   "source": [
+    "We can now create the steady state agent using the general equilibrium discount factor."
    ]
   },
   {
@@ -399,8 +430,8 @@
    "source": [
     "# Create a new agent\n",
     "HANK_Dict[\"DiscFac\"] = DiscFac\n",
-    "Agent_GE = NewKeynesianConsumerType(**HANK_Dict, verbose=False)\n",
-    "A_ss, C_ss = Agent_GE.compute_steady_state()"
+    "steadyHANK = NewKeynesianConsumerType(**HANK_Dict, verbose=False)\n",
+    "A_ss, C_ss = steadyHANK.compute_steady_state()"
    ]
   },
   {
@@ -426,9 +457,12 @@
     }
    ],
    "source": [
-    "# to make sure goods and asset markets clear\n",
-    "print(\"goods_clearing\", Y_ss - C_ss - calibration[\"delta\"] * K_ss)\n",
-    "print(\"asset_clearing\", A_ss - K_ss)"
+    "# To make sure goods and asset markets clear\n",
+    "print(\n",
+    "    \"Final goods clearing:\",\n",
+    "    calibration[\"Y_ss\"] - C_ss - calibration[\"delta\"] * calibration[\"K_ss\"],\n",
+    ")\n",
+    "print(\"Asset clearing:\", A_ss - calibration[\"K_ss\"])"
    ]
   },
   {
@@ -436,7 +470,15 @@
    "id": "a8167383",
    "metadata": {},
    "source": [
-    "# Computing Heterogenous Agent Jacobians"
+    "# Computing Jacobians using SSJ\n",
+    "\n",
+    "With the steady state agent in hand, we can compute the Jacobians of the steady state HANK agent's policy with respect to the aggregate state variables.\n",
+    "\n",
+    "The `calc_jacobian` method of our `NewKeynesianConsumerType` agent computes the Jacobians of the steady state aggregate consumption (`CJACW`) and assets (`AJACW`) with respect to a pertubation of a specified variable.\n",
+    "\n",
+    "Recall `CJACW[s,t]` is the time $t$ response to a shock at time $s$. \n",
+    "\n",
+    "Here is an example where we shock the wage rate. "
    ]
   },
   {
@@ -463,9 +505,9 @@
    "source": [
     "start = time.time()\n",
     "\n",
-    "CJACW, AJACW = Agent_GE.calc_jacobian(\"wage\", 300)  # Wage jacobians\n",
+    "CJACW, AJACW = steadyHANK.calc_jacobian(\"wage\", 300)  # Wage jacobians\n",
     "\n",
-    "print(\"Time taken to compute jacobians\", time.time() - start)"
+    "print(f\"Time taken to compute wage Jacobians: {time.time() - start:.3f} seconds\")"
    ]
   },
   {
@@ -493,15 +535,16 @@
     }
    ],
    "source": [
-    "plt.plot(CJACW.T[0])\n",
-    "plt.plot(CJACW.T[20])\n",
-    "plt.plot(CJACW.T[50])\n",
-    "plt.plot(CJACW.T[100])\n",
+    "plt.plot(CJACW.T[0], label=\"s=0\")\n",
+    "plt.plot(CJACW.T[20], label=\"s=20\")\n",
+    "plt.plot(CJACW.T[50], label=\"s=50\")\n",
+    "plt.plot(CJACW.T[100], label=\"s=100\")\n",
     "plt.xlim(-2, 300)\n",
     "plt.plot(np.arange(300), np.zeros(300), color=\"k\")\n",
-    "plt.title(\"Consumption Jacobian Wage\")\n",
-    "plt.xlabel(\"quarters\")\n",
-    "plt.ylabel(\"C response\")\n",
+    "plt.title(\"Consumption Jacobian (Wage) \")\n",
+    "plt.xlabel(\"Quarters\")\n",
+    "plt.ylabel(\"Consumption response\")\n",
+    "plt.legend()\n",
     "plt.show()"
    ]
   },
@@ -529,9 +572,19 @@
    "source": [
     "start = time.time()\n",
     "\n",
-    "CJACR, AJACR = Agent_GE.calc_jacobian(\"Rfree\", 300)  # Rfree jacobians\n",
+    "CJACR, AJACR = steadyHANK.calc_jacobian(\"Rfree\", 300)  # Rfree jacobians\n",
     "\n",
-    "print(\"Time taken to compute jacobians\", time.time() - start)"
+    "print(\n",
+    "    f\"Time taken to compute return factor Jacobians: {time.time() - start:.3f} seconds\"\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "311a3a60",
+   "metadata": {},
+   "source": [
+    "Here is an example where we shock the return factor. "
    ]
   },
   {
@@ -559,16 +612,25 @@
     }
    ],
    "source": [
-    "plt.plot(CJACR.T[0])\n",
-    "plt.plot(CJACR.T[30])\n",
-    "plt.plot(CJACR.T[50])\n",
+    "plt.plot(CJACR.T[0], label=\"s =0\")\n",
+    "plt.plot(CJACR.T[20], label=\"s=20\")\n",
+    "plt.plot(CJACR.T[50], label=\"s=50\")\n",
+    "plt.plot(CJACR.T[100], label=\"s=100\")\n",
     "plt.plot(np.arange(300), np.zeros(300), color=\"k\")\n",
-    "plt.title(\"Consumption Jacobian interest rate\")\n",
-    "plt.xlabel(\"quarters\")\n",
-    "plt.ylabel(\"C response\")\n",
+    "plt.title(\"Consumption Jacobian (Return Factor)\")\n",
+    "plt.xlabel(\"Quarters\")\n",
+    "plt.ylabel(\"Consumption response\")\n",
     "plt.show()"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "c5d0b70e",
+   "metadata": {},
+   "source": [
+    "Let's store the Jacobians in a dictionary for later use."
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": 18,
@@ -621,7 +683,7 @@
    "id": "d3669fab",
    "metadata": {},
    "source": [
-    "## Other Blocks of the Model"
+    "## Impulse Response Functions and Simulations "
    ]
   },
   {
@@ -675,7 +737,7 @@
    "id": "f3d569b5",
    "metadata": {},
    "source": [
-    "# Solving for Impulse Responses"
+    "### Solving for Impulse Responses"
    ]
   },
   {
@@ -692,8 +754,8 @@
    },
    "outputs": [],
    "source": [
-    "T = 300  # <-- the length of the IRF\n",
-    "rho_Z = 0.8  # persistence of IRF shock\n",
+    "T = 300  # Length of the IRF\n",
+    "rho_Z = 0.8  # Persistence of IRF shock\n",
     "dZ = 0.001 * Z_ss * rho_Z ** np.arange(T)\n",
     "shocks = {\"Z\": dZ}\n",
     "\n",
@@ -701,7 +763,6 @@
     "unknowns = [\"K\"]\n",
     "targets = [\"asset_mkt\"]\n",
     "\n",
-    "\n",
     "irfs_Z = ks.solve_impulse_linear(SteadyState_Dict, unknowns, targets, shocks)"
    ]
   },
@@ -723,7 +784,7 @@
     "    irfs_list,\n",
     "    variables,\n",
     "    labels=[\" \"],\n",
-    "    ylabel=r\"Percentage points (dev. from ss)\",\n",
+    "    ylabel=r\"Percentage point (dev. from ss)\",\n",
     "    T_plot=50,\n",
     "    figsize=(18, 6),\n",
     "):\n",
@@ -734,13 +795,17 @@
     "    for i in range(n_var):\n",
     "        # plot all irfs\n",
     "        for j, irf in enumerate(irfs_list):\n",
-    "            ax[i].plot(irf[variables[i]][:50], label=labels[j])\n",
+    "            ax[i].plot(irf[variables[i]][:50] * 1e2, label=labels[j])\n",
     "        ax[i].set_title(variables[i])\n",
-    "        ax[i].set_xlabel(r\"$t$\")\n",
+    "        ax[i].set_xlabel(r\"Quarters\")\n",
     "        if i == 0:\n",
     "            ax[i].set_ylabel(ylabel)\n",
-    "        ax[i].legend()\n",
-    "    plt.show()"
+    "        # ax[i].legend()\n",
+    "        # ax[i].x\n",
+    "\n",
+    "    plt.tight_layout()\n",
+    "    plt.show()\n",
+    "    # plt.savefig(\"irf.png\")"
    ]
   },
   {
@@ -777,7 +842,7 @@
    "id": "a0130fc3",
    "metadata": {},
    "source": [
-    "# Simulating the model"
+    "### Simulating the model"
    ]
   },
   {
@@ -829,7 +894,6 @@
     "rhos = {\"Z\": 0.8}\n",
     "impulses = {}\n",
     "\n",
-    "\n",
     "for i in inputs:\n",
     "    own_shock = {i: sigmas[i] * rhos[i] ** np.arange(T)}\n",
     "    impulses[i] = ks.solve_impulse_linear(\n",
@@ -844,14 +908,6 @@
     "data_simul = simulate(list(impulses.values()), outputs, T_sim)\n",
     "plot_timeseries(data_simul, (1, 3), figsize=(12, 4))"
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "27732edd",
-   "metadata": {},
-   "outputs": [],
-   "source": []
   }
  ],
  "metadata": {
diff --git a/examples/Journeys/Journey-Policymaker.ipynb b/examples/Journeys/Journey-Policymaker.ipynb
index 8a8fe51c9..d3ce854c2 100644
--- a/examples/Journeys/Journey-Policymaker.ipynb
+++ b/examples/Journeys/Journey-Policymaker.ipynb
@@ -325,7 +325,7 @@
    "metadata": {},
    "source": [
     "### 2.3 Simulation\n",
-    "After we solved the model backwards, we can simulate agents forward using Monte Carlo or [Transition Matrices](https://github.com/econ-ark/HARK/tree/master/examples/ConsIndShockModel/IndShockConsumerType_Transition_Matrix_Example.ipynb). These results can be used for in-model regressions or plotting distributions of assets or consumption."
+    "After we solved the model backwards, we can simulate agents forward using Monte Carlo or [Transition Matrices](https://github.com/econ-ark/HARK/tree/master/examples/ConsNewKeynesianModel/Transition_Matrix_Example.ipynb). These results can be used for in-model regressions or plotting distributions of assets or consumption."
    ]
   },
   {