Merge branch 'master' into i583-bl

Author: llorracc

HARK/ConsumptionSaving/ConsAggShockModel.py

@@ -10,12 +10,12 @@
 from builtins import range
 import numpy as np
 import scipy.stats as stats
+from HARK.distribution import DiscreteDistribution, combineIndepDstns, approxMeanOneLognormal
 from HARK.interpolation import LinearInterp, LinearInterpOnInterp1D, ConstantFunction, IdentityFunction,\
                                VariableLowerBoundFunc2D, BilinearInterp, LowerEnvelope2D, UpperEnvelope
 from HARK.utilities import CRRAutility, CRRAutilityP, CRRAutilityPP, CRRAutilityP_inv,\
-                           CRRAutility_invP, CRRAutility_inv, combineIndepDstns,\
-                           approxMeanOneLognormal
-from HARK.simulation import drawDiscrete, drawUniform
+                           CRRAutility_invP, CRRAutility_inv 
+from HARK.simulation import drawUniform
 from HARK.ConsumptionSaving.ConsIndShockModel import ConsumerSolution, IndShockConsumerType, init_idiosyncratic_shocks
 from HARK import HARKobject, Market, AgentType
 from copy import deepcopy
@@ -416,15 +416,18 @@ def addAggShkDstn(self, AggShkDstn):
         state. AggShkDstn is a list of aggregate productivity shock distributions
         for each Markov state.
         '''
-        if len(self.IncomeDstn[0][0]) > 3:
+        if len(self.IncomeDstn[0][0].X) > 2:
             self.IncomeDstn = self.IncomeDstnWithoutAggShocks
         else:
             self.IncomeDstnWithoutAggShocks = self.IncomeDstn
 
         IncomeDstnOut = []
         N = self.MrkvArray.shape[0]
         for t in range(self.T_cycle):
-            IncomeDstnOut.append([combineIndepDstns(self.IncomeDstn[t][n], AggShkDstn[n]) for n in range(N)])
+            IncomeDstnOut.append(
+                [combineIndepDstns(self.IncomeDstn[t][n],
+                                   AggShkDstn[n])
+                 for n in range(N)])
         self.IncomeDstn = IncomeDstnOut
 
     def updateSolutionTerminal(self):
@@ -475,8 +478,12 @@ def getShocks(self):
                 PermGroFacNow = self.PermGroFac[t-1]                # and permanent growth factor
                 Indices = np.arange(IncomeDstnNow[0].size)    # just a list of integers
                 # Get random draws of income shocks from the discrete distribution
-                EventDraws = drawDiscrete(N, X=Indices, P=IncomeDstnNow[0],
-                                          exact_match=True, seed=self.RNG.randint(0, 2**31-1))
+                EventDraws = DiscreteDistribution(
+                    IncomeDstnNow[0],
+                    Indices
+                ).drawDiscrete(N,
+                               exact_match=True,
+                               seed=self.RNG.randint(0, 2**31-1))
                 # permanent "shock" includes expected growth
                 PermShkNow[these] = IncomeDstnNow[1][EventDraws]*PermGroFacNow
                 TranShkNow[these] = IncomeDstnNow[2][EventDraws]
@@ -490,8 +497,10 @@ def getShocks(self):
             PermGroFacNow = self.PermGroFac[0]                # and permanent growth factor
             Indices = np.arange(IncomeDstnNow[0].size)        # just a list of integers
             # Get random draws of income shocks from the discrete distribution
-            EventDraws = drawDiscrete(N, X=Indices, P=IncomeDstnNow[0],
-                                      exact_match=False, seed=self.RNG.randint(0, 2**31-1))
+            EventDraws = DiscreteDistribution(
+                IncomeDstnNow[0], Indices 
+            ).drawDiscrete(N,
+                           seed=self.RNG.randint(0, 2**31-1))
             # permanent "shock" includes expected growth
             PermShkNow[these] = IncomeDstnNow[1][EventDraws]*PermGroFacNow
             TranShkNow[these] = IncomeDstnNow[2][EventDraws]
@@ -1121,10 +1130,13 @@ def makeAggShkHist(self):
         None
         '''
         sim_periods = self.act_T
-        Events = np.arange(self.AggShkDstn[0].size)  # just a list of integers
-        EventDraws = drawDiscrete(N=sim_periods, P=self.AggShkDstn[0], X=Events, seed=0)
-        PermShkAggHist = self.AggShkDstn[1][EventDraws]
-        TranShkAggHist = self.AggShkDstn[2][EventDraws]
+        Events = np.arange(self.AggShkDstn.pmf.size)  # just a list of integers
+        EventDraws = self.AggShkDstn.drawDiscrete(
+            N=sim_periods,
+            X=Events,
+            seed=0)
+        PermShkAggHist = self.AggShkDstn.X[0][EventDraws]
+        TranShkAggHist = self.AggShkDstn.X[1][EventDraws]
 
         # Store the histories
         self.PermShkAggHist = PermShkAggHist*self.PermGroFacAgg
@@ -1352,10 +1364,13 @@ def makeAggShkHist(self):
         None
         '''
         sim_periods = self.act_T
-        Events = np.arange(self.AggShkDstn[0].size)  # just a list of integers
-        EventDraws = drawDiscrete(N=sim_periods, P=self.AggShkDstn[0], X=Events, seed=0)
-        PermShkAggHist = self.AggShkDstn[1][EventDraws]
-        TranShkAggHist = self.AggShkDstn[2][EventDraws]
+        Events = np.arange(self.AggShkDstn.pmf.size)  # just a list of integers
+        EventDraws = self.AggShkDstn.drawDiscrete(
+            N=sim_periods,
+            X=Events,
+            seed=0)
+        PermShkAggHist = self.AggShkDstn.X[0][EventDraws]
+        TranShkAggHist = self.AggShkDstn.X[1][EventDraws]
 
         # Store the histories
         self.PermShkAggHist = PermShkAggHist
@@ -1543,10 +1558,13 @@ def makeAggShkHist(self):
         PermShkAggHistAll = np.zeros((StateCount, sim_periods))
         TranShkAggHistAll = np.zeros((StateCount, sim_periods))
         for i in range(StateCount):
-            Events = np.arange(self.AggShkDstn[i][0].size)  # just a list of integers
-            EventDraws = drawDiscrete(N=sim_periods, P=self.AggShkDstn[i][0], X=Events, seed=0)
-            PermShkAggHistAll[i, :] = self.AggShkDstn[i][1][EventDraws]
-            TranShkAggHistAll[i, :] = self.AggShkDstn[i][2][EventDraws]
+            Events = np.arange(self.AggShkDstn[i].pmf.size)  # just a list of integers
+            EventDraws = self.AggShkDstn[i].drawDiscrete(
+                N=sim_periods,
+                X=Events,
+                seed=0)
+            PermShkAggHistAll[i, :] = self.AggShkDstn[i].X[0][EventDraws]
+            TranShkAggHistAll[i, :] = self.AggShkDstn[i].X[1][EventDraws]
 
         # Select the actual history of aggregate shocks based on the sequence
         # of Markov states that the economy experiences

HARK/ConsumptionSaving/ConsGenIncProcessModel.py

@@ -11,12 +11,13 @@
 from copy import deepcopy
 import numpy as np
 from HARK import AgentType, HARKobject
+from HARK.distribution import DiscreteDistribution
 from HARK.interpolation import LowerEnvelope2D, BilinearInterp, VariableLowerBoundFunc2D, \
                                LinearInterpOnInterp1D, LinearInterp, CubicInterp, UpperEnvelope
 from HARK.utilities import CRRAutility, CRRAutilityP, CRRAutilityPP, CRRAutilityP_inv, \
                            CRRAutility_invP, CRRAutility_inv, CRRAutilityP_invP,\
                            getPercentiles
-from HARK.simulation import drawLognormal, drawDiscrete, drawUniform
+from HARK.simulation import drawLognormal, drawUniform
 from HARK.ConsumptionSaving.ConsIndShockModel import ConsIndShockSetup, ConsumerSolution, IndShockConsumerType, init_idiosyncratic_shocks
 
 __all__ = ['ValueFunc2D', 'MargValueFunc2D', 'MargMargValueFunc2D', 'pLvlFuncAR1',
@@ -1119,8 +1120,9 @@ def updatepLvlGrid(self):
             # Calculate distribution of persistent income in each period of lifecycle
             for t in range(len(self.PermShkStd)):
                 if t > 0:
-                    PermShkNow = drawDiscrete(N=self.AgentCount, P=self.PermShkDstn[t-1][0],
-                                              X=self.PermShkDstn[t-1][1], exact_match=False, seed=t)
+                    PermShkNow = self.PermShkDstn[t-1].drawDiscrete(
+                        N=self.AgentCount,
+                        seed=t)
                     pLvlNow = self.pLvlNextFunc[t-1](pLvlNow)*PermShkNow
                 pLvlGrid.append(getPercentiles(pLvlNow, percentiles=self.pLvlPctiles))
 
@@ -1136,10 +1138,12 @@ def updatepLvlGrid(self):
                 pLvlNow[who_dies] = drawLognormal(np.sum(who_dies), mu=self.pLvlInitMean,
                                                   sigma=self.pLvlInitStd, seed=t+92615)
                 t_cycle[who_dies] = 0
+                
                 for j in range(self.T_cycle):  # Update persistent income
                     these = t_cycle == j
-                    PermShkTemp = drawDiscrete(N=np.sum(these), P=self.PermShkDstn[j][0],
-                                               X=self.PermShkDstn[j][1], exact_match=False, seed=t+13*j)
+                    PermShkTemp = self.PermShkDstn[j].drawDiscrete(
+                        N=np.sum(these),
+                        seed=t+13*j)
                     pLvlNow[these] = self.pLvlNextFunc[j](pLvlNow[these])*PermShkTemp
                 t_cycle = t_cycle + 1
                 t_cycle[t_cycle == self.T_cycle] = 0

HARK/ConsumptionSaving/ConsIndShockModel.py

@@ -24,9 +24,9 @@
 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.simulation import drawDiscrete, drawLognormal, drawUniform
-from HARK.utilities import approxMeanOneLognormal, addDiscreteOutcomeConstantMean,\
-                           combineIndepDstns, makeGridExpMult, CRRAutility, CRRAutilityP, \
+from HARK.simulation import drawLognormal, drawUniform
+from HARK.distribution import DiscreteDistribution, approxMeanOneLognormal, addDiscreteOutcomeConstantMean, combineIndepDstns 
+from HARK.utilities import makeGridExpMult, CRRAutility, CRRAutilityP, \
                            CRRAutilityPP, CRRAutilityP_inv, CRRAutility_invP, CRRAutility_inv, \
                            CRRAutilityP_invP
 
@@ -768,11 +768,10 @@ def setAndUpdateValues(self,solution_next,IncomeDstn,LivPrb,DiscFac):
         ----------
         solution_next : ConsumerSolution
             The solution to next period's one period problem.
-        IncomeDstn : [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
-            probabilities, permanent shocks, transitory shocks.
+        IncomeDstn : distribution.DiscreteDistribution
+            A DiscreteDistribution with a pmf
+            and two point value arrays in X, order:
+            permanent shocks, transitory shocks.
         LivPrb : float
             Survival probability; likelihood of being alive at the beginning of
             the succeeding period.
@@ -783,10 +782,13 @@ def setAndUpdateValues(self,solution_next,IncomeDstn,LivPrb,DiscFac):
         -------
         None
         '''
+        if isinstance(IncomeDstn, list):
+            import pdb; pdb.set_trace()
+
         self.DiscFacEff       = DiscFac*LivPrb # "effective" discount factor
-        self.ShkPrbsNext      = IncomeDstn[0]
-        self.PermShkValsNext  = IncomeDstn[1]
-        self.TranShkValsNext  = IncomeDstn[2]
+        self.ShkPrbsNext      = IncomeDstn.pmf
+        self.PermShkValsNext  = IncomeDstn.X[0]
+        self.TranShkValsNext  = IncomeDstn.X[1]
         self.PermShkMinNext   = np.min(self.PermShkValsNext)
         self.TranShkMinNext   = np.min(self.TranShkValsNext)
         self.vPfuncNext       = solution_next.vPfunc
@@ -2157,11 +2159,12 @@ def getShocks(self):
             if N > 0:
                 IncomeDstnNow    = self.IncomeDstn[t-1] # set current income distribution
                 PermGroFacNow    = self.PermGroFac[t-1] # and permanent growth factor
-                Indices          = np.arange(IncomeDstnNow[0].size) # just a list of integers
                 # Get random draws of income shocks from the discrete distribution
-                EventDraws       = drawDiscrete(N,X=Indices,P=IncomeDstnNow[0],exact_match=False,seed=self.RNG.randint(0,2**31-1))
-                PermShkNow[these] = IncomeDstnNow[1][EventDraws]*PermGroFacNow # permanent "shock" includes expected growth
-                TranShkNow[these] = IncomeDstnNow[2][EventDraws]
+                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]
 
         # That procedure used the *last* period in the sequence for newborns, but that's not right
         # Redraw shocks for newborns, using the *first* period in the sequence.  Approximation.
@@ -2170,11 +2173,12 @@ def getShocks(self):
             these = newborn
             IncomeDstnNow    = self.IncomeDstn[0] # set current income distribution
             PermGroFacNow    = self.PermGroFac[0] # and permanent growth factor
-            Indices          = np.arange(IncomeDstnNow[0].size) # just a list of integers
+           
             # Get random draws of income shocks from the discrete distribution
-            EventDraws       = drawDiscrete(N,X=Indices,P=IncomeDstnNow[0],exact_match=False,seed=self.RNG.randint(0,2**31-1))
-            PermShkNow[these] = IncomeDstnNow[1][EventDraws]*PermGroFacNow # permanent "shock" includes expected growth
-            TranShkNow[these] = IncomeDstnNow[2][EventDraws]
+            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]
 #        PermShkNow[newborn] = 1.0
         TranShkNow[newborn] = 1.0
 
@@ -2375,7 +2379,8 @@ def checkFVAC(self,verbose,public_call):
         Evaluate and report on the Finite Value of Autarky Condition
         Hyperlink to paper: [url]/#Autarky-Value
         '''
-        EpShkuInv = np.dot(self.PermShkDstn[0][0],self.PermShkDstn[0][1]**(1-self.CRRA))
+        EpShkuInv = np.dot(self.PermShkDstn[0].pmf,
+                           self.PermShkDstn[0].X**(1-self.CRRA))
         if self.CRRA != 1.0:
             uInvEpShkuInv = EpShkuInv**(1/(1-self.CRRA)) # The term that gives a utility-consequence-adjusted utility growth
         else:
@@ -2431,8 +2436,9 @@ def checkConditions(self,verbose=False,public_call=True):
         # [url]/#Uncertainty-Modified-Conditions
 
         self.InvPermShkDstn=deepcopy(self.PermShkDstn)
-        self.InvPermShkDstn[0][1]=1/self.PermShkDstn[0][1]
-        EPermShkInv=np.dot(self.InvPermShkDstn[0][0],1/self.PermShkDstn[0][1]) # $\Ex_{t}[\psi^{-1}_{t+1}]$ (in first eqn in sec)
+        self.InvPermShkDstn[0].X = 1/self.PermShkDstn[0].X
+        EPermShkInv=np.dot(self.InvPermShkDstn[0].pmf,
+                           1/self.PermShkDstn[0].X) # $\Ex_{t}[\psi^{-1}_{t+1}]$ (in first eqn in sec)
 
         InvEPermShkInv=(1/EPermShkInv)                          # $\underline{\psi}$ in the paper (\bar{\isp} in private version)
         PermGroFacAdj=self.PermGroFac[0]*InvEPermShkInv                # [url]/#PGroAdj
@@ -2782,7 +2788,9 @@ def constructLognormalIncomeProcessUnemployment(parameters):
             PermShkValsRet  = np.array([1.0])
             TranShkValsRet  = np.array([1.0])
             ShkPrbsRet      = np.array([1.0])
-        IncomeDstnRet = [ShkPrbsRet,PermShkValsRet,TranShkValsRet]
+        IncomeDstnRet = DiscreteDistribution(ShkPrbsRet,
+                                             [PermShkValsRet,
+                                              TranShkValsRet])
 
     # Loop to fill in the list of IncomeDstn random variables.
     for t in range(T_cycle): # Iterate over all periods, counting forward