Fix CRRA=1.0 division by zero errors across consumption models

HARK/ConsumptionSaving/ConsBequestModel.py

@@ -419,7 +419,19 @@ def calc_vPPnext(S, a, R):
         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))
+        # Handle CRRA=1.0 case to avoid division by zero
+
+        if CRRA == 1.0:
+
+            # When CRRA=1.0, use a small epsilon to avoid division by zero
+
+            CRRA_safe = 1.0 + 1e-8
+
+            MPCminNvrs = MPCminNow ** (-CRRA_safe / (1.0 - CRRA_safe))
+
+        else:
+
+            MPCminNvrs = MPCminNow ** (-CRRA / (1.0 - CRRA))
         vNvrsFuncNow = CubicInterp(
             mNrm_temp, vNvrs_temp, vNvrsP_temp, MPCminNvrs * hNrmNow, MPCminNvrs
         )

HARK/ConsumptionSaving/ConsGenIncProcessModel.py

@@ -488,7 +488,19 @@ def calc_vPP_next(S, a, p):
         )
 
         # Add data at the lower bound of p
-        MPCminNvrs = MPCminNow ** (-CRRA / (1.0 - CRRA))
+        # Handle CRRA=1.0 case to avoid division by zero
+
+        if CRRA == 1.0:
+
+            # When CRRA=1.0, use a small epsilon to avoid division by zero
+
+            CRRA_safe = 1.0 + 1e-8
+
+            MPCminNvrs = MPCminNow ** (-CRRA_safe / (1.0 - CRRA_safe))
+
+        else:
+
+            MPCminNvrs = MPCminNow ** (-CRRA / (1.0 - 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

@@ -373,7 +373,21 @@ 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))
+
+    # Handle CRRA=1.0 case to avoid division by zero
+    if CRRA == 1.0:
+        # When CRRA=1.0, use a small epsilon to avoid division by zero
+        # This is a temporary fix until proper limit calculation is implemented
+        import warnings
+        warnings.warn(
+            "CRRA=1.0 detected. Using CRRA=1.0+1e-8 to avoid division by zero. "
+            "This is a temporary workaround.",
+            UserWarning
+        )
+        CRRA_safe = 1.0 + 1e-8
+        vNvrsSlopeMin = MPCminNow ** (-CRRA_safe / (1.0 - CRRA_safe))
+    else:
+        vNvrsSlopeMin = MPCminNow ** (-CRRA / (1.0 - 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 +800,16 @@ 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))
+
+        # Handle CRRA=1.0 case to avoid division by zero
+        if CRRA == 1.0:
+            # When CRRA=1.0, use a small epsilon to avoid division by zero
+            CRRA_safe = 1.0 + 1e-8
+            vNvrsP_temp = np.insert(vNvrsP_temp, 0, MPCmaxNow ** (-CRRA_safe / (1.0 - CRRA_safe)))
+            MPCminNvrs = MPCminNow ** (-CRRA_safe / (1.0 - CRRA_safe))
+        else:
+            vNvrsP_temp = np.insert(vNvrsP_temp, 0, MPCmaxNow ** (-CRRA / (1.0 - CRRA)))
+            MPCminNvrs = MPCminNow ** (-CRRA / (1.0 - CRRA))
         vNvrsFuncNow = CubicInterp(
             mNrm_temp,
             vNvrs_temp,
@@ -1053,8 +1075,16 @@ 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))
+
+        # Handle CRRA=1.0 case to avoid division by zero
+        if CRRA == 1.0:
+            # When CRRA=1.0, use a small epsilon to avoid division by zero
+            CRRA_safe = 1.0 + 1e-8
+            vNvrsP_temp = np.insert(vNvrsP_temp, 0, MPCmaxNow ** (-CRRA_safe / (1.0 - CRRA_safe)))
+            MPCminNvrs = MPCminNow ** (-CRRA_safe / (1.0 - CRRA_safe))
+        else:
+            vNvrsP_temp = np.insert(vNvrsP_temp, 0, MPCmaxNow ** (-CRRA / (1.0 - CRRA)))
+            MPCminNvrs = MPCminNow ** (-CRRA / (1.0 - CRRA))
         vNvrsFuncNow = CubicInterp(
             mNrm_temp,
             vNvrs_temp,

HARK/ConsumptionSaving/ConsIndShockModelFast.py

@@ -341,7 +341,19 @@ def _solveConsPerfForesightNumba(
     mNrmMinNow = mNrmNow[0]
 
     # See the PerfForesightConsumerType.ipynb documentation notebook for the derivations
-    vFuncNvrsSlope = MPCmin ** (-CRRA / (1.0 - CRRA))
+    # Handle CRRA=1.0 case to avoid division by zero
+
+    if CRRA == 1.0:
+
+        # When CRRA=1.0, use a small epsilon to avoid division by zero
+
+        CRRA_safe = 1.0 + 1e-8
+
+        vFuncNvrsSlope = MPCmin ** (-CRRA_safe / (1.0 - CRRA_safe))
+
+    else:
+
+        vFuncNvrsSlope = MPCmin ** (-CRRA / (1.0 - CRRA))
 
     return (
         mNrmNow,
@@ -869,7 +881,19 @@ def _add_vFuncNumba(
     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))
+    # Handle CRRA=1.0 case to avoid division by zero
+
+    if CRRA == 1.0:
+
+        # When CRRA=1.0, use a small epsilon to avoid division by zero
+
+        CRRA_safe = 1.0 + 1e-8
+
+        MPCminNvrs = MPCminNow ** (-CRRA_safe / (1.0 - CRRA_safe))
+
+    else:
+
+        MPCminNvrs = MPCminNow ** (-CRRA / (1.0 - CRRA))
 
     return (
         mNrmGrid,

HARK/ConsumptionSaving/ConsPrefShockModel.py

@@ -345,9 +345,15 @@ def calc_vPPnext(S, a, R):
         vNvrs_temp = uFunc.inv(v_temp)
         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))
+        vNvrs_temp = np.insert(vNvrs_temp, 0, 0.0)        # Handle CRRA=1.0 case to avoid division by zero
+        if CRRA == 1.0:
+            # When CRRA=1.0, use a small epsilon to avoid division by zero
+            CRRA_safe = 1.0 + 1e-8
+            vNvrsP_temp = np.insert(vNvrsP_temp, 0, MPCmaxEff ** (-CRRA_safe / (1.0 - CRRA_safe)))
+            MPCminNvrs = MPCminNow ** (-CRRA_safe / (1.0 - CRRA_safe))
+        else:
+            vNvrsP_temp = np.insert(vNvrsP_temp, 0, MPCmaxEff ** (-CRRA / (1.0 - CRRA)))
+            MPCminNvrs = MPCminNow ** (-CRRA / (1.0 - CRRA))
         vNvrsFuncNow = CubicInterp(
             mNrm_temp, vNvrs_temp, vNvrsP_temp, MPCminNvrs * hNrmNow, MPCminNvrs
         )
@@ -656,9 +662,15 @@ def calc_vPPnext(S, a, R):
         vNvrs_temp = uFunc.inv(v_temp)
         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))
+        vNvrs_temp = np.insert(vNvrs_temp, 0, 0.0)        # Handle CRRA=1.0 case to avoid division by zero
+        if CRRA == 1.0:
+            # When CRRA=1.0, use a small epsilon to avoid division by zero
+            CRRA_safe = 1.0 + 1e-8
+            vNvrsP_temp = np.insert(vNvrsP_temp, 0, MPCmaxEff ** (-CRRA_safe / (1.0 - CRRA_safe)))
+            MPCminNvrs = MPCminNow ** (-CRRA_safe / (1.0 - CRRA_safe))
+        else:
+            vNvrsP_temp = np.insert(vNvrsP_temp, 0, MPCmaxEff ** (-CRRA / (1.0 - CRRA)))
+            MPCminNvrs = MPCminNow ** (-CRRA / (1.0 - CRRA))
         vNvrsFuncNow = CubicInterp(
             mNrm_temp, vNvrs_temp, vNvrsP_temp, MPCminNvrs * hNrmNow, MPCminNvrs
         )