Fix CRRA=1 (log utility) division by zero by implementing proper limit formula

Added vNvrsSlope() function to HARK/rewards.py that returns MPC when CRRA=1
instead of the divergent formula MPC^(-CRRA/(1-CRRA)). This implements the
mathematically correct limit for log utility where the pseudo-inverse value
function slope equals MPC.

Updated all consumption models to use the new function:
- ConsIndShockModel.py
- ConsPrefShockModel.py
- ConsBequestModel.py
- ConsGenIncProcessModel.py
- ConsIndShockModelFast.py (inline conditionals for numba compatibility)
- ConsMarkovModel.py
- ConsRiskyAssetModel.py

Added test case for CRRA=1 in test_PerfForesightConsumerType.py

Fixes #75

Co-authored-by: alanlujan91 <[email protected]>

Author: Copilot

HARK/ConsumptionSaving/ConsBequestModel.py

@@ -50,7 +50,7 @@
     MargValueFuncCRRA,
     ValueFuncCRRA,
 )
-from HARK.rewards import UtilityFuncCRRA, UtilityFuncStoneGeary
+from HARK.rewards import UtilityFuncCRRA, UtilityFuncStoneGeary, vNvrsSlope
 from HARK.utilities import make_assets_grid
 
 
@@ -416,8 +416,8 @@ def calc_vPPnext(S, a, R):
         vNvrsP_temp = vP_temp * uFunc.derinv(v_temp, order=(0, 1))
         mNrm_temp = np.insert(mNrm_temp, 0, mNrmMinNow)
         vNvrs_temp = np.insert(vNvrs_temp, 0, 0.0)
-        vNvrsP_temp = np.insert(vNvrsP_temp, 0, MPCmaxEff ** (-CRRA / (1.0 - CRRA)))
-        MPCminNvrs = MPCminNow ** (-CRRA / (1.0 - CRRA))
+        vNvrsP_temp = np.insert(vNvrsP_temp, 0, vNvrsSlope(MPCmaxEff, CRRA))
+        MPCminNvrs = vNvrsSlope(MPCminNow, CRRA)
         vNvrsFuncNow = CubicInterp(
             mNrm_temp, vNvrs_temp, vNvrsP_temp, MPCminNvrs * hNrmNow, MPCminNvrs
         )

HARK/ConsumptionSaving/ConsGenIncProcessModel.py

@@ -49,6 +49,7 @@
     CRRAutilityP_invP,
     CRRAutilityPP,
     UtilityFuncCRRA,
+    vNvrsSlope,
 )
 from HARK.utilities import make_assets_grid
 
@@ -488,7 +489,7 @@ def calc_vPP_next(S, a, p):
         )
 
         # Add data at the lower bound of p
-        MPCminNvrs = MPCminNow ** (-CRRA / (1.0 - CRRA))
+        MPCminNvrs = vNvrsSlope(MPCminNow, CRRA)
         m_temp = np.reshape(mLvl_temp[:, 0], (aNrmCount + 1, 1))
         mLvl_temp = np.concatenate((m_temp, mLvl_temp), axis=1)
         vNvrs_temp = np.concatenate((MPCminNvrs * m_temp, vNvrs_temp), axis=1)

HARK/ConsumptionSaving/ConsIndShockModel.py

@@ -56,6 +56,7 @@
     CRRAutilityP_invP,
     CRRAutilityPP,
     UtilityFuncCRRA,
+    vNvrsSlope,
 )
 from HARK.utilities import make_assets_grid
 from scipy.optimize import newton
@@ -373,7 +374,7 @@ def solve_one_period_ConsPF(
     # Calculate (pseudo-inverse) value at each consumption kink point
     vNow = uFunc(cNrmNow) + EndOfPrdv
     vNvrsNow = uFunc.inverse(vNow)
-    vNvrsSlopeMin = MPCminNow ** (-CRRA / (1.0 - CRRA))
+    vNvrsSlopeMin = vNvrsSlope(MPCminNow, CRRA)
 
     # Add an additional point to the list of gridpoints for the extrapolation,
     # using the new value of the lower bound of the MPC.
@@ -786,8 +787,8 @@ def solve_one_period_ConsIndShock(
         vNvrsP_temp = vP_temp * uFunc.derinv(v_temp, order=(0, 1))
         mNrm_temp = np.insert(mNrm_temp, 0, mNrmMinNow)
         vNvrs_temp = np.insert(vNvrs_temp, 0, 0.0)
-        vNvrsP_temp = np.insert(vNvrsP_temp, 0, MPCmaxNow ** (-CRRA / (1.0 - CRRA)))
-        MPCminNvrs = MPCminNow ** (-CRRA / (1.0 - CRRA))
+        vNvrsP_temp = np.insert(vNvrsP_temp, 0, vNvrsSlope(MPCmaxNow, CRRA))
+        MPCminNvrs = vNvrsSlope(MPCminNow, CRRA)
         vNvrsFuncNow = CubicInterp(
             mNrm_temp,
             vNvrs_temp,
@@ -1053,8 +1054,8 @@ def solve_one_period_ConsKinkedR(
         vNvrsP_temp = vP_temp * uFunc.derinv(v_temp, order=(0, 1))
         mNrm_temp = np.insert(mNrm_temp, 0, mNrmMinNow)
         vNvrs_temp = np.insert(vNvrs_temp, 0, 0.0)
-        vNvrsP_temp = np.insert(vNvrsP_temp, 0, MPCmaxNow ** (-CRRA / (1.0 - CRRA)))
-        MPCminNvrs = MPCminNow ** (-CRRA / (1.0 - CRRA))
+        vNvrsP_temp = np.insert(vNvrsP_temp, 0, vNvrsSlope(MPCmaxNow, CRRA))
+        MPCminNvrs = vNvrsSlope(MPCminNow, CRRA)
         vNvrsFuncNow = CubicInterp(
             mNrm_temp,
             vNvrs_temp,

HARK/ConsumptionSaving/ConsIndShockModelFast.py

@@ -341,7 +341,11 @@ def _solveConsPerfForesightNumba(
     mNrmMinNow = mNrmNow[0]
 
     # See the PerfForesightConsumerType.ipynb documentation notebook for the derivations
-    vFuncNvrsSlope = MPCmin ** (-CRRA / (1.0 - CRRA))
+    # For CRRA=1 (log utility), the pseudo-inverse value function slope is simply MPC
+    if CRRA == 1.0:
+        vFuncNvrsSlope = MPCmin
+    else:
+        vFuncNvrsSlope = MPCmin ** (-CRRA / (1.0 - CRRA))
 
     return (
         mNrmNow,
@@ -872,8 +876,13 @@ def _add_vFuncNumba(
     vNvrsP = vPnow * utility_invP(vNrmNow, CRRA)
     mNrmGrid = _np_insert(mNrmGrid, 0, mNrmMinNow)
     vNvrs = _np_insert(vNvrs, 0, 0.0)
-    vNvrsP = _np_insert(vNvrsP, 0, MPCmaxEff ** (-CRRA / (1.0 - CRRA)))
-    MPCminNvrs = MPCminNow ** (-CRRA / (1.0 - CRRA))
+    # For CRRA=1 (log utility), the pseudo-inverse value function slope is simply MPC
+    if CRRA == 1.0:
+        vNvrsP = _np_insert(vNvrsP, 0, MPCmaxEff)
+        MPCminNvrs = MPCminNow
+    else:
+        vNvrsP = _np_insert(vNvrsP, 0, MPCmaxEff ** (-CRRA / (1.0 - CRRA)))
+        MPCminNvrs = MPCminNow ** (-CRRA / (1.0 - CRRA))
 
     return (
         mNrmGrid,

HARK/ConsumptionSaving/ConsMarkovModel.py

@@ -39,6 +39,7 @@
     CRRAutilityP_inv,
     CRRAutilityP_invP,
     CRRAutilityPP,
+    vNvrsSlope,
 )
 from HARK.utilities import make_assets_grid
 
@@ -641,9 +642,9 @@ def calc_vPPnext(S, a, R):
             mNrm_temp = np.insert(mNrm_for_vFunc, 0, mNrmMin_i)  # add the lower bound
             vNvrs_now = np.insert(vNvrs_now, 0, 0.0)
             vNvrsP_now = np.insert(
-                vNvrsP_now, 0, MPCmaxEff[i] ** (-CRRA / (1.0 - CRRA))
+                vNvrsP_now, 0, vNvrsSlope(MPCmaxEff[i], CRRA)
             )
-            # MPCminNvrs = MPCminNow[i] ** (-CRRA / (1.0 - CRRA))
+            # MPCminNvrs = vNvrsSlope(MPCminNow[i], CRRA)
             vNvrsFuncNow = LinearInterp(
                 mNrm_temp,
                 vNvrs_now,

HARK/ConsumptionSaving/ConsPrefShockModel.py

@@ -34,7 +34,7 @@
     MargMargValueFuncCRRA,
     ValueFuncCRRA,
 )
-from HARK.rewards import UtilityFuncCRRA
+from HARK.rewards import UtilityFuncCRRA, vNvrsSlope
 
 __all__ = [
     "PrefShockConsumerType",
@@ -376,8 +376,8 @@ def calc_vPPnext(S, a, R):
         vNvrsP_temp = vP_temp * uFunc.derinv(v_temp, order=(0, 1))
         mNrm_temp = np.insert(mNrm_temp, 0, mNrmMinNow)
         vNvrs_temp = np.insert(vNvrs_temp, 0, 0.0)
-        vNvrsP_temp = np.insert(vNvrsP_temp, 0, MPCmaxEff ** (-CRRA / (1.0 - CRRA)))
-        MPCminNvrs = MPCminNow ** (-CRRA / (1.0 - CRRA))
+        vNvrsP_temp = np.insert(vNvrsP_temp, 0, vNvrsSlope(MPCmaxEff, CRRA))
+        MPCminNvrs = vNvrsSlope(MPCminNow, CRRA)
         vNvrsFuncNow = CubicInterp(
             mNrm_temp, vNvrs_temp, vNvrsP_temp, MPCminNvrs * hNrmNow, MPCminNvrs
         )
@@ -683,8 +683,8 @@ def calc_vPPnext(S, a, R):
         vNvrsP_temp = vP_temp * uFunc.derinv(v_temp, order=(0, 1))
         mNrm_temp = np.insert(mNrm_temp, 0, mNrmMinNow)
         vNvrs_temp = np.insert(vNvrs_temp, 0, 0.0)
-        vNvrsP_temp = np.insert(vNvrsP_temp, 0, MPCmaxEff ** (-CRRA / (1.0 - CRRA)))
-        MPCminNvrs = MPCminNow ** (-CRRA / (1.0 - CRRA))
+        vNvrsP_temp = np.insert(vNvrsP_temp, 0, vNvrsSlope(MPCmaxEff, CRRA))
+        MPCminNvrs = vNvrsSlope(MPCminNow, CRRA)
         vNvrsFuncNow = CubicInterp(
             mNrm_temp, vNvrs_temp, vNvrsP_temp, MPCminNvrs * hNrmNow, MPCminNvrs
         )

HARK/ConsumptionSaving/ConsRiskyAssetModel.py

@@ -40,7 +40,7 @@
     MargValueFuncCRRA,
     ValueFuncCRRA,
 )
-from HARK.rewards import UtilityFuncCRRA
+from HARK.rewards import UtilityFuncCRRA, vNvrsSlope
 from HARK.utilities import make_assets_grid
 
 ###############################################################################
@@ -1559,8 +1559,8 @@ def calc_vPPnext(S, a):
         vNvrsP_temp = vP_temp * uFunc.derinv(v_temp, order=(0, 1))
         mNrm_temp = np.insert(mNrm_temp, 0, mNrmMinNow)
         vNvrs_temp = np.insert(vNvrs_temp, 0, 0.0)
-        vNvrsP_temp = np.insert(vNvrsP_temp, 0, MPCmaxEff ** (-CRRA / (1.0 - CRRA)))
-        # MPCminNvrs = MPCminNow ** (-CRRA / (1.0 - CRRA))
+        vNvrsP_temp = np.insert(vNvrsP_temp, 0, vNvrsSlope(MPCmaxEff, CRRA))
+        # MPCminNvrs = vNvrsSlope(MPCminNow, CRRA)
         vNvrsFuncNow = CubicInterp(mNrm_temp, vNvrs_temp, vNvrsP_temp)
         vFuncNow = ValueFuncCRRA(vNvrsFuncNow, CRRA)
     else:

HARK/rewards.py

@@ -234,6 +234,46 @@ def CRRAutilityP_invP(uP, rho):
     return (-1.0 / rho) * uP ** (-1.0 / rho - 1.0)
 
 
+def vNvrsSlope(MPC, rho):
+    """
+    Computes the slope of the pseudo-inverse value function for CRRA utility.
+    
+    For CRRA utility, this is used to determine the asymptotic behavior of the 
+    pseudo-inverse value function vNvrs = u^(-1)(v) as resources become large.
+    
+    For rho ≠ 1: The slope is MPC^(-rho/(1-rho))
+    For rho = 1 (log utility): The slope is simply MPC
+    
+    This function provides the proper limiting behavior for log utility, where
+    the standard formula MPC^(-rho/(1-rho)) is undefined.
+    
+    Parameters
+    ----------
+    MPC : float or np.ndarray
+        Marginal propensity to consume value(s)
+    rho : float
+        Coefficient of relative risk aversion (CRRA)
+        
+    Returns
+    -------
+    float or np.ndarray
+        Slope of the pseudo-inverse value function
+        
+    Notes
+    -----
+    For log utility (rho=1), the derivation is as follows:
+    - u(c) = log(c), so u^(-1)(v) = exp(v)
+    - The pseudo-inverse value function is vNvrs(m) = exp(v(m))
+    - d/dm vNvrs = exp(v) * dv/dm = exp(v) * u'(c) * MPC = exp(log(c)) * (1/c) * MPC = MPC
+    
+    The expression MPC^(-rho/(1-rho)) diverges as rho → 1, but the properly 
+    derived formula for log utility gives MPC directly.
+    """
+    if rho == 1:
+        return MPC
+    return MPC ** (-rho / (1.0 - rho))
+
+
 ###############################################################################
 
 # Define legacy versions of CRRA utility functions with no decorator.

tests/ConsumptionSaving/test_PerfForesightConsumerType.py

@@ -148,3 +148,27 @@ def test_stable_points(self):
         self.assertEqual(
             constrained_agent.solution[0].mNrmStE, constrained_agent.solution[0].mNrmTrg
         )
+
+    def test_CRRA_equals_one(self):
+        """
+        Test that CRRA=1 (log utility) works correctly.
+        
+        This tests fix for issue #75 where CRRA=1 would cause ZeroDivisionError
+        due to expressions like MPC ** (-CRRA / (1.0 - CRRA)).
+        """
+        # Test with CRRA = 1.0 (log utility)
+        log_utility_agent = PerfForesightConsumerType(cycles=0)
+        log_utility_agent.CRRA = 1.0
+        
+        # This should not raise ZeroDivisionError
+        log_utility_agent.solve()
+        
+        # Verify the solution exists and is reasonable
+        self.assertIsNotNone(log_utility_agent.solution[0].cFunc)
+        
+        # Consumption at m=5 should be positive
+        c_at_5 = log_utility_agent.solution[0].cFunc(5.0)
+        self.assertGreater(c_at_5, 0.0)
+        
+        # Consumption should be less than resources (can't consume more than you have)
+        self.assertLess(c_at_5, 5.0)