diff --git a/docs/book/content/theory/equilibrium.md b/docs/book/content/theory/equilibrium.md
index a7465353c..2dd49c006 100644
--- a/docs/book/content/theory/equilibrium.md
+++ b/docs/book/content/theory/equilibrium.md
@@ -112,7 +112,7 @@ The computational algorithm for solving for the steady-state follows the steps b
     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}$
-        2.  The capital-output ratio can be determined from the FOC for the firms' choice of capital: $\frac{\bar{K}_m}{\bar{Y}_m} = \frac{\bar{r} + \bar{\delta}_{M} -\bar{\delta}^{\tau}\bar{\tau}^{corp}_{m}}{\bar{p}_m(1-\bar{\tau}^{corp}_{m})\bar{Z}_m^{\frac{\varepsilon_m -1}{\varepsilon_m}}}^{-\varepsilon_m} \gamma_{m}$
+        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}{\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:
            ```{math}
@@ -123,7 +123,7 @@ The computational algorithm for solving for the steady-state follows the steps b
 
            ```{math}
            :label: EqSS_MPKg
-             \bar{K}_m^{r^*} = \bar{L}_m\left(\frac{\bar{w}}{\frac{\bar{r} + \bar{\delta}_M - \bar{\tau}^b_m\bar{\delta}^{\tau}_m}{1 - \bar{\tau}_m^b}}\right)^{\varepsilon_m} \frac{\gamma_m}{(1 - \gamma_m - \gamma_{g,m})}
+             \bar{K}_m^{r^*} = \bar{L}_m\left(\frac{\bar{w}}{\frac{\bar{r} + \bar{\delta}_M - \bar{\tau}^{corp}_m\bar{\delta}^{\tau}_m - \bar{\tau}^{inv}_m}{1 - \bar{\tau}^{corp}_m}}\right)^{\varepsilon_m} \frac{\gamma_m}{(1 - \gamma_m - \gamma_{g,m})}
            ```
 
     16. Determine factor demands and output for industry $M$:
@@ -134,8 +134,8 @@ The computational algorithm for solving for the steady-state follows the steps b
             2.  Using $\bar{D}^{d,i}$ we can then find domestic investors' holdings of private capital as the residual from their total asset holdings: , $\bar{K}^{d,i} = \bar{B}^i - \bar{D}^{d,i}$
             3.  Aggregate capital supply is then determined as $\bar{K}^{i'} = \bar{K}^{d,i} + \bar{K}^{f,i}$.
         4.  $\bar{K}_M = \bar{K}^{i'} - \sum_{m=1}^{M-1}\bar{K}_{m}$
-        5.  Use the factor demands and $\bar{K}_g$ in the production function for industry $M$ to find $\bar{Y}_M
-    17. Find an updated value for GDP, $\bar{Y}^{i'} = \sum_{m=1}^{M} \bar{p}_m \bar{Y}_m
+        5.  Use the factor demands and $\bar{K}_g$ in the production function for industry $M$ to find $\bar{Y}_M$.
+    17. Find an updated value for GDP, $\bar{Y}^{i'} = \sum_{m=1}^{M} \bar{p}_m \bar{Y}_m$.
     18. Find a updated values for $\bar{I}_{g}$ and $\bar{K}_g$ using  $\bar{Y}^{i'}$, equations {eq}`EqStnrz_Igt` and {eq}`EqStnrz_Kgmt`
 3. Given updated inner-loop values based on initial guesses for outer-loop variables $\{\bar{r}_p^i, \bar{r}^i, \bar{w}^i, \boldsymbol{\bar{p}}, \overline{BQ}^i, \overline{TR}^i, factor^i\}$, solve for updated values of outer-loop variables $\{\bar{r}_p^{i'}, \bar{r}^{i'}, \bar{w}^{i'}, \boldsymbol{\bar{p}}^{i'}, \overline{BQ}^{i'}, \overline{TR}^{i'}, factor^{i'}\}$ using the remaining equations:
 
@@ -339,7 +339,7 @@ The stationary non-steady state (transition path) solution algorithm has followi
 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}$
-    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}} = \frac{r_t + \delta_{M,t} -\delta_t^{\tau}\tau^{corp}_{m,t}}{p_{m,t}(1-\tau^{corp}_{m,t}){Z}_{m,t}^{\frac{\varepsilon_m -1}{\varepsilon_m}}}^{-\varepsilon_m} \gamma_{m}$
+    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}}{(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:
       ```{math}
@@ -349,7 +349,7 @@ The stationary non-steady state (transition path) solution algorithm has followi
     5\. Use the interest rate $r_t^*$ and labor demand $\hat{L}_{m,t}$ to solve for private capital demand at the world interest rate $\hat{K}_{m,t}^{r^*}$ using {eq}`EqFirmsMPKg_opt`
       ```{math}
       :label: EqTP_MPKg
-        \hat{K}_{m,t}^{r^*} = \hat{L}_{m,t}\left(\frac{\hat{w}_t}{\frac{r_t + \delta_{M,t} - \tau^b_{m,t}\delta^{\tau}_{m,t}}{1 - \bar{\tau}_{m,t}^b}}\right)^{\varepsilon_m} \frac{\gamma_m}{(1 - \gamma_m - \gamma_{g,m})}
+        \hat{K}_{m,t}^{r^*} = \hat{L}_{m,t}\left(\frac{\hat{w}_t}{\frac{r_t + \delta_{M,t} - \tau^{corp}_{m,t}\delta^{\tau}_{m,t} - \tau^{inv}_{m,t}}{1 - \bar{\tau}_{m,t}^b}}\right)^{\varepsilon_m} \frac{\gamma_m}{(1 - \gamma_m - \gamma_{g,m})}
       ```
 
 18.  Determine factor demands and output for industry $M$:
diff --git a/docs/book/content/theory/firms.md b/docs/book/content/theory/firms.md
index dd20a2348..ce540ac72 100644
--- a/docs/book/content/theory/firms.md
+++ b/docs/book/content/theory/firms.md
@@ -34,14 +34,14 @@ Industry $M$ in the model is unique in two respects.  First, we will define indu
   ```{math}
   :label: EqFirmsProfit
     PR_{m,t} &= (1 - \tau^{corp}_{m,t})\Bigl[p_{m,t}F(K_{m,t},K_{g,m,t},L_{m,t}) - w_t L_{m,t}\Bigr] - \\
-    &\qquad\qquad\quad \bigl(r_t + \delta_{M,t}\bigr)K_{m,t} + \tau^{corp}_{m,t}\delta^\tau_{m,t} K_{m,t} \quad\forall m,t
+    &\qquad\qquad\quad \bigl(r_t + \delta_{M,t}\bigr)K_{m,t} + \tau^{corp}_{m,t}\delta^\tau_{m,t}K_{m,t} + \tau^{inv}_{m,t}I_{m,t} \quad\forall m,t
   ```
 
   Gross income for the firms is $p_{m,t}F(K_{m,t},K_{g,m,t},L_{m,t})$. Labor costs to the firm are $w_t L_{m,t}$, and capital costs are $(r_t +\delta_{M,t})K_{m,t}$. The government supplies public capital $K_{g,m,t}$ to the firms at no cost. The per-period interest rate (rental rate) of capital for firms is $r_t$. The per-period economic depreciation rate for private capital is $\delta_{M,t}\in[0,1]$.[^delta_M] The $\delta^\tau_{m,t}$ parameter in the last term of the profit function governs how much of capital depreciation can be deducted from the corporate income tax.
 
   Taxes enter the firm's profit function {eq}`EqFirmsProfit` in two places. The first is the corporate income tax rate $\tau^{corp}_{m,t}$, which is a flat tax on corporate income that can vary by industry $m$. Corporate income is defined as gross income minus labor costs. This will cause the corporate tax to only have a direct effect on the firms' capital demand decision.
 
-  The tax policy also enters the profit function {eq}`EqFirmsProfit` through depreciation deductions at rate $\delta^\tau_{m,t}$, which then lower corporate tax liability. When $\delta^\tau_{m,t}=0$, no depreciation expense is deducted from the firm's tax liability. When $\delta^\tau_{m,t}=\delta_{M,t}$, all economic depreciation is deducted from corporate income.
+  The tax policy also enters the profit function {eq}`EqFirmsProfit` through depreciation deductions at rate $\delta^\tau_{m,t}$, which then lower corporate tax liability. When $\delta^\tau_{m,t}=0$, no depreciation expense is deducted from the firm's tax liability. When $\delta^\tau_{m,t}=\delta_{M,t}$, all economic depreciation is deducted from corporate income. The investment tax credit is characterized by the parameter $\tau^{inv}_{m,t}$ multiplied by $I_{m,t}$ investment by industry $m$, which represents the percent of industry-specific capital that can be credited against profits.
 
   Firms take as given prices $p_{m,t}$, $w_t$, and $r_t$ and the level of public capital supply $K_{g,m,t}$. Taking the derivative of the profit function {eq}`EqFirmsProfit` with respect to labor $L_{m,t}$ and setting it equal to zero (using the general CES form of the production function {eq}`EqFirmsCESprodfun`) and taking the derivative of the profit function with respect to private capital $K_{m,t}$ and setting it equal to zero, respectively, characterizes the optimal labor and capital demands.
 
@@ -52,7 +52,7 @@ Industry $M$ in the model is unique in two respects.  First, we will define indu
 
   ```{math}
   :label: EqFirmFOC_K
-    r_t = (1 - \tau^{corp}_{m,t})p_{m,t}(Z_{m,t})^\frac{\varepsilon_m-1}{\varepsilon_m}\left[\gamma_m\frac{Y_{m,t}}{K_{m,t}}\right]^\frac{1}{\varepsilon_m} - \delta_{M,t} + \tau^{corp}_{m,t}\delta^\tau_{m,t} \quad\forall m,t
+    r_t = (1 - \tau^{corp}_{m,t})p_{m,t}(Z_{m,t})^\frac{\varepsilon_m-1}{\varepsilon_m}\left[\gamma_m\frac{Y_{m,t}}{K_{m,t}}\right]^\frac{1}{\varepsilon_m} - \delta_{M,t} + \tau^{corp}_{m,t}\delta^\tau_{m,t} + \tau^{inv}_{m,t} \quad\forall m,t
   ```
 
   Note that the presence of the public capital good creates economic rents. These rents will accrue to the owners of capital via the financial intermediary. See Section Chapter {ref}`Chap_FinInt` for more details on the determination of the return to the household's portfolio. Because public capital is exogenous to the firm's decisions, the optimality condition for capital demand {eq}`EqFirmFOC_K` is only affected by public capital $K_{g,m,t}$ through the $Y_{m,t}$ term.
@@ -71,7 +71,7 @@ Industry $M$ in the model is unique in two respects.  First, we will define indu
 
   ```{math}
   :label: EqFirmsMPK_opt
-    MPK_{m,t} =  \frac{r_t + \delta_{M,t} - \tau^{corp}_{m,t}\delta^\tau_{m,t}}{p_{m,t}(1 - \tau^{corp}_{m,t})} \quad\forall m,t
+    MPK_{m,t} =  \frac{r_t + \delta_{M,t} - \tau^{corp}_{m,t}\delta^\tau_{m,t} - \tau^{inv}_{m,t}}{p_{m,t}(1 - \tau^{corp}_{m,t})} \quad\forall m,t
   ```
 
   Firm profit maximization for labor demand from equation {eq}`EqFirmFOC_L` implies that the marginal product of labor is equal to the real wage rate:
@@ -93,9 +93,9 @@ Industry $M$ in the model is unique in two respects.  First, we will define indu
   ```{math}
   :label: EqFirmsProfit_Kg
     \begin{split}
-      PR_{m,t} &= (1 - \tau^{corp}_{m,t})\Bigl[p_{m,t}Y_{m,t} - w_t L_{m,t}\Bigr] - \bigl(r_t + \delta_{M,t}\bigr)K_{m,t} + \tau^{corp}_{m,t}\delta^\tau_{m,t} K_{m,t} \\
-      &= (1 - \tau^{corp}_{m,t})\Biggl[\biggl(\frac{r_t + \delta_{M,t} - \tau^{corp}_{m,t}\delta^{\tau}_{m,t}}{1 - \tau^{corp}_{m,t}}\biggr)K_{m,t} + p_{m,t}MPK_{g,m,t}K_{g,m,t} + w_t L_{m,t}\Biggr] ... \\
-      &\quad\quad - (1 - \tau^{corp}_{m,t})w_t L_{m,t} - (r_t + \delta_{M,t})K_{m,t} + \tau^{corp}_{m,t}\delta^{\tau}_{m,t} K_{m,t} \\
+      PR_{m,t} &= (1 - \tau^{corp}_{m,t})\Bigl[p_{m,t}Y_{m,t} - w_t L_{m,t}\Bigr] - \bigl(r_t + \delta_{M,t}\bigr)K_{m,t} + \tau^{corp}_{m,t}\delta^\tau_{m,t}K_{m,t} + \tau^{inv}_{m,t}I_{m,t} \\
+      &= (1 - \tau^{corp}_{m,t})\Biggl[\biggl(\frac{r_t + \delta_{M,t} - \tau^{corp}_{m,t}\delta^{\tau}_{m,t} - \tau^{inv}_{m,t}}{1 - \tau^{corp}_{m,t}}\biggr)K_{m,t} + p_{m,t}MPK_{g,m,t}K_{g,m,t} + w_t L_{m,t}\Biggr] ... \\
+      &\quad\quad - (1 - \tau^{corp}_{m,t})w_t L_{m,t} - (r_t + \delta_{M,t})K_{m,t} + \tau^{corp}_{m,t}\delta^{\tau}_{m,t}K_{m,t} + \tau^{inv}_{m,t}I_{m,t} \\
       &= (1 - \tau^{corp}_{m,t})p_{m,t}MPK_{g,m,t}K_{g,m,t} \quad\forall m,t
     \end{split}
   ```
diff --git a/docs/book/content/theory/government.md b/docs/book/content/theory/government.md
index 704f0ad68..07bd78d9f 100644
--- a/docs/book/content/theory/government.md
+++ b/docs/book/content/theory/government.md
@@ -384,14 +384,14 @@ Businesses face a linear tax rate $\tau^{b}_{m,t}$, which can vary by industry a
   ```{math}
   :label: EqFirmsProfit2
     PR_{m,t} &= (1 - \tau^{corp}_{m,t})\Bigl[p_{m,t}F(K_{m,t},K_{g,m,t},L_{m,t}) - w_t L_{m,t}\Bigr] - \\
-    &\quad\quad\quad\quad\quad \bigl(r_t + \delta_{M,t}\bigr)K_{m,t} + \tau^{corp}_{m,t}\delta^\tau_{m,t} K_{m,t} \quad\forall m,t
+    &\quad\quad\quad\quad\quad \bigl(r_t + \delta_{M,t}\bigr)K_{m,t} + \left(\tau^{corp}_{m,t}\delta^\tau_{m,t} + \tau^{inv}_{m,t}\right)K_{m,t} \quad\forall m,t
   ```
 
   We define total government revenue from taxes in terms of the numeraire good as the following,
 
   ```{math}
   :label: EqUnbalGBCgovRev
-    Rev_t &= \underbrace{\sum_{m=1}^M\Bigl[\tau^{corp}_{m,t}\bigl(p_{m,t}Y_{m,t} - w_t L_t\bigr) - \tau^{corp}_{m,t}\delta^\tau_{m,t}K_{m,t}\Bigr]}_{\text{corporate tax revenue}} \\
+    Rev_t &= \underbrace{\sum_{m=1}^M\Bigl[\tau^{corp}_{m,t}\bigl(p_{m,t}Y_{m,t} - w_t L_t\bigr) - \tau^{corp}_{m,t}\delta^\tau_{m,t}K_{m,t} - \tau^{inv}_{m,t}I_{m,t}\Bigr]}_{\text{corporate tax revenue}} \\
     &\quad + \underbrace{\sum_{s=E+1}^{E+S}\sum_{j=1}^J\lambda_j\omega_{s,t}\tau^{etr}_{s,t}\left(x_{j,s,t},y_{j,s,t}\right)\bigl(x_{j,s,t} + y_{j,s,t}\bigr)}_{\text{household tax revenue}} \quad\forall t
   ```
 
diff --git a/docs/book/content/theory/stationarization.md b/docs/book/content/theory/stationarization.md
index 2d7657d46..bb4b923fc 100644
--- a/docs/book/content/theory/stationarization.md
+++ b/docs/book/content/theory/stationarization.md
@@ -133,7 +133,7 @@ The usual definition of equilibrium would be allocations and prices such that ho
   ```{math}
   :label: EqStnrzProfit
     \hat{PR}_{m,t} &= (1 - \tau^{corp}_{m,t})\Bigl[F(\hat{K}_{m,t},\hat{K}_{g,m,t},\hat{L}_{m,t}) - \hat{w}_t \hat{L}_{m,t}\Bigr] - ... \\
-    &\qquad\qquad\quad \bigl(r_t + \delta_{M,t}\bigr)\hat{K}_{m,t} + \tau^{corp}_{m,t}\delta^\tau_{m,t} \hat{K}_{m,t} \quad\forall m,t
+    &\qquad\qquad\quad \bigl(r_t + \delta_{M,t}\bigr)\hat{K}_{m,t} + \tau^{corp}_{m,t}\delta^\tau_{m,t}\hat{K}_{m,t} + \tau^{inv}_{m,t}\hat{I}_{m,t} \quad\forall m,t
   ```
 
   The firms' first order equation for labor demand {eq}`EqFirmFOC_L` is stationarized by dividing both sides by $e^{g_y t}$. This stationarizes the wage $\hat{w}_t$ on the left-hand-side and cancels out the $e^{g_y t}$ term in front of the right-hand-side. To complete the stationarization, we multiply and divide the $\frac{Y_{m,t}}{e^{g_y t}L_{m,t}}$ term on the right-hand-side by $\tilde{N}_t$.
@@ -147,7 +147,7 @@ The usual definition of equilibrium would be allocations and prices such that ho
 
   ```{math}
   :label: EqStnrzFOC_K
-    r_t = (1 - \tau^{corp}_{m,t})p_{m,t}(Z_{m,t})^\frac{\varepsilon_m-1}{\varepsilon_m}\left[\gamma_m\frac{\hat{Y}_{m,t}}{\hat{K}_{m,t}}\right]^\frac{1}{\varepsilon_m} - \delta_{M,t} + \tau^{corp}_{m,t}\delta^\tau_{m,t} \quad\forall m,t
+    r_t = (1 - \tau^{corp}_{m,t})p_{m,t}(Z_{m,t})^\frac{\varepsilon_m-1}{\varepsilon_m}\left[\gamma_m\frac{\hat{Y}_{m,t}}{\hat{K}_{m,t}}\right]^\frac{1}{\varepsilon_m} - \delta_{M,t} + \tau^{corp}_{m,t}\delta^\tau_{m,t} + \tau^{inv}_{m,t} \quad\forall m,t
   ```
 
   A stationary version of the firms' gross revenue attributed to each factor of production {eq}`EqFirmsMargRevEq` is found by dividing both sides of the equation by $e^{g_y t}\tilde{N}_t$.
@@ -205,7 +205,7 @@ The usual definition of equilibrium would be allocations and prices such that ho
   We can stationarize the expression for total government revenue $Rev_t$ in {eq}`EqUnbalGBCgovRev` by dividing both sides of the equation by $e^{g_y t}\tilde{N}_t$.
   ```{math}
   :label: EqStnrzGovRev
-    \hat{Rev}_t &= \underbrace{\sum_{m=1}^M\Bigl[\tau^{corp}_{m,t}\bigl(p_{m,t}\hat{Y}_{m,t} - \hat{w}_t\hat{L}_t\bigr) - \tau^{corp}_{m,t}\delta^\tau_{m,t}\hat{K}_{m,t}\Bigr]}_{\text{corporate tax revenue}} \\
+    \hat{Rev}_t &= \underbrace{\sum_{m=1}^M\Bigl[\tau^{corp}_{m,t}\bigl(p_{m,t}\hat{Y}_{m,t} - \hat{w}_t\hat{L}_t\bigr) - \tau^{corp}_{m,t}\delta^\tau_{m,t}\hat{K}_{m,t} - \tau^{inv}_{m,t}\hat{I}_{m,t}\Bigr]}_{\text{corporate tax revenue}} \\
     &\qquad + \underbrace{\sum_{s=E+1}^{E+S}\sum_{j=1}^J\lambda_j\hat{\omega}_{s,t}\tau^{etr}_{s,t}\left(\hat{x}_{j,s,t},\hat{y}_{j,s,t}\right)\bigl(\hat{x}_{j,s,t} + \hat{y}_{j,s,t}\bigr)}_{\text{household tax revenue}} \quad\forall t
   ```
 
diff --git a/ogcore/SS.py b/ogcore/SS.py
index 055e29d69..3c9ae0435 100644
--- a/ogcore/SS.py
+++ b/ogcore/SS.py
@@ -1152,42 +1152,45 @@ def run_SS(p, client=None):
         # Use the baseline solution to get starting values for the reform
         baseline_ss_dir = os.path.join(p.baseline_dir, "SS", "SS_vars.pkl")
         ss_solutions = utils.safe_read_pickle(baseline_ss_dir)
-        # use baseline solution as starting values if dimensions match
-        try:
-            if ss_solutions["bssmat_splus1"].shape == (
-                p.S,
-                p.J,
-            ) and ss_solutions["Y_vec_ss"].shape == (p.M):
-                print("Using previous solutions for SS")
-                (
-                    b_guess,
-                    n_guess,
-                    r_p_guess,
-                    rguess,
-                    wguess,
-                    p_m_guess,
-                    BQguess,
-                    TRguess,
-                    Yguess,
-                    factor,
-                ) = (
-                    ss_solutions["bssmat_splus1"],
-                    ss_solutions["nssmat"],
-                    float(ss_solutions["r_p_ss"]),
-                    float(ss_solutions["rss"]),
-                    float(ss_solutions["wss"]),
-                    ss_solutions[
-                        "p_m_ss"
-                    ],  # Not sure why need to index p_m,but otherwise its shape is off..
-                    ss_solutions["BQss"],
-                    float(ss_solutions["TR_ss"]),
-                    float(ss_solutions["Yss"]),
-                    ss_solutions["factor_ss"],
-                )
-                use_new_guesses = False
-            else:
+        if p.reform_use_baseline_solution:
+            # use baseline solution as starting values if dimensions match
+            try:
+                if ss_solutions["bssmat_splus1"].shape == (
+                    p.S,
+                    p.J,
+                ) and ss_solutions["Y_vec_ss"].shape == (p.M):
+                    print("Using previous solutions for SS")
+                    (
+                        b_guess,
+                        n_guess,
+                        r_p_guess,
+                        rguess,
+                        wguess,
+                        p_m_guess,
+                        BQguess,
+                        TRguess,
+                        Yguess,
+                        factor,
+                    ) = (
+                        ss_solutions["bssmat_splus1"],
+                        ss_solutions["nssmat"],
+                        float(ss_solutions["r_p_ss"]),
+                        float(ss_solutions["rss"]),
+                        float(ss_solutions["wss"]),
+                        ss_solutions[
+                            "p_m_ss"
+                        ],  # Not sure why need to index p_m,but otherwise its shape is off..
+                        ss_solutions["BQss"],
+                        float(ss_solutions["TR_ss"]),
+                        float(ss_solutions["Yss"]),
+                        ss_solutions["factor_ss"],
+                    )
+                    use_new_guesses = False
+                else:
+                    use_new_guesses = True
+            except KeyError:
                 use_new_guesses = True
-        except KeyError:
+        else:
             use_new_guesses = True
         if use_new_guesses:
             if p.use_zeta:
diff --git a/ogcore/default_parameters.json b/ogcore/default_parameters.json
index 87a9308df..5f2ae9c7c 100644
--- a/ogcore/default_parameters.json
+++ b/ogcore/default_parameters.json
@@ -879,6 +879,28 @@
             }
         }
     },
+    "inv_tax_credit": {
+        "title": "Investment tax credit rate",
+        "description": "Investment tax credit rate.  This credit reduced the cost of new investment by the specified rate.",
+        "section_1": "Fiscal Policy Parameters",
+        "section_2": "Taxes",
+        "notes": "This credit reduced the cost of new investment by the specified rate",
+        "type": "float",
+        "number_dims": 2,
+        "value": [
+            {
+                "value": [
+                    [0.0]
+                ]
+            }
+        ],
+        "validators": {
+            "range": {
+                "min": -1.0,
+                "max": 1.0
+            }
+        }
+    },
     "tau_c": {
         "title": "Consumption tax rate",
         "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.",
@@ -3215,6 +3237,19 @@
             }
         }
     },
+    "reform_use_baseline_solution":{
+        "title": "Whether or not the baseline SS solution is used for starting values when solving the reform",
+        "description": "Whether or not the baseline SS solution is used for starting values when solving the reform.",
+        "section_1": "Model Solution Parameters",
+        "notes": "",
+        "type": "bool",
+        "value": [
+            {
+                "value": true
+            }
+        ],
+        "validators": {}
+    },
     "initial_guess_r_SS": {
         "title": "Initial guess of r for the SS solution",
         "description": "Initial guess of r for the SS solution.",
diff --git a/ogcore/firm.py b/ogcore/firm.py
index 915a65998..1b6fb5cd1 100644
--- a/ogcore/firm.py
+++ b/ogcore/firm.py
@@ -179,7 +179,7 @@ def get_r(Y, K, p_m, p, method, m=-1):
     .. math::
         r_{t} = (1 - \tau^{corp}_t)Z_t^\frac{\varepsilon-1}{\varepsilon}
         \left[\gamma\frac{Y_t}{K_t}\right]^\frac{1}{\varepsilon} -
-        \delta + \tau^{corp}_t\delta^\tau_t
+        \delta + \tau^{corp}_t\delta^\tau_t + \tau^{inv}_{m,t}
 
     Args:
         Y (array_like): aggregate output
@@ -197,13 +197,15 @@ def get_r(Y, K, p_m, p, method, m=-1):
     if method == "SS":
         delta_tau = p.delta_tau[-1, m]
         tau_b = p.tau_b[-1, m]
+        tau_inv = p.inv_tax_credit[-1, m]
         p_mm = p_m[m]
     else:
         delta_tau = p.delta_tau[: p.T, m].reshape(p.T, 1)
         tau_b = p.tau_b[: p.T, m].reshape(p.T, 1)
+        tau_inv = p.inv_tax_credit[: p.T, m].reshape(p.T, 1)
         p_mm = p_m[:, m].reshape(p.T, 1)
     MPK = get_MPx(Y, K, p.gamma[m], p, method, m)
-    r = (1 - tau_b) * p_mm * MPK - p.delta + tau_b * delta_tau
+    r = (1 - tau_b) * p_mm * MPK - p.delta + tau_b * delta_tau + tau_inv
 
     return r
 
@@ -437,45 +439,6 @@ def get_K_KLonly(L, r, p, method, m=-1):
     return K
 
 
-def get_K_from_Y(Y, r, p, method):
-    r"""
-    Generates vector of aggregate capital. Use with the open economy
-    options.
-
-    .. math::
-        K_{t} = \frac{Y_{t}}{Y_{t}/K_{t}} \\
-        K_{t} = \frac{\gamma Z_t^{\varepsilon -1} Y_t}{
-            \left(\frac{r_t + \delta - \tau_t^{corp}\delta_t^\tau}
-            {1 - \tau_{t}^{corp}}\right)^\varepsilon}
-
-    Args:
-        Y (array_like): aggregate output
-        r (array_like): the real interest rate
-        p (OG-Core Specifications object): model parameters
-        method (str): adjusts calculation dimensions based on 'SS' or
-            'TPI'
-
-    Returns:
-        r (array_like): the real interest rate
-
-    """
-    if method == "SS":
-        Z = p.Z[-1]
-        tau_b = p.tau_b[-1]
-        delta_tau = p.delta_tau[-1]
-    else:
-        Z = p.Z[: p.T]
-        tau_b = p.tau_b[: p.T]
-        delta_tau = p.delta_tau[: p.T]
-    numerator = p.gamma * Z ** (p.epsilon - 1) * Y
-    denominator = (
-        (r + p.delta - tau_b * delta_tau) / (1 - tau_b)
-    ) ** p.epsilon
-    K = numerator / denominator
-
-    return K
-
-
 def get_L_from_Y(w, Y, p, method):
     r"""
     Find aggregate labor L from output Y and wages w
@@ -506,48 +469,6 @@ def get_L_from_Y(w, Y, p, method):
     return L
 
 
-def get_K_from_Y_and_L(Y, L, K_g, p, method):
-    r"""
-    Find aggregate private capital K from output Y, aggregate labor L,
-    and public capital K_g
-
-    .. math::
-        K_{t} = \left(\frac{\left(\frac{Y_t}{Z_t}\right)^{\frac{\varepsilon-1}
-        {\varepsilon}} -
-        (1-\gamma-\gamma_g)L_t^{\frac{\varepsilon-1}{\varepsilon}} -
-        \gamma_g^{\frac{1}{\varepsilon}}K_{g,t}^{\frac{\varepsilon-1}{\varepsilon}}}
-        {\gamma^{\frac{1}{\varepsilon}}}\right)^{\frac{\varepsilon}{\varepsilon-1}}
-
-    Args:
-        w (array_like): the wage rate
-        Y (array_like): aggregate output
-        L (array_like): aggregate labor
-        K_g (array_like): aggregate public capital
-        p (OG-Core Specifications object): model parameters
-        method (str): adjusts calculation dimensions based on 'SS' or
-            'TPI'
-
-    Returns:
-        K (array_like): firm capital demand
-
-    """
-    if method == "SS":
-        Z = p.Z[-1]
-    else:
-        Z = p.Z[: p.T]
-    K = (
-        (
-            (Y / Z) ** ((p.epsilon - 1) / p.epsilon)
-            - (1 - p.gamma - p.gamma_g) * L ** ((p.epsilon - 1) / p.epsilon)
-            - (p.gamma_g ** (1 / p.epsilon))
-            * (K_g ** ((p.epsilon - 1) / p.epsilon))
-        )
-        / (p.gamma ** (1 / p.epsilon))
-    ) ** (p.epsilon / (p.epsilon - 1))
-
-    return K
-
-
 def get_K(r, w, L, p, method, m=-1):
     r"""
     Get K from r, w, L.  For determining capital demand for open
@@ -580,7 +501,7 @@ def get_cost_of_capital(r, p, method, m=-1):
     Compute the cost of capital.
 
     .. math::
-        \rho_{m,t} = \frac{r_{t} + \delta_{M,t} - \tau^{b}_{m,t} \delta^{\tau}_{m,t}}{1 - \tau^{b}_{m,t}}
+        \rho_{m,t} = \frac{r_{t} + \delta_{M,t} - \tau^{b}_{m,t} \delta^{\tau}_{m,t} - \tau^{inv}_{m,t}}{1 - \tau^{b}_{m,t}}
 
     Args:
         r (array_like): the real interest rate
@@ -595,20 +516,24 @@ def get_cost_of_capital(r, p, method, m=-1):
         if method == "SS":
             tau_b = p.tau_b[-1, :]
             delta_tau = p.delta_tau[-1, :]
+            tau_inv = p.inv_tax_credit[-1, :]
         else:
             tau_b = p.tau_b[: p.T, :]
             delta_tau = p.delta_tau[: p.T, :]
+            tau_inv = p.inv_tax_credit[: p.T, :]
             r = r.reshape(p.T, 1)
     else:
         if method == "SS":
             tau_b = p.tau_b[-1, m]
             delta_tau = p.delta_tau[-1, m]
+            tau_inv = p.inv_tax_credit[-1, m]
         else:
             tau_b = p.tau_b[: p.T, m]
             delta_tau = p.delta_tau[: p.T, m]
+            tau_inv = p.inv_tax_credit[: p.T, m]
             r = r.reshape(p.T)
 
-    cost_of_capital = (r + p.delta - tau_b * delta_tau) / (1 - tau_b)
+    cost_of_capital = (r + p.delta - tau_b * delta_tau - tau_inv) / (1 - tau_b)
 
     return cost_of_capital
 
diff --git a/ogcore/parameters.py b/ogcore/parameters.py
index 6cfbc182f..1698c94bd 100644
--- a/ogcore/parameters.py
+++ b/ogcore/parameters.py
@@ -130,7 +130,6 @@ def compute_default_params(self):
             )
 
         # Extend parameters that may vary over the time path
-        # those that vary over m: 'Z', 'cit_rate',
         tp_param_list = [
             "alpha_G",
             "alpha_T",
@@ -166,7 +165,13 @@ def compute_default_params(self):
             )
             setattr(self, item, this_attr)
         # Deal with parameters that vary across industry and over time
-        tp_param_list2 = ["Z", "delta_tau_annual", "cit_rate", "tau_c"]
+        tp_param_list2 = [
+            "Z",
+            "delta_tau_annual",
+            "cit_rate",
+            "tau_c",
+            "inv_tax_credit",
+        ]
         for item in tp_param_list2:
             this_attr = getattr(self, item)
             if this_attr.ndim == 1:
diff --git a/ogcore/tax.py b/ogcore/tax.py
index 9fec991a6..fc5c0041b 100644
--- a/ogcore/tax.py
+++ b/ogcore/tax.py
@@ -99,7 +99,8 @@ def MTR_wealth(b, h_wealth, m_wealth, p_wealth):
 
     .. math::
         \frac{\partial T_{j,s,t}^{w}}{\partial b_{j,s,t}} =
-        \frac{h^{w}m^{w}p_{w}}{(b_{j,s,t}h^{w}m^{w})^{2}}
+        \frac{h^{w}p_{w}b_{j,s,t}}{(b_{j,s,t}h^{w}+m^{w})}\left[2 -
+        \frac{h^{w}p_{w}b_{j,s,t}}{(b_{j,s,t}h^{w}+m^{w})}\right]
 
     Args:
         b (Numpy array): savings
@@ -111,9 +112,9 @@ def MTR_wealth(b, h_wealth, m_wealth, p_wealth):
         tau_prime (Numpy array): marginal tax rate on wealth, size = SxJ
 
     """
-    tau_prime = (b * h_wealth * m_wealth * p_wealth) / (
-        (b * h_wealth + m_wealth) ** 2
-    ) + ETR_wealth(b, h_wealth, m_wealth, p_wealth)
+    tau_prime = ETR_wealth(b, h_wealth, m_wealth, p_wealth) * 2 - (
+        (h_wealth**2 * p_wealth * b**2) / ((b * h_wealth + m_wealth) ** 2)
+    )
     return tau_prime
 
 
@@ -204,7 +205,7 @@ def get_biz_tax(w, Y, L, K, p_m, p, m, method):
 
     .. math::
         R_{t}^{b} = \sum_{m=1}^{M}\tau_{m,t}^{b}(Y_{m,t} - w_{t}L_{m,t}) -
-        \tau_{m,t}^{b}\delta_{m,t}^{\tau}K_{m,t}^{\tau}
+        \tau_{m,t}^{b}\delta_{m,t}^{\tau}K_{m,t}^{\tau} - \tau^{inv}_{m,t}I_{m,t}
     Args:
         r (array_like): real interest rate
         Y (array_like): aggregate output for each industry
@@ -222,24 +223,44 @@ def get_biz_tax(w, Y, L, K, p_m, p, m, method):
         if method == "SS":
             delta_tau = p.delta_tau[-1, m]
             tau_b = p.tau_b[-1, m]
+            tau_inv = p.inv_tax_credit[-1, m]
             price = p_m[m]
+            Inv = p.delta * K  # compute gross investment
         else:
             delta_tau = p.delta_tau[: p.T, m].reshape(p.T)
             tau_b = p.tau_b[: p.T, m].reshape(p.T)
+            tau_inv = p.inv_tax_credit[: p.T, m].reshape(p.T)
             price = p_m[: p.T, m].reshape(p.T)
             w = w.reshape(p.T)
+            print("K = ", K.shape)
+            Inv = np.append(
+                K[1:] - K[:-1] + p.delta * K[:-1],
+                np.array([p.delta * K[-1]]),
+                axis=0,
+            )
     else:
         if method == "SS":
             delta_tau = p.delta_tau[-1, :]
             tau_b = p.tau_b[-1, :]
+            tau_inv = p.inv_tax_credit[-1, :]
             price = p_m
+            Inv = p.delta * K
         else:
             delta_tau = p.delta_tau[: p.T, :].reshape(p.T, p.M)
             tau_b = p.tau_b[: p.T, :].reshape(p.T, p.M)
+            tau_inv = p.inv_tax_credit[: p.T, :].reshape(p.T, p.M)
             price = p_m[: p.T, :].reshape(p.T, p.M)
             w = w.reshape(p.T, 1)
+            print("K = ", K.shape)
+            Inv = np.append(
+                K[1:, :] - K[:-1, :] + p.delta * K[:-1, :],
+                p.delta * K[-1, :].reshape(1, p.M),
+                axis=0,
+            )
 
-    business_revenue = tau_b * (price * Y - w * L) - tau_b * delta_tau * K
+    business_revenue = (
+        tau_b * (price * Y - w * L) - tau_b * delta_tau * K - tau_inv * Inv
+    )
     return business_revenue
 
 
diff --git a/tests/test_SS.py b/tests/test_SS.py
index 5e6e6c734..dde5c1133 100644
--- a/tests/test_SS.py
+++ b/tests/test_SS.py
@@ -457,14 +457,14 @@ def test_SS_solver_extra(baseline, param_updates, filename, dask_client):
 
 
 @pytest.mark.parametrize(
-    "baseline,param_updates,filename",
+    "baseline,r_p,param_updates,filename",
     [
-        (True, param_updates1, filename1),
-        (True, param_updates2, filename2),
-        (True, param_updates3, filename3),
-        (False, param_updates4, filename4),
-        (False, param_updates5, filename5),
-        (False, param_updates7, filename7),
+        (True, 0.03309231672773741, param_updates1, filename1),
+        (True, 0.05, param_updates2, filename2),
+        (True, 0.04260341179572245, param_updates3, filename3),
+        (False, 0.04260341179572245, param_updates4, filename4),
+        (False, 0.04260341179572245, param_updates5, filename5),
+        (False, 0.04759112768438152, param_updates7, filename7),
     ],
     ids=[
         "Baseline, Small Open",
@@ -475,7 +475,7 @@ def test_SS_solver_extra(baseline, param_updates, filename, dask_client):
         "Reform, M>1",
     ],
 )
-def test_inner_loop(baseline, param_updates, filename, dask_client):
+def test_inner_loop(baseline, r_p, param_updates, filename, dask_client):
     # Test SS.inner_loop function.  Provide inputs to function and
     # ensure that output returned matches what it has been before.
     p = Specifications(baseline=baseline, num_workers=NUM_WORKERS)
@@ -490,23 +490,7 @@ def test_inner_loop(baseline, param_updates, filename, dask_client):
     w = firm.get_w_from_r(r, p, "SS")
     TR = 0.12
     Y = 1.3
-
-    # Solve for r_p because of new sol'n algo
-    r_gov = fiscal.get_r_gov(r, p)
-    (
-        D,
-        D_d,
-        D_f,
-        new_borrowing,
-        debt_service,
-        new_borrowing_f,
-    ) = fiscal.get_D_ss(r_gov, Y, p)
-    I_g = fiscal.get_I_g(Y, p.alpha_I[-1])
-    K_g = fiscal.get_K_g(0, I_g, p, "SS")
-    MPKg = firm.get_MPx(Y, K_g, p.gamma_g, p, "SS")
-    K = firm.get_K_from_Y(Y, r, p, "SS")
     p_m = np.ones(p.M)
-    r_p = aggregates.get_r_p(r, r_gov, p_m, K, K_g, D, MPKg, p, "SS")
     factor = 100000
     BQ = np.ones(p.J) * 0.00019646295986015257
     if p.budget_balance:
@@ -607,22 +591,8 @@ def test_inner_loop_extra(baseline, param_updates, filename, dask_client):
     Y = 1.3
     factor = 100000
     BQ = np.ones(p.J) * 0.00019646295986015257
-    # Solve for r_p because of new sol'n algo
-    r_gov = fiscal.get_r_gov(r, p)
-    (
-        D,
-        D_d,
-        D_f,
-        new_borrowing,
-        debt_service,
-        new_borrowing_f,
-    ) = fiscal.get_D_ss(r_gov, Y, p)
-    I_g = fiscal.get_I_g(Y, p.alpha_I[-1])
-    K_g = fiscal.get_K_g(0, I_g, p, "SS")
-    MPKg = firm.get_MPx(Y, K_g, p.gamma_g, p, "SS")
-    K = firm.get_K_from_Y(Y, r, p, "SS")
     p_m = np.array([1.0])
-    r_p = aggregates.get_r_p(r, r_gov, p_m, K, K_g, D, MPKg, p, "SS")
+    r_p = 0.04260341179572245
     outer_loop_vars = (bssmat, nssmat, r_p, r, w, p_m, Y, BQ, TR, factor)
     test_tuple = SS.inner_loop(outer_loop_vars, p, dask_client)
     expected_tuple = utils.safe_read_pickle(
diff --git a/tests/test_firm.py b/tests/test_firm.py
index bd2d84af0..4850fcc07 100644
--- a/tests/test_firm.py
+++ b/tests/test_firm.py
@@ -275,13 +275,14 @@ def test_get_Y(K, K_g, L, p, method, m, expected):
     "delta_tau_annual": [[0.35], [0.2], [0.1]],
     "epsilon": [1.2],
     "delta_annual": 0.5,
+    "inv_tax_credit": [[0.07]],
     "T": 3,
     "S": 3,
     "eta": (np.ones((3, p5.J)) / (3 * p5.J)),
 }
 # update parameters instance with new values for test
 p5.update_specifications(new_param_values5)
-expected5 = np.array([-0.07814687, 0.48060411, 0.51451412])
+expected5 = np.array([-0.07814687, 0.48060411, 0.51451412]) + 0.07
 p_m = np.ones((p5.T, p5.M))
 
 
@@ -773,6 +774,7 @@ def test_get_MPx(Y, x, share, p, method, expected):
     "Z": [[2.0]],
     "delta_tau_annual": [[0.35]],
     "delta_annual": 0.25,
+    "inv_tax_credit": [[0.03]],
     "cit_rate": [[0.5]],
     "adjustment_factor_for_cit_receipts": [1.0],
     "c_corp_share_of_assets": 1.0,
@@ -784,9 +786,10 @@ def test_get_MPx(Y, x, share, p, method, expected):
 
 coc_expected1 = np.array([0.75])
 coc_expected2 = np.array([0.75, 0.75, 0.75])
-coc_expected3 = np.array([0.25, 0.25])
-coc_expected4 = np.array([[0.25, 0.25], [0.25, 0.25], [0.25, 0.25]])
-
+coc_expected3 = np.array([0.25, 0.25]) - (0.03 / (1 - 0.5))
+coc_expected4 = np.array([[0.25, 0.25], [0.25, 0.25], [0.25, 0.25]]) - (
+    0.03 / (1 - 0.5)
+)
 ky_expected1 = np.array([0.315478672])
 ky_expected2 = np.array([2.4, 2.4, 2.4])
 ky_expected3 = np.array([2.4])
diff --git a/tests/test_output_tables.py b/tests/test_output_tables.py
index d5cefc1ac..7c056ccfd 100644
--- a/tests/test_output_tables.py
+++ b/tests/test_output_tables.py
@@ -29,7 +29,13 @@
 reform_params = utils.safe_read_pickle(
     os.path.join(CUR_PATH, "test_io_data", "model_params_reform.pkl")
 )
-
+# add investment tax credit parameter that not in cached parameters
+base_params.inv_tax_credit = np.zeros(
+    (base_params.T + base_params.S, base_params.M)
+)
+reform_params.inv_tax_credit = np.zeros(
+    (reform_params.T + reform_params.S, reform_params.M)
+)
 
 test_data = [
     (base_tpi, base_params, reform_tpi, reform_params, "pct_diff"),
diff --git a/tests/test_tax.py b/tests/test_tax.py
index 7a679b8ec..b352e9887 100644
--- a/tests/test_tax.py
+++ b/tests/test_tax.py
@@ -796,6 +796,7 @@ def test_MTR_income(etr_params, mtr_params, params, mtr_capital, expected):
 new_param_values1 = {
     "cit_rate": [[0.20]],
     "delta_tau_annual": [[0.0023176377601205056]],
+    "inv_tax_credit": [[0.02]],
     "T": 3,
     "S": 3,
     "eta": (np.ones((3, p1.J)) / (3 * p1.J)),
@@ -807,7 +808,11 @@ def test_MTR_income(etr_params, mtr_params, params, mtr_capital, expected):
 p_m = np.ones((p1.T, 1))
 L = np.array([[2.0], [3.0], [2.0]])
 K = np.array([[5.0], [6.0], [5.0]])
-expected1 = np.array([[0.0102], [0.11356], [0.0102]])
+inv_tax_credit_amounts = (
+    np.array([[1 + p1.delta * 5.0], [-1 + p1.delta * 6.0], [p1.delta * 5.0]])
+    * 0.02
+)
+expected1 = np.array([[0.0102], [0.11356], [0.0102]]) - inv_tax_credit_amounts
 
 
 @pytest.mark.parametrize(