Merge branch 'master' into i569-timing

CHANGELOG.md

@@ -9,6 +9,22 @@ For more information on HARK, see [our Github organization](https://github.com/e
 
 ## Changes
 
+### 0.10.6
+
+Release Date: 17-04-2020
+
+#### Major Changes
+
+* Add Bellman equations for cyclical model example [#600](https://github.com/econ-ark/HARK/pull/600)
+
+* read_shocks now reads mortality as well [#613](https://github.com/econ-ark/HARK/pull/613)
+
+* Discrete probability distributions are now classes [#610](https://github.com/econ-ark/HARK/pull/610)
+
+#### Minor Changes 
+
+
+
 ### 0.10.5
 
 Release Date: 24-03-2020

Documentation/packaging.md

@@ -43,7 +43,7 @@ Follow the directions within https://packaging.python.org/tutorials/packaging-pr
 Deactivate the virtualenv, and make or use a different Python 2.7 virtualenv to test that you can `pip install` the `sdist` you've just made. One thing you can do to quickly verify that the package installs OK:
 
 `>>> import HARK.simulation`
-`>>> HARK.simulation.drawBernoulli(5)`
+`>>> HARK.simulation.Bernoulli().draw(5)`
 
 You should get something like:
 

HARK/ConsumptionSaving/ConsAggShockModel.py

@@ -15,7 +15,7 @@
                                VariableLowerBoundFunc2D, BilinearInterp, LowerEnvelope2D, UpperEnvelope
 from HARK.utilities import CRRAutility, CRRAutilityP, CRRAutilityPP, CRRAutilityP_inv,\
                            CRRAutility_invP, CRRAutility_inv 
-from HARK.simulation import drawUniform
+from HARK.distribution import Uniform
 from HARK.ConsumptionSaving.ConsIndShockModel import ConsumerSolution, IndShockConsumerType, init_idiosyncratic_shocks
 from HARK import HARKobject, Market, AgentType
 from copy import deepcopy
@@ -1629,7 +1629,7 @@ def makeMrkvHist(self):
 
         # Add histories until each state has been visited at least state_T_min times
         while go:
-            draws = drawUniform(N=self.act_T_orig, seed=loops)
+            draws = Uniform().draw(N=self.act_T_orig, seed=loops)
             for s in range(draws.size):  # Add act_T_orig more periods
                 MrkvNow_hist[t] = MrkvNow
                 MrkvNow = np.searchsorted(cutoffs[MrkvNow, :], draws[s])

HARK/ConsumptionSaving/ConsGenIncProcessModel.py

@@ -17,7 +17,7 @@
 from HARK.utilities import CRRAutility, CRRAutilityP, CRRAutilityPP, CRRAutilityP_inv, \
                            CRRAutility_invP, CRRAutility_inv, CRRAutilityP_invP,\
                            getPercentiles
-from HARK.simulation import drawLognormal, drawUniform
+from HARK.distribution import Lognormal, Uniform
 from HARK.ConsumptionSaving.ConsIndShockModel import ConsIndShockSetup, ConsumerSolution, IndShockConsumerType, init_idiosyncratic_shocks
 
 __all__ = ['ValueFunc2D', 'MargValueFunc2D', 'MargMargValueFunc2D', 'pLvlFuncAR1',
@@ -1111,7 +1111,8 @@ def updatepLvlGrid(self):
 
         # Simulate the distribution of persistent income levels by t_cycle in a lifecycle model
         if self.cycles == 1:
-            pLvlNow = drawLognormal(self.AgentCount, mu=self.pLvlInitMean, sigma=self.pLvlInitStd, seed=31382)
+            pLvlNow = Lognormal(self.pLvlInitMean,
+                                sigma=self.pLvlInitStd,).draw(self.AgentCount, seed=31382)
             pLvlGrid = []  # empty list of time-varying persistent income grids
             # Calculate distribution of persistent income in each period of lifecycle
             for t in range(len(self.PermShkStd)):
@@ -1125,14 +1126,15 @@ def updatepLvlGrid(self):
         # Calculate "stationary" distribution in infinite horizon (might vary across periods of cycle)
         elif self.cycles == 0:
             T_long = 1000  # Number of periods to simulate to get to "stationary" distribution
-            pLvlNow = drawLognormal(self.AgentCount, mu=self.pLvlInitMean, sigma=self.pLvlInitStd, seed=31382)
+            pLvlNow = Lognormal(mu=self.pLvlInitMean,
+                                sigma=self.pLvlInitStd).draw(self.AgentCount, seed=31382)
             t_cycle = np.zeros(self.AgentCount, dtype=int)
             for t in range(T_long):
                 LivPrb = LivPrbAll[t_cycle]  # Determine who dies and replace them with newborns
-                draws = drawUniform(self.AgentCount, seed=t)
+                draws = Uniform().draw(self.AgentCount, seed=t)
                 who_dies = draws > LivPrb
-                pLvlNow[who_dies] = drawLognormal(np.sum(who_dies), mu=self.pLvlInitMean,
-                                                  sigma=self.pLvlInitStd, seed=t+92615)
+                pLvlNow[who_dies] = Lognormal(self.pLvlInitMean,
+                                              self.pLvlInitStd).draw(np.sum(who_dies),  seed=t+92615)
                 t_cycle[who_dies] = 0
 
                 for j in range(self.T_cycle):  # Update persistent income
@@ -1175,10 +1177,13 @@ def simBirth(self, which_agents):
         '''
         # Get and store states for newly born agents
         N = np.sum(which_agents)  # Number of new consumers to make
-        aNrmNow_new = drawLognormal(N, mu=self.aNrmInitMean, sigma=self.aNrmInitStd,
+        aNrmNow_new = Lognormal(self.aNrmInitMean,
+                                self.aNrmInitStd).draw(
+                                    N,
                                     seed=self.RNG.randint(0, 2**31-1))
-        self.pLvlNow[which_agents] = drawLognormal(N, mu=self.pLvlInitMean, sigma=self.pLvlInitStd,
-                                                   seed=self.RNG.randint(0, 2**31-1))
+        self.pLvlNow[which_agents] = Lognormal(self.pLvlInitMean,
+                                               self.pLvlInitStd).draw(N, 
+                                                                       seed=self.RNG.randint(0, 2**31-1))
         self.aLvlNow[which_agents] = aNrmNow_new*self.pLvlNow[which_agents]
         self.t_age[which_agents] = 0  # How many periods since each agent was born
         self.t_cycle[which_agents] = 0  # Which period of the cycle each agent is currently in

HARK/ConsumptionSaving/ConsIndShockModel.py

@@ -24,7 +24,7 @@
 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 drawLognormal, drawUniform
+from HARK.distribution import Lognormal, Uniform
 from HARK.distribution import DiscreteDistribution, approxMeanOneLognormal, addDiscreteOutcomeConstantMean, combineIndepDstns 
 from HARK.utilities import makeGridExpMult, CRRAutility, CRRAutilityP, \
                            CRRAutilityPP, CRRAutilityP_inv, CRRAutility_invP, CRRAutility_inv, \
@@ -1743,9 +1743,11 @@ def simBirth(self,which_agents):
         '''
         # Get and store states for newly born agents
         N = np.sum(which_agents) # Number of new consumers to make
-        self.aNrmNow[which_agents] = drawLognormal(N,mu=self.aNrmInitMean,sigma=self.aNrmInitStd,seed=self.RNG.randint(0,2**31-1))
+        self.aNrmNow[which_agents] = Lognormal(mu=self.aNrmInitMean,
+                                               sigma=self.aNrmInitStd).draw(N,seed=self.RNG.randint(0,2**31-1))
         pLvlInitMeanNow = self.pLvlInitMean + np.log(self.PlvlAggNow) # Account for newer cohorts having higher permanent income
-        self.pLvlNow[which_agents] = drawLognormal(N,mu=pLvlInitMeanNow,sigma=self.pLvlInitStd,seed=self.RNG.randint(0,2**31-1))
+        self.pLvlNow[which_agents] = Lognormal(pLvlInitMeanNow,
+                                               self.pLvlInitStd,).draw(N,seed=self.RNG.randint(0,2**31-1))
         self.t_age[which_agents]   = 0 # How many periods since each agent was born
         self.t_cycle[which_agents] = 0 # Which period of the cycle each agent is currently in
         return None
@@ -1767,7 +1769,7 @@ def simDeath(self):
         # Determine who dies
         DiePrb_by_t_cycle = 1.0 - np.asarray(self.LivPrb)
         DiePrb = DiePrb_by_t_cycle[self.t_cycle-1] # Time has already advanced, so look back one
-        DeathShks = drawUniform(N=self.AgentCount,seed=self.RNG.randint(0,2**31-1))
+        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

HARK/ConsumptionSaving/ConsMarkovModel.py

@@ -12,7 +12,8 @@
 from HARK import AgentType
 from HARK.ConsumptionSaving.ConsIndShockModel import ConsIndShockSolver, ValueFunc, \
                              MargValueFunc, ConsumerSolution, IndShockConsumerType
-from HARK.simulation import drawUniform
+
+from HARK.distribution import DiscreteDistribution, Uniform
 from HARK.interpolation import CubicInterp, LowerEnvelope, LinearInterp
 from HARK.utilities import CRRAutility, CRRAutilityP, CRRAutilityPP, CRRAutilityP_inv, \
                            CRRAutility_invP, CRRAutility_inv, CRRAutilityP_invP
@@ -777,7 +778,7 @@ def updateSolutionTerminal(self):
     def initializeSim(self):
         IndShockConsumerType.initializeSim(self)
         if self.global_markov:  #Need to initialize markov state to be the same for all agents
-            base_draw = drawUniform(1,seed=self.RNG.randint(0,2**31-1))
+            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)
@@ -799,7 +800,7 @@ def simDeath(self):
         # 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 = drawUniform(N=self.AgentCount,seed=self.RNG.randint(0,2**31-1))
+        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
@@ -823,7 +824,7 @@ def simBirth(self,which_agents):
         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 = drawUniform(N,seed=self.RNG.randint(0,2**31-1))
+            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)
 
@@ -843,9 +844,9 @@ def getMarkovStates(self):
         '''
         # Draw random numbers that will be used to determine the next Markov state
         if self.global_markov:
-            base_draws = np.ones(self.AgentCount)*drawUniform(1,seed=self.RNG.randint(0,2**31-1))
+            base_draws = np.ones(self.AgentCount)*Uniform().draw(1,seed=self.RNG.randint(0,2**31-1))
         else:
-            base_draws = drawUniform(self.AgentCount,seed=self.RNG.randint(0,2**31-1))
+            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

HARK/ConsumptionSaving/ConsPortfolioModel.py

@@ -23,7 +23,7 @@
     MargValueFunc2D             # For representing 2D marginal value function
 )
 from HARK.distribution import approxLognormal, combineIndepDstns 
-from HARK.simulation import drawLognormal, drawBernoulli # Random draws for simulating agents
+from HARK.distribution import Lognormal, Bernoulli # Random draws for simulating agents
 from HARK.interpolation import(
         LinearInterp,           # Piecewise linear interpolation
         CubicInterp,            # Piecewise cubic interpolation
@@ -353,7 +353,7 @@ def getRisky(self):
 
         mu = np.log(RiskyAvg / (np.sqrt(1. + RiskyVar / RiskyAvgSqrd)))
         sigma = np.sqrt(np.log(1. + RiskyVar / RiskyAvgSqrd))
-        self.RiskyNow = drawLognormal(1, mu=mu, sigma=sigma, seed=self.RNG.randint(0, 2**31-1))
+        self.RiskyNow = Lognormal(mu, sigma).draw(1, seed=self.RNG.randint(0, 2**31-1))
 
 
     def getAdjust(self):
@@ -370,7 +370,7 @@ def getAdjust(self):
         -------
         None
         '''
-        self.AdjustNow = drawBernoulli(self.AgentCount, p=self.AdjustPrb, seed=self.RNG.randint(0, 2**31-1))
+        self.AdjustNow = Bernoulli(self.AdjustPrb).draw(self.AgentCount, seed=self.RNG.randint(0, 2**31-1))
 
 
     def getRfree(self):

HARK/ConsumptionSaving/ConsRepAgentModel.py

@@ -10,7 +10,7 @@
 from builtins import range
 import numpy as np
 from HARK.interpolation import LinearInterp
-from HARK.simulation import drawUniform
+from HARK.distribution import Uniform
 from HARK.ConsumptionSaving.ConsIndShockModel import IndShockConsumerType,\
           ConsumerSolution,MargValueFunc, init_idiosyncratic_shocks
 
@@ -304,7 +304,7 @@ def getShocks(self):
         None
         '''
         cutoffs = np.cumsum(self.MrkvArray[self.MrkvNow,:])
-        MrkvDraw = drawUniform(N=1,seed=self.RNG.randint(0,2**31-1))
+        MrkvDraw = Uniform().draw(N=1,seed=self.RNG.randint(0,2**31-1))
         self.MrkvNow = np.searchsorted(cutoffs,MrkvDraw)
 
         t = self.t_cycle[0]

HARK/ConsumptionSaving/TractableBufferStockModel.py

@@ -30,7 +30,7 @@
 from HARK.utilities import warnings  # Because of "patch" to warnings modules
 from HARK.utilities import CRRAutility, CRRAutilityP, CRRAutilityPP, CRRAutilityPPP, CRRAutilityPPPP, CRRAutilityP_inv, CRRAutility_invP, CRRAutility_inv
 from HARK.interpolation import CubicInterp
-from HARK.simulation import drawLognormal, drawBernoulli
+from HARK.distribution import Lognormal, Bernoulli
 from copy import copy
 from scipy.optimize import newton, brentq
 
@@ -378,7 +378,8 @@ def simBirth(self,which_agents):
         '''
         # Get and store states for newly born agents
         N = np.sum(which_agents) # Number of new consumers to make
-        self.aLvlNow[which_agents] = drawLognormal(N,mu=self.aLvlInitMean,sigma=self.aLvlInitStd,seed=self.RNG.randint(0,2**31-1))
+        self.aLvlNow[which_agents] = Lognormal(self.aLvlInitMean,
+                                               sigma=self.aLvlInitStd).draw(N,seed=self.RNG.randint(0,2**31-1))
         self.eStateNow[which_agents] = 1.0 # Agents are born employed
         self.t_age[which_agents]   = 0 # How many periods since each agent was born
         self.t_cycle[which_agents] = 0 # Which period of the cycle each agent is currently in
@@ -416,7 +417,8 @@ def getShocks(self):
         '''
         employed = self.eStateNow == 1.0
         N = int(np.sum(employed))
-        newly_unemployed = drawBernoulli(N,p=self.UnempPrb,seed=self.RNG.randint(0,2**31-1))
+        newly_unemployed = Bernoulli(self.UnempPrb).draw(N,
+                                                         seed=self.RNG.randint(0,2**31-1))
         self.eStateNow[employed] = 1.0 - newly_unemployed
 
     def getStates(self):

HARK/__init__.py

@@ -1,3 +1,3 @@
 from .core import *
 
-__version__ = "0.10.5"
+__version__ = "0.10.6"