Repository: econ-ark/HARK

Issue:
Allow discretization of continuous distributions
There are now classes for continuous distributions (Lognormal, etc.), a class for a generic discrete distribution, and functions for approximating discrete versions of continuous distributions.

The approximation functions should be turned into discetization functions on the continuous distribution classes.

Discussion:

Patch:
diff --git a/.gitignore b/.gitignore
index f7720865..8d6f08d5 100644
--- a/.gitignore
+++ b/.gitignore
@@ -263,8 +263,9 @@ Web
 # Mac system files
 .DS_Store
 
-# Emacs automatic backup files
+# Emacs automatic backup and hash files
 *~
+\#*#
 
 # Vim swap files
 *.swp
diff --git a/HARK/ConsumptionSaving/ConsAggShockModel.py b/HARK/ConsumptionSaving/ConsAggShockModel.py
index d5b99ea6..5a0ca601 100644
--- a/HARK/ConsumptionSaving/ConsAggShockModel.py
+++ b/HARK/ConsumptionSaving/ConsAggShockModel.py
@@ -10,7 +10,7 @@ from builtins import str
 from builtins import range
 import numpy as np
 import scipy.stats as stats
-from HARK.distribution import DiscreteDistribution, combineIndepDstns, approxMeanOneLognormal
+from HARK.distribution import DiscreteDistribution, combineIndepDstns, MeanOneLogNormal
 from HARK.interpolation import LinearInterp, LinearInterpOnInterp1D, ConstantFunction, IdentityFunction,\
                                VariableLowerBoundFunc2D, BilinearInterp, LowerEnvelope2D, UpperEnvelope
 from HARK.utilities import CRRAutility, CRRAutilityP, CRRAutilityPP, CRRAutilityP_inv,\
@@ -1095,8 +1095,12 @@ class CobbDouglasEconomy(Market):
         -------
         None
         '''
-        self.TranShkAggDstn = approxMeanOneLognormal(sigma=self.TranShkAggStd, N=self.TranShkAggCount)
-        self.PermShkAggDstn = approxMeanOneLognormal(sigma=self.PermShkAggStd, N=self.PermShkAggCount)
+        self.TranShkAggDstn = MeanOneLogNormal(
+            sigma=self.TranShkAggStd
+        ).approx(N=self.TranShkAggCount)
+        self.PermShkAggDstn = MeanOneLognormal(
+            sigma=self.PermShkAggStd
+        ).approx(N=self.PermShkAggCount)
         self.AggShkDstn = combineIndepDstns(self.PermShkAggDstn, self.TranShkAggDstn)
 
     def reset(self):
@@ -1312,8 +1316,12 @@ class SmallOpenEconomy(Market):
         -------
         None
         '''
-        self.TranShkAggDstn = approxMeanOneLognormal(sigma=self.TranShkAggStd, N=self.TranShkAggCount)
-        self.PermShkAggDstn = approxMeanOneLognormal(sigma=self.PermShkAggStd, N=self.PermShkAggCount)
+        self.TranShkAggDstn = MeanOneLogNormal(
+            sigma=self.TranShkAggStd
+        ).approx(N=self.TranShkAggCount)
+        self.PermShkAggDstn = MeanOneLogNormal(
+            sigma=self.PermShkAggStd
+        ).approx(N=self.PermShkAggCount)
         self.AggShkDstn = combineIndepDstns(self.PermShkAggDstn, self.TranShkAggDstn)
 
     def millRule(self):
@@ -1526,8 +1534,8 @@ class CobbDouglasMarkovEconomy(CobbDouglasEconomy):
         StateCount = self.MrkvArray.shape[0]
 
         for i in range(StateCount):
-            TranShkAggDstn.append(approxMeanOneLognormal(sigma=self.TranShkAggStd[i], N=self.TranShkAggCount))
-            PermShkAggDstn.append(approxMeanOneLognormal(sigma=self.PermShkAggStd[i], N=self.PermShkAggCount))
+            TranShkAggDstn.append(MeanOneLogNormal(sigma=self.TranShkAggStd[i]).approx(N=self.TranShkAggCount))
+            PermShkAggDstn.append(MeanOneLogNormal(sigma=self.PermShkAggStd[i]).approx(N=self.PermShkAggCount))
             AggShkDstn.append(combineIndepDstns(PermShkAggDstn[-1], TranShkAggDstn[-1]))
 
         self.TranShkAggDstn = TranShkAggDstn
diff --git a/HARK/ConsumptionSaving/ConsIndShockModel.py b/HARK/ConsumptionSaving/ConsIndShockModel.py
index 8a474f72..7c4e17ca 100644
--- a/HARK/ConsumptionSaving/ConsIndShockModel.py
+++ b/HARK/ConsumptionSaving/ConsIndShockModel.py
@@ -24,8 +24,8 @@ from scipy.optimize import newton
 from HARK import AgentType, Solution, NullFunc, HARKobject
 from HARK.utilities import warnings  # Because of "patch" to warnings modules
 from HARK.interpolation import CubicInterp, LowerEnvelope, LinearInterp
-from HARK.distribution import Lognormal, Uniform
-from HARK.distribution import DiscreteDistribution, approxMeanOneLognormal, addDiscreteOutcomeConstantMean, combineIndepDstns 
+from HARK.distribution import Lognormal, MeanOneLogNormal, Uniform
+from HARK.distribution import DiscreteDistribution, addDiscreteOutcomeConstantMean, combineIndepDstns 
 from HARK.utilities import makeGridExpMult, CRRAutility, CRRAutilityP, \
                            CRRAutilityPP, CRRAutilityP_inv, CRRAutility_invP, CRRAutility_inv, \
                            CRRAutilityP_invP
@@ -2257,11 +2257,19 @@ class IndShockConsumerType(PerfForesightConsumerType):
         if approx_inc_dstn:
             IncomeDstn = self.IncomeDstn[0]
         else:
-            TranShkDstn = approxMeanOneLognormal(N=200,sigma=self.TranShkStd[0],
-                                                 tail_N=50,tail_order=1.3, tail_bound=[0.05,0.95])
+            TranShkDstn = MeanOneLognormal(
+                sigma=self.TranShkStd[0]
+            ).approx(N=200,
+                     tail_N=50,
+                     tail_order=1.3,
+                     tail_bound=[0.05,0.95])
             TranShkDstn = addDiscreteOutcomeConstantMean(TranShkDstn,self.UnempPrb,self.IncUnemp)
-            PermShkDstn = approxMeanOneLognormal(N=200,sigma=self.PermShkStd[0],
-                                                 tail_N=50,tail_order=1.3, tail_bound=[0.05,0.95])
+            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)
 
         # Make a grid of market resources
@@ -2796,10 +2804,14 @@ def constructLognormalIncomeProcessUnemployment(parameters):
             TranShkDstn.append([ShkPrbsRet,TranShkValsRet])
         else:
             # We are in the "working life" periods.
-            TranShkDstn_t    = approxMeanOneLognormal(N=TranShkCount, sigma=TranShkStd[t], tail_N=0)
+            TranShkDstn_t    = MeanOneLogNormal(
+                sigma=TranShkStd[t] 
+            ).approx(TranShkCount, tail_N=0)
             if UnempPrb > 0:
                 TranShkDstn_t = addDiscreteOutcomeConstantMean(TranShkDstn_t, p=UnempPrb, x=IncUnemp)
-            PermShkDstn_t    = approxMeanOneLognormal(N=PermShkCount, sigma=PermShkStd[t], tail_N=0)
+            PermShkDstn_t    = MeanOneLogNormal(
+                sigma=PermShkStd[t]
+            ).approx(PermShkCount, tail_N=0)
             IncomeDstn.append(combineIndepDstns(PermShkDstn_t,TranShkDstn_t)) # mix the independent distributions
             PermShkDstn.append(PermShkDstn_t)
             TranShkDstn.append(TranShkDstn_t)
diff --git a/HARK/ConsumptionSaving/ConsMarkovModel.py.orig b/HARK/ConsumptionSaving/ConsMarkovModel.py.orig
new file mode 100644
index 00000000..f482770d
--- /dev/null
+++ b/HARK/ConsumptionSaving/ConsMarkovModel.py.orig
@@ -0,0 +1,1023 @@
+'''
+Classes to solve and simulate consumption-savings model with a discrete, exogenous,
+stochastic Markov state.  The only solver here extends ConsIndShockModel to
+include a Markov state; the interest factor, permanent growth factor, and income
+distribution can vary with the discrete state.
+'''
+from __future__ import division, print_function
+from __future__ import absolute_import
+from builtins import range
+from copy import deepcopy
+import numpy as np
+from HARK import AgentType
+from HARK.ConsumptionSaving.ConsIndShockModel import ConsIndShockSolver, ValueFunc, \
+                             MargValueFunc, ConsumerSolution, IndShockConsumerType
+<<<<<<< HEAD
+from HARK.distribution import Uniform
+=======
+
+
+from HARK.distribution import DiscreteDistribution
+from HARK.simulation import Uniform
+>>>>>>> 8fc1cd2b34dbb8b29e0c088182bc1d737c7f07ad
+from HARK.interpolation import CubicInterp, LowerEnvelope, LinearInterp
+from HARK.utilities import CRRAutility, CRRAutilityP, CRRAutilityPP, CRRAutilityP_inv, \
+                           CRRAutility_invP, CRRAutility_inv, CRRAutilityP_invP
+
+__all__ = ['ConsMarkovSolver', 'MarkovConsumerType']
+
+utility       = CRRAutility
+utilityP      = CRRAutilityP
+utilityPP     = CRRAutilityPP
+utilityP_inv  = CRRAutilityP_inv
+utility_invP  = CRRAutility_invP
+utility_inv   = CRRAutility_inv
+utilityP_invP = CRRAutilityP_invP
+
+class ConsMarkovSolver(ConsIndShockSolver):
+    '''
+    A class to solve a single period consumption-saving problem with risky income
+    and stochastic transitions between discrete states, in a Markov fashion.
+    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
+    now arrays or lists specifying those values in each (succeeding) Markov state.
+    '''
+    def __init__(self,solution_next,IncomeDstn_list,LivPrb,DiscFac,
+                      CRRA,Rfree_list,PermGroFac_list,MrkvArray,BoroCnstArt,
+                      aXtraGrid,vFuncBool,CubicBool):
+        '''
+        Constructor for a new solver for a one period problem with risky income
+        and transitions between discrete Markov states.  In the descriptions below,
+        N is the number of discrete states.
+
+        Parameters
+        ----------
+        solution_next : ConsumerSolution
+            The solution to next period's one period problem.
+        IncomeDstn_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
+            beginning of the succeeding period. Order: event probabilities,
+            permanent shocks, transitory shocks.
+        LivPrb : np.array
+            Survival probability; likelihood of being alive at the beginning of
+            the succeeding period for each Markov state.
+        DiscFac : float
+            Intertemporal discount factor for future utility.
+        CRRA : float
+            Coefficient of relative risk aversion.
+        Rfree_list : np.array
+            Risk free interest factor on end-of-period assets for each Markov
+            state in the succeeding period.
+        PermGroFac_list : np.array
+            Expected permanent income growth factor at the end of this period
+            for each Markov state in the succeeding period.
+        MrkvArray : np.array
+            An NxN array representing a Markov transition matrix between discrete
+            states.  The i,j-th element of MrkvArray is the probability of
+            moving from state i in period t to state j in period t+1.
+        BoroCnstArt: float or None
+            Borrowing constraint for the minimum allowable assets to end the
+            period with.  If it is less than the natural borrowing constraint,
+            then it is irrelevant; BoroCnstArt=None indicates no artificial bor-
+            rowing constraint.
+        aXtraGrid: np.array
+            Array of "extra" end-of-period asset values-- assets above the
+            absolute minimum acceptable level.
+        vFuncBool: boolean
+            An indicator for whether the value function should be computed and
+            included in the reported solution.
+        CubicBool: boolean
+            An indicator for whether the solver should use cubic or linear inter-
+            polation.
+
+        Returns
+        -------
+        None
+        '''
+        # Set basic attributes of the problem
+        ConsIndShockSolver.assignParameters(self,solution_next,np.nan,LivPrb,DiscFac,CRRA,np.nan,
+                                            np.nan,BoroCnstArt,aXtraGrid,vFuncBool,CubicBool)
+        self.defUtilityFuncs()
+
+        # Set additional attributes specific to the Markov model
+        self.IncomeDstn_list      = IncomeDstn_list
+        self.Rfree_list           = Rfree_list
+        self.PermGroFac_list      = PermGroFac_list
+        self.MrkvArray            = MrkvArray
+        self.StateCount           = MrkvArray.shape[0]
+
+    def solve(self):
+        '''
+        Solve the one period problem of the consumption-saving model with a Markov state.
+
+        Parameters
+        ----------
+        none
+
+        Returns
+        -------
+        solution : ConsumerSolution
+            The solution to the single period consumption-saving problem. Includes
+            a consumption function cFunc (using cubic or linear splines), a marg-
+            inal value function vPfunc, a minimum acceptable level of normalized
+            market resources mNrmMin, normalized human wealth hNrm, and bounding
+            MPCs MPCmin and MPCmax.  It might also have a value function vFunc
+            and marginal marginal value function vPPfunc.  All of these attributes
+            are lists or arrays, with elements corresponding to the current
+            Markov state.  E.g. solution.cFunc[0] is the consumption function
+            when in the i=0 Markov state this period.
+        '''
+        # Find the natural borrowing constraint in each current state
+        self.defBoundary()
+
+        # Initialize end-of-period (marginal) value functions
+        self.EndOfPrdvFunc_list  = []
+        self.EndOfPrdvPfunc_list = []
+        self.ExIncNextAll        = np.zeros(self.StateCount) + np.nan # expected income conditional on the next state
+        self.WorstIncPrbAll      = np.zeros(self.StateCount) + np.nan # probability of getting the worst income shock in each next period state
+
+        # Loop through each next-period-state and calculate the end-of-period
+        # (marginal) value function
+        for j in range(self.StateCount):
+            # Condition values on next period's state (and record a couple for later use)
+            self.conditionOnState(j)
+            self.ExIncNextAll[j]   = np.dot(self.ShkPrbsNext,self.PermShkValsNext*self.TranShkValsNext)
+            self.WorstIncPrbAll[j] = self.WorstIncPrb
+
+            # Construct the end-of-period marginal value function conditional
+            # on next period's state and add it to the list of value functions
+            EndOfPrdvPfunc_cond = self.makeEndOfPrdvPfuncCond()
+            self.EndOfPrdvPfunc_list.append(EndOfPrdvPfunc_cond)
+
+            # Construct the end-of-period value functional conditional on next
+            # period's state and add it to the list of value functions
+            if self.vFuncBool:
+                EndOfPrdvFunc_cond = self.makeEndOfPrdvFuncCond()
+                self.EndOfPrdvFunc_list.append(EndOfPrdvFunc_cond)
+
+        # EndOfPrdvP_cond is EndOfPrdvP conditional on *next* period's state.
+        # Take expectations to get EndOfPrdvP conditional on *this* period's state.
+        self.calcEndOfPrdvP()
+
+        # Calculate the bounding MPCs and PDV of human wealth for each state
+        self.calcHumWealthAndBoundingMPCs()
+
+        # Find consumption and market resources corresponding to each end-of-period
+        # assets point for each state (and add an additional point at the lower bound)
+        aNrm = np.asarray(self.aXtraGrid)[np.newaxis,:] + np.array(self.BoroCnstNat_list)[:,np.newaxis]
+        self.getPointsForInterpolation(self.EndOfPrdvP,aNrm)
+        cNrm = np.hstack((np.zeros((self.StateCount,1)),self.cNrmNow))
+        mNrm = np.hstack((np.reshape(self.mNrmMin_list,(self.StateCount,1)),self.mNrmNow))
+
+        # Package and return the solution for this period
+        self.BoroCnstNat = self.BoroCnstNat_list
+        solution = self.makeSolution(cNrm,mNrm)
+        return solution
+
+    def defBoundary(self):
+        '''
+        Find the borrowing constraint for each current state and save it as an
+        attribute of self for use by other methods.
+
+        Parameters
+        ----------
+        none
+
+        Returns
+        -------
+        none
+        '''
+        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])
+            self.BoroCnstNatAll[j] = (self.solution_next.mNrmMin[j] - TranShkMinNext)*\
+                                     (self.PermGroFac_list[j]*PermShkMinNext)/self.Rfree_list[j]
+
+        self.BoroCnstNat_list   = np.zeros(self.StateCount) + np.nan
+        self.mNrmMin_list       = np.zeros(self.StateCount) + np.nan
+        self.BoroCnstDependency = np.zeros((self.StateCount,self.StateCount)) + np.nan
+        # The natural borrowing constraint in each current state is the *highest*
+        # among next-state-conditional natural borrowing constraints that could
+        # occur from this current state.
+        for i in range(self.StateCount):
+            possible_next_states         = self.MrkvArray[i,:] > 0
+            self.BoroCnstNat_list[i]     = np.max(self.BoroCnstNatAll[possible_next_states])
+
+            # Explicitly handle the "None" case:
+            if self.BoroCnstArt is None:
+                self.mNrmMin_list[i]         = self.BoroCnstNat_list[i]
+            else:
+                self.mNrmMin_list[i]         = np.max([self.BoroCnstNat_list[i],self.BoroCnstArt])
+            self.BoroCnstDependency[i,:] = self.BoroCnstNat_list[i] == self.BoroCnstNatAll
+        # Also creates a Boolean array indicating whether the natural borrowing
+        # constraint *could* be hit when transitioning from i to j.
+
+    def conditionOnState(self,state_index):
+        '''
+        Temporarily assume that a particular Markov state will occur in the
+        succeeding period, and condition solver attributes on this assumption.
+        Allows the solver to construct the future-state-conditional marginal
+        value function (etc) for that future state.
+
+        Parameters
+        ----------
+        state_index : int
+            Index of the future Markov state to condition on.
+
+        Returns
+        -------
+        none
+        '''
+        # Set future-state-conditional values as attributes of self
+        self.IncomeDstn     = self.IncomeDstn_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.DiscFacEff     = self.DiscFac # survival probability LivPrb represents probability from
+                                           # *current* state, so DiscFacEff is just DiscFac for now
+
+        # These lines have to come after setAndUpdateValues to override the definitions there
+        self.vPfuncNext = self.solution_next.vPfunc[state_index]
+        if self.CubicBool:
+            self.vPPfuncNext= self.solution_next.vPPfunc[state_index]
+        if self.vFuncBool:
+            self.vFuncNext  = self.solution_next.vFunc[state_index]
+
+    def calcEndOfPrdvPP(self):
+        '''
+        Calculates end-of-period marginal marginal value using a pre-defined
+        array of next period market resources in self.mNrmNext.
+
+        Parameters
+        ----------
+        none
+
+        Returns
+        -------
+        EndOfPrdvPP : np.array
+            End-of-period marginal marginal value of assets at each value in
+            the grid of assets.
+        '''
+        EndOfPrdvPP = self.DiscFacEff*self.Rfree*self.Rfree*self.PermGroFac**(-self.CRRA-1.0)*\
+                      np.sum(self.PermShkVals_temp**(-self.CRRA-1.0)*self.vPPfuncNext(self.mNrmNext)
+                      *self.ShkPrbs_temp,axis=0)
+        return EndOfPrdvPP
+
+    def makeEndOfPrdvFuncCond(self):
+        '''
+        Construct the end-of-period value function conditional on next period's
+        state.  NOTE: It might be possible to eliminate this method and replace
+        it with ConsIndShockSolver.makeEndOfPrdvFunc, but the self.X_cond
+        variables must be renamed.
+
+        Parameters
+        ----------
+        none
+
+        Returns
+        -------
+        EndofPrdvFunc_cond : ValueFunc
+            The end-of-period value function conditional on a particular state
+            occuring in the next period.
+        '''
+        VLvlNext               = (self.PermShkVals_temp**(1.0-self.CRRA)*
+                                  self.PermGroFac**(1.0-self.CRRA))*self.vFuncNext(self.mNrmNext)
+        EndOfPrdv_cond         = self.DiscFacEff*np.sum(VLvlNext*self.ShkPrbs_temp,axis=0)
+        EndOfPrdvNvrs_cond     = self.uinv(EndOfPrdv_cond)
+        EndOfPrdvNvrsP_cond    = self.EndOfPrdvP_cond*self.uinvP(EndOfPrdv_cond)
+        EndOfPrdvNvrs_cond     = np.insert(EndOfPrdvNvrs_cond,0,0.0)
+        EndOfPrdvNvrsP_cond    = np.insert(EndOfPrdvNvrsP_cond,0,EndOfPrdvNvrsP_cond[0])
+        aNrm_temp              = np.insert(self.aNrm_cond,0,self.BoroCnstNat)
+        EndOfPrdvNvrsFunc_cond = CubicInterp(aNrm_temp,EndOfPrdvNvrs_cond,EndOfPrdvNvrsP_cond)
+        EndofPrdvFunc_cond     = ValueFunc(EndOfPrdvNvrsFunc_cond,self.CRRA)
+        return EndofPrdvFunc_cond
+
+
+    def calcEndOfPrdvPcond(self):
+        '''
+        Calculate end-of-period marginal value of assets at each point in aNrmNow
+        conditional on a particular state occuring in the next period.
+
+        Parameters
+        ----------
+        None
+
+        Returns
+        -------
+        EndOfPrdvP : np.array
+            A 1D array of end-of-period marginal value of assets.
+        '''
+        EndOfPrdvPcond = ConsIndShockSolver.calcEndOfPrdvP(self)
+        return EndOfPrdvPcond
+
+
+    def makeEndOfPrdvPfuncCond(self):
+        '''
+        Construct the end-of-period marginal value function conditional on next
+        period's state.
+
+        Parameters
+        ----------
+        None
+
+        Returns
+        -------
+        EndofPrdvPfunc_cond : MargValueFunc
+            The end-of-period marginal value function conditional on a particular
+            state occuring in the succeeding period.
+        '''
+        # Get data to construct the end-of-period marginal value function (conditional on next state)
+        self.aNrm_cond      = self.prepareToCalcEndOfPrdvP()
+        self.EndOfPrdvP_cond= self.calcEndOfPrdvPcond()
+        EndOfPrdvPnvrs_cond = self.uPinv(self.EndOfPrdvP_cond) # "decurved" marginal value
+        if self.CubicBool:
+            EndOfPrdvPP_cond = self.calcEndOfPrdvPP()
+            EndOfPrdvPnvrsP_cond = EndOfPrdvPP_cond*self.uPinvP(self.EndOfPrdvP_cond) # "decurved" marginal marginal value
+
+        # Construct the end-of-period marginal value function conditional on the next state.
+        if self.CubicBool:
+            EndOfPrdvPnvrsFunc_cond = CubicInterp(self.aNrm_cond,EndOfPrdvPnvrs_cond,
+                                                  EndOfPrdvPnvrsP_cond,lower_extrap=True)
+        else:
+            EndOfPrdvPnvrsFunc_cond = LinearInterp(self.aNrm_cond,EndOfPrdvPnvrs_cond,
+                                                   lower_extrap=True)
+        EndofPrdvPfunc_cond = MargValueFunc(EndOfPrdvPnvrsFunc_cond,self.CRRA) # "recurve" the interpolated marginal value function
+        return EndofPrdvPfunc_cond
+
+    def calcEndOfPrdvP(self):
+        '''
+        Calculates end of period marginal value (and marginal marginal) value
+        at each aXtra gridpoint for each current state, unconditional on the
+        future Markov state (i.e. weighting conditional end-of-period marginal
+        value by transition probabilities).
+
+        Parameters
+        ----------
+        none
+
+        Returns
+        -------
+        none
+        '''
+        # Find unique values of minimum acceptable end-of-period assets (and the
+        # current period states for which they apply).
+        aNrmMin_unique, state_inverse = np.unique(self.BoroCnstNat_list,return_inverse=True)
+        self.possible_transitions     = self.MrkvArray > 0
+
+        # Calculate end-of-period marginal value (and marg marg value) at each
+        # asset gridpoint for each current period state
+        EndOfPrdvP                    = np.zeros((self.StateCount,self.aXtraGrid.size))
+        EndOfPrdvPP                   = np.zeros((self.StateCount,self.aXtraGrid.size))
+        for k in range(aNrmMin_unique.size):
+            aNrmMin       = aNrmMin_unique[k]   # minimum assets for this pass
+            which_states  = state_inverse == k  # the states for which this minimum applies
+            aGrid         = aNrmMin + self.aXtraGrid # assets grid for this pass
+            EndOfPrdvP_all  = np.zeros((self.StateCount,self.aXtraGrid.size))
+            EndOfPrdvPP_all = np.zeros((self.StateCount,self.aXtraGrid.size))
+            for j in range(self.StateCount):
+                if np.any(np.logical_and(self.possible_transitions[:,j],which_states)): # only consider a future state if one of the relevant states could transition to it
+                    EndOfPrdvP_all[j,:] = self.EndOfPrdvPfunc_list[j](aGrid)
+                    if self.CubicBool: # Add conditional end-of-period (marginal) marginal value to the arrays
+                        EndOfPrdvPP_all[j,:] = self.EndOfPrdvPfunc_list[j].derivative(aGrid)
+            # Weight conditional marginal (marginal) values by transition probs
+            # to get unconditional marginal (marginal) value at each gridpoint.
+            EndOfPrdvP_temp = np.dot(self.MrkvArray,EndOfPrdvP_all)
+            EndOfPrdvP[which_states,:] = EndOfPrdvP_temp[which_states,:] # only take the states for which this asset minimum applies
+            if self.CubicBool:
+                EndOfPrdvPP_temp = np.dot(self.MrkvArray,EndOfPrdvPP_all)
+                EndOfPrdvPP[which_states,:] = EndOfPrdvPP_temp[which_states,:]
+
+        # Store the results as attributes of self, scaling end of period marginal value by survival probability from each current state
+        LivPrb_tiled = np.tile(np.reshape(self.LivPrb,(self.StateCount,1)),(1,self.aXtraGrid.size))
+        self.EndOfPrdvP = LivPrb_tiled*EndOfPrdvP
+        if self.CubicBool:
+            self.EndOfPrdvPP = LivPrb_tiled*EndOfPrdvPP
+
+    def calcHumWealthAndBoundingMPCs(self):
+        '''
+        Calculates human wealth and the maximum and minimum MPC for each current
+        period state, then stores them as attributes of self for use by other methods.
+
+        Parameters
+        ----------
+        none
+
+        Returns
+        -------
+        none
+        '''
+        # Upper bound on MPC at lower m-bound
+        WorstIncPrb_array = self.BoroCnstDependency*np.tile(np.reshape(self.WorstIncPrbAll,
+                            (1,self.StateCount)),(self.StateCount,1))
+        temp_array        = self.MrkvArray*WorstIncPrb_array
+        WorstIncPrbNow    = np.sum(temp_array,axis=1) # Probability of getting the "worst" income shock and transition from each current state
+        ExMPCmaxNext      = (np.dot(temp_array,self.Rfree_list**(1.0-self.CRRA)*
+                            self.solution_next.MPCmax**(-self.CRRA))/WorstIncPrbNow)**\
+                            (-1.0/self.CRRA)
+        DiscFacEff_temp   = self.DiscFac*self.LivPrb
+        self.MPCmaxNow    = 1.0/(1.0 + ((DiscFacEff_temp*WorstIncPrbNow)**
+                            (1.0/self.CRRA))/ExMPCmaxNext)
+        self.MPCmaxEff    = self.MPCmaxNow
+        self.MPCmaxEff[self.BoroCnstNat_list < self.mNrmMin_list] = 1.0
+        # State-conditional PDV of human wealth
+        hNrmPlusIncNext   = self.ExIncNextAll + self.solution_next.hNrm
+        self.hNrmNow      = np.dot(self.MrkvArray,(self.PermGroFac_list/self.Rfree_list)*
+                            hNrmPlusIncNext)
+        # Lower bound on MPC as m gets arbitrarily large
+        temp              = (DiscFacEff_temp*np.dot(self.MrkvArray,self.solution_next.MPCmin**
+                            (-self.CRRA)*self.Rfree_list**(1.0-self.CRRA)))**(1.0/self.CRRA)
+        self.MPCminNow    = 1.0/(1.0 + temp)
+
+    def makeSolution(self,cNrm,mNrm):
+        '''
+        Construct an object representing the solution to this period's problem.
+
+        Parameters
+        ----------
+        cNrm : np.array
+            Array of normalized consumption values for interpolation.  Each row
+            corresponds to a Markov state for this period.
+        mNrm : np.array
+            Array of normalized market resource values for interpolation.  Each
+            row corresponds to a Markov state for this period.
+
+        Returns
+        -------
+        solution : ConsumerSolution
+            The solution to the single period consumption-saving problem. Includes
+            a consumption function cFunc (using cubic or linear splines), a marg-
+            inal value function vPfunc, a minimum acceptable level of normalized
+            market resources mNrmMin, normalized human wealth hNrm, and bounding
+            MPCs MPCmin and MPCmax.  It might also have a value function vFunc
+            and marginal marginal value function vPPfunc.  All of these attributes
+            are lists or arrays, with elements corresponding to the current
+            Markov state.  E.g. solution.cFunc[0] is the consumption function
+            when in the i=0 Markov state this period.
+        '''
+        solution = ConsumerSolution() # An empty solution to which we'll add state-conditional solutions
+        # Calculate the MPC at each market resource gridpoint in each state (if desired)
+        if self.CubicBool:
+            dcda          = self.EndOfPrdvPP/self.uPP(np.array(self.cNrmNow))
+            MPC           = dcda/(dcda+1.0)
+            self.MPC_temp = np.hstack((np.reshape(self.MPCmaxNow,(self.StateCount,1)),MPC))
+            interpfunc    = self.makeCubiccFunc
+        else:
+            interpfunc    = self.makeLinearcFunc
+
+        # Loop through each current period state and add its solution to the overall solution
+        for i in range(self.StateCount):
+            # Set current-period-conditional human wealth and MPC bounds
+            self.hNrmNow_j   = self.hNrmNow[i]
+            self.MPCminNow_j = self.MPCminNow[i]
+            if self.CubicBool:
+                self.MPC_temp_j  = self.MPC_temp[i,:]
+
+            # Construct the consumption function by combining the constrained and unconstrained portions
+            self.cFuncNowCnst = LinearInterp([self.mNrmMin_list[i], self.mNrmMin_list[i]+1.0],
+                                             [0.0,1.0])
+            cFuncNowUnc       = interpfunc(mNrm[i,:],cNrm[i,:])
+            cFuncNow          = LowerEnvelope(cFuncNowUnc,self.cFuncNowCnst)
+
+            # Make the marginal value function and pack up the current-state-conditional solution
+            vPfuncNow     = MargValueFunc(cFuncNow,self.CRRA)
+            solution_cond = ConsumerSolution(cFunc=cFuncNow, vPfunc=vPfuncNow,
+                                             mNrmMin=self.mNrmMinNow)
+            if self.CubicBool: # Add the state-conditional marginal marginal value function (if desired)
+                solution_cond = self.addvPPfunc(solution_cond)
+
+            # Add the current-state-conditional solution to the overall period solution
+            solution.appendSolution(solution_cond)
+
+        # Add the lower bounds of market resources, MPC limits, human resources,
+        # and the value functions to the overall solution
+        solution.mNrmMin = self.mNrmMin_list
+        solution         = self.addMPCandHumanWealth(solution)
+        if self.vFuncBool:
+            vFuncNow = self.makevFunc(solution)
+            solution.vFunc = vFuncNow
+
+        # Return the overall solution to this period
+        return solution
+
+
+    def makeLinearcFunc(self,mNrm,cNrm):
+        '''
+        Make a linear interpolation to represent the (unconstrained) consumption
+        function conditional on the current period state.
+
+        Parameters
+        ----------
+        mNrm : np.array
+            Array of normalized market resource values for interpolation.
+        cNrm : np.array
+            Array of normalized consumption values for interpolation.
+
+        Returns
+        -------
+        cFuncUnc: an instance of HARK.interpolation.LinearInterp
+        '''
+        cFuncUnc = LinearInterp(mNrm,cNrm,self.MPCminNow_j*self.hNrmNow_j,self.MPCminNow_j)
+        return cFuncUnc
+
+
+    def makeCubiccFunc(self,mNrm,cNrm):
+        '''
+        Make a cubic interpolation to represent the (unconstrained) consumption
+        function conditional on the current period state.
+
+        Parameters
+        ----------
+        mNrm : np.array
+            Array of normalized market resource values for interpolation.
+        cNrm : np.array
+            Array of normalized consumption values for interpolation.
+
+        Returns
+        -------
+        cFuncUnc: an instance of HARK.interpolation.CubicInterp
+        '''
+        cFuncUnc = CubicInterp(mNrm,cNrm,self.MPC_temp_j,self.MPCminNow_j*self.hNrmNow_j,
+                               self.MPCminNow_j)
+        return cFuncUnc
+
+    def makevFunc(self,solution):
+        '''
+        Construct the value function for each current state.
+
+        Parameters
+        ----------
+        solution : ConsumerSolution
+            The solution to the single period consumption-saving problem. Must
+            have a consumption function cFunc (using cubic or linear splines) as
+            a list with elements corresponding to the current Markov state.  E.g.
+            solution.cFunc[0] is the consumption function when in the i=0 Markov
+            state this period.
+
+        Returns
+        -------
+        vFuncNow : [ValueFunc]
+            A list of value functions (defined over normalized market resources
+            m) for each current period Markov state.
+        '''
+        vFuncNow = [] # Initialize an empty list of value functions
+        # Loop over each current period state and construct the value function
+        for i in range(self.StateCount):
+            # Make state-conditional grids of market resources and consumption
+            mNrmMin       = self.mNrmMin_list[i]
+            mGrid         = mNrmMin + self.aXtraGrid
+            cGrid         = solution.cFunc[i](mGrid)
+            aGrid         = mGrid - cGrid
+
+            # Calculate end-of-period value at each gridpoint
+            EndOfPrdv_all   = np.zeros((self.StateCount,self.aXtraGrid.size))
+            for j in range(self.StateCount):
+                if self.possible_transitions[i,j]:
+                    EndOfPrdv_all[j,:] = self.EndOfPrdvFunc_list[j](aGrid)
+            EndOfPrdv     = np.dot(self.MrkvArray[i,:],EndOfPrdv_all)
+
+            # Calculate (normalized) value and marginal value at each gridpoint
+            vNrmNow       = self.u(cGrid) + EndOfPrdv
+            vPnow         = self.uP(cGrid)
+
+            # Make a "decurved" value function with the inverse utility function
+            vNvrs        = self.uinv(vNrmNow) # value transformed through inverse utility
+            vNvrsP       = vPnow*self.uinvP(vNrmNow)
+            mNrm_temp    = np.insert(mGrid,0,mNrmMin) # add the lower bound
+            vNvrs        = np.insert(vNvrs,0,0.0)
+            vNvrsP       = np.insert(vNvrsP,0,self.MPCmaxEff[i]**(-self.CRRA/(1.0-self.CRRA)))
+            MPCminNvrs   = self.MPCminNow[i]**(-self.CRRA/(1.0-self.CRRA))
+            vNvrsFunc_i  = CubicInterp(mNrm_temp,vNvrs,vNvrsP,MPCminNvrs*self.hNrmNow[i],MPCminNvrs)
+
+            # "Recurve" the decurved value function and add it to the list
+            vFunc_i     = ValueFunc(vNvrsFunc_i,self.CRRA)
+            vFuncNow.append(vFunc_i)
+        return vFuncNow
+
+
+def _solveConsMarkov(solution_next,IncomeDstn,LivPrb,DiscFac,CRRA,Rfree,PermGroFac,
+                                 MrkvArray,BoroCnstArt,aXtraGrid,vFuncBool,CubicBool):
+    '''
+    Solves a single period consumption-saving problem with risky income and
+    stochastic transitions between discrete states, in a Markov fashion.  Has
+    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
+    (succeeding) Markov state.
+
+    Parameters
+    ----------
+    solution_next : ConsumerSolution
+        The solution to next period's one period problem.
+    IncomeDstn_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
+        beginning of the succeeding period. Order: event probabilities,
+        permanent shocks, transitory shocks.
+    LivPrb : float
+        Survival probability; likelihood of being alive at the beginning of
+        the succeeding period.
+    DiscFac : float
+        Intertemporal discount factor for future utility.
+    CRRA : float
+        Coefficient of relative risk aversion.
+    Rfree_list : np.array
+        Risk free interest factor on end-of-period assets for each Markov
+        state in the succeeding period.
+    PermGroGac_list : float
+        Expected permanent income growth factor at the end of this period
+        for each Markov state in the succeeding period.
+    MrkvArray : numpy.array
+        An NxN array representing a Markov transition matrix between discrete
+        states.  The i,j-th element of MrkvArray is the probability of
+        moving from state i in period t to state j in period t+1.
+    BoroCnstArt: float or None
+        Borrowing constraint for the minimum allowable assets to end the
+        period with.  If it is less than the natural borrowing constraint,
+        then it is irrelevant; BoroCnstArt=None indicates no artificial bor-
+        rowing constraint.
+    aXtraGrid: np.array
+        Array of "extra" end-of-period asset values-- assets above the
+        absolute minimum acceptable level.
+    vFuncBool: boolean
+        An indicator for whether the value function should be computed and
+        included in the reported solution.
+    CubicBool: boolean
+        An indicator for whether the solver should use cubic or linear inter-
+        polation.
+
+    Returns
+    -------
+    solution : ConsumerSolution
+        The solution to the single period consumption-saving problem. Includes
+        a consumption function cFunc (using cubic or linear splines), a marg-
+        inal value function vPfunc, a minimum acceptable level of normalized
+        market resources mNrmMin, normalized human wealth hNrm, and bounding
+        MPCs MPCmin and MPCmax.  It might also have a value function vFunc
+        and marginal marginal value function vPPfunc.  All of these attributes
+        are lists or arrays, with elements corresponding to the current
+        Markov state.  E.g. solution.cFunc[0] is the consumption function
+        when in the i=0 Markov state this period.
+    '''
+    solver = ConsMarkovSolver(solution_next,IncomeDstn,LivPrb,DiscFac,CRRA,Rfree,
+                              PermGroFac,MrkvArray,BoroCnstArt,aXtraGrid,vFuncBool,CubicBool)
+    solution_now = solver.solve()
+    return solution_now
+
+
+
+####################################################################################################
+####################################################################################################
+
+class MarkovConsumerType(IndShockConsumerType):
+    '''
+    An agent in the Markov consumption-saving model.  His problem is defined by a sequence
+    of income distributions, survival probabilities, discount factors, and permanent
+    income growth rates, as well as time invariant values for risk aversion, the
+    interest rate, the grid of end-of-period assets, and how he is borrowing constrained.
+    '''
+    time_vary_ = IndShockConsumerType.time_vary_ + ['MrkvArray']
+    shock_vars_ = IndShockConsumerType.shock_vars_ + ['MrkvNow']
+
+    def __init__(self,cycles=1,time_flow=True,**kwds):
+        IndShockConsumerType.__init__(self,cycles=1,time_flow=True,**kwds)
+        self.solveOnePeriod = _solveConsMarkov
+        self.poststate_vars += ['MrkvNow']
+        if not hasattr(self, 'global_markov'):
+            self.global_markov = False
+
+    def checkMarkovInputs(self):
+        '''
+        Many parameters used by MarkovConsumerType are arrays.  Make sure those arrays are the
+        right shape.
+
+        Parameters
+        ----------
+        None
+
+        Returns
+        -------
+        None
+        '''
+        StateCount = self.MrkvArray[0].shape[0]
+
+        # Check that arrays are the right shape
+        if not isinstance(self.Rfree, np.ndarray) or self.Rfree.shape != (StateCount, ):
+            raise ValueError('Rfree not the right shape, it should an array of Rfree of all the states.')
+
+        # Check that arrays in lists are the right shape
+        for MrkvArray_t in self.MrkvArray:
+            if not isinstance(MrkvArray_t, np.ndarray) or MrkvArray_t.shape != (StateCount, StateCount):
+                raise ValueError('MrkvArray not the right shape, it should be of the size states*statres.')
+        for LivPrb_t in self.LivPrb:
+            if not isinstance(LivPrb_t, np.ndarray) or LivPrb_t.shape != (StateCount, ):
+                raise ValueError('Array in LivPrb is not the right shape, it should be an array of length equal to number of states')
+        for PermGroFac_t in self.PermGroFac:
+            if not isinstance(PermGroFac_t, np.ndarray) or PermGroFac_t.shape != (StateCount, ):
+                raise ValueError('Array in PermGroFac is not the right shape, it should be an array of length equal to number of states')
+
+        # Now check the income distribution.
+        # Note IncomeDstn 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:
+                raise ValueError('List in IncomeDstn is not the right length, it should be length equal to number of states')
+
+    def preSolve(self):
+        """
+        Check to make sure that the inputs that are specific to MarkovConsumerType
+        are of the right shape (if arrays) or length (if lists).
+
+        Parameters
+        ----------
+        None
+
+        Returns
+        -------
+        None
+        """
+        AgentType.preSolve(self)
+        self.checkMarkovInputs()
+
+    def updateSolutionTerminal(self):
+        '''
+        Update the terminal period solution.  This method should be run when a
+        new AgentType is created or when CRRA changes.
+
+        Parameters
+        ----------
+        none
+
+        Returns
+        -------
+        none
+        '''
+        IndShockConsumerType.updateSolutionTerminal(self)
+
+        # Make replicated terminal period solution: consume all resources, no human wealth, minimum m is 0
+        StateCount = self.MrkvArray[0].shape[0]
+        self.solution_terminal.cFunc   = StateCount*[self.cFunc_terminal_]
+        self.solution_terminal.vFunc   = StateCount*[self.solution_terminal.vFunc]
+        self.solution_terminal.vPfunc  = StateCount*[self.solution_terminal.vPfunc]
+        self.solution_terminal.vPPfunc = StateCount*[self.solution_terminal.vPPfunc]
+        self.solution_terminal.mNrmMin = np.zeros(StateCount)
+        self.solution_terminal.hRto    = np.zeros(StateCount)
+        self.solution_terminal.MPCmax  = np.ones(StateCount)
+        self.solution_terminal.MPCmin  = np.ones(StateCount)
+
+    def initializeSim(self):
+        IndShockConsumerType.initializeSim(self)
+        if self.global_markov:  #Need to initialize markov state to be the same for all agents
+            base_draw = Uniform().draw(1,seed=self.RNG.randint(0,2**31-1))
+            Cutoffs = np.cumsum(np.array(self.MrkvPrbsInit))
+            self.MrkvNow = np.ones(self.AgentCount)*np.searchsorted(Cutoffs,base_draw).astype(int)
+        self.MrkvNow = self.MrkvNow.astype(int)
+
+    def simDeath(self):
+        '''
+        Determines which agents die this period and must be replaced.  Uses the sequence in LivPrb
+        to determine survival probabilities for each agent.
+
+        Parameters
+        ----------
+        None
+
+        Returns
+        -------
+        which_agents : np.array(bool)
+            Boolean array of size AgentCount indicating which agents die.
+        '''
+        # Determine who dies
+        LivPrb = np.array(self.LivPrb)[self.t_cycle-1,self.MrkvNow] # Time has already advanced, so look back one
+        DiePrb = 1.0 - LivPrb
+        DeathShks = Uniform().draw(N=self.AgentCount,seed=self.RNG.randint(0,2**31-1))
+        which_agents = DeathShks < DiePrb
+        if self.T_age is not None: # Kill agents that have lived for too many periods
+            too_old = self.t_age >= self.T_age
+            which_agents = np.logical_or(which_agents,too_old)
+        return which_agents
+
+    def simBirth(self,which_agents):
+        '''
+        Makes new Markov consumer by drawing initial normalized assets, permanent income levels, and
+        discrete states. Calls IndShockConsumerType.simBirth, then draws from initial Markov distribution.
+
+        Parameters
+        ----------
+        which_agents : np.array(Bool)
+            Boolean array of size self.AgentCount indicating which agents should be "born".
+
+        Returns
+        -------
+        None
+        '''
+        IndShockConsumerType.simBirth(self,which_agents) # Get initial assets and permanent income
+        if not self.global_markov:  #Markov state is not changed if it is set at the global level
+            N = np.sum(which_agents)
+            base_draws = Uniform().draw(N,seed=self.RNG.randint(0,2**31-1))
+            Cutoffs = np.cumsum(np.array(self.MrkvPrbsInit))
+            self.MrkvNow[which_agents] = np.searchsorted(Cutoffs,base_draws).astype(int)
+            
+            
+    def getMarkovStates(self):
+        '''
+        Draw new Markov states for each agent in the simulated population, using
+        the attribute MrkvArray to determine transition probabilities.
+        
+        Parameters
+        ----------
+        None
+
+        Returns
+        -------
+        None
+        '''
+        # Draw random numbers that will be used to determine the next Markov state
+        if self.global_markov:
+            base_draws = np.ones(self.AgentCount)*Uniform().draw(1,seed=self.RNG.randint(0,2**31-1))
+        else:
+            base_draws = Uniform().draw(self.AgentCount,seed=self.RNG.randint(0,2**31-1))
+        dont_change = self.t_age == 0 # Don't change Markov state for those who were just born (unless global_markov)
+        if self.t_sim == 0: # Respect initial distribution of Markov states
+            dont_change[:] = True
+        
+        # Determine which agents are in which states right now
+        J = self.MrkvArray[0].shape[0]
+        MrkvPrev = self.MrkvNow
+        MrkvNow = np.zeros(self.AgentCount,dtype=int)
+        MrkvBoolArray = np.zeros((J,self.AgentCount))
+        for j in range(J):
+            MrkvBoolArray[j,:] = MrkvPrev == j
+        
+        # Draw new Markov states for each agent
+        for t in range(self.T_cycle):
+            Cutoffs = np.cumsum(self.MrkvArray[t],axis=1)
+            right_age = self.t_cycle == t
+            for j in range(J):
+                these = np.logical_and(right_age, MrkvBoolArray[j,:])
+                MrkvNow[these] = np.searchsorted(Cutoffs[j,:],base_draws[these]).astype(int)
+        if not self.global_markov:
+            MrkvNow[dont_change] = MrkvPrev[dont_change]
+        self.MrkvNow = MrkvNow.astype(int)
+    
+
+    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.
+
+        Parameters
+        ----------
+        None
+
+        Returns
+        -------
+        None
+        '''
+        self.getMarkovStates()
+        MrkvNow = self.MrkvNow
+
+        # Now get income shocks for each consumer, by cycle-time and discrete state
+        PermShkNow = np.zeros(self.AgentCount) # Initialize shock arrays
+        TranShkNow = np.zeros(self.AgentCount)
+        for t in range(self.T_cycle):
+            for j in range(self.MrkvArray[t].shape[0]):
+                these = np.logical_and(t == self.t_cycle, j == MrkvNow)
+                N = np.sum(these)
+                if N > 0:
+                    IncomeDstnNow    = self.IncomeDstn[t-1][j] # set current income distribution
+                    PermGroFacNow    = self.PermGroFac[t-1][j] # and permanent growth factor
+                    Indices          = np.arange(IncomeDstnNow.pmf.size) # just a list of integers
+                    # Get random draws of income shocks from the discrete distribution
+                    EventDraws       = IncomeDstnNow.draw_events(
+                        N,
+                        seed=self.RNG.randint(0,2**31-1))
+                    PermShkNow[these] = IncomeDstnNow.X[0][EventDraws]*PermGroFacNow # permanent "shock" includes expected growth
+                    TranShkNow[these] = IncomeDstnNow.X[1][EventDraws]
+        newborn = self.t_age == 0
+        PermShkNow[newborn] = 1.0
+        TranShkNow[newborn] = 1.0
+        self.PermShkNow = PermShkNow
+        self.TranShkNow = TranShkNow
+
+    def readShocks(self):
+        '''
+        A slight modification of AgentType.readShocks that makes sure that MrkvNow is int, not float.
+
+        Parameters
+        ----------
+        None
+
+        Returns
+        -------
+        None
+        '''
+        IndShockConsumerType.readShocks(self)
+        self.MrkvNow = self.MrkvNow.astype(int)
+
+    def getRfree(self):
+        '''
+        Returns an array of size self.AgentCount with interest factor that varies with discrete state.
+
+        Parameters
+        ----------
+        None
+
+        Returns
+        -------
+        RfreeNow : np.array
+             Array of size self.AgentCount with risk free interest rate for each agent.
+        '''
+        RfreeNow = self.Rfree[self.MrkvNow]
+        return RfreeNow
+
+    def getControls(self):
+        '''
+        Calculates consumption for each consumer of this type using the consumption functions.
+
+        Parameters
+        ----------
+        None
+
+        Returns
+        -------
+        None
+        '''
+        cNrmNow = np.zeros(self.AgentCount) + np.nan
+        MPCnow = np.zeros(self.AgentCount) + np.nan
+        J = self.MrkvArray[0].shape[0]
+        
+        MrkvBoolArray = np.zeros((J,self.AgentCount), dtype=bool)
+        for j in range(J):
+            MrkvBoolArray[j,:] = j == self.MrkvNow
+        
+        for t in range(self.T_cycle):
+            right_t = t == self.t_cycle
+            for j in range(J):
+                these = np.logical_and(right_t, MrkvBoolArray[j,:])
+                cNrmNow[these], MPCnow[these] = self.solution[t].cFunc[j].eval_with_derivative(self.mNrmNow[these])
+        self.cNrmNow = cNrmNow
+        self.MPCnow  = MPCnow
+        
+
+    def calcBoundingValues(self):
+        '''
+        Calculate human wealth plus minimum and maximum MPC in an infinite
+        horizon model with only one period repeated indefinitely.  Store results
+        as attributes of self.  Human wealth is the present discounted value of
+        expected future income after receiving income this period, ignoring mort-
+        ality.  The maximum MPC is the limit of the MPC as m --> mNrmMin.  The
+        minimum MPC is the limit of the MPC as m --> infty.  Results are all
+        np.array with elements corresponding to each Markov state.
+
+        NOT YET IMPLEMENTED FOR THIS CLASS
+
+        Parameters
+        ----------
+        None
+
+        Returns
+        -------
+        None
+        '''
+        raise NotImplementedError()
+
+    def makeEulerErrorFunc(self,mMax=100,approx_inc_dstn=True):
+        '''
+        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
+        or to use a (temporary) very dense approximation.
+
+        NOT YET IMPLEMENTED FOR THIS CLASS
+
+        Parameters
+        ----------
+        mMax : float
+            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
+            discrete approximation instead.  When True, uses approximation in
+            IncomeDstn; when False, makes and uses a very dense approximation.
+
+        Returns
+        -------
+        None
+        '''
+        raise NotImplementedError()
diff --git a/HARK/ConsumptionSaving/ConsPortfolioModel.py b/HARK/ConsumptionSaving/ConsPortfolioModel.py
index 7bdbf42d..60b926b9 100644
--- a/HARK/ConsumptionSaving/ConsPortfolioModel.py
+++ b/HARK/ConsumptionSaving/ConsPortfolioModel.py
@@ -22,7 +22,7 @@ from HARK.ConsumptionSaving.ConsGenIncProcessModel import(
     ValueFunc2D,                # For representing 2D value function
     MargValueFunc2D             # For representing 2D marginal value function
 )
-from HARK.distribution import approxLognormal, combineIndepDstns 
+from HARK.distribution import combineIndepDstns 
 from HARK.distribution import Lognormal, Bernoulli # Random draws for simulating agents
 from HARK.interpolation import(
         LinearInterp,           # Piecewise linear interpolation
@@ -242,7 +242,7 @@ class PortfolioConsumerType(IndShockConsumerType):
                 RiskyVar = self.RiskyStd[t] ** 2
                 mu = np.log(self.RiskyAvg[t] / (np.sqrt(1. + RiskyVar / RiskyAvgSqrd)))
                 sigma = np.sqrt(np.log(1. + RiskyVar / RiskyAvgSqrd))
-                RiskyDstn.append(approxLognormal(self.RiskyCount, mu=mu, sigma=sigma))
+                RiskyDstn.append(Lognormal(mu=mu, sigma=sigma).approx(self.RiskyCount))
             self.RiskyDstn = RiskyDstn
             self.addToTimeVary('RiskyDstn')
                 
@@ -253,7 +253,7 @@ class PortfolioConsumerType(IndShockConsumerType):
             RiskyVar = self.RiskyStd ** 2
             mu = np.log(self.RiskyAvg / (np.sqrt(1. + RiskyVar / RiskyAvgSqrd)))
             sigma = np.sqrt(np.log(1. + RiskyVar / RiskyAvgSqrd))
-            self.RiskyDstn = approxLognormal(self.RiskyCount, mu=mu, sigma=sigma)
+            self.RiskyDstn = Lognormal(mu=mu, sigma=sigma).approx(self.RiskyCount)
             self.addToTimeInv('RiskyDstn')
             
             
diff --git a/HARK/ConsumptionSaving/ConsPrefShockModel.py b/HARK/ConsumptionSaving/ConsPrefShockModel.py
index a8e1b392..62c67c3f 100644
--- a/HARK/ConsumptionSaving/ConsPrefShockModel.py
+++ b/HARK/ConsumptionSaving/ConsPrefShockModel.py
@@ -11,7 +11,7 @@ from __future__ import absolute_import
 from builtins import str
 from builtins import range
 import numpy as np
-from HARK.distribution import approxMeanOneLognormal
+from HARK.distribution import MeanOneLogNormal
 from HARK.ConsumptionSaving.ConsIndShockModel import IndShockConsumerType, ConsumerSolution, ConsIndShockSolver, \
                                    ValueFunc, MargValueFunc, KinkedRconsumerType, ConsKinkedRsolver, \
                                    init_idiosyncratic_shocks, init_kinked_R
@@ -111,8 +111,11 @@ class PrefShockConsumerType(IndShockConsumerType):
         PrefShkDstn = [] # discrete distributions of preference shocks
         for t in range(len(self.PrefShkStd)):
             PrefShkStd = self.PrefShkStd[t]
-            PrefShkDstn.append(approxMeanOneLognormal(N=self.PrefShkCount,
-                                                      sigma=PrefShkStd,tail_N=self.PrefShk_tail_N))
+            PrefShkDstn.append(
+                MeanOneLogNormal(
+                    sigma=PrefShkStd
+                ).approx(N=self.PrefShkCount,
+                         tail_N=self.PrefShk_tail_N))
 
         # Store the preference shocks in self (time-varying) and restore time flow
         self.PrefShkDstn = PrefShkDstn
@@ -137,8 +140,9 @@ class PrefShockConsumerType(IndShockConsumerType):
             N = np.sum(these)
             if N > 0:
                 PrefShkNow[these] = self.RNG.permutation(
-                    approxMeanOneLognormal(N,
-                                           sigma=self.PrefShkStd[t]).X)
+                    MeanOneLogNormal(
+                        sigma=self.PrefShkStd[t]
+                    ).approx(N).X)
         self.PrefShkNow = PrefShkNow
 
     def getControls(self):
diff --git a/HARK/distribution.py b/HARK/distribution.py
index 4f0a431a..cf694cb4 100644
--- a/HARK/distribution.py
+++ b/HARK/distribution.py
@@ -72,6 +72,90 @@ class Lognormal():
                                                size=N))
         return draws
 
+    def approx(self, N, tail_N=0, tail_bound=[0.02,0.98], tail_order=np.e):
+        '''
+        Construct a discrete approximation to a lognormal distribution with underlying
+        normal distribution N(mu,sigma).  Makes an equiprobable distribution by
+        default, but user can optionally request augmented tails with exponentially
+        sized point masses.  This can improve solution accuracy in some models.
+
+        Parameters
+        ----------
+        N: int
+            Number of discrete points in the "main part" of the approximation.
+        tail_N: int
+            Number of points in each "tail part" of the approximation; 0 = no tail.
+        tail_bound: [float]
+            CDF boundaries of the tails vs main portion; tail_bound[0] is the lower
+            tail bound, tail_bound[1] is the upper tail bound.  Inoperative when
+            tail_N = 0.  Can make "one tailed" approximations with 0.0 or 1.0.
+        tail_order: float
+            Factor by which consecutive point masses in a "tail part" differ in
+            probability.  Should be >= 1 for sensible spacing.
+
+        Returns
+        -------
+        d : DiscreteDistribution
+            Probability associated with each point in array of discrete
+            points for discrete probability mass function.
+        '''
+        # Find the CDF boundaries of each segment
+        if self.sigma > 0.0:
+            if tail_N > 0:
+                lo_cut     = tail_bound[0]
+                hi_cut     = tail_bound[1]
+            else:
+                lo_cut     = 0.0
+                hi_cut     = 1.0
+            inner_size     = hi_cut - lo_cut
+            inner_CDF_vals = [lo_cut + x*N**(-1.0)*inner_size for x in range(1, N)]
+            if inner_size < 1.0:
+                scale      = 1.0/tail_order
+                mag        = (1.0-scale**tail_N)/(1.0-scale)
+            lower_CDF_vals = [0.0]
+            if lo_cut > 0.0:
+                for x in range(tail_N-1,-1,-1):
+                    lower_CDF_vals.append(lower_CDF_vals[-1] + lo_cut*scale**x/mag)
+            upper_CDF_vals  = [hi_cut]
+            if hi_cut < 1.0:
+                for x in range(tail_N):
+                    upper_CDF_vals.append(upper_CDF_vals[-1] + (1.0-hi_cut)*scale**x/mag)
+            CDF_vals       = lower_CDF_vals + inner_CDF_vals + upper_CDF_vals
+            temp_cutoffs   = list(stats.lognorm.ppf(CDF_vals[1:-1], s=self.sigma, loc=0, scale=np.exp(self.mu)))
+            cutoffs        = [0] + temp_cutoffs + [np.inf]
+            CDF_vals       = np.array(CDF_vals)
+
+            # Construct the discrete approximation by finding the average value within each segment.
+            # This codeblock ignores warnings because it throws a "divide by zero encountered in log"
+            # warning due to computing erf(infty) at the tail boundary.  This is irrelevant and
+            # apparently freaks new users out.
+            with warnings.catch_warnings():
+                warnings.simplefilter("ignore")
+                K              = CDF_vals.size-1 # number of points in approximation
+                pmf            = CDF_vals[1:(K+1)] - CDF_vals[0:K]
+                X              = np.zeros(K)
+                for i in range(K):
+                    zBot  = cutoffs[i]
+                    zTop = cutoffs[i+1]
+                    tempBot = (self.mu+self.sigma**2-np.log(zBot))/(np.sqrt(2)*self.sigma)
+                    tempTop = (self.mu+self.sigma**2-np.log(zTop))/(np.sqrt(2)*self.sigma)
+                    if tempBot <= 4:
+                        X[i] = -0.5*np.exp(self.mu+(self.sigma**2)*0.5)*(erf(tempTop) - erf(tempBot))/pmf[i]
+                    else:
+                        X[i] = -0.5*np.exp(self.mu+(self.sigma**2)*0.5)*(erfc(tempBot) - erfc(tempTop))/pmf[i]
+
+        else:
+            pmf = np.ones(N)/N
+            X   = np.exp(self.mu)*np.ones(N)
+        return DiscreteDistribution(pmf, X)
+
+class MeanOneLogNormal(Lognormal):
+
+    def __init__(self, sigma=1.0):
+        mu = -0.5*sigma**2
+        super().__init__(mu=mu, sigma=sigma)
+
+
 class Normal():
     """
     A Normal distribution.
@@ -92,8 +176,8 @@ class Normal():
             One or more standard deviations. Number of elements T in sigma
             determines number of rows of output.
         '''
-        self.mu = 0.0
-        self.sigma = 1.0
+        self.mu = mu
+        self.sigma = sigma
 
     def draw(self, N, seed=0):
         '''
@@ -127,6 +211,17 @@ class Normal():
                 draws.append(self.sigma[t]*RNG.randn(N) + self.mu[t])
         return draws
 
+    def approx(self, N):
+        """
+        Returns a discrete approximation of this distribution.
+        """
+        x, w = np.polynomial.hermite.hermgauss(N)
+        # normalize w
+        pmf = w*np.pi**-0.5
+        # correct x
+        X = math.sqrt(2.0)*self.sigma*x + self.mu
+        return DiscreteDistribution(pmf, X)
+
 class Weibull():
     '''
     A Weibull distribution
@@ -245,34 +340,26 @@ class Uniform():
                 draws.append(self.bot[t] + (self.top[t] - self.bot[t])*RNG.rand(N))
         return draws
 
+    def approx(self, N):
+        '''
+        Makes a discrete approximation to this uniform distribution.
 
-def drawMeanOneLognormal(N, sigma=1.0, seed=0):
-    '''
-    Generate arrays of mean one lognormal draws. The sigma input can be a number
-    or list-like.  If a number, output is a length N array of draws from the
-    lognormal distribution with standard deviation sigma. If a list, output is
-    a length T list whose t-th entry is a length N array of draws from the
-    lognormal with standard deviation sigma[t].
-
-    Parameters
-    ----------
-    N : int
-        Number of draws in each row.
-    sigma : float or [float]
-        One or more standard deviations. Number of elements T in sigma
-        determines number of rows of output.
-    seed : int
-        Seed for random number generator.
-
-    Returns:
-    ------------
-    draws : np.array or [np.array]
-        T-length list of arrays of mean one lognormal draws each of size N, or
-        a single array of size N (if sigma is a scalar).
-    '''
-    mu = -0.5*sigma**2
+        Parameters
+        ----------
+        N : int
+            The number of points in the discrete approximation
 
-    return Lognormal(mu,sigma).draw(N,seed=seed)
+        Returns
+        -------
+        d : DiscreteDistribution
+            Probability associated with each point in array of discrete
+            points for discrete probability mass function.
+        '''
+        pmf = np.ones(N)/float(N)
+        center = (self.top+self.bot)/2.0
+        width = (self.top-self.bot)/2.0
+        X = center + width*np.linspace(-(N-1.0)/2.0,(N-1.0)/2.0,N)/(N/2.0)
+        return DiscreteDistribution(pmf,X)
 
 
 ### DISCRETE DISTRIBUTIONS
@@ -432,121 +519,8 @@ class DiscreteDistribution():
 
         return draws
 
-def approxLognormal(N, mu=0.0, sigma=1.0, tail_N=0, tail_bound=[0.02,0.98], tail_order=np.e):
-    '''
-    Construct a discrete approximation to a lognormal distribution with underlying
-    normal distribution N(mu,sigma).  Makes an equiprobable distribution by
-    default, but user can optionally request augmented tails with exponentially
-    sized point masses.  This can improve solution accuracy in some models.
-
-    Parameters
-    ----------
-    N: int
-        Number of discrete points in the "main part" of the approximation.
-    mu: float
-        Mean of underlying normal distribution.
-    sigma: float
-        Standard deviation of underlying normal distribution.
-    tail_N: int
-        Number of points in each "tail part" of the approximation; 0 = no tail.
-    tail_bound: [float]
-        CDF boundaries of the tails vs main portion; tail_bound[0] is the lower
-        tail bound, tail_bound[1] is the upper tail bound.  Inoperative when
-        tail_N = 0.  Can make "one tailed" approximations with 0.0 or 1.0.
-    tail_order: float
-        Factor by which consecutive point masses in a "tail part" differ in
-        probability.  Should be >= 1 for sensible spacing.
-
-    Returns
-    -------
-    d : DiscreteDistribution
-        Probability associated with each point in array of discrete
-        points for discrete probability mass function.
-    '''
-    # Find the CDF boundaries of each segment
-    if sigma > 0.0:
-        if tail_N > 0:
-            lo_cut     = tail_bound[0]
-            hi_cut     = tail_bound[1]
-        else:
-            lo_cut     = 0.0
-            hi_cut     = 1.0
-        inner_size     = hi_cut - lo_cut
-        inner_CDF_vals = [lo_cut + x*N**(-1.0)*inner_size for x in range(1, N)]
-        if inner_size < 1.0:
-            scale      = 1.0/tail_order
-            mag        = (1.0-scale**tail_N)/(1.0-scale)
-        lower_CDF_vals = [0.0]
-        if lo_cut > 0.0:
-            for x in range(tail_N-1,-1,-1):
-                lower_CDF_vals.append(lower_CDF_vals[-1] + lo_cut*scale**x/mag)
-        upper_CDF_vals  = [hi_cut]
-        if hi_cut < 1.0:
-            for x in range(tail_N):
-                upper_CDF_vals.append(upper_CDF_vals[-1] + (1.0-hi_cut)*scale**x/mag)
-        CDF_vals       = lower_CDF_vals + inner_CDF_vals + upper_CDF_vals
-        temp_cutoffs   = list(stats.lognorm.ppf(CDF_vals[1:-1], s=sigma, loc=0, scale=np.exp(mu)))
-        cutoffs        = [0] + temp_cutoffs + [np.inf]
-        CDF_vals       = np.array(CDF_vals)
-
-        # Construct the discrete approximation by finding the average value within each segment.
-        # This codeblock ignores warnings because it throws a "divide by zero encountered in log"
-        # warning due to computing erf(infty) at the tail boundary.  This is irrelevant and
-        # apparently freaks new users out.
-        with warnings.catch_warnings():
-            warnings.simplefilter("ignore")
-            K              = CDF_vals.size-1 # number of points in approximation
-            pmf            = CDF_vals[1:(K+1)] - CDF_vals[0:K]
-            X              = np.zeros(K)
-            for i in range(K):
-                zBot  = cutoffs[i]
-                zTop = cutoffs[i+1]
-                tempBot = (mu+sigma**2-np.log(zBot))/(np.sqrt(2)*sigma)
-                tempTop = (mu+sigma**2-np.log(zTop))/(np.sqrt(2)*sigma)
-                if tempBot <= 4:
-                    X[i] = -0.5*np.exp(mu+(sigma**2)*0.5)*(erf(tempTop) - erf(tempBot))/pmf[i]
-                else:
-                    X[i] = -0.5*np.exp(mu+(sigma**2)*0.5)*(erfc(tempBot) - erfc(tempTop))/pmf[i]
-
-    else:
-        pmf = np.ones(N)/N
-        X   = np.exp(mu)*np.ones(N)
-    return DiscreteDistribution(pmf, X)
-
-@memoize
-def approxMeanOneLognormal(N, sigma=1.0, **kwargs):
-    '''
-    Calculate a discrete approximation to a mean one lognormal distribution.
-    Based on function approxLognormal; see that function's documentation for
-    further notes.
-
-    Parameters
-    ----------
-    N : int
-        Size of discrete space vector to be returned.
-    sigma : float
-        standard deviation associated with underlying normal probability distribution.
-
-    Returns
-    -------
-    d : DiscreteDistribution
-        Probability associated with each point in array of discrete
-        points for discrete probability mass function.
-    '''
-    mu_adj = - 0.5*sigma**2;
-    dist = approxLognormal(N=N, mu=mu_adj, sigma=sigma, **kwargs)
-    return dist
-
-def approxNormal(N, mu=0.0, sigma=1.0):
-    x, w = np.polynomial.hermite.hermgauss(N)
-    # normalize w
-    pmf = w*np.pi**-0.5
-    # correct x
-    X = math.sqrt(2.0)*sigma*x + mu
-    return DiscreteDistribution(pmf, X)
-
 def approxLognormalGaussHermite(N, mu=0.0, sigma=1.0):
-    d = approxNormal(N, mu, sigma)
+    d = Normal(mu, sigma).approx(N)
     return DiscreteDistribution(d.pmf, np.exp(d.X))
 
 
@@ -591,32 +565,6 @@ def approxBeta(N,a=1.0,b=1.0):
     pmf  = np.ones(N)/float(N)
     return DiscreteDistribution(pmf, X)
 
-def approxUniform(N,bot=0.0,top=1.0):
-    '''
-    Makes a discrete approximation to a uniform distribution, given its bottom
-    and top limits and number of points.
-
-    Parameters
-    ----------
-    N : int
-        The number of points in the discrete approximation
-    bot : float
-        The bottom of the uniform distribution
-    top : float
-        The top of the uniform distribution
-
-    Returns
-    -------
-    d : DiscreteDistribution
-        Probability associated with each point in array of discrete
-        points for discrete probability mass function.
-    '''
-    pmf = np.ones(N)/float(N)
-    center = (top+bot)/2.0
-    width = (top-bot)/2.0
-    X = center + width*np.linspace(-(N-1.0)/2.0,(N-1.0)/2.0,N)/(N/2.0)
-    return DiscreteDistribution(pmf,X)
-
 
 def makeMarkovApproxToNormal(x_grid,mu,sigma,K=351,bound=3.5):
     '''
diff --git a/HARK/lifecycle.py b/HARK/lifecycle.py
new file mode 100644
index 00000000..73016ad0
--- /dev/null
+++ b/HARK/lifecycle.py
@@ -0,0 +1,29 @@
+from HARK.ConsumptionSaving.ConsIndShockModel import \
+    IndShockConsumerType, ConsIndShockSolverBasic
+import HARK.ConsumptionSaving.ConsumerParameters as Params
+
+LifecycleExample = IndShockConsumerType(
+    **Params.init_lifecycle)
+LifecycleExample.cycles = 1
+LifecycleExample.solve()
+
+# test the solution_terminal
+assert(
+    LifecycleExample.solution[10].cFunc(2).tolist() == 2)
+
+print(LifecycleExample.solution[9].cFunc(1))
+assert(LifecycleExample.solution[9].cFunc(1) == 0.97769632)
+self.assertAlmostEqual(LifecycleExample.solution[8].cFunc(1),
+                       0.96624445)
+self.assertAlmostEqual(LifecycleExample.solution[7].cFunc(1),
+                       0.95691449)
+
+self.assertAlmostEqual(
+    LifecycleExample.solution[0].cFunc(1).tolist(),
+    0.87362789)
+self.assertAlmostEqual(
+    LifecycleExample.solution[1].cFunc(1).tolist(),
+    0.9081621)
+self.assertAlmostEqual(
+    LifecycleExample.solution[2].cFunc(1).tolist(),
+    0.9563899)
diff --git a/HARK/policy.py b/HARK/policy.py
new file mode 100644
index 00000000..2c2fe720
--- /dev/null
+++ b/HARK/policy.py
@@ -0,0 +1,11 @@
+from core import HARKobject
+
+
+class Solution(HARKobject):
+    '''
+    A superclass for representing the "solution" to a single period problem in a
+    dynamic microeconomic model.
+
+    NOTE: This can be deprecated now that HARKobject exists, but this requires
+    replacing each instance of Solution with HARKobject in the other modules.
+    '''
diff --git a/HARK/solver.py b/HARK/solver.py
new file mode 100644
index 00000000..1f8223f9
--- /dev/null
+++ b/HARK/solver.py
@@ -0,0 +1,10 @@
+from core import HARKobject
+
+class Solution(HARKobject):
+    '''
+    A superclass for representing the "solution" to a single period problem in a
+    dynamic microeconomic model.
+
+    NOTE: This can be deprecated now that HARKobject exists, but this requires
+    replacing each instance of Solution with HARKobject in the other modules.
+    '''


Test patch:
diff --git a/HARK/tests/test_approxDstns.py b/HARK/tests/test_approxDstns.py
index 0a24f632..fbe1adf9 100644
--- a/HARK/tests/test_approxDstns.py
+++ b/HARK/tests/test_approxDstns.py
@@ -16,7 +16,7 @@ class testsForDCEGM(unittest.TestCase):
     def test_mu_normal(self):
         for muNormal in self.muNormals:
             for stdNormal in self.stdNormals:
-                d = distribution.approxNormal(40, muNormal)
+                d = distribution.Normal(muNormal).approx(40)
                 self.assertTrue(sum(d.pmf*d.X)-muNormal<1e-12)
 
     def test_mu_lognormal_from_normal(self):
diff --git a/HARK/tests/test_distribution.py b/HARK/tests/test_distribution.py
index c8b27aa5..b0fecef2 100644
--- a/HARK/tests/test_distribution.py
+++ b/HARK/tests/test_distribution.py
@@ -26,7 +26,7 @@ class DistributionClassTests(unittest.TestCase):
 
     def test_drawMeanOneLognormal(self):
         self.assertEqual(
-            distribution.drawMeanOneLognormal(1)[0],
+            distribution.MeanOneLogNormal().draw(1)[0],
             3.5397367004222002)
 
     def test_Lognormal(self):