diff --git a/Documentation/CHANGELOG.md b/Documentation/CHANGELOG.md index 634ef4eef..0382530be 100644 --- a/Documentation/CHANGELOG.md +++ b/Documentation/CHANGELOG.md @@ -25,6 +25,7 @@ Release Date: TBD - Fixes bug in `AgentPopulation` that caused discretization of distributions to not work. [1275](https://github.com/econ-ark/HARK/pull/1275) - Adds support for distributions, booleans, and callables as parameters in the `Parameters` class. [1387](https://github.com/econ-ark/HARK/pull/1387) - Removes a specific way of accounting for ``employment'' in the idiosyncratic-shocks income process. [1473](https://github.com/econ-ark/HARK/pull/1473) +- Improved tracking of the bounds of support for distributions, and (some) solvers now respect those bounds when computing the "worst outcome". [1479](https://github.com/econ-ark/HARK/pull/1479) ### 0.15.1 diff --git a/HARK/Calibration/Income/IncomeProcesses.py b/HARK/Calibration/Income/IncomeProcesses.py index 1221207b4..3095916f9 100644 --- a/HARK/Calibration/Income/IncomeProcesses.py +++ b/HARK/Calibration/Income/IncomeProcesses.py @@ -75,14 +75,27 @@ class LognormPermIncShk(DiscreteDistribution): def __init__(self, sigma, n_approx, neutral_measure=False, seed=0): # Construct an auxiliary discretized normal - logn_approx = MeanOneLogNormal(sigma).discretize( + lognormal_dstn = MeanOneLogNormal(sigma) + logn_approx = lognormal_dstn.discretize( n_approx if sigma > 0.0 else 1, method="equiprobable", tail_N=0 ) + + limit = { + "dist": lognormal_dstn, + "method": "equiprobable", + "N": n_approx, + "endpoints": False, + "infimum": logn_approx.limit["infimum"], + "supremum": logn_approx.limit["supremum"], + } + # Change the pmv if necessary if neutral_measure: logn_approx.pmv = (logn_approx.atoms * logn_approx.pmv).flatten() - super().__init__(pmv=logn_approx.pmv, atoms=logn_approx.atoms, seed=seed) + super().__init__( + pmv=logn_approx.pmv, atoms=logn_approx.atoms, limit=limit, seed=seed + ) class MixtureTranIncShk(DiscreteDistribution): @@ -111,15 +124,22 @@ class MixtureTranIncShk(DiscreteDistribution): """ def __init__(self, sigma, UnempPrb, IncUnemp, n_approx, seed=0): - dstn_approx = MeanOneLogNormal(sigma).discretize( + dstn_orig = MeanOneLogNormal(sigma) + dstn_approx = dstn_orig.discretize( n_approx if sigma > 0.0 else 1, method="equiprobable", tail_N=0 ) + if UnempPrb > 0.0: dstn_approx = add_discrete_outcome_constant_mean( dstn_approx, p=UnempPrb, x=IncUnemp ) - super().__init__(pmv=dstn_approx.pmv, atoms=dstn_approx.atoms, seed=seed) + super().__init__( + pmv=dstn_approx.pmv, + atoms=dstn_approx.atoms, + limit=dstn_approx.limit, + seed=seed, + ) class MixtureTranIncShk_HANK(DiscreteDistribution): @@ -174,7 +194,12 @@ def __init__( # Rescale the transitory shock values to account for new features TranShkMean_temp = (1.0 - tax_rate) * labor * wage dstn_approx.atoms *= TranShkMean_temp - super().__init__(pmv=dstn_approx.pmv, atoms=dstn_approx.atoms, seed=seed) + super().__init__( + pmv=dstn_approx.pmv, + atoms=dstn_approx.atoms, + limit=dstn_approx.limit, + seed=seed, + ) class BufferStockIncShkDstn(DiscreteDistributionLabeled): @@ -240,6 +265,7 @@ def __init__( var_names=["PermShk", "TranShk"], pmv=joint_dstn.pmv, atoms=joint_dstn.atoms, + limit=joint_dstn.limit, seed=seed, ) @@ -319,6 +345,7 @@ def __init__( var_names=["PermShk", "TranShk"], pmv=joint_dstn.pmv, atoms=joint_dstn.atoms, + limit=joint_dstn.limit, seed=seed, ) diff --git a/HARK/ConsumptionSaving/ConsIndShockModel.py b/HARK/ConsumptionSaving/ConsIndShockModel.py index 61b7725cc..d864d8dd4 100644 --- a/HARK/ConsumptionSaving/ConsIndShockModel.py +++ b/HARK/ConsumptionSaving/ConsIndShockModel.py @@ -395,20 +395,26 @@ def solve_one_period_ConsPF( return solution_now -def calc_worst_inc_prob(inc_shk_dstn): +def calc_worst_inc_prob(inc_shk_dstn, use_infimum=True): """Calculate the probability of the worst income shock. Args: inc_shk_dstn (DiscreteDistribution): Distribution of shocks to income. + use_infimum (bool): Indicator for whether to try to use the infimum of the limiting (true) income distribution. """ probs = inc_shk_dstn.pmv perm, tran = inc_shk_dstn.atoms income = perm * tran - worst_inc = np.min(income) + if use_infimum: + worst_inc = np.prod(inc_shk_dstn.limit["infimum"]) + else: + worst_inc = np.min(income) return np.sum(probs[income == worst_inc]) -def calc_boro_const_nat(m_nrm_min_next, inc_shk_dstn, rfree, perm_gro_fac): +def calc_boro_const_nat( + m_nrm_min_next, inc_shk_dstn, rfree, perm_gro_fac, use_infimum=True +): """Calculate the natural borrowing constraint. Args: @@ -416,10 +422,18 @@ def calc_boro_const_nat(m_nrm_min_next, inc_shk_dstn, rfree, perm_gro_fac): inc_shk_dstn (DiscreteDstn): Distribution of shocks to income. rfree (float): Risk free interest factor. perm_gro_fac (float): Permanent income growth factor. + use_infimum (bool): Indicator for whether to use the infimum of the limiting (true) income distribution """ - perm, tran = inc_shk_dstn.atoms - temp_fac = (perm_gro_fac * np.min(perm)) / rfree - return (m_nrm_min_next - np.min(tran)) * temp_fac + if use_infimum: + perm_min, tran_min = inc_shk_dstn.limit["infimum"] + else: + perm, tran = inc_shk_dstn.atoms + perm_min = np.min(perm) + tran_min = np.min(tran) + + temp_fac = (perm_gro_fac * perm_min) / rfree + boro_cnst_nat = (m_nrm_min_next - tran_min) * temp_fac + return boro_cnst_nat def calc_m_nrm_min(boro_const_art, boro_const_nat): @@ -823,7 +837,7 @@ def solve_one_period_ConsKinkedR( DiscFacEff = DiscFac * LivPrb # "effective" discount factor # Calculate the probability that we get the worst possible income draw - WorstIncPrb = calc_worst_inc_prob(IncShkDstn) + WorstIncPrb = calc_worst_inc_prob(IncShkDstn, use_infimum=False) # WorstIncPrb is the "Weierstrass p" concept: the odds we get the WORST thing Ex_IncNext = expected(lambda x: x["PermShk"] * x["TranShk"], IncShkDstn) hNrmNow = calc_human_wealth(solution_next.hNrm, PermGroFac, Rsave, Ex_IncNext) @@ -835,7 +849,11 @@ def solve_one_period_ConsKinkedR( # Calculate the minimum allowable value of money resources in this period BoroCnstNat = calc_boro_const_nat( - solution_next.mNrmMin, IncShkDstn, Rboro, PermGroFac + solution_next.mNrmMin, + IncShkDstn, + Rboro, + PermGroFac, + use_infimum=False, ) # Set the minimum allowable (normalized) market resources based on the natural # and artificial borrowing constraints diff --git a/HARK/ConsumptionSaving/tests/test_ConsNewKeynesianModel.py b/HARK/ConsumptionSaving/tests/test_ConsNewKeynesianModel.py index c2b654882..ec15814b3 100644 --- a/HARK/ConsumptionSaving/tests/test_ConsNewKeynesianModel.py +++ b/HARK/ConsumptionSaving/tests/test_ConsNewKeynesianModel.py @@ -39,7 +39,7 @@ def test_calc_tran_matrix(self): AggC = np.dot(gridc.flatten(), vecDstn) # Aggregate Consumption AggA = np.dot(grida.flatten(), vecDstn) # Aggregate Assets - self.assertAlmostEqual(AggA[0], 0.82984, places=4) + self.assertAlmostEqual(AggA[0], 0.82960, places=4) self.assertAlmostEqual(AggC[0], 1.00780, places=4) @@ -52,6 +52,6 @@ def test_calc_jacobian(self): Agent.compute_steady_state() CJAC_Perm, AJAC_Perm = Agent.calc_jacobian("PermShkStd", 50) - self.assertAlmostEqual(CJAC_Perm.T[30][29], -0.10503, places=HARK_PRECISION) + self.assertAlmostEqual(CJAC_Perm.T[30][29], -0.10508, places=HARK_PRECISION) self.assertAlmostEqual(CJAC_Perm.T[30][30], 0.10316, places=HARK_PRECISION) self.assertAlmostEqual(CJAC_Perm.T[30][31], 0.09059, places=HARK_PRECISION) diff --git a/HARK/ConsumptionSaving/tests/test_IndShockConsumerType.py b/HARK/ConsumptionSaving/tests/test_IndShockConsumerType.py index 7f8e0a421..6178791d8 100644 --- a/HARK/ConsumptionSaving/tests/test_IndShockConsumerType.py +++ b/HARK/ConsumptionSaving/tests/test_IndShockConsumerType.py @@ -50,17 +50,17 @@ def test_ConsIndShockSolverBasic(self): self.assertAlmostEqual( LifecycleExample.solution[0].cFunc(1).tolist(), - 0.75062, + 0.75074, places=HARK_PRECISION, ) self.assertAlmostEqual( LifecycleExample.solution[1].cFunc(1).tolist(), - 0.75864, + 0.75876, places=HARK_PRECISION, ) self.assertAlmostEqual( LifecycleExample.solution[2].cFunc(1).tolist(), - 0.76812, + 0.76824, places=HARK_PRECISION, ) @@ -223,13 +223,15 @@ def test_infinite_horizon(self): IndShockExample.solve() self.assertAlmostEqual( - IndShockExample.solution[0].mNrmStE, 1.54882, places=HARK_PRECISION - ) - self.assertAlmostEqual( - IndShockExample.solution[0].cFunc.functions[0].x_list[0], - -0.25018, - places=HARK_PRECISION, + IndShockExample.solution[0].mNrmStE, 1.54765, places=HARK_PRECISION ) + # self.assertAlmostEqual( + # IndShockExample.solution[0].cFunc.functions[0].x_list[0], + # -0.25018, + # places=HARK_PRECISION, + # ) + # This test is commented out because it was trivialized by revisions to the "worst income shock" code. + # The bottom x value of the unconstrained consumption function will definitely be zero, so this is pointless. IndShockExample.track_vars = ["aNrm", "mNrm", "cNrm", "pLvl"] IndShockExample.initialize_sim() @@ -297,7 +299,7 @@ def test_lifecyle(self): self.assertAlmostEqual( LifecycleExample.solution[5].cFunc(3).tolist(), - 2.12998, + 2.13004, places=HARK_PRECISION, ) @@ -371,7 +373,7 @@ def test_cyclical(self): self.assertAlmostEqual( CyclicalExample.solution[3].cFunc(3).tolist(), - 1.59584, + 1.59597, places=HARK_PRECISION, ) @@ -379,7 +381,7 @@ def test_cyclical(self): CyclicalExample.simulate() self.assertAlmostEqual( - CyclicalExample.state_now["aLvl"][1], 3.32431, places=HARK_PRECISION + CyclicalExample.state_now["aLvl"][1], 3.31950, places=HARK_PRECISION ) diff --git a/HARK/ConsumptionSaving/tests/test_modelcomparisons.py b/HARK/ConsumptionSaving/tests/test_modelcomparisons.py index 0b4669e16..3cbc81a5c 100644 --- a/HARK/ConsumptionSaving/tests/test_modelcomparisons.py +++ b/HARK/ConsumptionSaving/tests/test_modelcomparisons.py @@ -47,7 +47,7 @@ def setUp(self): test_dictionary["UnempPrb"] = 0.0 test_dictionary["T_cycle"] = 1 test_dictionary["T_retire"] = 0 - test_dictionary["BoroCnstArt"] = None + test_dictionary["BoroCnstArt"] = 0.0 InfiniteType = IndShockConsumerType(**test_dictionary) InfiniteType.cycles = 0 @@ -79,7 +79,7 @@ def diffFunc(m): m ) - self.InfiniteType.cFunc[0](m) - points = np.arange(0.5, mNrmMinInf + aXtraMin, mNrmMinInf + aXtraMax) + points = np.arange(0.5, mNrmMinInf + aXtraMax, mNrmMinInf + aXtraMin) difference = diffFunc(points) max_difference = np.max(np.abs(difference)) diff --git a/HARK/distribution.py b/HARK/distribution.py index c53ea166d..cec18956c 100644 --- a/HARK/distribution.py +++ b/HARK/distribution.py @@ -53,6 +53,10 @@ def __init__(self, seed: Optional[int] = 0) -> None: self._seed: int = _seed self._rng: random.Generator = random.default_rng(self._seed) + # Bounds of distribution support should be overwritten by subclasses + self.infimum = np.array([]) + self.supremum = np.array([]) + @property def seed(self) -> int: """ @@ -139,7 +143,7 @@ def discretize( Returns ------- - DiscreteDistribution + discretized_dstn : DiscreteDistribution Discretized distribution. Raises @@ -158,7 +162,10 @@ def discretize( ) approx = getattr(self, approx_method) - return approx(N, endpoints, **kwds) + discretized_dstn = approx(N, endpoints, **kwds) + discretized_dstn.limit["infimum"] = self.infimum.copy() + discretized_dstn.limit["supremum"] = self.supremum.copy() + return discretized_dstn # CONTINUOUS DISTRIBUTIONS @@ -214,6 +221,9 @@ def __init__(self, mu=0.0, sigma=1.0, seed=0): super().__init__(stats.norm, loc=mu, scale=sigma, seed=seed) + self.infimum = -np.inf * np.ones(self.mu.size) + self.supremum = np.inf * np.ones(self.mu.size) + def discretize(self, N, method="hermite", endpoints=False): """ For normal distributions, the Gauss-Hermite quadrature rule is @@ -227,8 +237,6 @@ def _approx_hermite(self, N, endpoints=False): Returns a discrete approximation of this distribution using the Gauss-Hermite quadrature rule. - TODO: add endpoints option - Parameters ---------- N : int @@ -259,8 +267,6 @@ def _approx_equiprobable(self, N, endpoints=False): """ Returns a discrete equiprobable approximation of this distribution. - TODO: add endpoints option - Parameters ---------- N : int @@ -371,6 +377,9 @@ def __init__( stats.lognorm, s=self.sigma, scale=np.exp(self.mu), loc=0, seed=seed ) + self.infimum = np.array([0.0]) + self.supremum = np.array([np.inf]) + def _approx_equiprobable( self, N, endpoints=False, tail_N=0, tail_bound=None, tail_order=np.e ): @@ -594,6 +603,9 @@ def __init__(self, bot=0.0, top=1.0, seed=0): stats.uniform, loc=self.bot, scale=self.top - self.bot, seed=seed ) + self.infimum = np.array([self.bot]) + self.supremum = np.array([self.top]) + def _approx_equiprobable(self, N, endpoints=False): """ Makes a discrete approximation to this uniform distribution. @@ -623,7 +635,12 @@ def _approx_equiprobable(self, N, endpoints=False): atoms = np.concatenate(([self.bot], atoms, [self.top])) pmv = np.concatenate(([0.0], pmv, [0.0])) - limit = {"dist": self, "method": "equiprobable", "N": N, "endpoints": endpoints} + limit = { + "dist": self, + "method": "equiprobable", + "N": N, + "endpoints": endpoints, + } return DiscreteDistribution( pmv, @@ -653,9 +670,13 @@ class Weibull(ContinuousFrozenDistribution): def __init__(self, scale=1.0, shape=1.0, seed=0): self.scale = np.asarray(scale) self.shape = np.asarray(shape) + # Set up the RNG super().__init__(stats.weibull_min, c=shape, scale=scale, seed=seed) + self.infimum = np.array([0.0]) + self.supremum = np.array([np.inf]) + # MULTIVARIATE DISTRIBUTIONS @@ -682,6 +703,9 @@ def __init__(self, mu=[1, 1], Sigma=[[1, 0], [0, 1]], seed=0): multivariate_normal_frozen.__init__(self, mean=self.mu, cov=self.Sigma) Distribution.__init__(self, seed=seed) + self.infimum = -np.inf * np.ones(self.M) + self.supremum = np.inf * np.ones(self.M) + def discretize(self, N, method="hermite", endpoints=False): """ For multivariate normal distributions, the Gauss-Hermite @@ -786,6 +810,8 @@ def __init__(self, p=0.5, seed=0): self.pmv = [1 - self.p, self.p] self.atoms = [0, 1] self.limit = {"dist": self} + self.infimum = np.array([0.0]) + self.supremum = np.array([1.0]) class DiscreteDistribution(Distribution): @@ -801,6 +827,10 @@ class DiscreteDistribution(Distribution): For multivariate distributions, the last dimension of atoms must index "atom" or the random realization. For instance, if atoms.shape == (2,6,4), the random variable has 4 possible realizations and each of them has shape (2,6). + limit : dict + Dictionary with information about the continuous distribution from which + this distribution was generated. The reference distribution is in the entry + called 'dist'. seed : int Seed for random number generator. """ @@ -816,6 +846,11 @@ def __init__( self.pmv = np.asarray(pmv) self.atoms = np.atleast_2d(atoms) + if limit is None: + limit = { + "infimum": np.min(self.atoms, axis=-1), + "supremum": np.max(self.atoms, axis=-1), + } self.limit = limit # Check that pmv and atoms have compatible dimensions. @@ -1095,7 +1130,12 @@ def __init__( ) attrs = {} if attrs is None else attrs - limit = {} if limit is None else limit + if limit is None: + limit = { + "infimum": np.min(self.atoms, axis=-1), + "supremum": np.max(self.atoms, axis=-1), + } + self.limit = limit attrs.update(limit) attrs["name"] = name @@ -1526,6 +1566,7 @@ def draw(self, condition): return draws +# TODO: This function does not generate the limit attribute def approx_lognormal_gauss_hermite(N, mu=0.0, sigma=1.0, seed=0): d = Normal(mu, sigma).discretize(N, method="hermite") return DiscreteDistribution(d.pmv, np.exp(d.atoms), seed=seed) @@ -1763,8 +1804,17 @@ def add_discrete_outcome_constant_mean(distribution, x, p, sort=False): Probability associated with each point in array of discrete points for discrete probability mass function. """ + if type(distribution) == TimeVaryingDiscreteDistribution: + # apply recursively on all the internal distributions + return TimeVaryingDiscreteDistribution( + [ + add_discrete_outcome_constant_mean(d, x, p) + for d in distribution.distributions + ], + seed=distribution.seed, + ) - if type(distribution) != TimeVaryingDiscreteDistribution: + else: atoms = np.append(x, distribution.atoms * (1 - p * x) / (1 - p)) pmv = np.append(p, distribution.pmv * (1 - p)) @@ -1773,16 +1823,37 @@ def add_discrete_outcome_constant_mean(distribution, x, p, sort=False): atoms = atoms[indices] pmv = pmv[indices] - return DiscreteDistribution(pmv, atoms) - elif type(distribution) == TimeVaryingDiscreteDistribution: - # apply recursively on all the internal distributions - return TimeVaryingDiscreteDistribution( - [ - add_discrete_outcome_constant_mean(d, x, p) - for d in distribution.distributions - ], - seed=distribution.seed, - ) + # Update infimum and supremum + temp_x = np.array(x, ndmin=1) + try: + infimum = np.array( + [ + np.minimum(temp_x[i], distribution.limit["infimum"][i]) + for i in range(temp_x.size) + ] + ) + except: + infimum = np.min(atoms, axis=-1, keepdims=True) + try: + supremum = np.array( + [ + np.maximum(temp_x[i], distribution.limit["supremum"][i]) + for i in range(temp_x.size) + ] + ) + except: + supremum = np.max(atoms, axis=-1, keepdims=True) + + limit = { + "dist": distribution, + "method": "add_discrete_outcome_constant_mean", + "x": x, + "p": p, + "infimum": infimum, + "supremum": supremum, + } + + return DiscreteDistribution(pmv, atoms, seed=distribution.seed, limit=limit) def add_discrete_outcome(distribution, x, p, sort=False): @@ -1814,7 +1885,38 @@ def add_discrete_outcome(distribution, x, p, sort=False): atoms = atoms[indices] pmv = pmv[indices] - return DiscreteDistribution(pmv, atoms) + # Update infimum and supremum + # Update infimum and supremum + temp_x = np.array(x, ndmin=1) + try: + infimum = np.array( + [ + np.minimum(temp_x[i], distribution.limit["infimum"][i]) + for i in range(temp_x.size) + ] + ) + except: + infimum = np.min(atoms, axis=-1, keepdims=True) + try: + supremum = np.array( + [ + np.maximum(temp_x[i], distribution.limit["supremum"][i]) + for i in range(temp_x.size) + ] + ) + except: + supremum = np.max(atoms, axis=-1, keepdims=True) + + limit = { + "dist": distribution, + "method": "add_discrete_outcome", + "x": x, + "p": p, + "infimum": infimum, + "supremum": supremum, + } + + return DiscreteDistribution(pmv, atoms, seed=distribution.seed, limit=limit) def combine_indep_dstns(*distributions, seed=0): @@ -1827,6 +1929,8 @@ def combine_indep_dstns(*distributions, seed=0): distributions : DiscreteDistribution Arbitrary number of discrete distributionss to combine. Their realizations must be vector-valued (for each D in distributions, it must be the case that len(D.dim())==1). + seed : int, optional + Value to use as the RNG seed for the combined distribution, default is 0. Returns ------- @@ -1854,7 +1958,7 @@ def combine_indep_dstns(*distributions, seed=0): else: var_labels += tuple([""] * dist.dim()[0]) - number_of_distributions = len(distributions) + dstn_count = len(distributions) all_labeled = all(dist_is_labeled) labels_are_unique = len(var_labels) == len(set(var_labels)) @@ -1878,11 +1982,28 @@ def combine_indep_dstns(*distributions, seed=0): assert np.isclose(np.sum(P_out), 1), "Probabilities do not sum to 1!" + # Make the limit dictionary + infimum = np.concatenate( + [distributions[i].limit["infimum"] for i in range(dstn_count)] + ) + supremum = np.concatenate( + [distributions[i].limit["supremum"] for i in range(dstn_count)] + ) + limit = { + "dist": distributions, + "method": "combine_indep_dstns", + "infimum": infimum, + "supremum": supremum, + } + # except: + # limit=None + if all_labeled and labels_are_unique: combined_dstn = DiscreteDistributionLabeled( pmv=P_out, atoms=atoms_out, var_names=var_labels, + limit=limit, seed=seed, ) else: @@ -1890,7 +2011,7 @@ def combine_indep_dstns(*distributions, seed=0): warn( "There are duplicated labels in the provided distributions. Returning a non-labeled combination" ) - combined_dstn = DiscreteDistribution(P_out, atoms_out, seed=seed) + combined_dstn = DiscreteDistribution(P_out, atoms_out, limit=limit, seed=seed) return combined_dstn