From 52cc680de47cde941de2352889a03f7f53e2e83f Mon Sep 17 00:00:00 2001 From: "Matthew N. White" Date: Thu, 25 Jul 2024 16:03:39 -0400 Subject: [PATCH 1/8] Add infimum and supremum information to limit dicts The limit attribute of a DiscreteDistribution has data about the continuous distribution from which it was derived. This dictionary now contains infimum and supremum fields, which are also derived from the original distribution. It's present for *at least* the income distributions used by ConsIndShock. --- HARK/Calibration/Income/IncomeProcesses.py | 29 +++- HARK/distribution.py | 146 ++++++++++++++++++--- 2 files changed, 151 insertions(+), 24 deletions(-) diff --git a/HARK/Calibration/Income/IncomeProcesses.py b/HARK/Calibration/Income/IncomeProcesses.py index 1221207b4..3542a53e6 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): @@ -240,6 +260,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/distribution.py b/HARK/distribution.py index c53ea166d..0a5ac7afe 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 = self.bot + self.supremum = self.top + def _approx_equiprobable(self, N, endpoints=False): """ Makes a discrete approximation to this uniform distribution. @@ -653,9 +665,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 +698,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 +805,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 +822,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. """ @@ -1526,6 +1551,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 +1789,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 +1808,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 +1870,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 +1914,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 +1943,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 +1967,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 +1996,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 From 7bba4d74f1cf7d97c48b4ed85bbb4729e7214520 Mon Sep 17 00:00:00 2001 From: "Matthew N. White" Date: Tue, 13 Aug 2024 15:05:28 -0400 Subject: [PATCH 2/8] ConsIndShock solver draws on infimum of IncShkDstn I made these changes almost three weeks ago but forgot to commit. The solver for ConsIndShock now uses information in the limit dictionary of IncShkDstn when determining the "worst income shock", under the assumption that the "joint infimum" exists. As a default backup, it uses the lowest discretized shock, like before. This does not really change the behavior of the solution unless you both add points *extremely* close to zero into aXtraGrid *and* actually include the infimum (as a tiny point mass) in the discretization itself. --- HARK/ConsumptionSaving/ConsIndShockModel.py | 18 ++++++++++++++---- 1 file changed, 14 insertions(+), 4 deletions(-) diff --git a/HARK/ConsumptionSaving/ConsIndShockModel.py b/HARK/ConsumptionSaving/ConsIndShockModel.py index 61b7725cc..4f63c2ff2 100644 --- a/HARK/ConsumptionSaving/ConsIndShockModel.py +++ b/HARK/ConsumptionSaving/ConsIndShockModel.py @@ -404,7 +404,10 @@ def calc_worst_inc_prob(inc_shk_dstn): probs = inc_shk_dstn.pmv perm, tran = inc_shk_dstn.atoms income = perm * tran - worst_inc = np.min(income) + try: + worst_inc = np.prod(inc_shk_dstn.limit['infimum']) + except: + worst_inc = np.min(income) return np.sum(probs[income == worst_inc]) @@ -417,9 +420,16 @@ def calc_boro_const_nat(m_nrm_min_next, inc_shk_dstn, rfree, perm_gro_fac): rfree (float): Risk free interest factor. perm_gro_fac (float): Permanent income growth factor. """ - 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 + try: + perm_min, tran_min = inc_shk_dstn.limit['infimum'] + except: + 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): From d18bc4ef48c8b5a64675a856f2f39b9c11fc3402 Mon Sep 17 00:00:00 2001 From: "Matthew N. White" Date: Tue, 13 Aug 2024 15:17:31 -0400 Subject: [PATCH 3/8] Forgot to run pre-commit --- HARK/ConsumptionSaving/ConsIndShockModel.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/HARK/ConsumptionSaving/ConsIndShockModel.py b/HARK/ConsumptionSaving/ConsIndShockModel.py index 4f63c2ff2..22bdeca65 100644 --- a/HARK/ConsumptionSaving/ConsIndShockModel.py +++ b/HARK/ConsumptionSaving/ConsIndShockModel.py @@ -405,7 +405,7 @@ def calc_worst_inc_prob(inc_shk_dstn): perm, tran = inc_shk_dstn.atoms income = perm * tran try: - worst_inc = np.prod(inc_shk_dstn.limit['infimum']) + worst_inc = np.prod(inc_shk_dstn.limit["infimum"]) except: worst_inc = np.min(income) return np.sum(probs[income == worst_inc]) @@ -421,12 +421,12 @@ def calc_boro_const_nat(m_nrm_min_next, inc_shk_dstn, rfree, perm_gro_fac): perm_gro_fac (float): Permanent income growth factor. """ try: - perm_min, tran_min = inc_shk_dstn.limit['infimum'] + perm_min, tran_min = inc_shk_dstn.limit["infimum"] except: 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 From db51e9ead201fd06281dbc241871f0ad291dfe93 Mon Sep 17 00:00:00 2001 From: "Matthew N. White" Date: Tue, 13 Aug 2024 16:06:08 -0400 Subject: [PATCH 4/8] Adjust test targets, fix HANK income distribution Some consumption function test targets moved in the 5th digit because of the code changes. Targets have been adjusted. Probably still some tests failing. --- HARK/Calibration/Income/IncomeProcesses.py | 8 +++++- .../tests/test_ConsNewKeynesianModel.py | 4 +-- .../tests/test_IndShockConsumerType.py | 26 ++++++++++--------- 3 files changed, 23 insertions(+), 15 deletions(-) diff --git a/HARK/Calibration/Income/IncomeProcesses.py b/HARK/Calibration/Income/IncomeProcesses.py index 3542a53e6..3095916f9 100644 --- a/HARK/Calibration/Income/IncomeProcesses.py +++ b/HARK/Calibration/Income/IncomeProcesses.py @@ -194,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): @@ -340,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/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 ) From 739400c6c8701eceac80d6aa075effd5b334ee3d Mon Sep 17 00:00:00 2001 From: "Matthew N. White" Date: Tue, 13 Aug 2024 16:38:50 -0400 Subject: [PATCH 5/8] Add inf and sup to DiscreteDistribution(Labeled) Also fixed inf/sup for Uniform class. One or two tests might still fail. --- HARK/distribution.py | 23 +++++++++++++++++++---- 1 file changed, 19 insertions(+), 4 deletions(-) diff --git a/HARK/distribution.py b/HARK/distribution.py index 0a5ac7afe..4ec611b87 100644 --- a/HARK/distribution.py +++ b/HARK/distribution.py @@ -603,8 +603,8 @@ 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 = self.bot - self.supremum = self.top + self.infimum = np.array([self.bot]) + self.supremum = np.array([self.top]) def _approx_equiprobable(self, N, endpoints=False): """ @@ -635,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, @@ -841,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, keepdims=True), + "supremum": np.max(self.atoms, axis=-1, keepdims=True), + } self.limit = limit # Check that pmv and atoms have compatible dimensions. @@ -1120,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, keepdims=True), + "supremum": np.max(self.atoms, axis=-1, keepdims=True), + } + self.limit = limit attrs.update(limit) attrs["name"] = name From ea28e3a540de264292c8a7e1fdbf47922003128b Mon Sep 17 00:00:00 2001 From: "Matthew N. White" Date: Tue, 13 Aug 2024 17:13:08 -0400 Subject: [PATCH 6/8] Fix two more tests Order of arguments was wrong on np.arange, and BoroCnstArt needs to be zero to compare to riskless ConsIndShock after changes. Also fixed issue with dimensionality of supremum and infimum. --- HARK/ConsumptionSaving/tests/test_modelcomparisons.py | 4 ++-- HARK/distribution.py | 8 ++++---- 2 files changed, 6 insertions(+), 6 deletions(-) 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 4ec611b87..cec18956c 100644 --- a/HARK/distribution.py +++ b/HARK/distribution.py @@ -848,8 +848,8 @@ def __init__( self.atoms = np.atleast_2d(atoms) if limit is None: limit = { - "infimum": np.min(self.atoms, axis=-1, keepdims=True), - "supremum": np.max(self.atoms, axis=-1, keepdims=True), + "infimum": np.min(self.atoms, axis=-1), + "supremum": np.max(self.atoms, axis=-1), } self.limit = limit @@ -1132,8 +1132,8 @@ def __init__( attrs = {} if attrs is None else attrs if limit is None: limit = { - "infimum": np.min(self.atoms, axis=-1, keepdims=True), - "supremum": np.max(self.atoms, axis=-1, keepdims=True), + "infimum": np.min(self.atoms, axis=-1), + "supremum": np.max(self.atoms, axis=-1), } self.limit = limit attrs.update(limit) From fcffb80938591dedecc856cd078f2898727d02f8 Mon Sep 17 00:00:00 2001 From: "Matthew N. White" Date: Tue, 13 Aug 2024 17:43:01 -0400 Subject: [PATCH 7/8] Fixed last test In order for the KinkedR model to actually do anything, the infimum of the true income distribution has to be ignored. --- HARK/ConsumptionSaving/ConsIndShockModel.py | 24 ++++++++++++++------- 1 file changed, 16 insertions(+), 8 deletions(-) diff --git a/HARK/ConsumptionSaving/ConsIndShockModel.py b/HARK/ConsumptionSaving/ConsIndShockModel.py index 22bdeca65..d864d8dd4 100644 --- a/HARK/ConsumptionSaving/ConsIndShockModel.py +++ b/HARK/ConsumptionSaving/ConsIndShockModel.py @@ -395,23 +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 - try: + if use_infimum: worst_inc = np.prod(inc_shk_dstn.limit["infimum"]) - except: + 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: @@ -419,10 +422,11 @@ 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 """ - try: + if use_infimum: perm_min, tran_min = inc_shk_dstn.limit["infimum"] - except: + else: perm, tran = inc_shk_dstn.atoms perm_min = np.min(perm) tran_min = np.min(tran) @@ -833,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) @@ -845,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 From ca1f734c55a6bd9887d25df34a665e732db26f24 Mon Sep 17 00:00:00 2001 From: "Matthew N. White" Date: Tue, 20 Aug 2024 13:23:12 -0400 Subject: [PATCH 8/8] Update CHANGELOG.md --- Documentation/CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) 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