Repository: econ-ark/HARK

Issue:
"Income" should be replaced with IncShkDstn or something like that
As it stands currently, most of the components of HARK that are keeping track of the income shocks use the tag "Income".  That is confusing to users, who might not realize that the object being referenced is a distribution of income *shocks*, rather than the level of income.  

At the next overhaul, we should do a global search-and-replace.  I'd propose using IncShk instead of Income

Discussion:

Patch:
diff --git a/Documentation/CHANGELOG.md b/Documentation/CHANGELOG.md
index 7e73b40f..d828e39d 100644
--- a/Documentation/CHANGELOG.md
+++ b/Documentation/CHANGELOG.md
@@ -28,10 +28,10 @@ interpolation for problems with CRRA utility. See [#888](https://github.com/econ
 
 * Move AgentType constructor parameters docs to class docstring so it is rendered by Sphinx.
 * Remove uses of deprecated time.clock [#887](https://github.com/econ-ark/HARK/pull/887)
+* Rename IncomeDstn to IncShkDstn
 * AgentType simulate() method now returns history. [#916](https://github.com/econ-ark/HARK/pull/916)
 * Rename DiscreteDistribution.drawDiscrete() to draw()
 
-
 ### 0.10.8
 
 Release Data: Nov. 05 2020
diff --git a/HARK/ConsumptionSaving/ConsAggShockModel.py b/HARK/ConsumptionSaving/ConsAggShockModel.py
index 535d9777..2252255c 100644
--- a/HARK/ConsumptionSaving/ConsAggShockModel.py
+++ b/HARK/ConsumptionSaving/ConsAggShockModel.py
@@ -210,7 +210,7 @@ class AggShockConsumerType(IndShockConsumerType):
 
     def addAggShkDstn(self, AggShkDstn):
         """
-        Updates attribute IncomeDstn by combining idiosyncratic shocks with aggregate shocks.
+        Updates attribute IncShkDstn by combining idiosyncratic shocks with aggregate shocks.
 
         Parameters
         ----------
@@ -223,12 +223,12 @@ class AggShockConsumerType(IndShockConsumerType):
         -------
         None
         """
-        if len(self.IncomeDstn[0].X) > 2:
-            self.IncomeDstn = self.IncomeDstnWithoutAggShocks
+        if len(self.IncShkDstn[0].X) > 2:
+            self.IncShkDstn = self.IncShkDstnWithoutAggShocks
         else:
-            self.IncomeDstnWithoutAggShocks = self.IncomeDstn
-        self.IncomeDstn = [
-            combineIndepDstns(self.IncomeDstn[t], AggShkDstn)
+            self.IncShkDstnWithoutAggShocks = self.IncShkDstn
+        self.IncShkDstn = [
+            combineIndepDstns(self.IncShkDstn[t], AggShkDstn)
             for t in range(self.T_cycle)
         ]
 
@@ -400,7 +400,7 @@ class AggShockConsumerType(IndShockConsumerType):
         Creates a "normalized Euler error" function for this instance, mapping
         from market resources to "consumption error per dollar of consumption."
         Stores result in attribute eulerErrorFunc as an interpolated function.
-        Has option to use approximate income distribution stored in self.IncomeDstn
+        Has option to use approximate income distribution stored in self.IncShkDstn
         or to use a (temporary) very dense approximation.
 
         NOT YET IMPLEMENTED FOR THIS CLASS
@@ -411,9 +411,9 @@ class AggShockConsumerType(IndShockConsumerType):
             Maximum normalized market resources for the Euler error function.
         approx_inc_dstn : Boolean
             Indicator for whether to use the approximate discrete income distri-
-            bution stored in self.IncomeDstn[0], or to use a very accurate
+            bution stored in self.IncShkDstn[0], or to use a very accurate
             discrete approximation instead.  When True, uses approximation in
-            IncomeDstn; when False, makes and uses a very dense approximation.
+            IncShkDstn; when False, makes and uses a very dense approximation.
 
         Returns
         -------
@@ -467,21 +467,21 @@ class AggShockMarkovConsumerType(AggShockConsumerType):
         state. AggShkDstn is a list of aggregate productivity shock distributions
         for each Markov state.
         """
-        if len(self.IncomeDstn[0][0].X) > 2:
-            self.IncomeDstn = self.IncomeDstnWithoutAggShocks
+        if len(self.IncShkDstn[0][0].X) > 2:
+            self.IncShkDstn = self.IncShkDstnWithoutAggShocks
         else:
-            self.IncomeDstnWithoutAggShocks = self.IncomeDstn
+            self.IncShkDstnWithoutAggShocks = self.IncShkDstn
 
-        IncomeDstnOut = []
+        IncShkDstnOut = []
         N = self.MrkvArray.shape[0]
         for t in range(self.T_cycle):
-            IncomeDstnOut.append(
+            IncShkDstnOut.append(
                 [
-                    combineIndepDstns(self.IncomeDstn[t][n], AggShkDstn[n])
+                    combineIndepDstns(self.IncShkDstn[t][n], AggShkDstn[n])
                     for n in range(N)
                 ]
             )
-        self.IncomeDstn = IncomeDstnOut
+        self.IncShkDstn = IncShkDstnOut
 
     def updateSolutionTerminal(self):
         """
@@ -509,7 +509,7 @@ class AggShockMarkovConsumerType(AggShockConsumerType):
 
     def getShocks(self):
         """
-        Gets permanent and transitory income shocks for this period.  Samples from IncomeDstn for
+        Gets permanent and transitory income shocks for this period.  Samples from IncShkDstn for
         each period in the cycle.  This is a copy-paste from IndShockConsumerType, with the
         addition of the Markov macroeconomic state.  Unfortunately, the getShocks method for
         MarkovConsumerType cannot be used, as that method assumes that MrkvNow is a vector
@@ -530,13 +530,13 @@ class AggShockMarkovConsumerType(AggShockConsumerType):
             these = t == self.t_cycle
             N = np.sum(these)
             if N > 0:
-                IncomeDstnNow = self.IncomeDstn[t - 1][
+                IncShkDstnNow = self.IncShkDstn[t - 1][
                     self.MrkvNow
                 ]  # set current income distribution
                 PermGroFacNow = self.PermGroFac[t - 1]  # and permanent growth factor
 
                 # Get random draws of income shocks from the discrete distribution
-                ShockDraws = IncomeDstnNow.draw(N, exact_match=True)
+                ShockDraws = IncShkDstnNow.draw(N, exact_match=True)
                 # Permanent "shock" includes expected growth
                 PermShkNow[these] = ShockDraws[0] * PermGroFacNow
                 TranShkNow[these] = ShockDraws[1]
@@ -546,13 +546,14 @@ class AggShockMarkovConsumerType(AggShockConsumerType):
         N = np.sum(newborn)
         if N > 0:
             these = newborn
-            IncomeDstnNow = self.IncomeDstn[0][
+            IncShkDstnNow = self.IncShkDstn[0][
                 self.MrkvNow
             ]  # set current income distribution
             PermGroFacNow = self.PermGroFac[0]  # and permanent growth factor
 
             # Get random draws of income shocks from the discrete distribution
-            ShockDraws = IncomeDstnNow.draw(N, exact_match=True)
+            ShockDraws = IncShkDstnNow.draw(N, exact_match=True)
+
             # Permanent "shock" includes expected growth
             PermShkNow[these] = ShockDraws[0] * PermGroFacNow
             TranShkNow[these] = ShockDraws[1]
@@ -1032,7 +1033,7 @@ class KrusellSmithType(AgentType):
 
 def solveConsAggShock(
     solution_next,
-    IncomeDstn,
+    IncShkDstn,
     LivPrb,
     DiscFac,
     CRRA,
@@ -1054,7 +1055,7 @@ def solveConsAggShock(
     ----------
     solution_next : ConsumerSolution
         The solution to the succeeding one period problem.
-    IncomeDstn : [np.array]
+    IncShkDstn : [np.array]
         A list containing five arrays of floats, representing a discrete
         approximation to the income process between the period being solved
         and the one immediately following (in solution_next). Order: event
@@ -1100,11 +1101,11 @@ def solveConsAggShock(
     mNrmMinNext = solution_next.mNrmMin
 
     # Unpack the income shocks
-    ShkPrbsNext = IncomeDstn.pmf
-    PermShkValsNext = IncomeDstn.X[0]
-    TranShkValsNext = IncomeDstn.X[1]
-    PermShkAggValsNext = IncomeDstn.X[2]
-    TranShkAggValsNext = IncomeDstn.X[3]
+    ShkPrbsNext = IncShkDstn.pmf
+    PermShkValsNext = IncShkDstn.X[0]
+    TranShkValsNext = IncShkDstn.X[1]
+    PermShkAggValsNext = IncShkDstn.X[2]
+    TranShkAggValsNext = IncShkDstn.X[3]
     ShkCount = ShkPrbsNext.size
 
     # Make the grid of end-of-period asset values, and a tiled version
@@ -1232,7 +1233,7 @@ def solveConsAggShock(
     return solution_now
 
 
-def solveConsAggShockNEW(solution_next, IncomeDstn, LivPrb, DiscFac, CRRA, PermGroFac,
+def solveConsAggShockNEW(solution_next, IncShkDstn, LivPrb, DiscFac, CRRA, PermGroFac,
                       PermGroFacAgg, aXtraGrid, BoroCnstArt, Mgrid, AFunc, Rfunc, wFunc, DeprFac):
     '''
     Solve one period of a consumption-saving problem with idiosyncratic and
@@ -1244,7 +1245,7 @@ def solveConsAggShockNEW(solution_next, IncomeDstn, LivPrb, DiscFac, CRRA, PermG
     ----------
     solution_next : ConsumerSolution
         The solution to the succeeding one period problem.
-    IncomeDstn : [np.array]
+    IncShkDstn : [np.array]
         A list containing five arrays of floats, representing a discrete
         approximation to the income process between the period being solved
         and the one immediately following (in solution_next). Order: event
@@ -1286,10 +1287,10 @@ def solveConsAggShockNEW(solution_next, IncomeDstn, LivPrb, DiscFac, CRRA, PermG
         tions) and marginal value function vPfunc.
     '''
     # Unpack the income shocks and get grid sizes
-    PermShkValsNext = IncomeDstn.X[0]
-    TranShkValsNext = IncomeDstn.X[1]
-    PermShkAggValsNext = IncomeDstn.X[2]
-    TranShkAggValsNext = IncomeDstn.X[3]
+    PermShkValsNext = IncShkDstn.X[0]
+    TranShkValsNext = IncShkDstn.X[1]
+    PermShkAggValsNext = IncShkDstn.X[2]
+    TranShkAggValsNext = IncShkDstn.X[3]
     aCount = aXtraGrid.size
     Mcount = Mgrid.size
     
@@ -1326,7 +1327,7 @@ def solveConsAggShockNEW(solution_next, IncomeDstn, LivPrb, DiscFac, CRRA, PermG
     
     # Compute end-of-period marginal value of assets
     MaggNow = np.tile(np.reshape(Mgrid,(1,Mcount)),(aCount,1)) # Tiled Mgrid
-    EndOfPrdvP = DiscFac*LivPrb*calcExpectation(IncomeDstn,vPnextFunc,[aNrmNow,MaggNow])
+    EndOfPrdvP = DiscFac*LivPrb*calcExpectation(IncShkDstn,vPnextFunc,[aNrmNow,MaggNow])
 
     # Calculate optimal consumption from each asset gridpoint and endogenous m_t gridpoint
     cNrmNow = EndOfPrdvP**(-1.0/CRRA)
@@ -1365,7 +1366,7 @@ def solveConsAggShockNEW(solution_next, IncomeDstn, LivPrb, DiscFac, CRRA, PermG
 
 def solveConsAggMarkov(
     solution_next,
-    IncomeDstn,
+    IncShkDstn,
     LivPrb,
     DiscFac,
     CRRA,
@@ -1390,7 +1391,7 @@ def solveConsAggMarkov(
     ----------
     solution_next : ConsumerSolution
         The solution to the succeeding one period problem.
-    IncomeDstn : [[np.array]]
+    IncShkDstn : [[np.array]]
         A list of lists, each containing five arrays of floats, representing a
         discrete approximation to the income process between the period being
         solved and the one immediately following (in solution_next). Order: event
@@ -1452,11 +1453,11 @@ def solveConsAggMarkov(
         mNrmMinNext = solution_next.mNrmMin[j]
 
         # Unpack the income shocks
-        ShkPrbsNext = IncomeDstn[j].pmf
-        PermShkValsNext = IncomeDstn[j].X[0]
-        TranShkValsNext = IncomeDstn[j].X[1]
-        PermShkAggValsNext = IncomeDstn[j].X[2]
-        TranShkAggValsNext = IncomeDstn[j].X[3]
+        ShkPrbsNext = IncShkDstn[j].pmf
+        PermShkValsNext = IncShkDstn[j].X[0]
+        TranShkValsNext = IncShkDstn[j].X[1]
+        PermShkAggValsNext = IncShkDstn[j].X[2]
+        TranShkAggValsNext = IncShkDstn[j].X[3]
         ShkCount = ShkPrbsNext.size
         aXtra_tiled = np.tile(
             np.reshape(aXtraGrid, (1, aCount, 1)), (Mcount, 1, ShkCount)
diff --git a/HARK/ConsumptionSaving/ConsGenIncProcessModel.py b/HARK/ConsumptionSaving/ConsGenIncProcessModel.py
index e10f49ad..42f64fa1 100644
--- a/HARK/ConsumptionSaving/ConsGenIncProcessModel.py
+++ b/HARK/ConsumptionSaving/ConsGenIncProcessModel.py
@@ -119,7 +119,7 @@ class ConsGenIncProcessSolver(ConsIndShockSetup):
     ----------
     solution_next : ConsumerSolution
         The solution to next period's one period problem.
-    IncomeDstn : [np.array]
+    IncShkDstn : [np.array]
         A list containing three arrays of floats, representing a discrete
         approximation to the income process between the period being solved
         and the one immediately following (in solution_next). Order: event
@@ -153,7 +153,7 @@ class ConsGenIncProcessSolver(ConsIndShockSetup):
     def __init__(
         self,
         solution_next,
-        IncomeDstn,
+        IncShkDstn,
         LivPrb,
         DiscFac,
         CRRA,
@@ -167,7 +167,7 @@ class ConsGenIncProcessSolver(ConsIndShockSetup):
     ):
         self.assignParameters(
             solution_next=solution_next,
-            IncomeDstn=IncomeDstn,
+            IncShkDstn=IncShkDstn,
             LivPrb=LivPrb,
             DiscFac=DiscFac,
             CRRA=CRRA,
@@ -182,7 +182,7 @@ class ConsGenIncProcessSolver(ConsIndShockSetup):
         )  # dummy 0.0 variable why PermGroFac?
         self.defUtilityFuncs()
 
-    def setAndUpdateValues(self, solution_next, IncomeDstn, LivPrb, DiscFac):
+    def setAndUpdateValues(self, solution_next, IncShkDstn, LivPrb, DiscFac):
         """
         Unpacks some of the inputs (and calculates simple objects based on them),
         storing the results in self for use by other methods.  These include:
@@ -195,7 +195,7 @@ class ConsGenIncProcessSolver(ConsIndShockSetup):
         ----------
         solution_next : ConsumerSolution
             The solution to next period's one period problem.
-        IncomeDstn : [np.array]
+        IncShkDstn : [np.array]
             A list containing three arrays of floats, representing a discrete
             approximation to the income process between the period being solved
             and the one immediately following (in solution_next). Order: event
@@ -212,7 +212,7 @@ class ConsGenIncProcessSolver(ConsIndShockSetup):
         """
         # Run basic version of this method
         ConsIndShockSetup.setAndUpdateValues(
-            self, solution_next, IncomeDstn, LivPrb, DiscFac
+            self, solution_next, IncShkDstn, LivPrb, DiscFac
         )
         self.mLvlMinNext = solution_next.mLvlMin
 
diff --git a/HARK/ConsumptionSaving/ConsIndShockModel.py b/HARK/ConsumptionSaving/ConsIndShockModel.py
index df1abb23..4bb60b85 100644
--- a/HARK/ConsumptionSaving/ConsIndShockModel.py
+++ b/HARK/ConsumptionSaving/ConsIndShockModel.py
@@ -468,7 +468,7 @@ class ConsIndShockSetup(ConsPerfForesightSolver):
     ----------
     solution_next : ConsumerSolution
         The solution to next period's one period problem.
-    IncomeDstn : [np.array]
+    IncShkDstn : [np.array]
         A list containing three arrays of floats, representing a discrete
         approximation to the income process between the period being solved
         and the one immediately following (in solution_next). Order: event
@@ -503,7 +503,7 @@ class ConsIndShockSetup(ConsPerfForesightSolver):
     def __init__(
         self,
         solution_next,
-        IncomeDstn,
+        IncShkDstn,
         LivPrb,
         DiscFac,
         CRRA,
@@ -516,7 +516,7 @@ class ConsIndShockSetup(ConsPerfForesightSolver):
     ):
         self.assignParameters(
             solution_next=solution_next,
-            IncomeDstn=IncomeDstn,
+            IncShkDstn=IncShkDstn,
             LivPrb=LivPrb,
             DiscFac=DiscFac,
             CRRA=CRRA,
@@ -550,7 +550,7 @@ class ConsIndShockSetup(ConsPerfForesightSolver):
         if self.vFuncBool:
             self.uinv = lambda u: utility_inv(u, gam=self.CRRA)
 
-    def setAndUpdateValues(self, solution_next, IncomeDstn, LivPrb, DiscFac):
+    def setAndUpdateValues(self, solution_next, IncShkDstn, LivPrb, DiscFac):
         """
         Unpacks some of the inputs (and calculates simple objects based on them),
         storing the results in self for use by other methods.  These include:
@@ -562,7 +562,7 @@ class ConsIndShockSetup(ConsPerfForesightSolver):
         ----------
         solution_next : ConsumerSolution
             The solution to next period's one period problem.
-        IncomeDstn : distribution.DiscreteDistribution
+        IncShkDstn : distribution.DiscreteDistribution
             A DiscreteDistribution with a pmf
             and two point value arrays in X, order:
             permanent shocks, transitory shocks.
@@ -577,10 +577,10 @@ class ConsIndShockSetup(ConsPerfForesightSolver):
         None
         """
         self.DiscFacEff = DiscFac * LivPrb  # "effective" discount factor
-        self.IncomeDstn = IncomeDstn
-        self.ShkPrbsNext = IncomeDstn.pmf
-        self.PermShkValsNext = IncomeDstn.X[0]
-        self.TranShkValsNext = IncomeDstn.X[1]
+        self.IncShkDstn = IncShkDstn
+        self.ShkPrbsNext = IncShkDstn.pmf
+        self.PermShkValsNext = IncShkDstn.X[0]
+        self.TranShkValsNext = IncShkDstn.X[1]
         self.PermShkMinNext = np.min(self.PermShkValsNext)
         self.TranShkMinNext = np.min(self.TranShkValsNext)
         self.vPfuncNext = solution_next.vPfunc
@@ -672,7 +672,7 @@ class ConsIndShockSetup(ConsPerfForesightSolver):
         none
         """
         self.setAndUpdateValues(
-            self.solution_next, self.IncomeDstn, self.LivPrb, self.DiscFac
+            self.solution_next, self.IncShkDstn, self.LivPrb, self.DiscFac
         )
         self.defBoroCnst(self.BoroCnstArt)
 
@@ -763,7 +763,7 @@ class ConsIndShockSolverBasic(ConsIndShockSetup):
             * self.Rfree
             * self.PermGroFac ** (-self.CRRA)
             * calcExpectation(
-                self.IncomeDstn,
+                self.IncShkDstn,
                 vp_next,
                 self.aNrmNow
             )
@@ -964,7 +964,7 @@ class ConsIndShockSolver(ConsIndShockSolverBasic):
             * self.Rfree
             * self.PermGroFac ** (-self.CRRA - 1.0)
             * calcExpectation(
-                self.IncomeDstn,
+                self.IncShkDstn,
                 vpp_next,
                 self.aNrmNow
             )
@@ -999,7 +999,7 @@ class ConsIndShockSolver(ConsIndShockSolverBasic):
             * self.PermGroFac ** (1.0 - self.CRRA)
             ) * self.vFuncNext(self.m_nrm_next(shocks, a_nrm))
         EndOfPrdv = self.DiscFacEff * calcExpectation(
-            self.IncomeDstn, v_lvl_next, self.aNrmNow
+            self.IncShkDstn, v_lvl_next, self.aNrmNow
         )
         EndOfPrdvNvrs = self.uinv(
             EndOfPrdv
@@ -1157,7 +1157,7 @@ class ConsKinkedRsolver(ConsIndShockSolver):
     ----------
     solution_next : ConsumerSolution
         The solution to next period's one period problem.
-    IncomeDstn : [np.array]
+    IncShkDstn : [np.array]
         A list containing three arrays of floats, representing a discrete
         approximation to the income process between the period being solved
         and the one immediately following (in solution_next). Order: event
@@ -1196,7 +1196,7 @@ class ConsKinkedRsolver(ConsIndShockSolver):
     def __init__(
         self,
         solution_next,
-        IncomeDstn,
+        IncShkDstn,
         LivPrb,
         DiscFac,
         CRRA,
@@ -1217,7 +1217,7 @@ class ConsKinkedRsolver(ConsIndShockSolver):
         ConsIndShockSolver.__init__(
             self,
             solution_next,
-            IncomeDstn,
+            IncShkDstn,
             LivPrb,
             DiscFac,
             CRRA,
@@ -1903,14 +1903,14 @@ class IndShockConsumerType(PerfForesightConsumerType):
         none
         """
         (
-            IncomeDstn,
+            IncShkDstn,
             PermShkDstn,
             TranShkDstn,
         ) = self.constructLognormalIncomeProcessUnemployment()
-        self.IncomeDstn = IncomeDstn
+        self.IncShkDstn = IncShkDstn
         self.PermShkDstn = PermShkDstn
         self.TranShkDstn = TranShkDstn
-        self.addToTimeVary("IncomeDstn", "PermShkDstn", "TranShkDstn")
+        self.addToTimeVary("IncShkDstn", "PermShkDstn", "TranShkDstn")
 
     def updateAssetsGrid(self):
         """
@@ -1950,7 +1950,7 @@ class IndShockConsumerType(PerfForesightConsumerType):
         Reset the RNG behavior of this type.  This method is called automatically
         by initializeSim(), ensuring that each simulation run uses the same sequence
         of random shocks; this is necessary for structural estimation to work.
-        This method extends AgentType.resetRNG() to also reset elements of IncomeDstn.
+        This method extends AgentType.resetRNG() to also reset elements of IncShkDstn.
 
         Parameters
         ----------
@@ -1962,14 +1962,14 @@ class IndShockConsumerType(PerfForesightConsumerType):
         """
         PerfForesightConsumerType.resetRNG(self)
 
-        # Reset IncomeDstn if it exists (it might not because resetRNG is called at init)
-        if hasattr(self, "IncomeDstn"):
-            for dstn in self.IncomeDstn:
+        # Reset IncShkDstn if it exists (it might not because resetRNG is called at init)
+        if hasattr(self, "IncShkDstn"):
+            for dstn in self.IncShkDstn:
                 dstn.reset()
 
     def getShocks(self):
         """
-        Gets permanent and transitory income shocks for this period.  Samples from IncomeDstn for
+        Gets permanent and transitory income shocks for this period.  Samples from IncShkDstn for
         each period in the cycle.
 
         Parameters
@@ -1987,12 +1987,13 @@ class IndShockConsumerType(PerfForesightConsumerType):
             these = t == self.t_cycle
             N = np.sum(these)
             if N > 0:
-                IncomeDstnNow = self.IncomeDstn[
+                IncShkDstnNow = self.IncShkDstn[
                     t - 1
                 ]  # set current income distribution
                 PermGroFacNow = self.PermGroFac[t - 1]  # and permanent growth factor
                 # Get random draws of income shocks from the discrete distribution
-                IncShks = IncomeDstnNow.draw(N)
+                IncShks = IncShkDstnNow.draw(N)
+
                 PermShkNow[these] = (
                     IncShks[0, :] * PermGroFacNow
                 )  # permanent "shock" includes expected growth
@@ -2003,15 +2004,15 @@ class IndShockConsumerType(PerfForesightConsumerType):
         N = np.sum(newborn)
         if N > 0:
             these = newborn
-            IncomeDstnNow = self.IncomeDstn[0]  # set current income distribution
+            IncShkDstnNow = self.IncShkDstn[0]  # set current income distribution
             PermGroFacNow = self.PermGroFac[0]  # and permanent growth factor
 
             # Get random draws of income shocks from the discrete distribution
-            EventDraws = IncomeDstnNow.draw_events(N)
+            EventDraws = IncShkDstnNow.draw_events(N)
             PermShkNow[these] = (
-                IncomeDstnNow.X[0][EventDraws] * PermGroFacNow
+                IncShkDstnNow.X[0][EventDraws] * PermGroFacNow
             )  # permanent "shock" includes expected growth
-            TranShkNow[these] = IncomeDstnNow.X[1][EventDraws]
+            TranShkNow[these] = IncShkDstnNow.X[1][EventDraws]
         #        PermShkNow[newborn] = 1.0
         TranShkNow[newborn] = 1.0
 
@@ -2040,9 +2041,9 @@ class IndShockConsumerType(PerfForesightConsumerType):
         None
         """
         # Unpack the income distribution and get average and worst outcomes
-        PermShkValsNext = self.IncomeDstn[0][1]
-        TranShkValsNext = self.IncomeDstn[0][2]
-        ShkPrbsNext = self.IncomeDstn[0][0]
+        PermShkValsNext = self.IncShkDstn[0][1]
+        TranShkValsNext = self.IncShkDstn[0][2]
+        ShkPrbsNext = self.IncShkDstn[0][0]
         ExIncNext = np.dot(ShkPrbsNext, PermShkValsNext * TranShkValsNext)
         PermShkMinNext = np.min(PermShkValsNext)
         TranShkMinNext = np.min(TranShkValsNext)
@@ -2077,7 +2078,7 @@ class IndShockConsumerType(PerfForesightConsumerType):
         Creates a "normalized Euler error" function for this instance, mapping
         from market resources to "consumption error per dollar of consumption."
         Stores result in attribute eulerErrorFunc as an interpolated function.
-        Has option to use approximate income distribution stored in self.IncomeDstn
+        Has option to use approximate income distribution stored in self.IncShkDstn
         or to use a (temporary) very dense approximation.
 
         Only works on (one period) infinite horizon models at this time, will
@@ -2089,9 +2090,9 @@ class IndShockConsumerType(PerfForesightConsumerType):
             Maximum normalized market resources for the Euler error function.
         approx_inc_dstn : Boolean
             Indicator for whether to use the approximate discrete income distri-
-            bution stored in self.IncomeDstn[0], or to use a very accurate
+            bution stored in self.IncShkDstn[0], or to use a very accurate
             discrete approximation instead.  When True, uses approximation in
-            IncomeDstn; when False, makes and uses a very dense approximation.
+            IncShkDstn; when False, makes and uses a very dense approximation.
 
         Returns
         -------
@@ -2104,7 +2105,7 @@ class IndShockConsumerType(PerfForesightConsumerType):
         """
         # Get the income distribution (or make a very dense one)
         if approx_inc_dstn:
-            IncomeDstn = self.IncomeDstn[0]
+            IncShkDstn = self.IncShkDstn[0]
         else:
             TranShkDstn = MeanOneLogNormal(sigma=self.TranShkStd[0]).approx(
                 N=200, tail_N=50, tail_order=1.3, tail_bound=[0.05, 0.95]
@@ -2115,7 +2116,7 @@ class IndShockConsumerType(PerfForesightConsumerType):
             PermShkDstn = MeanOneLogNormal(sigma=self.PermShkStd[0]).approx(
                 N=200, tail_N=50, tail_order=1.3, tail_bound=[0.05, 0.95]
             )
-            IncomeDstn = combineIndepDstns(PermShkDstn, TranShkDstn)
+            IncShkDstn = combineIndepDstns(PermShkDstn, TranShkDstn)
 
         # Make a grid of market resources
         mNowMin = self.solution[0].mNrmMin + 10 ** (
@@ -2134,12 +2135,12 @@ class IndShockConsumerType(PerfForesightConsumerType):
         aNowGrid = mNowGrid - cNowGrid
 
         # Tile the grids for fast computation
-        ShkCount = IncomeDstn[0].size
+        ShkCount = IncShkDstn[0].size
         aCount = aNowGrid.size
         aNowGrid_tiled = np.tile(aNowGrid, (ShkCount, 1))
-        PermShkVals_tiled = (np.tile(IncomeDstn[1], (aCount, 1))).transpose()
-        TranShkVals_tiled = (np.tile(IncomeDstn[2], (aCount, 1))).transpose()
-        ShkPrbs_tiled = (np.tile(IncomeDstn[0], (aCount, 1))).transpose()
+        PermShkVals_tiled = (np.tile(IncShkDstn[1], (aCount, 1))).transpose()
+        TranShkVals_tiled = (np.tile(IncShkDstn[2], (aCount, 1))).transpose()
+        ShkPrbs_tiled = (np.tile(IncShkDstn[0], (aCount, 1))).transpose()
 
         # Calculate marginal value next period for each gridpoint and each shock
         mNextArray = (
@@ -2471,7 +2472,7 @@ class IndShockConsumerType(PerfForesightConsumerType):
 
         Returns
         -------
-        IncomeDstn :  [[np.array]]
+        IncShkDstn :  [[np.array]]
             A list with T_cycle elements, each of which is a list of three arrays
             representing a discrete approximation to the income process in a period.
             Order: probabilities, permanent shocks, transitory shocks.
@@ -2494,7 +2495,7 @@ class IndShockConsumerType(PerfForesightConsumerType):
         UnempPrbRet = self.UnempPrbRet
         IncUnempRet = self.IncUnempRet
 
-        IncomeDstn = []  # Discrete approximations to income process in each period
+        IncShkDstn = []  # Discrete approximations to income process in each period
         PermShkDstn = []  # Discrete approximations to permanent income shocks
         TranShkDstn = []  # Discrete approximations to transitory income shocks
 
@@ -2516,18 +2517,18 @@ class IndShockConsumerType(PerfForesightConsumerType):
                 PermShkValsRet = np.array([1.0])
                 TranShkValsRet = np.array([1.0])
                 ShkPrbsRet = np.array([1.0])
-            IncomeDstnRet = DiscreteDistribution(
+            IncShkDstnRet = DiscreteDistribution(
                 ShkPrbsRet,
                 [PermShkValsRet, TranShkValsRet],
                 seed=self.RNG.randint(0, 2 ** 31 - 1),
             )
 
-        # Loop to fill in the list of IncomeDstn random variables.
+        # Loop to fill in the list of IncShkDstn random variables.
         for t in range(T_cycle):  # Iterate over all periods, counting forward
 
             if T_retire > 0 and t >= T_retire:
                 # Then we are in the "retirement period" and add a retirement income object.
-                IncomeDstn.append(deepcopy(IncomeDstnRet))
+                IncShkDstn.append(deepcopy(IncShkDstnRet))
                 PermShkDstn.append([np.array([1.0]), np.array([1.0])])
                 TranShkDstn.append([ShkPrbsRet, TranShkValsRet])
             else:
@@ -2544,7 +2545,7 @@ class IndShockConsumerType(PerfForesightConsumerType):
                 )
                 ### REPLACE
                 ###REPLACE
-                IncomeDstn.append(
+                IncShkDstn.append(
                     combineIndepDstns(
                         PermShkDstn_t,
                         TranShkDstn_t,
@@ -2553,7 +2554,7 @@ class IndShockConsumerType(PerfForesightConsumerType):
                 )  # mix the independent distributions
                 PermShkDstn.append(PermShkDstn_t)
                 TranShkDstn.append(TranShkDstn_t)
-        return IncomeDstn, PermShkDstn, TranShkDstn
+        return IncShkDstn, PermShkDstn, TranShkDstn
 
 
 # Make a dictionary to specify a "kinked R" idiosyncratic shock consumer
@@ -2621,11 +2622,11 @@ class KinkedRconsumerType(IndShockConsumerType):
         None
         """
         # Unpack the income distribution and get average and worst outcomes
-        PermShkValsNext = self.IncomeDstn[0][1]
-        TranShkValsNext = self.IncomeDstn[0][2]
-        ShkPrbsNext = self.IncomeDstn[0][0]
+        PermShkValsNext = self.IncShkDstn[0][1]
+        TranShkValsNext = self.IncShkDstn[0][2]
+        ShkPrbsNext = self.IncShkDstn[0][0]
         ExIncNext = calcExpectation(
-            IncomeDstn,
+            IncShkDstn,
             lambda trans, perm : trans * perm
         )
         PermShkMinNext = np.min(PermShkValsNext)
@@ -2664,7 +2665,7 @@ class KinkedRconsumerType(IndShockConsumerType):
         Creates a "normalized Euler error" function for this instance, mapping
         from market resources to "consumption error per dollar of consumption."
         Stores result in attribute eulerErrorFunc as an interpolated function.
-        Has option to use approximate income distribution stored in self.IncomeDstn
+        Has option to use approximate income distribution stored in self.IncShkDstn
         or to use a (temporary) very dense approximation.
 
         SHOULD BE INHERITED FROM ConsIndShockModel
@@ -2675,9 +2676,9 @@ class KinkedRconsumerType(IndShockConsumerType):
             Maximum normalized market resources for the Euler error function.
         approx_inc_dstn : Boolean
             Indicator for whether to use the approximate discrete income distri-
-            bution stored in self.IncomeDstn[0], or to use a very accurate
+            bution stored in self.IncShkDstn[0], or to use a very accurate
             discrete approximation instead.  When True, uses approximation in
-            IncomeDstn; when False, makes and uses a very dense approximation.
+            IncShkDstn; when False, makes and uses a very dense approximation.
 
         Returns
         -------
@@ -2728,7 +2729,7 @@ class KinkedRconsumerType(IndShockConsumerType):
 
 
 def applyFlatIncomeTax(
-    IncomeDstn, tax_rate, T_retire, unemployed_indices=None, transitory_index=2
+    IncShkDstn, tax_rate, T_retire, unemployed_indices=None, transitory_index=2
 ):
     """
     Applies a flat income tax rate to all employed income states during the working
@@ -2736,7 +2737,7 @@ def applyFlatIncomeTax(
 
     Parameters
     ----------
-    IncomeDstn : [income distributions]
+    IncShkDstn : [income distributions]
         The discrete approximation to the income distribution in each time period.
     tax_rate : float
         A flat income tax rate to be applied to all employed income.
@@ -2745,24 +2746,24 @@ def applyFlatIncomeTax(
     unemployed_indices : [int]
         Indices of transitory shocks that represent unemployment states (no tax).
     transitory_index : int
-        The index of each element of IncomeDstn representing transitory shocks.
+        The index of each element of IncShkDstn representing transitory shocks.
 
     Returns
     -------
-    IncomeDstn_new : [income distributions]
+    IncShkDstn_new : [income distributions]
         The updated income distributions, after applying the tax.
     """
     unemployed_indices = (
         unemployed_indices if unemployed_indices is not None else list()
     )
-    IncomeDstn_new = deepcopy(IncomeDstn)
+    IncShkDstn_new = deepcopy(IncShkDstn)
     i = transitory_index
-    for t in range(len(IncomeDstn)):
+    for t in range(len(IncShkDstn)):
         if t < T_retire:
-            for j in range((IncomeDstn[t][i]).size):
+            for j in range((IncShkDstn[t][i]).size):
                 if j not in unemployed_indices:
-                    IncomeDstn_new[t][i][j] = IncomeDstn[t][i][j] * (1 - tax_rate)
-    return IncomeDstn_new
+                    IncShkDstn_new[t][i][j] = IncShkDstn[t][i][j] * (1 - tax_rate)
+    return IncShkDstn_new
 
 
 # =======================================================
diff --git a/HARK/ConsumptionSaving/ConsIndShockModelFast.py b/HARK/ConsumptionSaving/ConsIndShockModelFast.py
index 2fab53b7..a7a4d6e5 100644
--- a/HARK/ConsumptionSaving/ConsIndShockModelFast.py
+++ b/HARK/ConsumptionSaving/ConsIndShockModelFast.py
@@ -585,9 +585,9 @@ class ConsIndShockSolverBasicFast(ConsIndShockSolverBasic):
         none
         """
 
-        self.ShkPrbsNext = self.IncomeDstn.pmf
-        self.PermShkValsNext = self.IncomeDstn.X[0]
-        self.TranShkValsNext = self.IncomeDstn.X[1]
+        self.ShkPrbsNext = self.IncShkDstn.pmf
+        self.PermShkValsNext = self.IncShkDstn.X[0]
+        self.TranShkValsNext = self.IncShkDstn.X[1]
 
         (
             self.DiscFacEff,
diff --git a/HARK/ConsumptionSaving/ConsLaborModel.py b/HARK/ConsumptionSaving/ConsLaborModel.py
index ef676771..6804e0a7 100644
--- a/HARK/ConsumptionSaving/ConsLaborModel.py
+++ b/HARK/ConsumptionSaving/ConsLaborModel.py
@@ -437,7 +437,7 @@ class LaborIntMargConsumerType(IndShockConsumerType):
         Creates a "normalized Euler error" function for this instance, mapping
         from market resources to "consumption error per dollar of consumption."
         Stores result in attribute eulerErrorFunc as an interpolated function.
-        Has option to use approximate income distribution stored in self.IncomeDstn
+        Has option to use approximate income distribution stored in self.IncShkDstn
         or to use a (temporary) very dense approximation.
 
         NOT YET IMPLEMENTED FOR THIS CLASS
@@ -448,9 +448,9 @@ class LaborIntMargConsumerType(IndShockConsumerType):
             Maximum normalized market resources for the Euler error function.
         approx_inc_dstn : Boolean
             Indicator for whether to use the approximate discrete income distri-
-            bution stored in self.IncomeDstn[0], or to use a very accurate
+            bution stored in self.IncShkDstn[0], or to use a very accurate
             discrete approximation instead.  When True, uses approximation in
-            IncomeDstn; when False, makes and uses a very dense approximation.
+            IncShkDstn; when False, makes and uses a very dense approximation.
 
         Returns
         -------
diff --git a/HARK/ConsumptionSaving/ConsMarkovModel.py b/HARK/ConsumptionSaving/ConsMarkovModel.py
index 3d91399d..69479469 100644
--- a/HARK/ConsumptionSaving/ConsMarkovModel.py
+++ b/HARK/ConsumptionSaving/ConsMarkovModel.py
@@ -57,14 +57,14 @@ class ConsMarkovSolver(ConsIndShockSolver):
     Extends ConsIndShockSolver, with identical inputs but for a discrete
     Markov state, whose transition rule is summarized in MrkvArray.  Markov
     states can differ in their interest factor, permanent growth factor, live probability, and
-    income distribution, so the inputs Rfree, PermGroFac, IncomeDstn, and LivPrb are
+    income distribution, so the inputs Rfree, PermGroFac, IncShkDstn, and LivPrb are
     now arrays or lists specifying those values in each (succeeding) Markov state.
     """
 
     def __init__(
         self,
         solution_next,
-        IncomeDstn_list,
+        IncShkDstn_list,
         LivPrb,
         DiscFac,
         CRRA,
@@ -85,7 +85,7 @@ class ConsMarkovSolver(ConsIndShockSolver):
         ----------
         solution_next : ConsumerSolution
             The solution to next period's one period problem.
-        IncomeDstn_list : [[np.array]]
+        IncShkDstn_list : [[np.array]]
             A length N list of income distributions in each succeeding Markov
             state.  Each income distribution contains three arrays of floats,
             representing a discrete approximation to the income process at the
@@ -131,7 +131,7 @@ class ConsMarkovSolver(ConsIndShockSolver):
 
         self.assignParameters(
             solution_next=solution_next,
-            IncomeDstn_list=IncomeDstn_list,
+            IncShkDstn_list=IncShkDstn_list,
             LivPrb=LivPrb,
             DiscFac=DiscFac,
             CRRA=CRRA,
@@ -241,8 +241,8 @@ class ConsMarkovSolver(ConsIndShockSolver):
         self.BoroCnstNatAll = np.zeros(self.StateCount) + np.nan
         # Find the natural borrowing constraint conditional on next period's state
         for j in range(self.StateCount):
-            PermShkMinNext = np.min(self.IncomeDstn_list[j].X[0])
-            TranShkMinNext = np.min(self.IncomeDstn_list[j].X[1])
+            PermShkMinNext = np.min(self.IncShkDstn_list[j].X[0])
+            TranShkMinNext = np.min(self.IncShkDstn_list[j].X[1])
             self.BoroCnstNatAll[j] = (
                 (self.solution_next.mNrmMin[j] - TranShkMinNext)
                 * (self.PermGroFac_list[j] * PermShkMinNext)
@@ -289,14 +289,14 @@ class ConsMarkovSolver(ConsIndShockSolver):
         none
         """
         # Set future-state-conditional values as attributes of self
-        self.IncomeDstn = self.IncomeDstn_list[state_index]
+        self.IncShkDstn = self.IncShkDstn_list[state_index]
         self.Rfree = self.Rfree_list[state_index]
         self.PermGroFac = self.PermGroFac_list[state_index]
         self.vPfuncNext = self.solution_next.vPfunc[state_index]
         self.mNrmMinNow = self.mNrmMin_list[state_index]
         self.BoroCnstNat = self.BoroCnstNatAll[state_index]
         self.setAndUpdateValues(
-            self.solution_next, self.IncomeDstn, self.LivPrb, self.DiscFac
+            self.solution_next, self.IncShkDstn, self.LivPrb, self.DiscFac
         )
         self.DiscFacEff = (
             self.DiscFac
@@ -335,7 +335,7 @@ class ConsMarkovSolver(ConsIndShockSolver):
             * self.Rfree
             * self.PermGroFac ** (-self.CRRA - 1.0)
             * calcExpectation(
-                self.IncomeDstn,
+                self.IncShkDstn,
                 vpp_next,
                 self.aNrmNow
             )
@@ -736,7 +736,7 @@ class ConsMarkovSolver(ConsIndShockSolver):
 
 def _solveConsMarkov(
     solution_next,
-    IncomeDstn,
+    IncShkDstn,
     LivPrb,
     DiscFac,
     CRRA,
@@ -754,14 +754,14 @@ def _solveConsMarkov(
     identical inputs as solveConsIndShock, except for a discrete
     Markov transitionrule MrkvArray.  Markov states can differ in their interest
     factor, permanent growth factor, and income distribution, so the inputs Rfree,
-    PermGroFac, and IncomeDstn are arrays or lists specifying those values in each
+    PermGroFac, and IncShkDstn are arrays or lists specifying those values in each
     (succeeding) Markov state.
 
     Parameters
     ----------
     solution_next : ConsumerSolution
         The solution to next period's one period problem.
-    IncomeDstn_list : [[np.array]]
+    IncShkDstn_list : [[np.array]]
         A length N list of income distributions in each succeeding Markov
         state.  Each income distribution contains three arrays of floats,
         representing a discrete approximation to the income process at the
@@ -814,7 +814,7 @@ def _solveConsMarkov(
     """
     solver = ConsMarkovSolver(
         solution_next,
-        IncomeDstn,
+        IncShkDstn,
         LivPrb,
         DiscFac,
         CRRA,
@@ -899,16 +899,16 @@ class MarkovConsumerType(IndShockConsumerType):
                 )
 
         # Now check the income distribution.
-        # Note IncomeDstn is (potentially) time-varying, so it is in time_vary.
+        # Note IncShkDstn is (potentially) time-varying, so it is in time_vary.
         # Therefore it is a list, and each element of that list responds to the income distribution
         # at a particular point in time.  Each income distribution at a point in time should itself
         # be a list, with each element corresponding to the income distribution
         # conditional on a particular Markov state.
         # TODO: should this be a numpy array too?
-        for IncomeDstn_t in self.IncomeDstn:
-            if not isinstance(IncomeDstn_t, list) or len(IncomeDstn_t) != StateCount:
+        for IncShkDstn_t in self.IncShkDstn:
+            if not isinstance(IncShkDstn_t, list) or len(IncShkDstn_t) != StateCount:
                 raise ValueError(
-                    "List in IncomeDstn is not the right length, it should be length equal to number of states"
+                    "List in IncShkDstn is not the right length, it should be length equal to number of states"
                 )
 
     def preSolve(self):
@@ -982,11 +982,11 @@ class MarkovConsumerType(IndShockConsumerType):
         """
         PerfForesightConsumerType.resetRNG(self)
 
-        # Reset IncomeDstn if it exists (it might not because resetRNG is called at init)
-        if hasattr(self, "IncomeDstn"):
-            T = len(self.IncomeDstn)
+        # Reset IncShkDstn if it exists (it might not because resetRNG is called at init)
+        if hasattr(self, "IncShkDstn"):
+            T = len(self.IncShkDstn)
             for t in range(T):
-                for dstn in self.IncomeDstn[t]:
+                for dstn in self.IncShkDstn[t]:
                     dstn.reset()
 
     def simDeath(self):
@@ -1098,7 +1098,7 @@ class MarkovConsumerType(IndShockConsumerType):
     def getShocks(self):
         """
         Gets new Markov states and permanent and transitory income shocks for this period.  Samples
-        from IncomeDstn for each period-state in the cycle.
+        from IncShkDstn for each period-state in the cycle.
 
         Parameters
         ----------
@@ -1119,7 +1119,7 @@ class MarkovConsumerType(IndShockConsumerType):
                 these = np.logical_and(t == self.t_cycle, j == MrkvNow)
                 N = np.sum(these)
                 if N > 0:
-                    IncomeDstnNow = self.IncomeDstn[t - 1][
+                    IncShkDstnNow = self.IncShkDstn[t - 1][
                         j
                     ]  # set current income distribution
                     PermGroFacNow = self.PermGroFac[t - 1][
@@ -1127,11 +1127,11 @@ class MarkovConsumerType(IndShockConsumerType):
                     ]  # and permanent growth factor
 
                     # Get random draws of income shocks from the discrete distribution
-                    EventDraws = IncomeDstnNow.draw_events(N)
+                    EventDraws = IncShkDstnNow.draw_events(N)
                     PermShkNow[these] = (
-                        IncomeDstnNow.X[0][EventDraws] * PermGroFacNow
+                        IncShkDstnNow.X[0][EventDraws] * PermGroFacNow
                     )  # permanent "shock" includes expected growth
-                    TranShkNow[these] = IncomeDstnNow.X[1][EventDraws]
+                    TranShkNow[these] = IncShkDstnNow.X[1][EventDraws]
         newborn = self.t_age == 0
         PermShkNow[newborn] = 1.0
         TranShkNow[newborn] = 1.0
@@ -1226,7 +1226,7 @@ class MarkovConsumerType(IndShockConsumerType):
         Creates a "normalized Euler error" function for this instance, mapping
         from market resources to "consumption error per dollar of consumption."
         Stores result in attribute eulerErrorFunc as an interpolated function.
-        Has option to use approximate income distribution stored in self.IncomeDstn
+        Has option to use approximate income distribution stored in self.IncShkDstn
         or to use a (temporary) very dense approximation.
 
         NOT YET IMPLEMENTED FOR THIS CLASS
@@ -1237,9 +1237,9 @@ class MarkovConsumerType(IndShockConsumerType):
             Maximum normalized market resources for the Euler error function.
         approx_inc_dstn : Boolean
             Indicator for whether to use the approximate discrete income distri-
-            bution stored in self.IncomeDstn[0], or to use a very accurate
+            bution stored in self.IncShkDstn[0], or to use a very accurate
             discrete approximation instead.  When True, uses approximation in
-            IncomeDstn; when False, makes and uses a very dense approximation.
+            IncShkDstn; when False, makes and uses a very dense approximation.
 
         Returns
         -------
diff --git a/HARK/ConsumptionSaving/ConsMedModel.py b/HARK/ConsumptionSaving/ConsMedModel.py
index d77f8152..d8b83800 100644
--- a/HARK/ConsumptionSaving/ConsMedModel.py
+++ b/HARK/ConsumptionSaving/ConsMedModel.py
@@ -885,7 +885,7 @@ class ConsMedShockSolver(ConsGenIncProcessSolver):
     ----------
     solution_next : ConsumerSolution
         The solution to next period's one period problem.
-    IncomeDstn : [np.array]
+    IncShkDstn : [np.array]
         A list containing three arrays of floats, representing a discrete
         approximation to the income process between the period being solved
         and the one immediately following (in solution_next). Order: event
@@ -927,7 +927,7 @@ class ConsMedShockSolver(ConsGenIncProcessSolver):
     def __init__(
         self,
         solution_next,
-        IncomeDstn,
+        IncShkDstn,
         MedShkDstn,
         LivPrb,
         DiscFac,
@@ -944,7 +944,7 @@ class ConsMedShockSolver(ConsGenIncProcessSolver):
     ):
         self.assignParameters(
             solution_next=solution_next,
-            IncomeDstn=IncomeDstn,
+            IncShkDstn=IncShkDstn,
             MedShkDstn=MedShkDstn,
             LivPrb=LivPrb,
             DiscFac=DiscFac,
@@ -962,7 +962,7 @@ class ConsMedShockSolver(ConsGenIncProcessSolver):
         )  # dummy value required?
         self.defUtilityFuncs()
 
-    def setAndUpdateValues(self, solution_next, IncomeDstn, LivPrb, DiscFac):
+    def setAndUpdateValues(self, solution_next, IncShkDstn, LivPrb, DiscFac):
         """
         Unpacks some of the inputs (and calculates simple objects based on them),
         storing the results in self for use by other methods.  These include:
@@ -975,7 +975,7 @@ class ConsMedShockSolver(ConsGenIncProcessSolver):
         ----------
         solution_next : ConsumerSolution
             The solution to next period's one period problem.
-        IncomeDstn : [np.array]
+        IncShkDstn : [np.array]
             A list containing three arrays of floats, representing a discrete
             approximation to the income process between the period being solved
             and the one immediately following (in solution_next). Order: event
@@ -992,7 +992,7 @@ class ConsMedShockSolver(ConsGenIncProcessSolver):
         """
         # Run basic version of this method
         ConsGenIncProcessSolver.setAndUpdateValues(
-            self, self.solution_next, self.IncomeDstn, self.LivPrb, self.DiscFac
+            self, self.solution_next, self.IncShkDstn, self.LivPrb, self.DiscFac
         )
 
         # Also unpack the medical shock distribution
diff --git a/HARK/ConsumptionSaving/ConsPortfolioModel.py b/HARK/ConsumptionSaving/ConsPortfolioModel.py
index 09c1cbfb..e68b3922 100644
--- a/HARK/ConsumptionSaving/ConsPortfolioModel.py
+++ b/HARK/ConsumptionSaving/ConsPortfolioModel.py
@@ -272,12 +272,12 @@ class PortfolioConsumerType(IndShockConsumerType):
         """
         if "RiskyDstn" in self.time_vary:
             self.ShockDstn = [
-                combineIndepDstns(self.IncomeDstn[t], self.RiskyDstn[t])
+                combineIndepDstns(self.IncShkDstn[t], self.RiskyDstn[t])
                 for t in range(self.T_cycle)
             ]
         else:
             self.ShockDstn = [
-                combineIndepDstns(self.IncomeDstn[t], self.RiskyDstn)
+                combineIndepDstns(self.IncShkDstn[t], self.RiskyDstn)
                 for t in range(self.T_cycle)
             ]
         self.addToTimeVary("ShockDstn")
@@ -508,7 +508,7 @@ class PortfolioConsumerType(IndShockConsumerType):
 def solveConsPortfolio(
     solution_next,
     ShockDstn,
-    IncomeDstn,
+    IncShkDstn,
     RiskyDstn,
     LivPrb,
     DiscFac,
@@ -536,7 +536,7 @@ def solveConsPortfolio(
         transitory income shocks, and risky returns.  This is only used if the
         input IndepDstnBool is False, indicating that income and return distributions
         can't be assumed to be independent.
-    IncomeDstn : [np.array]
+    IncShkDstn : [np.array]
         List with three arrays: discrete probabilities, permanent income shocks,
         and transitory income shocks.  This is only used if the input IndepDsntBool
         is True, indicating that income and return distributions are independent.
@@ -616,9 +616,9 @@ def solveConsPortfolio(
     # Major method fork: (in)dependent risky asset return and income distributions
     if IndepDstnBool:  # If the distributions ARE independent...
         # Unpack the shock distribution
-        IncPrbs_next = IncomeDstn.pmf
-        PermShks_next = IncomeDstn.X[0]
-        TranShks_next = IncomeDstn.X[1]
+        IncPrbs_next = IncShkDstn.pmf
+        PermShks_next = IncShkDstn.X[0]
+        TranShks_next = IncShkDstn.X[1]
         Rprbs_next = RiskyDstn.pmf
         Risky_next = RiskyDstn.X
         zero_bound = (
@@ -646,7 +646,7 @@ def solveConsPortfolio(
         Income_N = IncPrbs_next.size
         Risky_N = Rprbs_next.size
 
-        # Make tiled arrays to calculate future realizations of mNrm and Share when integrating over IncomeDstn
+        # Make tiled arrays to calculate future realizations of mNrm and Share when integrating over IncShkDstn
         bNrm_tiled = np.tile(
             np.reshape(bNrmGrid, (bNrm_N, 1, 1)), (1, Share_N, Income_N)
         )
diff --git a/HARK/ConsumptionSaving/ConsPrefShockModel.py b/HARK/ConsumptionSaving/ConsPrefShockModel.py
index 964d443b..f5ead3d1 100644
--- a/HARK/ConsumptionSaving/ConsPrefShockModel.py
+++ b/HARK/ConsumptionSaving/ConsPrefShockModel.py
@@ -229,7 +229,7 @@ class PrefShockConsumerType(IndShockConsumerType):
         Creates a "normalized Euler error" function for this instance, mapping
         from market resources to "consumption error per dollar of consumption."
         Stores result in attribute eulerErrorFunc as an interpolated function.
-        Has option to use approximate income distribution stored in self.IncomeDstn
+        Has option to use approximate income distribution stored in self.IncShkDstn
         or to use a (temporary) very dense approximation.
 
         NOT YET IMPLEMENTED FOR THIS CLASS
@@ -240,9 +240,9 @@ class PrefShockConsumerType(IndShockConsumerType):
             Maximum normalized market resources for the Euler error function.
         approx_inc_dstn : Boolean
             Indicator for whether to use the approximate discrete income distri-
-            bution stored in self.IncomeDstn[0], or to use a very accurate
+            bution stored in self.IncShkDstn[0], or to use a very accurate
             discrete approximation instead.  When True, uses approximation in
-            IncomeDstn; when False, makes and uses a very dense approximation.
+            IncShkDstn; when False, makes and uses a very dense approximation.
 
         Returns
         -------
@@ -301,7 +301,7 @@ class ConsPrefShockSolver(ConsIndShockSolver):
     ----------
     solution_next : ConsumerSolution
         The solution to the succeeding one period problem.
-    IncomeDstn : [np.array]
+    IncShkDstn : [np.array]
         A list containing three arrays of floats, representing a discrete
         approximation to the income process between the period being solved
         and the one immediately following (in solution_next). Order: event
@@ -339,7 +339,7 @@ class ConsPrefShockSolver(ConsIndShockSolver):
     def __init__(
         self,
         solution_next,
-        IncomeDstn,
+        IncShkDstn,
         PrefShkDstn,
         LivPrb,
         DiscFac,
@@ -363,7 +363,7 @@ class ConsPrefShockSolver(ConsIndShockSolver):
         ConsIndShockSolver.__init__(
             self,
             solution_next,
-            IncomeDstn,
+            IncShkDstn,
             LivPrb,
             DiscFac,
             CRRA,
@@ -522,7 +522,7 @@ class ConsPrefShockSolver(ConsIndShockSolver):
 
 def solveConsPrefShock(
     solution_next,
-    IncomeDstn,
+    IncShkDstn,
     PrefShkDstn,
     LivPrb,
     DiscFac,
@@ -542,7 +542,7 @@ def solveConsPrefShock(
     ----------
     solution_next : ConsumerSolution
         The solution to the succeeding one period problem.
-    IncomeDstn : [np.array]
+    IncShkDstn : [np.array]
         A list containing three arrays of floats, representing a discrete
         approximation to the income process between the period being solved
         and the one immediately following (in solution_next). Order: event
@@ -590,7 +590,7 @@ def solveConsPrefShock(
     """
     solver = (
         solution_next,
-        IncomeDstn,
+        IncShkDstn,
         PrefShkDstn,
         LivPrb,
         DiscFac,
@@ -620,7 +620,7 @@ class ConsKinkyPrefSolver(ConsPrefShockSolver, ConsKinkedRsolver):
     ----------
     solution_next : ConsumerSolution
         The solution to the succeeding one period problem.
-    IncomeDstn : [np.array]
+    IncShkDstn : [np.array]
         A list containing three arrays of floats, representing a discrete
         approximation to the income process between the period being solved
         and the one immediately following (in solution_next). Order: event
@@ -662,7 +662,7 @@ class ConsKinkyPrefSolver(ConsPrefShockSolver, ConsKinkedRsolver):
     def __init__(
         self,
         solution_next,
-        IncomeDstn,
+        IncShkDstn,
         PrefShkDstn,
         LivPrb,
         DiscFac,
@@ -678,7 +678,7 @@ class ConsKinkyPrefSolver(ConsPrefShockSolver, ConsKinkedRsolver):
         ConsKinkedRsolver.__init__(
             self,
             solution_next,
-            IncomeDstn,
+            IncShkDstn,
             LivPrb,
             DiscFac,
             CRRA,
diff --git a/HARK/ConsumptionSaving/ConsRepAgentModel.py b/HARK/ConsumptionSaving/ConsRepAgentModel.py
index ba491dbc..b2c8730a 100644
--- a/HARK/ConsumptionSaving/ConsRepAgentModel.py
+++ b/HARK/ConsumptionSaving/ConsRepAgentModel.py
@@ -22,7 +22,7 @@ __all__ = ["RepAgentConsumerType", "RepAgentMarkovConsumerType"]
 
 
 def solveConsRepAgent(
-    solution_next, DiscFac, CRRA, IncomeDstn, CapShare, DeprFac, PermGroFac, aXtraGrid
+    solution_next, DiscFac, CRRA, IncShkDstn, CapShare, DeprFac, PermGroFac, aXtraGrid
 ):
     """
     Solve one period of the simple representative agent consumption-saving model.
@@ -35,7 +35,7 @@ def solveConsRepAgent(
         Intertemporal discount factor for future utility.
     CRRA : float
         Coefficient of relative risk aversion.
-    IncomeDstn : [np.array]
+    IncShkDstn : [np.array]
         A list containing three arrays of floats, representing a discrete
         approximation to the income process between the period being solved
         and the one immediately following (in solution_next). Order: event
@@ -58,9 +58,9 @@ def solveConsRepAgent(
     """
     # Unpack next period's solution and the income distribution
     vPfuncNext = solution_next.vPfunc
-    ShkPrbsNext = IncomeDstn.pmf
-    PermShkValsNext = IncomeDstn.X[0]
-    TranShkValsNext = IncomeDstn.X[1]
+    ShkPrbsNext = IncShkDstn.pmf
+    PermShkValsNext = IncShkDstn.X[0]
+    TranShkValsNext = IncShkDstn.X[1]
 
     # Make tiled versions of end-of-period assets, shocks, and probabilities
     aNrmNow = aXtraGrid
@@ -115,7 +115,7 @@ def solveConsRepAgentMarkov(
     MrkvArray,
     DiscFac,
     CRRA,
-    IncomeDstn,
+    IncShkDstn,
     CapShare,
     DeprFac,
     PermGroFac,
@@ -135,7 +135,7 @@ def solveConsRepAgentMarkov(
         Intertemporal discount factor for future utility.
     CRRA : float
         Coefficient of relative risk aversion.
-    IncomeDstn : [[np.array]]
+    IncShkDstn : [[np.array]]
         A list of lists containing three arrays of floats, representing a discrete
         approximation to the income process between the period being solved
         and the one immediately following (in solution_next). Order: event
@@ -167,9 +167,9 @@ def solveConsRepAgentMarkov(
     for j in range(StateCount):
         # Define next-period-state conditional objects
         vPfuncNext = solution_next.vPfunc[j]
-        ShkPrbsNext = IncomeDstn[j].pmf
-        PermShkValsNext = IncomeDstn[j].X[0]
-        TranShkValsNext = IncomeDstn[j].X[1]
+        ShkPrbsNext = IncShkDstn[j].pmf
+        PermShkValsNext = IncShkDstn[j].X[0]
+        TranShkValsNext = IncShkDstn[j].X[1]
 
         # Make tiled versions of end-of-period assets, shocks, and probabilities
         ShkCount = ShkPrbsNext.size
@@ -355,14 +355,14 @@ class RepAgentMarkovConsumerType(RepAgentConsumerType):
 
         t = self.t_cycle[0]
         i = self.MrkvNow[0]
-        IncomeDstnNow = self.IncomeDstn[t - 1][i]  # set current income distribution
+        IncShkDstnNow = self.IncShkDstn[t - 1][i]  # set current income distribution
         PermGroFacNow = self.PermGroFac[t - 1][i]  # and permanent growth factor
         # Get random draws of income shocks from the discrete distribution
-        EventDraw = IncomeDstnNow.draw_events(1)
+        EventDraw = IncShkDstnNow.draw_events(1)
         PermShkNow = (
-            IncomeDstnNow.X[0][EventDraw] * PermGroFacNow
+            IncShkDstnNow.X[0][EventDraw] * PermGroFacNow
         )  # permanent "shock" includes expected growth
-        TranShkNow = IncomeDstnNow.X[1][EventDraw]
+        TranShkNow = IncShkDstnNow.X[1][EventDraw]
         self.shocks["PermShkNow"] = np.array(PermShkNow)
         self.shocks["TranShkNow"] = np.array(TranShkNow)
 


Test patch:
diff --git a/HARK/ConsumptionSaving/tests/test_ConsAggShockModel.py b/HARK/ConsumptionSaving/tests/test_ConsAggShockModel.py
index 8f767ae6..a24bfa4f 100644
--- a/HARK/ConsumptionSaving/tests/test_ConsAggShockModel.py
+++ b/HARK/ConsumptionSaving/tests/test_ConsAggShockModel.py
@@ -63,7 +63,7 @@ class testAggShockMarkovConsumerType(unittest.TestCase):
         self.agent = AggShockMarkovConsumerType()
         self.agent.cycles = 0
         self.agent.AgentCount = 1000  # Low number of simulated agents
-        self.agent.IncomeDstn[0] = 2 * [self.agent.IncomeDstn[0]]  ## see #557
+        self.agent.IncShkDstn[0] = 2 * [self.agent.IncShkDstn[0]]  ## see #557
         self.economy = CobbDouglasMarkovEconomy(agents=[self.agent])
 
     def test_agent(self):
diff --git a/HARK/ConsumptionSaving/tests/test_ConsMarkovModel.py b/HARK/ConsumptionSaving/tests/test_ConsMarkovModel.py
index 8a46d5cd..7e4af5c9 100644
--- a/HARK/ConsumptionSaving/tests/test_ConsMarkovModel.py
+++ b/HARK/ConsumptionSaving/tests/test_ConsMarkovModel.py
@@ -65,7 +65,7 @@ class test_ConsMarkovSolver(unittest.TestCase):
         unemployed_income_dist = DiscreteDistribution(
             np.ones(1), [np.ones(1), np.zeros(1)]
         )  # Definitely don't
-        self.model.IncomeDstn = [
+        self.model.IncShkDstn = [
             [
                 employed_income_dist,
                 unemployed_income_dist,
diff --git a/HARK/ConsumptionSaving/tests/test_ConsRepAgentModel.py b/HARK/ConsumptionSaving/tests/test_ConsRepAgentModel.py
index 15345788..a1c42c97 100644
--- a/HARK/ConsumptionSaving/tests/test_ConsRepAgentModel.py
+++ b/HARK/ConsumptionSaving/tests/test_ConsRepAgentModel.py
@@ -27,7 +27,7 @@ class testRepAgentConsumerType(unittest.TestCase):
 class testRepAgentMarkovConsumerType(unittest.TestCase):
     def setUp(self):
         self.agent = RepAgentMarkovConsumerType()
-        self.agent.IncomeDstn[0] = 2 * [self.agent.IncomeDstn[0]]
+        self.agent.IncShkDstn[0] = 2 * [self.agent.IncShkDstn[0]]
         self.agent.solve()
 
     def test_solution(self):
diff --git a/HARK/ConsumptionSaving/tests/test_IndShockConsumerType.py b/HARK/ConsumptionSaving/tests/test_IndShockConsumerType.py
index 05bbe56e..f920e2f0 100644
--- a/HARK/ConsumptionSaving/tests/test_IndShockConsumerType.py
+++ b/HARK/ConsumptionSaving/tests/test_IndShockConsumerType.py
@@ -51,7 +51,7 @@ class testIndShockConsumerType(unittest.TestCase):
 
         solver = ConsIndShockSolverBasic(
             LifecycleExample.solution[1],
-            LifecycleExample.IncomeDstn[0],
+            LifecycleExample.IncShkDstn[0],
             LifecycleExample.LivPrb[0],
             LifecycleExample.DiscFac,
             LifecycleExample.CRRA,
diff --git a/HARK/ConsumptionSaving/tests/test_modelcomparisons.py b/HARK/ConsumptionSaving/tests/test_modelcomparisons.py
index 2879d5a2..34b4b226 100644
--- a/HARK/ConsumptionSaving/tests/test_modelcomparisons.py
+++ b/HARK/ConsumptionSaving/tests/test_modelcomparisons.py
@@ -161,7 +161,7 @@ class Compare_TBS_and_Markov(unittest.TestCase):
         unemployed_income_dist = DiscreteDistribution(
             np.ones(1), [np.ones(1), np.zeros(1)]
         )
-        MarkovType.IncomeDstn = [[employed_income_dist, unemployed_income_dist]]
+        MarkovType.IncShkDstn = [[employed_income_dist, unemployed_income_dist]]
 
         MarkovType.solve()
         MarkovType.unpack("cFunc")
diff --git a/HARK/tests/test_distribution.py b/HARK/tests/test_distribution.py
index 05ec6bfe..bc38da8d 100644
--- a/HARK/tests/test_distribution.py
+++ b/HARK/tests/test_distribution.py
@@ -77,17 +77,17 @@ class DiscreteDistributionTests(unittest.TestCase):
 
         PermShkDstn = MeanOneLogNormal().approx(200)
         TranShkDstn = MeanOneLogNormal().approx(200)
-        IncomeDstn = combineIndepDstns(PermShkDstn, TranShkDstn)
+        IncShkDstn = combineIndepDstns(PermShkDstn, TranShkDstn)
 
         ce8 = calcExpectation(
-            IncomeDstn,
+            IncShkDstn,
             lambda X: X[0] + X[1]
         )
 
         self.assertAlmostEqual(ce8, 2.0)
 
         ce9 = calcExpectation(
-            IncomeDstn,
+            IncShkDstn,
             lambda X, a, r: r / X[0] * a + X[1],
             np.array([0,1,2,3,4,5]), # an aNrmNow grid?
             1.05 # an interest rate?