diff --git a/docs/book/content/theory/derivations.md b/docs/book/content/theory/derivations.md
index b9ac1022a..444aeb529 100644
--- a/docs/book/content/theory/derivations.md
+++ b/docs/book/content/theory/derivations.md
@@ -5,15 +5,15 @@
 This appendix contains derivations from the theory in the body of this book.
 
 
-(SecAppDerivIndSpecCons)=
+(SecAppDerivHHcons)=
 ## Household first order condition for industry-specific consumption demand
 
-  The derivation for the household first order condition for industry-specific consumption demand {eq}`EqHHFOCcm` is the following:
+  The derivation for the household first order condition for industry-specific consumption demand {eq}`EqHHFOCci` is the following:
   ```{math}
-  :label: EqAppDerivHHIndSpecConsFOC
-    \tilde{p}_{m,t} = \tilde{p}_{j,s,t}\alpha_m(c_{j,m,s,t} - c_{min,m})^{\alpha_m-1}\prod_{u\neq m}^M\left(c_{j,u,s,t} - c_{min,u}\right)^{\alpha_u} \\
-    \tilde{p}_{m,t}(c_{j,m,s,t} - c_{min,m}) = \tilde{p}_{j,s,t}\alpha_m(c_{j,m,s,t} - c_{min,m})^{\alpha_m}\prod_{u\neq m}^M\left(c_{j,u,s,t} - c_{min,u}\right)^{\alpha_u} \\
-    \tilde{p}_{m,t}(c_{j,m,s,t} - c_{min,m}) = \tilde{p}_{j,s,t}\alpha_m\prod_{m=1}^M\left(c_{j,m,s,t} - c_{min,m}\right)^{\alpha_m} = \alpha_m \tilde{p}_{j,s,t}c_{j,s,t}
+  :label: EqAppDerivHHconsFOC
+    \tilde{p}_{i,t} = \tilde{p}_{j,s,t}\alpha_m(c_{i,j,s,t} - c_{min,i})^{\alpha_i-1}\prod_{u\neq i}^I\left(c_{u,j,s,t} - c_{min,u}\right)^{\alpha_u} \\
+    \tilde{p}_{i,t}(c_{i,j,s,t} - c_{min,i}) = \tilde{p}_{j,s,t}\alpha_i(c_{i,j,s,t} - c_{min,i})^{\alpha_i}\prod_{u\neq i}^I\left(c_{u,j,s,t} - c_{min,u}\right)^{\alpha_u} \\
+    \tilde{p}_{i,t}(c_{i,j,s,t} - c_{min,i}) = \tilde{p}_{j,s,t}\alpha_i\prod_{i=1}^I\left(c_{i,j,s,t} - c_{min,i}\right)^{\alpha_i} = \alpha_i \tilde{p}_{j,s,t}c_{j,s,t}
   ```
 
 
diff --git a/docs/book/content/theory/equilibrium.md b/docs/book/content/theory/equilibrium.md
index a366c6e02..57ad0f5e9 100644
--- a/docs/book/content/theory/equilibrium.md
+++ b/docs/book/content/theory/equilibrium.md
@@ -60,10 +60,11 @@ The computational algorithm for solving for the steady-state follows the steps b
 
 2. Choose an initial guess for the values of the steady-state interest rate (the after-tax marginal product of capital) $\bar{r}^i$, wage rate $\bar{w}^i$,  portfolio rate of return $\bar{r}_p^i$, output prices $\boldsymbol{\bar{p}}^i$ (note that $\bar{p}_M =1$ since it's the numeraire good), total bequests $\overline{BQ}^{\,i}$, total household transfers $\overline{TR}^{\,i}$, and income multiplier $factor^i$, where superscript $i$ is the index of the iteration number of the guess.
 
-    1. Given $\boldsymbol{\bar{p}}^i$ find the price of the composite good, $\bar{p}$ using equation {eq}`EqCompPnorm2`
-    2. Using {eq}`Eq_tr` with $\overline{TR}^{\,i}$, find transfers to each household, $\overline{tr}_{j,s}^i$
-    3. Using the bequest transfer process, {eq}`Eq_bq` and aggregate bequests, $\overline{BQ}^{\,i}$, find $bq_{j,s}^i$
-    4. Given values $\bar{p}$, $\bar{r}_{p}^i$, $\bar{w}^i$ $\overline{bq}_{j,s}^i$, $\overline{tr}_{j,s}^i$, and $factor^i$, solve for the steady-state household labor supply $\bar{n}_{j,s}$ and savings $\bar{b}_{j,s+1}$ decisions for all $j$ and $E+1\leq s\leq E+S$.
+    1. Given $\boldsymbol{\bar{p}}^i$ find the price of consumption goods using {eq}`EqHH_pi2`
+    2. From price of consumption goods, determine the price of the composite consmpution good, $\bar{p}$ using equation {eq}`EqCompPnorm2`
+    3. Using {eq}`Eq_tr` with $\overline{TR}^{\,i}$, find transfers to each household, $\overline{tr}_{j,s}^i$
+    4. Using the bequest transfer process, {eq}`Eq_bq` and aggregate bequests, $\overline{BQ}^{\,i}$, find $bq_{j,s}^i$
+    5. Given values $\bar{p}$, $\bar{r}_{p}^i$, $\bar{w}^i$ $\overline{bq}_{j,s}^i$, $\overline{tr}_{j,s}^i$, and $factor^i$, solve for the steady-state household labor supply $\bar{n}_{j,s}$ and savings $\bar{b}_{j,s+1}$ decisions for all $j$ and $E+1\leq s\leq E+S$.
 
         1. Each of the $j\in 1,2,...J$ sets of $2S$ steady-state Euler equations can be solved separately. `OG-Core` parallelizes this process using the maximum number of processors possible (up to $J$ processors). Solve each system of Euler equations using a multivariate root-finder to solve the $2S$ necessary conditions of the household given by the following steady-state versions of stationarized household Euler equations {eq}`EqStnrz_eul_n`, {eq}`EqStnrz_eul_b`, and {eq}`EqStnrz_eul_bS` simultaneously for each $j$.
 
@@ -89,9 +90,9 @@ The computational algorithm for solving for the steady-state follows the steps b
         :label: EqSS_HHeul_bS
           (\bar{c}_{j,E+S})^{-\sigma} = e^{-\sigma g_y}\chi^b_j(\bar{b}_{j,E+S+1})^{-\sigma} \quad\forall j
         ```
-    5. Determine from the quantity of the composite consumption good consumed by each household, $\bar{c}_{j,s}$, use equation {eq}`EqHH_cmDem` to determine consumption of each output good, $\bar{c}_{m,j,s}$
-    6. Using $\bar{c}_{m,j,s}$ in {eq}`EqCmt`, solve for aggregate consumption of each output good, $\bar{C}_{m}$
-    7. Given values for $\bar{n}_{j,s}$ and $\bar{b}_{j,s+1}$ for all $j$ and $s$, solve for steady-state labor supply, $\bar{L}$, savings, $\bar{B}$
+    6. Determine from the quantity of the composite consumption good consumed by each household, $\bar{c}_{j,s}$, use equation {eq}`EqHH_cmDem` to determine consumption of each output good, $\bar{c}_{m,j,s}$
+    7. Using $\bar{c}_{m,j,s}$ in {eq}`EqCmt`, solve for aggregate consumption of each output good, $\bar{C}_{m}$
+    8. Given values for $\bar{n}_{j,s}$ and $\bar{b}_{j,s+1}$ for all $j$ and $s$, solve for steady-state labor supply, $\bar{L}$, savings, $\bar{B}$
         1. Use $\bar{n}_{j,s}$ and the steady-state version of the stationarized labor market clearing equation {eq}`EqStnrzMarkClrLab` to get a value for $\bar{L}^{i}$.
 
            ```{math}
@@ -104,14 +105,14 @@ The computational algorithm for solving for the steady-state follows the steps b
            :label: EqSS_Bt
              \bar{B} \equiv \frac{1}{1 + \bar{g}_{n}}\sum_{s=E+2}^{E+S+1}\sum_{j=1}^{J}\Bigl(\bar{\omega}_{s-1}\lambda_j\bar{b}_{j,s} + i_s\bar{\omega}_{s}\lambda_j\bar{b}_{j,s}\Bigr)
            ```
-    8. Solve for the exogenous government interest rate $\bar{r}_{gov}^{i}$ using equation {eq}`EqUnbalGBC_rate_wedge`.
-    9. Use {eq}`EqStnrzTfer` to find $\bar{Y}^i$ from the guess of $\overline{TR}^i$
-    10. Use {eq}`EqStnrz_DY` to find $\bar{D}^i$ from $\bar{Y}^i$
-    11. Using $\bar{D}^i$, we can find foreign investor holdings of debt, $\bar{D}^{f,i}$ from {eq}`EqMarkClr_zetaD2` and then solve for domestic debt holdings through the debt market clearing condition: $\bar{D}^{d,i} = \bar{D}^i - \bar{D}^{f,i}$
-    12. Using $\bar{Y}^i$, find government infrastructure investment, $\bar{I}_{g}$ from {eq}`EqStnrz_Igt`
-    13. Using the law of motion of the stock of infrastructure, {eq}`EqStnrz_Kgmt`, and $\bar{I}_{g}$, solve for $\bar{K}_{g}^{i}$
-    14. Find output and factor demands for M-1 industries:
-        1. By {eq}`EqMarkClrGoods_Mm1`, $\bar{Y}_{m}=\bar{C}_{m}$
+    9. Solve for the exogenous government interest rate $\bar{r}_{gov}^{i}$ using equation {eq}`EqUnbalGBC_rate_wedge`.
+    10. Use {eq}`EqStnrzTfer` to find $\bar{Y}^i$ from the guess of $\overline{TR}^i$
+    11. Use {eq}`EqStnrz_DY` to find $\bar{D}^i$ from $\bar{Y}^i$
+    12. Using $\bar{D}^i$, we can find foreign investor holdings of debt, $\bar{D}^{f,i}$ from {eq}`EqMarkClr_zetaD2` and then solve for domestic debt holdings through the debt market clearing condition: $\bar{D}^{d,i} = \bar{D}^i - \bar{D}^{f,i}$
+    13. Using $\bar{Y}^i$, find government infrastructure investment, $\bar{I}_{g}$ from {eq}`EqStnrz_Igt`
+    14. Using the law of motion of the stock of infrastructure, {eq}`EqStnrz_Kgmt`, and $\bar{I}_{g}$, solve for $\bar{K}_{g}^{i}$
+    15. Find output and factor demands for M-1 industries:
+        1. By {eq}`EqMarkClrGoods_Mm1`, $\hat{Y}_{m,t}=\hat{C}_{m,t}$, where $\hat{C}_{m,t}$ is determined by {eq}`EqStnrzEqCmt`
         2. The capital-output ratio can be determined from the FOC for the firms' choice of capital: $\frac{\bar{K}_m}{\bar{Y}_m} = \gamma_m\left[\frac{\bar{r} +\bar{\delta}_M - \bar{\tau}^{corp}_m\bar{\delta}^{\tau}_m - \bar{\tau}^{inv}_m\bar{\delta}_M}{\left(1 - \bar{\tau}^{corp}_m\right)\bar{p}_m(\bar{Z}_m)^\frac{\varepsilon_m-1}{\varepsilon_m}}\right]^{-\varepsilon_m}$
         3. Capital demand can thus be found: $\bar{K}_{m} = \frac{\bar{K}_m}{\bar{Y}_m} * \bar{Y}_m$
         4. Labor demand can be found by inverting the production function:
@@ -323,8 +324,10 @@ The stationary non-steady state (transition path) solution algorithm has followi
 4. Using {eq}`Eq_tr` with $\boldsymbol{\hat{TR}}^{\,i}$, find transfers to each household, $\boldsymbol{\hat{tr}}_{j,s}^i$
 5. Using the bequest transfer process, {eq}`Eq_bq` and aggregate bequests, $\boldsymbol{\hat{BQ}}^{\,i}$, find $\boldsymbol{\hat{bq}}_{j,s}^i$
 6. Given time path guesses $\{\boldsymbol{r}_p^i, \boldsymbol{\hat{w}}^i, \boldsymbol{p}^i, \boldsymbol{\hat{bq}}^i, \boldsymbol{\hat{tr}}^i\}$, we can solve for each household's lifetime decisions $\{n_{j,s,t},\hat{b}_{j,s+1,t+1}\}_{s=E+1}^{E+S}$ for all $j$, $E+1\leq s \leq E+S$, and $1\leq t\leq T_2+S-1$.
-   1. The household problem can be solved with a multivariate root finder solving the $2S$ equations and unknowns at once for each $j$ and $1\leq t\leq T+S-1$. The root finder uses $2S$ household Euler equations  {eq}`EqStnrz_eul_n`, {eq}`EqStnrz_eul_b`, and {eq}`EqStnrz_eul_bS` to solve for each household's $2S$ lifetime decisions. The household decision rules for each type and birth cohort are solved separately.
-   2. After solving the first iteration of time path iteration, subsequent initial values for the $J$, $2S$ root finding problems are based on the solution in the prior iteration. This speeds up computation further and makes the initial guess for the highly nonlinear system of equations start closer to the solution value.
+   1. Given $\boldsymbol{p}^i$ find the price of consumption goods using {eq}`EqHH_pi2`
+   2. From price of consumption goods, determine the price of the composite consmpution good, $\bar{p}$ using equation {eq}`EqCompPnorm2`
+   3. The household problem can be solved with a multivariate root finder solving the $2S$ equations and unknowns at once for each $j$ and $1\leq t\leq T+S-1$. The root finder uses $2S$ household Euler equations  {eq}`EqStnrz_eul_n`, {eq}`EqStnrz_eul_b`, and {eq}`EqStnrz_eul_bS` to solve for each household's $2S$ lifetime decisions. The household decision rules for each type and birth cohort are solved separately.
+   4. After solving the first iteration of time path iteration, subsequent initial values for the $J$, $2S$ root finding problems are based on the solution in the prior iteration. This speeds up computation further and makes the initial guess for the highly nonlinear system of equations start closer to the solution value.
 7. Determine from the quantity of the composite consumption good consumed by each household, $\hat{c}_{j,s,t}$, use equation {eq}`EqHH_cmDem` to determine consumption of each output good, $\hat{c}_{m,j,s,t}$
 8. Using $\hat{c}_{m,j,s,t}$ in {eq}`EqCmt`, solve for aggregate consumption of each output good, $\hat{C}_{m,t}$
 9.  Given values for $n_{j,s,t}$ and $\hat{b}_{j,s+1,t+1}$ for all $j$, $s$, and $t$, solve for aggregate labor supply, $\hat{L}_t$, and savings, $B_t$ in each period
@@ -338,7 +341,7 @@ The stationary non-steady state (transition path) solution algorithm has followi
 15. Using $\hat{Y}_t^i$, find government infrastructure investment, $\hat{I}_{g,t}$ from {eq}`EqStnrz_Igt`
 16. Using the law of motion of the stock of infrastructure, {eq}`EqStnrz_Kgmt`, and $\hat{I}_{g,t}$, solve for $\hat{K}_{g,t}^{i}$
 17. Find output and factor demands for M-1 industries:
-    1. By {eq}`EqMarkClrGoods_Mm1`, $\hat{Y}_{m,t}=\hat{C}_{m,t}$
+    1. By {eq}`EqMarkClrGoods_Mm1`, $\hat{Y}_{m,t}=\hat{C}_{m,t}$, where $\hat{C}_{m,t}$ is determined by {eq}`EqStnrzEqCmt`
     2. The capital-output ratio can be determined from the FOC for the firms' choice of capital: $\frac{\hat{K}_{m,t}}{\hat{Y}_{m,t}} = \gamma_m\left[\frac{r_t + \delta_{M,t} - \tau^{corp}_{m,t}\delta^{\tau}_{m,t} - \tau^{inv}_{m,t}\delta_{M,t}}{(1-\tau^{corp}_{m,t})p_{m,t}({Z}_{m,t})^\frac{\varepsilon_m -1}{\varepsilon_m}}\right]^{-\varepsilon_m}$
     3. Capital demand can thus be found: $\hat{K}_{m,t} = \frac{\hat{K}_{m,t}}{\hat{Y}_{m,t}} * \hat{Y}_{m,t}$
     4. Labor demand can be found by inverting the production function:
diff --git a/docs/book/content/theory/government.md b/docs/book/content/theory/government.md
index d49a7ca0e..2f1936978 100644
--- a/docs/book/content/theory/government.md
+++ b/docs/book/content/theory/government.md
@@ -17,9 +17,10 @@ Government levies taxes on households and firms, funds public pensions, and make
 Income taxes are modeled through the total tax liability function $T_{s,t}$, which can be decomposed into the effective tax rate times total income {eq}`EqTaxCalcLiabETR2`. In this chapter, we detail the household tax component of government activity $T_{s,t}$ in `OG-Core`, along with our method of incorporating detailed microsimulation data into a dynamic general equilibrium model.
 
 ```{math}
-:label: EqHHBC
-  c_{j,s,t} + b_{j,s+1,t+1} &= (1 + r_{hh,t})b_{j,s,t} + w_t e_{j,s} n_{j,s,t} + \\
-  &\quad\quad\zeta_{j,s}\frac{BQ_t}{\lambda_j\omega_{s,t}} + \eta_{j,s,t}\frac{TR_{t}}{\lambda_j\omega_{s,t}} + ubi_{j,s,t} - T_{s,t}  \\
+:label: EqHHBC2
+  p_t c_{j,s,t} + &\sum_{i=1}^I (1 + \tau^{c}_{i,t})p_{i,t}c_{min,i} + b_{j,s+1,t+1} = \\
+  &(1 + r_{p,t})b_{j,s,t} + w_t e_{j,s} n_{j,s,t} + \\
+  &\quad\quad\zeta_{j,s}\frac{BQ_t}{\lambda_j\omega_{s,t}} + \eta_{j,s,t}\frac{TR_{t}}{\lambda_j\omega_{s,t}} + ubi_{j,s,t} - T_{j,s,t}  \\
   &\quad\forall j,t\quad\text{and}\quad s\geq E+1 \quad\text{where}\quad b_{j,E+1,t}=0\quad\forall j,t
 ```
 
@@ -280,7 +281,7 @@ When computing the role of compliance on the effective tax rate, we take a weigh
 
 #### Consumption taxes
 
-Linear consumption taxes, $\tau^c_{m,t}$ can vary over time and by output good.
+Linear consumption taxes, $\tau^c_{i,t}$ can vary over time and by consumption good.
 
 #### Wealth taxes
 
@@ -389,8 +390,8 @@ Businesses face a linear tax rate $\tau^{b}_{m,t}$, which can vary by industry a
   We see from the household's budget constraint that taxes $T_{j,s,t}$ and transfers $TR_{t}$ enter into the household's decision,
 
   ```{math}
-  :label: EqHHBC2
-    p_t c_{j,s,t} + &\sum_{m=1}^M p_{m,t}c_{min,m} + b_{j,s+1,t+1} = \\
+  :label: EqHHBC3
+    p_t c_{j,s,t} + &\sum_{i=1}^I (1 + \tau^{c}_{i,t})p_{i,t}c_{min,i} + b_{j,s+1,t+1} = \\
     &(1 + r_{p,t})b_{j,s,t} + w_t e_{j,s} n_{j,s,t} + \\
     &\quad\quad\zeta_{j,s}\frac{BQ_t}{\lambda_j\omega_{s,t}} + \eta_{j,s,t}\frac{TR_{t}}{\lambda_j\omega_{s,t}} + ubi_{j,s,t} - T_{j,s,t}  \\
     &\quad\forall j,t\quad\text{and}\quad s\geq E+1 \quad\text{where}\quad b_{j,E+1,t}=0\quad\forall j,t
diff --git a/docs/book/content/theory/households.md b/docs/book/content/theory/households.md
index da8102e01..deb8559c1 100644
--- a/docs/book/content/theory/households.md
+++ b/docs/book/content/theory/households.md
@@ -20,86 +20,96 @@ In this section, we describe what is arguably the most important economic agent
 
 <!-- This last sentence will not be true once we start carefully modeling government benefit programs, which focus on individual members of the household. -->
 
-(SecHH_IndSpecCons)=
-## Household Industry-specific Consumption
+(SecHH_Cons)=
+## Household Consumption
 
   We describe the derivation and dynamics of the population distribution in the {ref}`Chap_Demog` chapter in this documentation and in more detail in the calibration chapter on demographics in the country-specific repository documentation. A measure $\omega_{1,t}$ of households is born each period, become economically relevant at age $s=E+1$ if they survive to that age, and live for up to $E+S$ periods ($S$ economically active periods), with the population of age-$s$ individuals in period $t$ being $\omega_{s,t}$. Let the age of a household be indexed by $s = \{1,2,...E+S\}$.
 
   At birth, each household age $s=1$ is randomly assigned one of $J$ ability groups, indexed by $j$. Let $\lambda_j$ represent the fraction of individuals in each ability group, such that $\sum_j\lambda_j=1$. Note that this implies that the distribution across ability types in each age is given by $\boldsymbol{\lambda}=[\lambda_1,\lambda_2,...\lambda_J]$. Once an household is born and assigned to an ability type, it remains that ability type for its entire lifetime. This is deterministic ability heterogeneity as described in the calibration chapter on the lifetime earnings process in the country-specific repository documentation. Let $e_{j,s}>0$ be a matrix of ability-levels such that an individual of ability type $j$ will have lifetime abilities of $[e_{j,1},e_{j,2},...e_{j,E+S}]$.
 
-  Individuals in this economy choose how much to work each period $n_{j,s,t}$ and how much to consume among $M$ different industry-specific consumption goods $c_{m,j,s,t}$. We assume that households aggregate these industry-specific consumption goods in their preferences into a composite consumption good $c_{j,s,t}$ every period in every individual's preferences according to the following Stone-Geary version of a Cobb-Douglas consumption aggregator,
+  Individuals in this economy choose how much to work each period $n_{j,s,t}$ and how much to consume among $I$ different consumption goods $c_{i,j,s,t}$. We assume that households aggregate these industry-specific consumption goods in their preferences into a composite consumption good $c_{j,s,t}$ every period in every individual's preferences according to the following Stone-Geary version of a Cobb-Douglas consumption aggregator,
   ```{math}
   :label: EqHHCompCons
-    c_{j,s,t} \equiv \prod_{m=1}^M \left(c_{m,j,s,t} - c_{min,m}\right)^{\alpha_m} \quad\forall j,s,t \quad\text{with}\quad \sum_{m=1}^M\alpha_m=1
+    c_{j,s,t} \equiv \prod_{i=1}^I \left(c_{i,j,s,t} - c_{min,i}\right)^{\alpha_i} \quad\forall j,s,t \quad\text{with}\quad \sum_{i=1}^I\alpha_i=1
   ```
-  where $c_{min,m}$ is the minimum consumption of good $m$ allowed.[^StoneGeary]
+  where $c_{min,i}$ is the minimum consumption of good $i$ allowed.[^StoneGeary]
 
-  Assume that the non-normalized price of each individual consumption good is $\tilde{p}_{m,t}$. We can solve for the optimal good-$m$ consumption demands $c_{m,j,s,t}$ as a function of composite consumption $c_{j,s,t}$ by minimizing the total after-tax expenditure on consumption given that individual consumption adds up to composite consumption according to {eq}`EqHHCompCons`. Letting $\tau^{c}_{m,t}$ represent the consumption tax rate on goods of type $m$, the Lagrangian for this expenditure minimization problem is the following.
+  Assume that the non-normalized price of each individual consumption good is $\tilde{p}_{i,t}$. We can solve for the optimal good-$i$ consumption demands $c_{i,j,s,t}$ as a function of composite consumption $c_{j,s,t}$ by minimizing the total after-tax expenditure on consumption given that individual consumption adds up to composite consumption according to {eq}`EqHHCompCons`. Letting $\tau^{c}_{i,t}$ represent the consumption tax rate on goods of type $i$, the Lagrangian for this expenditure minimization problem is the following.
   ```{math}
   :label: EqHHCostMinLagr}
-    \mathcal{L} = \sum_{m=1}^M (1 + \tau^{c}_{m,t})\tilde{p}_{m,t}c_{m,j,s,t} + \lambda_{j,s,t}\Bigl[c_{j,s,t} - \prod_{m=1}^M \left(c_{m,j,s,t} - c_{min,m}\right)^{\alpha_m}\Bigr] \quad\forall j,s,t
+    \mathcal{L} = \sum_{i=1}^I (1 + \tau^{c}_{i,t})\tilde{p}_{i,t}c_{i,j,s,t} + \lambda_{j,s,t}\Bigl[c_{j,s,t} - \prod_{i=1}^I \left(c_{i,j,s,t} - c_{min,i}\right)^{\alpha_i}\Bigr] \quad\forall j,s,t
   ```
   Because the Lagrangian multiplier on the constraint $\lambda_{j,s,t}$ represents the shadow price of an extra unit of composite consumption, we can relabel it as the price of composite consumption $\tilde{p}_{j,s,t}$.
   ```{math}
   :label: EqHHCostMinLagr2
-    \mathcal{L} = \sum_{m=1}^M(1 + \tau^{c}_{m,t}) \tilde{p}_{m,t}c_{m,j,s,t} + \tilde{p}_{j,s,t}\Bigl[c_{j,s,t} - \prod_{m=1}^M \left(c_{m,j,s,t} - c_{min,m}\right)^{\alpha_m}\Bigr] \quad\forall j,s,t
+    \mathcal{L} = \sum_{i=1}^I(1 + \tau^{c}_{i,t}) \tilde{p}_{i,t}c_{i,j,s,t} + \tilde{p}_{j,s,t}\Bigl[c_{j,s,t} - \prod_{i=1}^I \left(c_{i,j,s,t} - c_{min,i}\right)^{\alpha_i}\Bigr] \quad\forall j,s,t
   ```
   Note that the price of composite consumption in period $t$ can be different for each ability-$j$ and age-$s$ individual at this point.
 
-  The $M+1$ first order conditions of this constrained minimization problem are the following $M$ first order conditions {eq}`EqHHFOCcm` plus the composite consumption aggregator {eq}`EqHHCompCons`.[^IndSpecConsDeriv]
+  The $I+1$ first order conditions of this constrained minimization problem are the following $i$ first order conditions {eq}`EqHHFOCci` plus the composite consumption aggregator {eq}`EqHHCompCons`.[^ConsDeriv]
   ```{math}
-  :label: EqHHFOCcm
-    (1 + \tau^{c}_{m,t})\tilde{p}_{m,t} = \alpha_m \tilde{p}_{j,s,t}\left(\frac{c_{j,s,t}}{c_{m,j,s,t} - c_{min,m}}\right) \quad\forall m,j,s,t
+  :label: EqHHFOCci
+    (1 + \tau^{c}_{i,t})\tilde{p}_{i,t} = \alpha_i \tilde{p}_{j,s,t}\left(\frac{c_{j,s,t}}{c_{i,j,s,t} - c_{min,i}}\right) \quad\forall i,j,s,t
   ```
-  Solving {eq}`EqHHFOCcm` for $c_{m,j,s,t}$ gives the optimal demand function for consumption of good $m$ by ability-$j$ and age-$s$ individual in period $t$.
+  Solving {eq}`EqHHFOCci` for $c_{i,j,s,t}$ gives the optimal demand function for consumption of good $i$ by ability-$j$ and age-$s$ individual in period $t$.
   ```{math}
-  :label: EqHH_cmDem
-    c_{m,j,s,t} = \alpha_m\left(\frac{(1 + \tau^{c}_{m,t})\tilde{p}_{m,t}}{\tilde{p}_{j,s,t}}\right)^{-1}c_{j,s,t} + c_{min,m} \quad\forall m,j,s,t
+  :label: EqHH_ciDem
+    c_{i,j,s,t} = \alpha_i\left(\frac{(1 + \tau^{c}_{i,t})\tilde{p}_{i,t}}{\tilde{p}_{j,s,t}}\right)^{-1}c_{j,s,t} + c_{min,i} \quad\forall i,j,s,t
   ```
-  This household demand function for good-$m$ shows that $c_{m,j,s,t}$ is a fraction of total composite consumption $c_{j,s,t}$, and that fraction is negatively correlated with the relative price of good-$m$ to the composite good price.
+  This household demand function for good-$i$ shows that $c_{i,j,s,t}$ is a fraction of total composite consumption $c_{j,s,t}$, and that fraction is negatively correlated with the relative price of good-$i$ to the composite good price.
 
-  Substituting the demand equations {eq}`EqHH_cmDem` back into the composite consumption definition {eq}`EqHHCompCons` gives us the expression for the non-normalized composite price $\tilde{p}_{j,s,t}$ as a function of each non-normalized industry-$m$ good price $\tilde{p}_{m,t}$.
+  Substituting the demand equations {eq}`EqHH_ciDem` back into the composite consumption definition {eq}`EqHHCompCons` gives us the expression for the non-normalized composite price $\tilde{p}_{j,s,t}$ as a function of each non-normalized industry-$i$ good price $\tilde{p}_{i,t}$.
   ```{math}
   :label: EqCompPnonnorm
-    \tilde{p}_{j,s,t} = \prod_{m=1}^M\left(\frac{(1 + \tau^{c}_{m,t})\tilde{p}_{m,t}}{\alpha_m}\right)^{\alpha_m} \quad\forall j,s,t
+    \tilde{p}_{j,s,t} = \prod_{i=1}^I\left(\frac{(1 + \tau^{c}_{i,t})\tilde{p}_{i,t}}{\alpha_i}\right)^{\alpha_i} \quad\forall j,s,t
   ```
   Because nothing on the right-hand-side of {eq}`EqCompPnonnorm` is a function of $j$ or $s$, then $\tilde{p}_{j,s,t}=\tilde{p}_t$ for all $j$ and $s$.
   ```{math}
   :label: EqCompPnonnorm2
-    \tilde{p}_{t} = \prod_{m=1}^M\left(\frac{(1 + \tau^{c}_{m,t})\tilde{p}_{m,t}}{\alpha_m}\right)^{\alpha_m} \quad\forall t
+    \tilde{p}_{t} = \prod_{i=1}^I\left(\frac{(1 + \tau^{c}_{i,t})\tilde{p}_{i,t}}{\alpha_i}\right)^{\alpha_i} \quad\forall t
   ```
 
-  Finally, we assume that the consumption good in industry $M$ is the numeraire.[^Numeraire] We can normalize the composite consumption price $\tilde{p}_t$ and the remaining $M-1$ prices $\tilde{p}_{m,t}$ for $m=1,2,...M-1$ in every period $t$ by dividing all the equations with prices by the industry-$M$ price $\tilde{p}_{M,t}$. Then we can rewrite the optimal consumption demand {eq}`EqHH_cmDem` and composite price index {eq}`EqCompPnonnorm2` equations as the following functions of normalized prices,
+  Consumption goods are determined from production goods through a fixed $I\times M$ coefficient matrix, $\Pi^I$. Each element, $\pi_{i,m}$, of $\Pi^I$ represents the fraction of good $i$ that is made up of output from industry $m$. It follows that the prices of consumption goods can be expressed as a function of the prices of production goods and the fixed coefficient matrix $\Pi^I$.
   ```{math}
-  :label: EqHH_cmDem2
-    c_{m,j,s,t} = \alpha_m\left(\frac{(1 + \tau^{c}_{m,t})p_{m,t}}{p_t}\right)^{-1}c_{j,s,t} + c_{min,m} \quad\forall m,j,s,t
+  :label: EqHH_pi
+    \tilde{p}_{i,t} = \sum_{m=1}^M \pi_{i,m}\tilde{p}_{m,t} \quad\forall i,t
+  ```
+
+ We assume that the production good in industry $M$ is the numeraire.[^Numeraire] We can normalize the composite consumption price $\tilde{p}_t$ and the remaining $M-1$ output goods prices $\tilde{p}_{m,t}$ for $m=1,2,...M-1$ in every period $t$ by dividing all the equations with prices by the industry-$m$ price $\tilde{p}_{m,t}$. Simlarly, we can divide through all the consumption good prices in every period, $\tilde{p}_{i,t}$ for $i=1,2,...I$ and rewrite the optimal consumption demand {eq}`EqHH_ciDem` and composite price index {eq}`EqCompPnonnorm2` equations as the following functions of normalized prices,
+   ```{math}
+    :label: EqHH_pi2
+    p_{i,t} = \sum_{m=1}^M \pi_{i,m}p_{m,t} \quad\forall i,t
+  ```
+  ```{math}
+  :label: EqHH_ciDem2
+    c_{i,j,s,t} = \alpha_i\left(\frac{(1 + \tau^{c}_{i,t})p_{i,t}}{p_t}\right)^{-1}c_{j,s,t} + c_{min,i} \quad\forall m,j,s,t
   ```
   ```{math}
   :label: EqCompPnorm2
-    p_t = \prod_{m=1}^M\left(\frac{(1 + \tau^{c}_{m,t})p_{m,t}}{\alpha_m}\right)^{\alpha_m} \quad\forall t
+    p_t = \prod_{i=1}^I\left(\frac{(1 + \tau^{c}_{i,t})p_{i,t}}{\alpha_i}\right)^{\alpha_i} \quad\forall t
   ```
   ```{math}
   :label: EqPmPcompNormDef
     \text{where}\quad &p_{m,t} \equiv \frac{\tilde{p}_{m,t}}{\tilde{p}_{M,t}} \quad\forall m, t \quad\Rightarrow\quad p_{M,t} = 1 \quad\forall t \\
-    &\text{and}\quad p_t \equiv\frac{\tilde{p}_t}{\tilde{p}_{M,t}} \quad\forall t
+    \quad\text{and}\quad  &p_{i,t} \equiv \frac{\tilde{p}_{i,t}}{\tilde{p}_{M,t}} \quad\forall i, t \quad\text{and}\quad p_t \equiv\frac{\tilde{p}_t}{\tilde{p}_{M,t}} \quad\forall t
   ```
-  where $p_{m,t}$ and $p_t$ defined in {eq}`EqPmPcompNormDef` are normalized industry prices and normalized composite goods price, respectively, with the $M$th industry good being the numeraire.
+  where $p_{i,t}$ and $p_t$ defined in {eq}`EqPmPcompNormDef` are normalized consumption goods prices and normalized composite goods price, respectively, with the $M$th industry good being the numeraire.
 
 
 (SecHHBC)=
 ## Budget Constraint
 
-  Because the household's industry-specific demand problem from Section {ref}`SecHH_IndSpecCons` is characterized by equations {eq}`EqHHCompCons`, {eq}`EqHH_cmDem2`, and {eq}`EqCompPnorm2` is determined by functions of composite consumption $c_{j,s,t}$ and normalized industry prices $p_t$ and $p_{m,t}$, we can write the individual's utility maximization in terms of composite consumption $c_{j,s,t}$. An ability-$j$ and age-$s$ individual faces the following per-period budget constraint.
+  Because the household's industry-specific demand problem from Section {ref}`SecHH_Cons` is characterized by equations {eq}`EqHHCompCons`, {eq}`EqHH_ciDem2`, and {eq}`EqCompPnorm2` is determined by functions of composite consumption $c_{j,s,t}$ and normalized industry prices $p_t$ and $p_{i,t}$, we can write the individual's utility maximization in terms of composite consumption $c_{j,s,t}$. An ability-$j$ and age-$s$ individual faces the following per-period budget constraint.
 
   ```{math}
   :label: EqHHBC
-    p_t c_{j,s,t} + &\sum_{m=1}^M (1 + \tau^{c}_{m,t})p_{m,t}c_{min,m} + b_{j,s+1,t+1} = \\
+    p_t c_{j,s,t} + &\sum_{i=1}^I (1 + \tau^{c}_{i,t})p_{i,t}c_{min,i} + b_{j,s+1,t+1} = \\
     &(1 + r_{p,t})b_{j,s,t} + w_t e_{j,s} n_{j,s,t} + \\
     &\quad\quad\zeta_{j,s}\frac{BQ_t}{\lambda_j\omega_{s,t}} + \eta_{j,s,t}\frac{TR_{t}}{\lambda_j\omega_{s,t}} + ubi_{j,s,t} - T_{j,s,t}  \\
     &\quad\forall j,t\quad\text{and}\quad s\geq E+1 \quad\text{where}\quad b_{j,E+1,t}=0\quad\forall j,t
   ```
 
-  where $c_{j,s,t}$ is consumption, $b_{j,s+1,t+1}$ is savings for the next period, $r_{p,t}$ is the normalized interest rate (return) on household savings invested in the financial intermediary, $b_{j,s,t}$ is current period wealth (savings from last period), $w_t$ is the normalized wage, and $n_{j,s,t}$ is labor supply. Equations {eq}`eq_rK` and {eq}`eq_portfolio_return` of Chapter {ref}`Chap_FinInt` show how the rate of return from the financial intermediary $r_{p,t}$ might differ from the marginal product of capital $r_t$ and from the interest rate the government pays $r_{gov,t}$. Note that we must add in the cost of minimum consumption $c_{min,m}$ for all $m$ because that amount is subtracted out of composite consumption in {eq}`EqHHCompCons`.
+  where $c_{j,s,t}$ is consumption, $b_{j,s+1,t+1}$ is savings for the next period, $r_{p,t}$ is the normalized interest rate (return) on household savings invested in the financial intermediary, $b_{j,s,t}$ is current period wealth (savings from last period), $w_t$ is the normalized wage, and $n_{j,s,t}$ is labor supply. Equations {eq}`eq_rK` and {eq}`eq_portfolio_return` of Chapter {ref}`Chap_FinInt` show how the rate of return from the financial intermediary $r_{p,t}$ might differ from the marginal product of capital $r_t$ and from the interest rate the government pays $r_{gov,t}$. Note that we must add in the cost of minimum consumption $c_{min,i}$ for all $i$ because that amount is subtracted out of composite consumption in {eq}`EqHHCompCons`.
 
   The third term on the right-hand-side of the budget constraint {eq}`EqHHBC` represents the portion of total bequests $BQ_t$ that go to the age-$s$, income-group-$j$ household. Let $\zeta_{j,s}$ be the fraction of total bequests $BQ_t$ that go to the age-$s$, income-group-$j$ household, such that $\sum_{s=E+1}^{E+S}\sum_{j=1}^J\zeta_{j,s}=1$. We must divide that amount by the population of $(j,s)$ households $\lambda_j\omega_{s,t}$. The calibration chapter on beqests in the country-specific repository documentation details how to calibrate the $\zeta_{j,s}$ values from consumer finance data.
 
@@ -183,7 +193,7 @@ In this section, we describe what is arguably the most important economic agent
 
   It is necessary to multiply the disutility of labor in {eq}`EqHHPerUtil` by $e^{g_y(1-\sigma)}$ because labor supply $n_{j,s,t}$ is stationary, but both consumption $c_{j,s,t}$ and savings $b_{j,s+1,t+1}$ are growing at the rate of technological progress (see Chapter {ref}`Chap_Stnrz`). The $e^{g_y(1-\sigma)}$ term keeps the relative utility values of consumption, labor supply, and savings in the same units.
 
-  The final term in the period utility function {eq}`EqHHPerUtil` is the "warm glow" bequest motive. It is a CRRA utility of savings, discounted by the mortality rate $\rho_s$.[^mort_rates_note]  Intuitively, it signifies the utility a household gets in the event that they don't live to the next period with probability $\rho_s$. It is a utility of savings beyond its usual benefit of allowing for more consumption in the next period. This utility of bequests also has constant $\chi^b_j$ which adjusts the utility of bequests relative to consumption and can vary by lifetime income group $j$. This is helpful for calibrating the model to match wealth distribution moments. See the calibration chapter on beqests in the country-specific repository documentation for a discussion of the calibration. Note that any bequest before age $E+S$ is unintentional as it was bequeathed due an event of death that was uncertain. Intentional bequests are all bequests given in the final period of life in which death is certain $b_{j,E+S+1,t}$.
+  The final term in the period utility function {eq}`EqHHPerUtil` is the "warm glow" bequest motive. It is a CRRA utility of savings, discounted by the mortality rate $\rho_s$.[^Iort_rates_note]  Intuitively, it signifies the utility a household gets in the event that they don't live to the next period with probability $\rho_s$. It is a utility of savings beyond its usual benefit of allowing for more consumption in the next period. This utility of bequests also has constant $\chi^b_j$ which adjusts the utility of bequests relative to consumption and can vary by lifetime income group $j$. This is helpful for calibrating the model to match wealth distribution moments. See the calibration chapter on beqests in the country-specific repository documentation for a discussion of the calibration. Note that any bequest before age $E+S$ is unintentional as it was bequeathed due an event of death that was uncertain. Intentional bequests are all bequests given in the final period of life in which death is certain $b_{j,E+S+1,t}$.
 
   The household lifetime optimization problem is to choose consumption $c_{j,s,t}$, labor supply $n_{j,s,t}$, and savings $b_{j,s+1,t+1}$ in every period of life to maximize expected discounted lifetime utility, subject to budget constraints and upper-bound and lower-bound constraints.
 
@@ -194,7 +204,7 @@ In this section, we describe what is arguably the most important economic agent
 
   ```{math}
   :label: EqHHBC2
-    \text{s.t.}\quad &p_t c_{j,s,t} + \sum_{m=1}^M (1 + \tau^{c}_{m,t})p_{m,t}c_{min,m} + b_{j,s+1,t+1} = \\
+    \text{s.t.}\quad &p_t c_{j,s,t} + \sum_{i=1}^I (1 + \tau^{c}_{i,t})p_{i,t}c_{min,i} + b_{j,s+1,t+1} = \\
     &\quad (1 + r_{p,t})b_{j,s,t} + w_t e_{j,s} n_{j,s,t} + \zeta_{j,s}\frac{BQ_t}{\lambda_j\omega_{s,t}} + \eta_{j,s,t}\frac{TR_{t}}{\lambda_j\omega_{s,t}} + ubi_{j,s,t} - T_{s,t} \\
     &\qquad\text{and}\quad c_{j,s,t}\geq 0,\: n_{j,s,t} \in[0,\tilde{l}],\:\text{and}\: b_{j,1,t}=0 \quad\forall j, t, \:\text{and}\: E+1\leq s\leq E+S \nonumber
   ```
@@ -238,7 +248,7 @@ In this section, we describe what is arguably the most important economic agent
 
   The tax functions $\tau^{etr}_{s,t}$, $\tau^{mtrx}_{s,t}$, and $\tau^{mtry}_{s,t}$ are estimated in each country calibration model based on the currency units of the corresponding income data. However, the consumption units of the `OG-Core` model or any of its country calibrations are not in the same units as income data. For this reason, we have to transform the model income units $x$ and $y$ by a $factor$ so that they are in the same units as the income data on which the tax functions were estimated.
 
-  The tax rate functions are each functions of capital income and labor income $\tau(x,y)$. In order to make the tax functions return accurate tax rates associated with the correct levels of income, we multiply the model income $x^m$ and $y^m$ by a $factor$ so that they are in the same units as the real-world income data $\tau(factor\times x^m, factor\times y^m)$. We define the $factor$ such that average steady-state household total income in the model times the $factor$ equals the U.S. data average total income.
+  The tax rate functions are each functions of capital income and labor income $\tau(x,y)$. In order to make the tax functions return accurate tax rates associated with the correct levels of income, we multiply the model income $x^I$ and $y^I$ by a $factor$ so that they are in the same units as the real-world income data $\tau(factor\times x^I, factor\times y^I)$. We define the $factor$ such that average steady-state household total income in the model times the $factor$ equals the U.S. data average total income.
 
   ```{math}
   :label: EqIncFactor
@@ -294,15 +304,15 @@ If `use_zeta=False`, then bequests from households of lifetime earnings type `j`
 
   [^StoneGeary]: This functional form was originally proposed as a utility function by in a short comment by {cite}`Geary:1950` that aggregates differentiated goods into a scalar utility value. It is differentiated from Cobb-Douglas utility by the subsistence consumption amount in each term of the product. This function was further developed and operationalized by {cite}`Stone:1954`.
 
-  [^IndSpecConsDeriv]: See section {ref}`SecAppDerivIndSpecCons` in the {ref}`Chap_Deriv` Chapter for the derivation of the household industry-specific consumption demand.
+  [^ConsDeriv]: See section {ref}`SecAppDerivHHcons` in the {ref}`Chap_Deriv` Chapter for the derivation of the household industry-specific consumption demand.
 
-  [^Numeraire]: We can normalize the model by any of the $M$ industry-specific consumption good prices $\tilde{p}_{m,t}$ or we could normalize the model by the composite good price $\tilde{p}_t$. We choose to normalize by the $M$th industry good price $\tilde{p}_{M,t}$ because that industry is the only one the output of which can be used as investment, government spending, or government debt. Furthermore, this nicely nests the case of one industry in which all the other industries share in consumption is set to zero $\alpha_m=0$ for $m=1,2,...M-1$.
+  [^Numeraire]: We can normalize the model by any of the $I$ consumption good prices $\tilde{p}_{i,t}$, any of the $M$ industry-specific princes $\tilde{p}_{m,t}$, or we could normalize the model by the composite good price $\tilde{p}_t$. We choose to normalize by the $M$th industry good price $\tilde{p}_{M,t}$ because that industry is the only one the output of which can be used as investment, government spending, or government debt. Furthermore, this nicely nests the case of one industry in which all the other industries share in consumption is set to zero $\alpha_m=0$ for $m=1,2,...M-1$.
 
   [^sav_util_note]: Savings enters the period utility function to provide a "warm glow" bequest motive.
 
   [^frisch_note]: {cite}`Peterman:2016` shows that in a U.S. macro-model that has only an intensive margin of labor supply and no extensive margin and represents a broad composition of individuals supplying labor---such as `OG-Core`---a Frisch elasticity of around 0.9 is probably appropriate. He tests the implied macro elasticity when the assumed micro elasticities are small on the intensive margin but only macro aggregates---which include both extensive and intensive margin agents---are observed.
 
-  [^mort_rates_note]: See Section the mortality rate section of the calibration chapter on demographics in the country-specific repository documentation for a detailed discussion of mortality rates for the specific country calibration interfacing with `OG-Core`.
+  [^Iort_rates_note]: See Section the mortality rate section of the calibration chapter on demographics in the country-specific repository documentation for a detailed discussion of mortality rates for the specific country calibration interfacing with `OG-Core`.
 
   [^constraint_note]: It is important to note that savings also has an implicit upper bound $b_{j,s,t}\leq k$ above which consumption would be negative in current period. However, this upper bound on savings in taken care of by the Inada condition on consumption.
 
diff --git a/docs/book/content/theory/intro.md b/docs/book/content/theory/intro.md
index 22c51cb0d..836276152 100644
--- a/docs/book/content/theory/intro.md
+++ b/docs/book/content/theory/intro.md
@@ -10,7 +10,7 @@ The main characteristic that differentiates the overlapping generations model fr
 * Households
     * overlapping generations of finitely lived households
     * households are forward looking and see to maximize their expected lifetime utility, which is a function of consumption, labor supply, and bequests
-    * households choose consumption of $M$ different consumption goods, composite consumption, savings, and labor supply every period.
+    * households choose consumption of $i\in\{1,2,...I\}$ different consumption goods, composite consumption, savings, and labor supply every period.
     * the only uncertainty households face is with respect to their mortality risk
     * realistic demographics: mortality rates, fertility rates, immigration rates, population growth, and population distribution dynamics
     * heterogeneous lifetime income groups within each age cohort, calibrated from U.S. tax data
diff --git a/docs/book/content/theory/market_clearing.md b/docs/book/content/theory/market_clearing.md
index e49fa7a22..5c5198412 100644
--- a/docs/book/content/theory/market_clearing.md
+++ b/docs/book/content/theory/market_clearing.md
@@ -93,15 +93,22 @@ We also characterize here the law of motion for total bequests $BQ_t$. Although
   (SecMarkClrMktClr_goods)=
   ### Goods market clearing
 
-  All $M$ industry goods markets must clear. We make a simplifying assumption that only the $M$th industry output can be used as investment, government spending, or government debt. This means that total consumption of good $m$ equals total output of good $m$ in the first $M-1$ industries.
+  All $M$ industry goods markets must clear. Total demand of production good $m$ for consumption can be written as a function of total household demand for each consumption good $i$
+
   ```{math}
-  :label: EqMarkClrGoods_Mm1
-    Y_{m,t} = C_{m,t} \quad\forall t \quad\text{and}\quad m=1,2,...M-1
+  :label: EqMarkConsDemand
+    C_{m,t} = \sum_{i=1}^{I} \pi_{i,m} C_{i,t} \quad\forall t \quad\text{and}\quad m=1,2,...M
   ```
   where
   ```{math}
   :label: EqCmt
-    C_{m,t} \equiv \sum_{s=E+1}^{E+S}\sum_{j=1}^{J}\omega_{s,t}\lambda_j c_{m,j,s,t} \quad\forall m,t
+    C_{i,t} \equiv \sum_{s=E+1}^{E+S}\sum_{j=1}^{J}\omega_{s,t}\lambda_j c_{i,j,s,t} \quad\forall i,t
+  ```
+
+  Because we make a simplifying assumption that only the $M$th industry output can be used as investment, government spending, or government debt, consumption demand equals total output of good $m$ in the first $M-1$ industries.
+  ```{math}
+  :label: EqMarkClrGoods_Mm1
+    Y_{m,t} = C_{m,t} \quad\forall t \quad\text{and}\quad m=1,2,...M-1
   ```
 
   The output of the $M$th industry can be used for private investment, infrastructure investment, government spending, and government debt.[^M_ind] As such, the market clearing condition in the $M$th industry will look more like the traditional $Y=C+I+G+NX$ expression.[^RCrates_note]
diff --git a/docs/book/content/theory/stationarization.md b/docs/book/content/theory/stationarization.md
index d2607bf1e..458d4d1b5 100644
--- a/docs/book/content/theory/stationarization.md
+++ b/docs/book/content/theory/stationarization.md
@@ -29,9 +29,9 @@ The previous chapters derive all the equations necessary to solve for the steady
   -
   - $\hat{BQ}_{j,t}\equiv\frac{BQ_{j,t}}{e^{g_y t}\tilde{N}_t}$
   - $r_{p,t}$
-* - $\hat{c}_{m,j,s,t}\equiv \frac{c_{m,j,s,t}}{e^{g_y t}}$
+* - $\hat{c}_{i,j,s,t}\equiv \frac{c_{i,j,s,t}}{e^{g_y t}}$
   -
-  - $\hat{C}_{m,t}\equiv\frac{C_{m,t}}{e^{g_y t}\tilde{N}_t}$
+  - $\hat{C}_{i,t}\equiv\frac{C_{i,t}}{e^{g_y t}\tilde{N}_t}$
   - $r_{gov,t}$
 * - $\hat{tr}_{j,s,t}\equiv \frac{tr_{j,s,t}}{e^{g_y t}}$
   -
@@ -40,18 +40,18 @@ The previous chapters derive all the equations necessary to solve for the steady
 * - $\hat{ubi}_{j,s,t}\equiv\frac{ubi_{j,s,t}}{e^{g_y t}}$
   -
   - $\hat{TR}_t\equiv\frac{TR_t}{e^{g_y t}\tilde{N}_t}$
-  - $p_t \equiv \frac{\tilde{p}_t}{\tilde{p}_{M,t}}$
+  - $p_{i,t} \equiv \frac{\tilde{p}_{i,t}}{\tilde{p}_{M,t}}$
 * - $\hat{T}_{j,s,t}\equiv \frac{T_{j,s,t}}{e^{g_y t}}$
   -
   - $\hat{UBI}_t\equiv\frac{UBI_t}{e^{g_y t}\tilde{N}_t}$
-  - $p_{m,t} \equiv \frac{\tilde{p}_{m,t}}{\tilde{p}_{M,t}}$
+  - $p_t \equiv \frac{\tilde{p}_t}{\tilde{p}_{M,t}}$
 * - $\hat{w}_t\equiv \frac{w_t}{e^{g_y t}}$
   -
   - $\hat{D}_t\equiv\frac{D_t}{e^{g_y t}\tilde{N}_t}$
-  -
+  - $p_{m,t} \equiv \frac{\tilde{p}_{m,t}}{\tilde{p}_{M,t}}$
 ```
 
-The usual definition of equilibrium would be allocations and prices such that households optimize {eq}`EqHH_cmDem2`, {eq}`EqHHeul_n`, {eq}`EqHHeul_b`, and {eq}`EqHHeul_bS`, firms optimize {eq}`EqFirmFOC_L` and {eq}`EqFirmFOC_K`, and markets clear {eq}`EqMarkClrLab`, {eq}`EqMarkClr_DtDdDf`, {eq}`EqMarkClr_KtKdKf`, {eq}`EqMarkClrGoods_Mm1`, {eq}`EqMarkClrGoods_M`, and {eq}`EqMarkClrBQ`. In this chapter, we show how to stationarize each of these characterizing equations so that we can use our fixed point methods described in Sections {ref}`SecEqlbSSsoln` and {ref}`SecEqlbNSSsoln` of Chapter {ref}`Chap_Eqm` to solve for the equilibria in the steady-state and transition path equilibrium definitions.
+The usual definition of equilibrium would be allocations and prices such that households optimize {eq}`EqHH_ciDem2`, {eq}`EqHHeul_n`, {eq}`EqHHeul_b`, and {eq}`EqHHeul_bS`, firms optimize {eq}`EqFirmFOC_L` and {eq}`EqFirmFOC_K`, and markets clear {eq}`EqMarkClrLab`, {eq}`EqMarkClr_DtDdDf`, {eq}`EqMarkClr_KtKdKf`, {eq}`EqMarkClrGoods_Mm1`, {eq}`EqMarkClrGoods_M`, and {eq}`EqMarkClrBQ`. In this chapter, we show how to stationarize each of these characterizing equations so that we can use our fixed point methods described in Sections {ref}`SecEqlbSSsoln` and {ref}`SecEqlbNSSsoln` of Chapter {ref}`Chap_Eqm` to solve for the equilibria in the steady-state and transition path equilibrium definitions.
 
 
 (SecStnrzHH)=
@@ -61,30 +61,30 @@ The usual definition of equilibrium would be allocations and prices such that ho
 
   ```{math}
   :label: EqStnrzCompCons
-    \hat{c}_{j,s,t} \equiv \prod_{m=1}^M \left(\hat{c}_{m,j,s,t} - \hat{c}_{min,m,t}\right)^{\alpha_m} \quad\forall j,s,t \quad\text{with}\quad \sum_{m=1}^M\alpha_m=1
+    \hat{c}_{j,s,t} \equiv \prod_{i=1}^I \left(\hat{c}_{i,j,s,t} - \hat{c}_{min,i,t}\right)^{\alpha_i} \quad\forall j,s,t \quad\text{with}\quad \sum_{i=1}^I\alpha_i=1
   ```
   ```{math}
   :label: EqStnrz_cmDem2
-    \hat{c}_{m,j,s,t} = \alpha_m\left(\frac{p_{m,t}}{p_t}\right)^{-1}\hat{c}_{j,s,t} + \hat{c}_{min,m,t} \quad\forall m,j,s,t
+    \hat{c}_{i,j,s,t} = \alpha_i\left(\frac{p_{i,t}}{p_t}\right)^{-1}\hat{c}_{j,s,t} + \hat{c}_{min,i,t} \quad\forall i,j,s,t
   ```
   ```{math}
   :label: EqStnrz_cmin
-    \hat{c}_{min,m,t} \equiv
+    \hat{c}_{min,i,t} \equiv
     \begin{cases}
-      \frac{c_{min,m}}{e^{g_y t}} \quad\text{for}\quad t < T \\
-      \frac{c_{min,m}}{e^{g_y T}} \quad\text{for}\quad t \geq T
-    \end{cases} \quad\forall m
+      \frac{c_{min,i}}{e^{g_y t}} \quad\text{for}\quad t < T \\
+      \frac{c_{min,i}}{e^{g_y T}} \quad\text{for}\quad t \geq T
+    \end{cases} \quad\forall i
   ```
 
-  where {eq}`EqStnrzCompCons` is the stationarized Stone-Geary industry-specific consumption aggregator for composite consumption and  {eq}`EqStnrzCompCons` is the stationarized household demand for industry-specific consumption. The composite price aggregation equation {eq}`EqCompPnorm2` is already stationary.
+  where {eq}`EqStnrzCompCons` is the stationarized Stone-Geary consumption aggregator for composite consumption and  {eq}`EqStnrzCompCons` is the stationarized household demand for the composite consumption good. The composite price aggregation equation {eq}`EqCompPnorm2` is already stationary.
 
-  Note that the only way to stationarize the consumption aggregator {eq}`EqStnrzCompCons` and consumption demand {eq}`EqStnrz_cmDem2` is to divide $c_{min,m}$ by the growth rate $e^{g_y t}$. However, $c_{min,m}$ is already stationary. It is constant for each $m$. Therefore, the version of $\hat{c}_{min,m,t}$ divided by $e^{g_y t}$ would be changing over time (nonstationary) for $g_y\neq 0$. For this reason, we define $\hat{c}_{min,m,t}$ in {eq}`EqStnrz_cmin` as being constant after the steady-state period $T$ at whatever value it reaches at that period. In most cases with $g_y>0$, that value will be close to zero. But we use $\bar{c}_{min,m} = c_{min,m}/e^{g_y T}$ from {eq}`EqStnrz_cmin` as the steady-state value of $c_{min,m}$.
+  Note that the only way to stationarize the consumption aggregator {eq}`EqStnrzCompCons` and consumption demand {eq}`EqStnrz_cmDem2` is to divide $c_{min,i}$ by the growth rate $e^{g_y t}$. However, $c_{min,i}$ is already stationary. It is constant for each $m$. Therefore, the version of $\hat{c}_{min,i,t}$ divided by $e^{g_y t}$ would be changing over time (nonstationary) for $g_y\neq 0$. For this reason, we define $\hat{c}_{min,i,t}$ in {eq}`EqStnrz_cmin` as being constant after the steady-state period $T$ at whatever value it reaches at that period. In most cases with $g_y>0$, that value will be close to zero. But we use $\bar{c}_{min,i} = c_{min,i}/e^{g_y T}$ from {eq}`EqStnrz_cmin` as the steady-state value of $c_{min,i}$.
 
   The stationary version of the household budget constraint {eq}`EqHHBC` is found by dividing both sides of the equation by $e^{g_y t}$. For the savings term $b_{j,s+1,t+1}$, we must multiply and divide by $e^{g_y(t+1)}$, which leaves an $e^{g_y} = \frac{e^{g_y(t+1)}}{e^{g_y t}}$ in front of the stationarized variable.
 
   ```{math}
   :label: EqStnrzHHBC
-    p_t\hat{c}_{j,s,t} + &\sum_{m=1}^M p_{m,t}\hat{c}_{min,m} + e^{g_y}\hat{b}_{j,s+1,t+1} = \\
+    p_t\hat{c}_{j,s,t} + &\sum_{i=1}^I (1 + \tau^{c}_{i,t})p_{i,t}\hat{c}_{min,i} + e^{g_y}\hat{b}_{j,s+1,t+1} = \\
     &(1 + r_{p,t})\hat{b}_{j,s,t} + \hat{w}_t e_{j,s} n_{j,s,t} + \\
     &\quad\quad\zeta_{j,s}\frac{\hat{BQ}_t}{\lambda_j\hat{\omega}_{s,t}} + \eta_{j,s,t}\frac{\hat{TR}_{t}}{\lambda_j\hat{\omega}_{s,t}} + \hat{ubi}_{j,s,t} - \hat{T}_{j,s,t}  \\
     &\quad\forall j,t\quad\text{and}\quad s\geq E+1 \quad\text{where}\quad \hat{b}_{j,E+1,t}=0\quad\forall j,t
@@ -368,8 +368,8 @@ The usual definition of equilibrium would be allocations and prices such that ho
   ```
   where
   ```{math}
-  :label: EqCmt
-    \hat{C}_{m,t} \equiv \sum_{s=E+1}^{E+S}\sum_{j=1}^{J}\hat{\omega}_{s,t}\lambda_j \hat{c}_{m,j,s,t} \quad\forall m,t
+  :label: EqStnrzEqCmt
+    \hat{C}_{m,t} \equiv \sum_{i=1}^{I}\sum_{s=E+1}^{E+S}\sum_{j=1}^{J}\hat{\omega}_{s,t}\lambda_j \pi_{i,m} \hat{c}_{i,j,s,t} \quad\forall m,t
   ```
   and
   ```{math}
diff --git a/ogcore/SS.py b/ogcore/SS.py
index 2530c4177..18d1b7fd1 100644
--- a/ogcore/SS.py
+++ b/ogcore/SS.py
@@ -176,7 +176,7 @@ def inner_loop(outer_loop_vars, p, client):
         r_p (scalar): return on household investment portfolio
         r (scalar): real interest rate
         w (scalar): real wage rate
-        p_m (array_like): good prices
+        p_m (array_like): production goods prices
         BQ (array_like): aggregate bequest amount(s)
         TR (scalar): lump sum transfer amount
         Y (scalar): real GDP
@@ -196,7 +196,7 @@ def inner_loop(outer_loop_vars, p, client):
             * new_r_p (scalar): real interest rate on household
                 portfolio
             * new_w (scalar): real wage rate
-            * new_p_m (array_like): good prices
+            * new_p_i (array_like): good prices
             * K_vec (array_like): capital demand for each industry
             * L_vec (array_like): labor demand for each industry
             * Y_vec (array_like): output from each industry
@@ -213,11 +213,12 @@ def inner_loop(outer_loop_vars, p, client):
     bssmat, nssmat, r_p, r, w, p_m, Y, BQ, TR, factor = outer_loop_vars
 
     p_m = np.array(p_m)  # TODO: why is this a list otherwise?
+    p_i = np.dot(p.io_matrix, p_m)
     BQ = np.array(BQ)
     # initialize array for euler errors
     euler_errors = np.zeros((2 * p.S, p.J))
 
-    p_tilde = aggr.get_ptilde(p_m, p.tau_c[-1, :], p.alpha_c)
+    p_tilde = aggr.get_ptilde(p_i, p.tau_c[-1, :], p.alpha_c)
     bq = household.get_bq(BQ, None, p, "SS")
     tr = household.get_tr(TR, None, p, "SS")
     ubi = p.ubi_nom_array[-1, :, :] / factor
@@ -304,7 +305,7 @@ def inner_loop(outer_loop_vars, p, client):
         p.e,
         p,
     )
-    c_m = household.get_cm(c_s, p_m, p_tilde, p.tau_c[-1, :], p.alpha_c)
+    c_i = household.get_ci(c_s, p_i, p_tilde, p.tau_c[-1, :], p.alpha_c)
     L = aggr.get_L(nssmat, p, "SS")
     B = aggr.get_B(bssmat, p, "SS", False)
 
@@ -324,14 +325,13 @@ def inner_loop(outer_loop_vars, p, client):
     # Find output, labor demand, capital demand for M-1 industries
     L_vec = np.zeros(p.M)
     K_vec = np.zeros(p.M)
-    Y_vec = np.zeros(p.M)
-    C_vec = np.zeros(p.M)
+    C_vec = np.zeros(p.I)
     K_demand_open_vec = np.zeros(p.M)
+    for i_ind in range(p.I):
+        C_vec[i_ind] = aggr.get_C(c_i[i_ind, :, :], p, "SS")
+    Y_vec = np.dot(p.io_matrix.T, C_vec)
     for m_ind in range(p.M - 1):
-        C_m = aggr.get_C(c_m[m_ind, :, :], p, "SS")
-        C_vec[m_ind] = C_m
         KYrat_m = firm.get_KY_ratio(r, p_m, p, "SS", m_ind)
-        Y_vec[m_ind] = C_m
         K_vec[m_ind] = KYrat_m * Y_vec[m_ind]
         L_vec[m_ind] = firm.solve_L(
             Y_vec[m_ind], K_vec[m_ind], K_g, p, "SS", m_ind
@@ -348,7 +348,6 @@ def inner_loop(outer_loop_vars, p, client):
         B, K_demand_open_vec.sum(), D_d, p.zeta_K[-1]
     )
     K_M = max(0.001, K - K_vec.sum())  # make sure K_M > 0
-    C_vec[-1] = np.squeeze(aggr.get_C(c_m[-1, :], p, "SS"))
     L_vec[-1] = L_M
     K_vec[-1] = K_M
     Y_vec[-1] = firm.get_Y(K_vec[-1], K_g, L_vec[-1], p, "SS", -1)
@@ -396,7 +395,8 @@ def inner_loop(outer_loop_vars, p, client):
     # Find updated goods prices
     new_p_m = firm.get_pm(new_w, Y_vec, L_vec, p, "SS")
     new_p_m = new_p_m / new_p_m[-1]  # normalize prices by industry M
-    new_p_tilde = aggr.get_ptilde(new_p_m, p.tau_c[-1, :], p.alpha_c)
+    new_p_i = np.dot(p.io_matrix, new_p_m)
+    new_p_tilde = aggr.get_ptilde(new_p_i, p.tau_c[-1, :], p.alpha_c)
 
     etr_params_3D = np.tile(
         np.reshape(p.etr_params[-1, :, :], (p.S, 1, p.etr_params.shape[2])),
@@ -449,7 +449,7 @@ def inner_loop(outer_loop_vars, p, client):
         b_s,
         nssmat,
         new_bq,
-        c_m,
+        c_i,
         Y_vec,
         L_vec,
         K_vec,
@@ -487,6 +487,7 @@ def inner_loop(outer_loop_vars, p, client):
 
     G_vec = np.zeros(p.M)
     G_vec[-1] = G
+    C_m_vec = np.dot(p.io_matrix.T, C_vec)
     I_d_vec = np.zeros(p.M)
     I_d = aggr.get_I(b_splus1, K_d, K_d, p, "SS")
     I_d_vec[-1] = I_d
@@ -499,7 +500,7 @@ def inner_loop(outer_loop_vars, p, client):
     net_capital_outflows_vec = np.zeros(p.M)
     net_capital_outflows_vec[-1] = net_capital_outflows
     rc_error = (
-        Y_vec - C_vec - G_vec - I_d_vec - I_g_vec - net_capital_outflows_vec
+        Y_vec - C_m_vec - G_vec - I_d_vec - I_g_vec - net_capital_outflows_vec
     )
 
     return (
@@ -659,7 +660,8 @@ def SS_solver(
     Y_vec_ss = new_Y_vec
     r_gov_ss = fiscal.get_r_gov(rss, p)
     p_m_ss = new_p_m
-    p_tilde_ss = aggr.get_ptilde(p_m_ss, p.tau_c[-1, :], p.alpha_c)
+    p_i_ss = np.dot(p.io_matrix, p_m_ss)
+    p_tilde_ss = aggr.get_ptilde(p_i_ss, p.tau_c[-1, :], p.alpha_c)
     TR_ss = new_TR
     Yss = new_Y
     I_g_ss = fiscal.get_I_g(Yss, p.alpha_I[-1])
@@ -681,7 +683,6 @@ def SS_solver(
         K_demand_open_ss[m] = firm.get_K(
             p.world_int_rate[-1], w_open, L_vec_ss[m], p, "SS", m
         )
-    print("K_demand_open_ss = ", K_demand_open_ss, L_vec_ss)
     Kss, K_d_ss, K_f_ss = aggr.get_K_splits(
         Bss, K_demand_open_ss.sum(), D_d_ss, p.zeta_K[-1]
     )
@@ -800,15 +801,15 @@ def SS_solver(
     )
     yss_before_tax_mat = household.get_y(r_p_ss, wss, bssmat_s, nssmat, p)
     Css = aggr.get_C(cssmat, p, "SS")
-    c_m_ss_mat = household.get_cm(
-        cssmat, p_m_ss, p_tilde_ss, p.tau_c[-1, :], p.alpha_c
+    c_i_ss_mat = household.get_ci(
+        cssmat, p_i_ss, p_tilde_ss, p.tau_c[-1, :], p.alpha_c
     )
-    C_vec_ss = np.zeros(p.M)
-    for m_ind in range(
-        p.M
-    ):  # TODO: update aggr.get_C to take full MxSxJ array
-        C_vec_ss[m_ind] = aggr.get_C(
-            c_m_ss_mat[m_ind, :, :],
+    C_vec_ss = np.zeros(p.I)
+    for i_ind in range(
+        p.I
+    ):  # TODO: update aggr.get_C to take full IxSxJ array
+        C_vec_ss[i_ind] = aggr.get_C(
+            c_i_ss_mat[i_ind, :, :],
             p,
             "SS",
         )
@@ -830,7 +831,7 @@ def SS_solver(
         bssmat_s,
         nssmat,
         bqssmat,
-        c_m_ss_mat,
+        c_i_ss_mat,
         Y_vec_ss,
         L_vec_ss,
         K_vec_ss,
@@ -866,6 +867,9 @@ def SS_solver(
     )
     # Fill in arrays, noting that M-1 industries only produce consumption goods
     G_vec_ss = np.zeros(p.M)
+    # Map consumption goods back to demands for production goods
+    print("IO: ", p.io_matrix.T.shape, ", C: ", C_vec_ss.shape)
+    C_m_vec_ss = np.dot(p.io_matrix.T, C_vec_ss)
     G_vec_ss[-1] = Gss
     I_d_vec_ss = np.zeros(p.M)
     I_d_vec_ss[-1] = I_d_ss
@@ -875,7 +879,7 @@ def SS_solver(
     net_capital_outflows_vec[-1] = net_capital_outflows
     RC = aggr.resource_constraint(
         Y_vec_ss,
-        C_vec_ss,
+        C_m_vec_ss,
         G_vec_ss,
         I_d_vec_ss,
         I_g_vec_ss,
diff --git a/ogcore/TPI.py b/ogcore/TPI.py
index b6173593d..edb7e7793 100644
--- a/ogcore/TPI.py
+++ b/ogcore/TPI.py
@@ -363,7 +363,11 @@ def inner_loop(guesses, outer_loop_vars, initial_values, ubi, j, ind, p):
     r_p, r, w, p_m, BQ, TR, theta = outer_loop_vars
 
     # compute composite good price
-    p_tilde = aggr.get_ptilde(p_m[:, :], p.tau_c[:, :], p.alpha_c, "TPI")
+    p_i = (
+        np.tile(p.io_matrix.reshape(1, p.I, p.M), (p.T + p.S, 1, 1))
+        * np.tile(p_m.reshape(p.T + p.S, 1, p.M), (1, p.I, 1))
+    ).sum(axis=2)
+    p_tilde = aggr.get_ptilde(p_i[:, :], p.tau_c[:, :], p.alpha_c, "TPI")
     # compute bq
     bq = household.get_bq(BQ, None, p, "TPI")
     # compute tr
@@ -600,7 +604,11 @@ def run_TPI(p, client=None):
     p_m = p_m / p_m[:, -1].reshape(
         p.T + p.S, 1
     )  # normalize prices by industry M
-    p_tilde = aggr.get_ptilde(p_m[:, :], p.tau_c[:, :], p.alpha_c, "TPI")
+    p_i = (
+        np.tile(p.io_matrix.reshape(1, p.I, p.M), (p.T + p.S, 1, 1))
+        * np.tile(p_m.reshape(p.T + p.S, 1, p.M), (1, p.I, 1))
+    ).sum(axis=2)
+    p_tilde = aggr.get_ptilde(p_i[:, :], p.tau_c[:, :], p.alpha_c, "TPI")
     if not any(p.zeta_K == 1):
         w[: p.T] = np.squeeze(
             firm.get_w(Y[: p.T], L[: p.T], p_m[: p.T, :], p, "TPI")
@@ -612,7 +620,11 @@ def run_TPI(p, client=None):
     p_m = p_m / p_m[:, -1].reshape(
         p.T + p.S, 1
     )  # normalize prices by industry M
-    p_tilde = aggr.get_ptilde(p_m[:, :], p.tau_c[:, :], p.alpha_c, "TPI")
+    p_i = (
+        np.tile(p.io_matrix.reshape(1, p.I, p.M), (p.T + p.S, 1, 1))
+        * np.tile(p_m.reshape(p.T + p.S, 1, p.M), (1, p.I, 1))
+    ).sum(axis=2)
+    p_tilde = aggr.get_ptilde(p_i[:, :], p.tau_c[:, :], p.alpha_c, "TPI")
     # path for interest rates
     r = np.zeros_like(Y)
     r[: p.T] = np.squeeze(
@@ -697,6 +709,12 @@ def run_TPI(p, client=None):
     while (TPIiter < p.maxiter) and (TPIdist >= p.mindist_TPI):
 
         outer_loop_vars = (r_p, r, w, p_m, BQ, TR, theta)
+        # compute composite good price
+        p_i = (
+            np.tile(p.io_matrix.reshape(1, p.I, p.M), (p.T + p.S, 1, 1))
+            * np.tile(p_m.reshape(p.T + p.S, 1, p.M), (1, p.I, 1))
+        ).sum(axis=2)
+        p_tilde = aggr.get_ptilde(p_i[:, :], p.tau_c[:, :], p.alpha_c, "TPI")
 
         euler_errors = np.zeros((p.T, 2 * p.S, p.J))
         lazy_values = []
@@ -778,9 +796,10 @@ def run_TPI(p, client=None):
             p,
         )
         C = aggr.get_C(c_mat, p, "TPI")
-        c_m = household.get_cm(
+
+        c_i = household.get_ci(
             c_mat[: p.T, :, :],
-            p_m[: p.T, :],
+            p_i[: p.T, :],
             p_tilde[: p.T],
             p.tau_c[: p.T, :],
             p.alpha_c,
@@ -801,16 +820,18 @@ def run_TPI(p, client=None):
         # Find output, labor demand, capital demand for M-1 industries
         L_vec = np.zeros((p.T, p.M))
         K_vec = np.zeros((p.T, p.M))
-        Y_vec = np.zeros((p.T, p.M))
-        C_vec = np.zeros((p.T, p.M))
+        C_vec = np.zeros((p.T, p.I))
         K_demand_open_vec = np.zeros((p.T, p.M))
+        for i_ind in range(p.I):
+            C_vec[:, i_ind] = aggr.get_C(c_i[: p.T, i_ind, :, :], p, "TPI")
+        Y_vec = (
+            np.tile(p.io_matrix.reshape(1, p.I, p.M), (p.T, 1, 1))
+            * np.tile(C_vec[: p.T, :].reshape(p.T, p.I, 1), (1, 1, p.M))
+        ).sum(axis=1)
         for m_ind in range(p.M - 1):
-            C_m = aggr.get_C(c_m[: p.T, m_ind, :, :], p, "TPI")
-            C_vec[:, m_ind] = C_m
             KYrat_m = firm.get_KY_ratio(
                 r[: p.T], p_m[: p.T, :], p, "TPI", m_ind
             )
-            Y_vec[:, m_ind] = C_m
             K_vec[:, m_ind] = KYrat_m * Y_vec[:, m_ind]
             L_vec[:, m_ind] = firm.solve_L(
                 Y_vec[:, m_ind], K_vec[:, m_ind], K_g, p, "TPI", m_ind
@@ -841,7 +862,6 @@ def run_TPI(p, client=None):
             np.ones(p.T) * 0.001, K[: p.T] - K_vec[: p.T, :].sum(-1)
         )  # make sure K_M > 0
 
-        C_vec[:, -1] = np.squeeze(aggr.get_C(c_m[: p.T, -1, :], p, "TPI"))
         L_vec[:, -1] = L_M
         K_vec[:, -1] = K_M
         Y_vec[:, -1] = firm.get_Y(
@@ -867,7 +887,7 @@ def run_TPI(p, client=None):
             bmat_s,
             n_mat[: p.T, :, :],
             bqmat[: p.T, :, :],
-            c_m[: p.T, :, :, :],
+            c_i[: p.T, :, :, :],
             Y_vec[: p.T, :],
             L_vec[: p.T, :],
             K_vec[: p.T, :],
@@ -973,7 +993,7 @@ def run_TPI(p, client=None):
             bmat_s,
             n_mat[: p.T, :, :],
             bqmat_new[: p.T, :, :],
-            c_m[: p.T, :, :, :],
+            c_i[: p.T, :, :, :],
             Y_vec[: p.T, :],
             L_vec[: p.T, :],
             K_vec[: p.T, :],
@@ -1151,6 +1171,11 @@ def run_TPI(p, client=None):
     # Fill in arrays, noting that M-1 industries only produce consumption goods
     G_vec = np.zeros((p.T, p.M))
     G_vec[:, -1] = G[: p.T]
+    # Map consumption goods back to demands for production goods
+    C_m_vec = (
+        np.tile(p.io_matrix.reshape(1, p.I, p.M), (p.T, 1, 1))
+        * np.tile(C_vec[: p.T, :].reshape(p.T, p.I, 1), (1, 1, p.M))
+    ).sum(axis=1)
     I_d_vec = np.zeros((p.T, p.M))
     I_d_vec[:, -1] = I_d[: p.T]
     I_g_vec = np.zeros((p.T, p.M))
@@ -1158,7 +1183,7 @@ def run_TPI(p, client=None):
     net_capital_outflows_vec = np.zeros((p.T, p.M))
     net_capital_outflows_vec[:, -1] = net_capital_outflows[: p.T]
     RC_error = aggr.resource_constraint(
-        Y_vec, C_vec, G_vec, I_d_vec, I_g_vec, net_capital_outflows_vec
+        Y_vec, C_m_vec, G_vec, I_d_vec, I_g_vec, net_capital_outflows_vec
     )
     # Compute total investment (not just domestic)
     I_total = aggr.get_I(None, K[1 : p.T + 1], K[: p.T], p, "total_tpi")
diff --git a/ogcore/__init__.py b/ogcore/__init__.py
index 77faa93a0..774c155d2 100644
--- a/ogcore/__init__.py
+++ b/ogcore/__init__.py
@@ -19,4 +19,4 @@
 from ogcore.txfunc import *
 from ogcore.utils import *
 
-__version__ = "0.9.2"
+__version__ = "0.10.0"
diff --git a/ogcore/aggregates.py b/ogcore/aggregates.py
index 5f9e846c4..81a44c6cb 100644
--- a/ogcore/aggregates.py
+++ b/ogcore/aggregates.py
@@ -341,6 +341,7 @@ def revenue(
     bq_tax_liab = tax.bequest_tax_liab(r, b, bq, 0, None, method, p)
     w_tax_liab = tax.wealth_tax_liab(r, b, 0, None, method, p)
     if method == "SS":
+        p_i = np.dot(p.io_matrix, p_m)
         pop_weights = np.transpose(p.omega_SS * p.lambdas)
         iit_payroll_tax_revenue = (inc_pay_tax_liab * pop_weights).sum()
         agg_pension_outlays = (pension_benefits * pop_weights).sum()
@@ -348,11 +349,15 @@ def revenue(
         wealth_tax_revenue = (w_tax_liab * pop_weights).sum()
         bequest_tax_revenue = (bq_tax_liab * pop_weights).sum()
         cons_tax_revenue = (
-            ((p.tau_c[-1, :] * p_m).reshape(p.M, 1, 1) * c).sum(axis=0)
+            ((p.tau_c[-1, :] * p_i).reshape(p.I, 1, 1) * c).sum(axis=0)
             * pop_weights
         ).sum()
         payroll_tax_revenue = p.frac_tax_payroll[-1] * iit_payroll_tax_revenue
     elif method == "TPI":
+        p_i = (
+            np.tile(p.io_matrix.reshape(1, p.I, p.M), (p.T, 1, 1))
+            * np.tile(p_m[: p.T, :].reshape(p.T, 1, p.M), (1, p.I, 1))
+        ).sum(axis=2)
         pop_weights = np.squeeze(p.lambdas) * np.tile(
             np.reshape(p.omega[: p.T, :], (p.T, p.S, 1)), (1, 1, p.J)
         )
@@ -365,7 +370,7 @@ def revenue(
         bequest_tax_revenue = (bq_tax_liab * pop_weights).sum(1).sum(1)
         cons_tax_revenue = (
             (
-                ((p.tau_c[: p.T, :] * p_m).reshape(p.T, p.M, 1, 1) * c).sum(
+                ((p.tau_c[: p.T, :] * p_i).reshape(p.T, p.I, 1, 1) * c).sum(
                     axis=1
                 )
                 * pop_weights
@@ -548,25 +553,26 @@ def get_K_splits(B, K_demand_open, D_d, zeta_K):
     return K, K_d, K_f
 
 
-def get_ptilde(p_m, tau_c, alpha_c, method="SS"):
+def get_ptilde(p_i, tau_c, alpha_c, method="SS"):
     r"""
     Calculate price of composite good.
 
     .. math::
-        \tilde{p}_{t} = \prod_{m=1}^{M} \left(\frac{(1 + \tau^{c}_{m,t})p_{m,j}}{\alpha_{m,j}}\right)^{\alpha_{m,j}}
+        \tilde{p}_{t} = \prod_{i=1}^{I} \left(\frac{(1 +
+        \tau^{c}_{i,t})p_{i,j}}{\alpha_{i,j}}\right)^{\alpha_{i,j}}
 
     Args:
-        p_m (array_like): prices for consumption good m
-        tau_c (array_like): consumption taxes on good m
+        p_i (array_like): prices for consumption good i
+        tau_c (array_like): consumption taxes on good i
         alpha_c (array_like): consumption share parameters
 
     Returns:
         p_tilde (array_like): tax-inclusive price of composite good
     """
     if method == "SS":
-        p_tilde = np.prod((((1 + tau_c) * p_m) / alpha_c) ** alpha_c)
+        p_tilde = np.prod((((1 + tau_c) * p_i) / alpha_c) ** alpha_c)
     else:  # TPI case
         alpha_c = alpha_c.reshape(1, alpha_c.shape[0])
-        p_tilde = np.prod((((1 + tau_c) * p_m) / alpha_c) ** alpha_c, axis=1)
+        p_tilde = np.prod((((1 + tau_c) * p_i) / alpha_c) ** alpha_c, axis=1)
 
     return p_tilde
diff --git a/ogcore/default_parameters.json b/ogcore/default_parameters.json
index ceeba470d..074635222 100644
--- a/ogcore/default_parameters.json
+++ b/ogcore/default_parameters.json
@@ -119,6 +119,25 @@
             }
         }
     },
+    "I": {
+        "title": "Number of different consumption goods",
+        "description": "Number of different consumption goods.",
+        "section_1": "Household Parameters",
+        "section_2": "Model Dimensions",
+        "notes": "",
+        "type": "int",
+        "value": [
+            {
+                "value": 1
+            }
+        ],
+        "validators": {
+            "range": {
+                "min": 1,
+                "max": 50
+            }
+        }
+    },
     "lambdas": {
         "title": "Fraction of population of each labor productivity type",
         "description": "Fraction of population of each labor productivity type.",
@@ -444,6 +463,28 @@
             }
         }
     },
+    "io_matrix": {
+        "title": "Input-output matrix",
+        "description": "Input-output matrix used to map production outputs into consumption goods using a fixed coefficient model.  This matrix has dimensions I x M, where I is the number of distinct consumption goods and M is the number of distinct production goods.  The sum each row of this matrix must be 1.",
+        "section_1": "Firm Parameters",
+        "section_2": "Production Function",
+        "notes": "",
+        "type": "float",
+        "number_dims": 2,
+        "value": [
+            {
+                "value": [[
+                    1.0
+                ]]
+            }
+        ],
+        "validators": {
+            "range": {
+                "min": 0.0,
+                "max": 1.0
+            }
+        }
+    },
     "Z": {
         "title": "Total factor productivity in firm production function",
         "description": "Total factor productivity in firm production function.  Set value for base year, click '+' to add value for next year.  All future years not specified are set to last value entered.",
@@ -906,7 +947,7 @@
         "description": "Consumption tax rate.  Set value for base year, click '+' to add value for next year.  All future years not specified are set to last value entered.",
         "section_1": "Fiscal Policy Parameters",
         "section_2": "Taxes",
-        "notes": "This policy parameter represents the effective consumption tax rate from sales taxes, VATs, and excise taxes.  To capture exemptions of certain goods, this is assumed to vary by age and ability type.  It is thus a SxJ matrix.",
+        "notes": "This policy parameter represents the effective consumption tax rate from sales taxes, VATs, and excise taxes by consumption good. Tax rates cab vary over time. It is thus a TxI array.",
         "type": "float",
         "number_dims": 2,
         "value": [
diff --git a/ogcore/household.py b/ogcore/household.py
index 6d7567969..ef2b45f4d 100644
--- a/ogcore/household.py
+++ b/ogcore/household.py
@@ -269,47 +269,47 @@ def get_cons(r, w, p_tilde, b, b_splus1, n, bq, net_tax, e, p):
     return cons
 
 
-def get_cm(c_s, p_m, p_tilde, tau_c, alpha_c, method="SS"):
+def get_ci(c_s, p_i, p_tilde, tau_c, alpha_c, method="SS"):
     r"""
-    Compute consumption of good m given amount of composite consumption
+    Compute consumption of good i given amount of composite consumption
     and prices.
 
     .. math::
-        c_{m,j,s,t} = \frac{c_{s,j,t}}{\alpha_{m,j}p_{m,j}}
+        c_{i,j,s,t} = \frac{c_{s,j,t}}{\alpha_{i,j}p_{i,j}}
 
     Args:
         c_s (array_like): composite consumption
-        p_m (array_like): prices for consumption good m
+        p_i (array_like): prices for consumption good i
         p_tilde (array_like): composite good price
         tau_c (array_like): consumption tax rate
         alpha_c (array_like): consumption share parameters
         method (str): adjusts calculation dimensions based on 'SS' or 'TPI'
 
     Returns:
-        c_sm (array_like): consumption of good m
+        c_si (array_like): consumption of good i
     """
     if method == "SS":
-        M = alpha_c.shape[0]
+        I = alpha_c.shape[0]
         S = c_s.shape[0]
         J = c_s.shape[1]
-        tau_c = tau_c.reshape(M, 1, 1)
-        alpha_c = alpha_c.reshape(M, 1, 1)
+        tau_c = tau_c.reshape(I, 1, 1)
+        alpha_c = alpha_c.reshape(I, 1, 1)
         p_tilde.reshape(1, 1, 1)
-        p_m = p_m.reshape(M, 1, 1)
+        p_i = p_i.reshape(I, 1, 1)
         c_s = c_s.reshape(1, S, J)
-        c_sm = alpha_c * (((1 + tau_c) * p_m) / p_tilde) ** (-1) * c_s
+        c_si = alpha_c * (((1 + tau_c) * p_i) / p_tilde) ** (-1) * c_s
     else:  # Time path case
-        M = alpha_c.shape[0]
-        T = p_m.shape[0]
+        I = alpha_c.shape[0]
+        T = p_i.shape[0]
         S = c_s.shape[1]
         J = c_s.shape[2]
-        tau_c = tau_c.reshape(T, M, 1, 1)
-        alpha_c = alpha_c.reshape(1, M, 1, 1)
+        tau_c = tau_c.reshape(T, I, 1, 1)
+        alpha_c = alpha_c.reshape(1, I, 1, 1)
         p_tilde = p_tilde.reshape(T, 1, 1, 1)
-        p_m = p_m.reshape(T, M, 1, 1)
+        p_i = p_i.reshape(T, I, 1, 1)
         c_s = c_s.reshape(T, 1, S, J)
-        c_sm = alpha_c * (((1 + tau_c) * p_m) / p_tilde) ** (-1) * c_s
-    return c_sm
+        c_si = alpha_c * (((1 + tau_c) * p_i) / p_tilde) ** (-1) * c_s
+    return c_si
 
 
 def FOC_savings(
diff --git a/ogcore/parameters.py b/ogcore/parameters.py
index a9c18f2c4..3354d336e 100644
--- a/ogcore/parameters.py
+++ b/ogcore/parameters.py
@@ -1,4 +1,3 @@
-from inspect import Parameter
 import os
 import numpy as np
 import scipy.interpolate as si
@@ -169,7 +168,6 @@ def compute_default_params(self):
             "Z",
             "delta_tau_annual",
             "cit_rate",
-            "tau_c",
             "inv_tax_credit",
         ]
         for item in tp_param_list2:
@@ -228,6 +226,64 @@ def compute_default_params(self):
                     )
                 )
             setattr(self, item, this_attr)
+        # Deal with parameters that vary across consumption good and over time
+        tp_param_list3 = ["tau_c"]
+        for item in tp_param_list3:
+            this_attr = getattr(self, item)
+            if this_attr.ndim == 1:
+                # case where enter single number, so assume constant
+                # across years and industries
+                if this_attr.shape[0] == 1:
+                    this_attr = (
+                        np.ones((self.T + self.S, self.I)) * this_attr[0]
+                    )
+                # case where user enters just one year for all industries
+                if this_attr.shape[0] == self.I:
+                    this_attr = np.tile(
+                        this_attr.reshape(1, self.I), (self.T + self.S, 1)
+                    )
+                else:
+                    # case where user enters multiple years for one industry
+                    # will assume they implied values the same across industries
+                    this_attr = np.concatenate(
+                        (
+                            this_attr,
+                            np.ones((self.T + self.S - this_attr.size))
+                            * this_attr[-1],
+                        )
+                    )
+                    this_attr = np.tile(
+                        this_attr.reshape(self.T + self.S, 1), (1, self.I)
+                    )
+                this_attr = np.squeeze(this_attr, axis=2)
+            elif this_attr.ndim == 2:
+                if this_attr.shape[1] > 1 and this_attr.shape[1] != self.I:
+                    print(
+                        "please provide values of "
+                        + item
+                        + " for each industry (or one if common across "
+                        + "industries"
+                    )
+                    assert False
+                if this_attr.shape[1] == 1:
+                    this_attr = np.tile(
+                        this_attr.reshape(this_attr.shape[0], 1), (1, self.I)
+                    )
+                if this_attr.shape[0] > self.T + self.S:
+                    this_attr = this_attr[: self.T + self.S, :]
+                this_attr = np.concatenate(
+                    (
+                        this_attr,
+                        np.ones(
+                            (
+                                self.T + self.S - this_attr.shape[0],
+                                this_attr.shape[1],
+                            )
+                        )
+                        * this_attr[-1, :],
+                    )
+                )
+            setattr(self, item, this_attr)
         # Deal with parameters that vary across J and over time
         tp_param_list3 = [
             "labor_income_tax_noncompliance_rate",
diff --git a/setup.py b/setup.py
index 0c6d2d550..6f31d12eb 100755
--- a/setup.py
+++ b/setup.py
@@ -5,7 +5,7 @@
 
 setuptools.setup(
     name="ogcore",
-    version="0.9.2",
+    version="0.10.0",
     author="Jason DeBacker and Richard W. Evans",
     license="CC0 1.0 Universal (CC0 1.0) Public Domain Dedication",
     description="A general equilibribum overlapping generations model for fiscal policy analysis",
diff --git a/tests/test_SS.py b/tests/test_SS.py
index 473e55f2c..139d0ec34 100644
--- a/tests/test_SS.py
+++ b/tests/test_SS.py
@@ -185,6 +185,8 @@ def dask_client():
 p7.update_specifications(
     {
         "M": 4,
+        "I": 4,
+        "io_matrix": np.eye(4),
         "alpha_c": [0.1, 0.5, 0.3, 0.1],
         "epsilon": [1.0, 1.0, 1.0, 1.0],
         "gamma": [0.3, 0.4, 0.35, 0.45],
@@ -448,12 +450,32 @@ def test_SS_solver_extra(baseline, param_updates, filename, dask_client):
 filename5 = "inner_loop_outputs_reform_baselinespending.pkl"
 param_updates7 = {
     "M": 4,
+    "I": 4,
+    "io_matrix": np.eye(4),
     "alpha_c": [0.1, 0.5, 0.3, 0.1],
     "epsilon": [1.0, 1.0, 1.0, 1.0],
     "gamma": [0.3, 0.4, 0.35, 0.45],
     "gamma_g": [0.0, 0.0, 0.0, 0.0],
 }
 filename7 = "inner_loop_outputs_reform_M4.pkl"
+param_updates8 = {
+    "M": 4,
+    "I": 5,
+    "io_matrix": np.array(
+        [
+            [0.2, 0.2, 0.2, 0.4],
+            [0.3, 0.1, 0.4, 0.2],
+            [0.25, 0.25, 0.25, 0.25],
+            [0.1, 0.7, 0.0, 0.2],
+            [0.0, 0.0, 1.0, 0.0],
+        ]
+    ),
+    "alpha_c": [0.1, 0.4, 0.3, 0.1, 0.1],
+    "epsilon": [1.0, 1.0, 1.0, 1.0],
+    "gamma": [0.3, 0.4, 0.35, 0.45],
+    "gamma_g": [0.0, 0.0, 0.0, 0.0],
+}
+filename8 = "inner_loop_outputs_reform_MneI.pkl"
 
 
 @pytest.mark.parametrize(
@@ -465,6 +487,7 @@ def test_SS_solver_extra(baseline, param_updates, filename, dask_client):
         (False, 0.04260341179572245, param_updates4, filename4),
         (False, 0.04260341179572245, param_updates5, filename5),
         (False, 0.04759112768438152, param_updates7, filename7),
+        (False, 0.04759112768438152, param_updates8, filename8),
     ],
     ids=[
         "Baseline, Small Open",
@@ -473,6 +496,7 @@ def test_SS_solver_extra(baseline, param_updates, filename, dask_client):
         "Reform",
         "Reform, baseline spending",
         "Reform, M>1",
+        "Reform, I!=>M",
     ],
 )
 def test_inner_loop(baseline, r_p, param_updates, filename, dask_client):
@@ -1032,6 +1056,8 @@ def test_euler_equation_solver(input_tuple, ubi_j, p, expected):
     "frisch": 0.41,
     "cit_rate": [[0.21, 0.25, 0.35]],
     "M": 3,
+    "I": 3,
+    "io_matrix": np.eye(3),
     "epsilon": [1.0, 1.0, 1.0],
     "gamma": [0.3, 0.35, 0.4],
     "gamma_g": [0.1, 0.05, 0.15],
@@ -1049,6 +1075,8 @@ def test_euler_equation_solver(input_tuple, ubi_j, p, expected):
     "frisch": 0.41,
     "cit_rate": [[0.21, 0.25, 0.35]],
     "M": 3,
+    "I": 3,
+    "io_matrix": np.eye(3),
     "epsilon": [1.0, 1.0, 1.0],
     "gamma": [0.3, 0.35, 0.4],
     "gamma_g": [0.0, 0.0, 0.0],
diff --git a/tests/test_TPI.py b/tests/test_TPI.py
index f627c3f9d..a9d41013b 100644
--- a/tests/test_TPI.py
+++ b/tests/test_TPI.py
@@ -266,6 +266,8 @@ def test_inner_loop():
     "frisch": 0.41,
     "cit_rate": [[0.21, 0.25, 0.35]],
     "M": 3,
+    "I": 3,
+    "io_matrix": np.eye(3),
     "epsilon": [1.0, 1.0, 1.0],
     "gamma": [0.3, 0.35, 0.4],
     "gamma_g": [0.1, 0.05, 0.15],
@@ -294,6 +296,8 @@ def test_inner_loop():
     "frisch": 0.41,
     "cit_rate": [[0.21, 0.25, 0.35]],
     "M": 3,
+    "I": 3,
+    "io_matrix": np.eye(3),
     "epsilon": [1.0, 1.0, 1.0],
     "gamma": [0.3, 0.35, 0.4],
     "gamma_g": [0.0, 0.0, 0.0],
@@ -307,6 +311,34 @@ def test_inner_loop():
 filename10 = os.path.join(
     CUR_PATH, "test_io_data", "run_TPI_baseline_M3_Kg_zero.pkl"
 )
+param_updates11 = {
+    "start_year": 2023,
+    "budget_balance": True,
+    "frisch": 0.41,
+    "cit_rate": [[0.21, 0.25, 0.35]],
+    "M": 3,
+    "I": 4,
+    "io_matrix": np.array(
+        [
+            [0.3, 0.3, 0.4],
+            [0.6, 0.1, 0.3],
+            [0.25, 0.5, 0.25],
+            [0.0, 1.0, 0.0],
+        ]
+    ),
+    "epsilon": [1.0, 1.0, 1.0],
+    "gamma": [0.3, 0.35, 0.4],
+    "gamma_g": [0.0, 0.0, 0.0],
+    "alpha_c": [0.2, 0.4, 0.3, 0.1],
+    "initial_guess_r_SS": 0.11,
+    "initial_guess_TR_SS": 0.07,
+    "debt_ratio_ss": 1.5,
+    "alpha_T": alpha_T.tolist(),
+    "alpha_G": alpha_G.tolist(),
+}
+filename11 = os.path.join(
+    CUR_PATH, "test_io_data", "run_TPI_baseline_MneI.pkl"
+)
 
 
 @pytest.mark.local
@@ -323,6 +355,7 @@ def test_inner_loop():
         (True, param_updates8, filename8),
         (True, param_updates9, filename9),
         (True, param_updates10, filename10),
+        (True, param_updates11, filename11),
     ],
     ids=[
         "Baseline, balanced budget",
@@ -335,6 +368,7 @@ def test_inner_loop():
         "Baseline, Kg > 0",
         "Baseline, M=3 non-zero Kg",
         "Baseline, M=3 zero Kg",
+        "Baseline, M!=I",
     ],
 )
 def test_run_TPI_full_run(
diff --git a/tests/test_household.py b/tests/test_household.py
index e8c8fdc9f..f237b79c9 100644
--- a/tests/test_household.py
+++ b/tests/test_household.py
@@ -1008,16 +1008,16 @@ def test_constraint_checker_TPI(bssmat, nssmat, cssmat, ltilde):
     assert True
 
 
-def test_get_cm():
+def test_ci():
     """
-    Test of the get_cm function
+    Test of the get_ci function
     """
     c_s = np.array([2.0, 3.0, 5.0, 7.0]).reshape(4, 1)
-    p_m = np.array([1.1, 0.8, 1.0])
+    p_i = np.array([1.1, 0.8, 1.0])
     p_tilde = np.array([2.3])
     tau_c = np.array([0.2, 0.3, 0.5])
     alpha_c = np.array([0.5, 0.3, 0.2])
-    expected_cm = np.array(
+    expected_ci = np.array(
         [
             [1.742424242, 2.613636364, 4.356060606, 6.098484848],
             [1.326923077, 1.990384615, 3.317307692, 4.644230769],
@@ -1025,6 +1025,6 @@ def test_get_cm():
         ]
     ).reshape(3, 4, 1)
 
-    test_cm = household.get_cm(c_s, p_m, p_tilde, tau_c, alpha_c)
+    test_ci = household.get_ci(c_s, p_i, p_tilde, tau_c, alpha_c)
 
-    assert np.allclose(test_cm, expected_cm)
+    assert np.allclose(test_ci, expected_ci)
diff --git a/tests/test_io_data/inner_loop_outputs_reform_MneI.pkl b/tests/test_io_data/inner_loop_outputs_reform_MneI.pkl
new file mode 100644
index 000000000..f3b3e277e
Binary files /dev/null and b/tests/test_io_data/inner_loop_outputs_reform_MneI.pkl differ
diff --git a/tests/test_io_data/run_TPI_baseline_M3_Kg_nonzero.pkl b/tests/test_io_data/run_TPI_baseline_M3_Kg_nonzero.pkl
index 8502a966a..7664d6f43 100644
Binary files a/tests/test_io_data/run_TPI_baseline_M3_Kg_nonzero.pkl and b/tests/test_io_data/run_TPI_baseline_M3_Kg_nonzero.pkl differ
diff --git a/tests/test_io_data/run_TPI_baseline_M3_Kg_zero.pkl b/tests/test_io_data/run_TPI_baseline_M3_Kg_zero.pkl
index 5ff087493..7c398aff2 100644
Binary files a/tests/test_io_data/run_TPI_baseline_M3_Kg_zero.pkl and b/tests/test_io_data/run_TPI_baseline_M3_Kg_zero.pkl differ
diff --git a/tests/test_io_data/run_TPI_baseline_MneI.pkl b/tests/test_io_data/run_TPI_baseline_MneI.pkl
new file mode 100644
index 000000000..403782898
Binary files /dev/null and b/tests/test_io_data/run_TPI_baseline_MneI.pkl differ