Reproduce PR #1479

Alan rejiggered the contents of distribution.py, so I have recreated the changes from my prior PR into this one. Committing and then will run tests.

Author: mnwhite

Documentation/CHANGELOG.md

@@ -34,6 +34,7 @@ Release Date: TBD
 - Fixes notation errors in IndShockConsumerType.make_euler_error_func from prior changes. [1495](https://github.com/econ-ark/HARK/pull/1495)
 - Fixes typos in IdentityFunction interpolator class. [1492](https://github.com/econ-ark/HARK/pull/1492)
 - Expands functionality of Cobb-Douglas aggregator for CRRA utility. [1363](https://github.com/econ-ark/HARK/pull/1363)
+- Improved tracking of the bounds of support for distributions, and (some) solvers now respect those bounds when computing the "worst outcome". [XXXX](https://github.com/econ-ark/HARK/pull/1479)
 
 ### 0.15.1
 

HARK/Calibration/Income/IncomeProcesses.py

@@ -74,14 +74,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):
@@ -110,15 +123,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):
@@ -173,7 +193,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):
@@ -239,6 +264,7 @@ def __init__(
             var_names=["PermShk", "TranShk"],
             pmv=joint_dstn.pmv,
             atoms=joint_dstn.atoms,
+            limit=joint_dstn.limit,
             seed=seed,
         )
 

HARK/ConsumptionSaving/ConsIndShockModel.py

@@ -395,31 +395,45 @@ 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:
         m_nrm_min_next (float): Minimum normalized market resources next period.
         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

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)

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,
         )
 
@@ -227,13 +227,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()
@@ -309,7 +311,7 @@ def test_lifecyle(self):
 
         self.assertAlmostEqual(
             LifecycleExample.solution[5].cFunc(3).tolist(),
-            2.12998,
+            2.13004,
             places=HARK_PRECISION,
         )
 
@@ -383,15 +385,15 @@ def test_cyclical(self):
 
         self.assertAlmostEqual(
             CyclicalExample.solution[3].cFunc(3).tolist(),
-            1.59584,
+            1.59597,
             places=HARK_PRECISION,
         )
 
         CyclicalExample.initialize_sim()
         CyclicalExample.simulate()
 
         self.assertAlmostEqual(
-            CyclicalExample.state_now["aLvl"][1], 3.32431, places=HARK_PRECISION
+            CyclicalExample.state_now["aLvl"][1], 3.31950, places=HARK_PRECISION
         )
 
 

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))
 

HARK/distributions/base.py

@@ -43,6 +43,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:
@@ -130,7 +134,7 @@ def discretize(
 
         Returns
         -------
-        DiscreteDistribution
+        discretized_dstn : DiscreteDistribution
             Discretized distribution.
 
         Raises
@@ -149,7 +153,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
 
 
 class MarkovProcess(Distribution):

HARK/distributions/continuous.py

@@ -61,6 +61,8 @@ 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):
         """
@@ -75,8 +77,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
@@ -107,8 +107,6 @@ def _approx_equiprobable(self, N, endpoints=False):
         """
         Returns a discrete equiprobable approximation of this distribution.
 
-        TODO: add endpoints option
-
         Parameters
         ----------
         N : int
@@ -223,6 +221,8 @@ def __init__(
         super().__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
@@ -349,8 +349,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
@@ -446,6 +444,8 @@ def __init__(self, bot=0.0, top=1.0, seed=0):
         super().__init__(
             stats.uniform, loc=self.bot, scale=self.top - self.bot, seed=seed
         )
+        self.infimum = np.array([0.0])
+        self.supremum = np.array([np.inf])
 
     def _approx_equiprobable(self, N, endpoints=False):
         """
@@ -476,7 +476,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,
@@ -506,5 +511,8 @@ 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])