diff --git a/ogcore/firm.py b/ogcore/firm.py
index 5d454f679..dfddeb4a3 100644
--- a/ogcore/firm.py
+++ b/ogcore/firm.py
@@ -8,6 +8,10 @@
 ------------------------------------------------------------------------
 """
 
+# Floor for quantities raised to potentially negative powers in CES.
+# Prevents 0^(negative) → inf/NaN when epsilon != 1.
+_FLOOR = 1e-12
+
 """
 ------------------------------------------------------------------------
     Functions
@@ -60,18 +64,21 @@ def get_Y(K, K_g, L, p, method, m=-1):
                     * (L ** (1 - gamma - gamma_g))
                 )
             else:
+                K_f = np.maximum(K, _FLOOR)
+                K_g_f = np.maximum(K_g, _FLOOR)
+                L_f = np.maximum(L, _FLOOR)
                 Y = Z * (
                     (
                         (gamma ** (1 / epsilon))
-                        * (K ** ((epsilon - 1) / epsilon))
+                        * (K_f ** ((epsilon - 1) / epsilon))
                     )
                     + (
                         (gamma_g ** (1 / epsilon))
-                        * (K_g ** ((epsilon - 1) / epsilon))
+                        * (K_g_f ** ((epsilon - 1) / epsilon))
                     )
                     + (
                         ((1 - gamma - gamma_g) ** (1 / epsilon))
-                        * (L ** ((epsilon - 1) / epsilon))
+                        * (L_f ** ((epsilon - 1) / epsilon))
                     )
                 ) ** (epsilon / (epsilon - 1))
         else:
@@ -85,15 +92,18 @@ def get_Y(K, K_g, L, p, method, m=-1):
             gamma = p.gamma
             epsilon = p.epsilon
             Z = p.Z[-1, :]
+            K_f = np.maximum(K, _FLOOR)
+            K_g_f = np.maximum(K_g, _FLOOR)
+            L_f = np.maximum(L, _FLOOR)
             Y = Z * (
-                ((gamma ** (1 / epsilon)) * (K ** ((epsilon - 1) / epsilon)))
+                ((gamma ** (1 / epsilon)) * (K_f ** ((epsilon - 1) / epsilon)))
                 + (
                     (gamma_g ** (1 / epsilon))
-                    * (K_g ** ((epsilon - 1) / epsilon))
+                    * (K_g_f ** ((epsilon - 1) / epsilon))
                 )
                 + (
                     ((1 - gamma - gamma_g) ** (1 / epsilon))
-                    * (L ** ((epsilon - 1) / epsilon))
+                    * (L_f ** ((epsilon - 1) / epsilon))
                 )
             ) ** (epsilon / (epsilon - 1))
             Y2 = Z * (K**gamma) * (K_g**gamma_g) * (L ** (1 - gamma - gamma_g))
@@ -117,18 +127,21 @@ def get_Y(K, K_g, L, p, method, m=-1):
                     * (L ** (1 - gamma - gamma_g))
                 )
             else:
+                K_f = np.maximum(K, _FLOOR)
+                K_g_f = np.maximum(K_g, _FLOOR)
+                L_f = np.maximum(L, _FLOOR)
                 Y = Z * (
                     (
                         (gamma ** (1 / epsilon))
-                        * (K ** ((epsilon - 1) / epsilon))
+                        * (K_f ** ((epsilon - 1) / epsilon))
                     )
                     + (
                         (gamma_g ** (1 / epsilon))
-                        * (K_g ** ((epsilon - 1) / epsilon))
+                        * (K_g_f ** ((epsilon - 1) / epsilon))
                     )
                     + (
                         ((1 - gamma - gamma_g) ** (1 / epsilon))
-                        * (L ** ((epsilon - 1) / epsilon))
+                        * (L_f ** ((epsilon - 1) / epsilon))
                     )
                 ) ** (epsilon / (epsilon - 1))
         else:
@@ -141,15 +154,18 @@ def get_Y(K, K_g, L, p, method, m=-1):
             gamma = p.gamma
             epsilon = p.epsilon
             Z = p.Z[: p.T, :]
+            K_f = np.maximum(K, _FLOOR)
+            K_g_f = np.maximum(K_g, _FLOOR)
+            L_f = np.maximum(L, _FLOOR)
             Y = Z * (
-                ((gamma ** (1 / epsilon)) * (K ** ((epsilon - 1) / epsilon)))
+                ((gamma ** (1 / epsilon)) * (K_f ** ((epsilon - 1) / epsilon)))
                 + (
                     (gamma_g ** (1 / epsilon))
-                    * (K_g ** ((epsilon - 1) / epsilon))
+                    * (K_g_f ** ((epsilon - 1) / epsilon))
                 )
                 + (
                     ((1 - gamma - gamma_g) ** (1 / epsilon))
-                    * (L ** ((epsilon - 1) / epsilon))
+                    * (L_f ** ((epsilon - 1) / epsilon))
                 )
             ) ** (epsilon / (epsilon - 1))
             Y2 = Z * (K**gamma) * (K_g**gamma_g) * (L ** (1 - gamma - gamma_g))
@@ -273,9 +289,10 @@ def get_KLratio_KLonly(r, p, method, m=-1):
         # General CES case
         cost_of_capital = get_cost_of_capital(r, p, method, m)
         bracket = cost_of_capital * (Z * (gamma ** (1 / epsilon))) ** -1
+        denom = (bracket ** (epsilon - 1)) - (gamma ** (1 / epsilon))
+        denom = np.maximum(denom, _FLOOR)
         KLratio = (
-            ((1 - gamma) ** (1 / epsilon))
-            / ((bracket ** (epsilon - 1)) - (gamma ** (1 / epsilon)))
+            ((1 - gamma) ** (1 / epsilon)) / denom
         ) ** (epsilon / (epsilon - 1))
 
     return KLratio
@@ -337,12 +354,13 @@ def get_MPx(Y, x, share, p, method, m=-1):
         Z = p.Z[: p.T, m].reshape(p.T, 1)
         Y = Y[: p.T].reshape(p.T, 1)
         x = x[: p.T].reshape(p.T, 1)
-    if np.any(x) == 0:
+    if np.any(x == 0):
         MPx = np.zeros_like(Y)
     else:
-        MPx = Z ** ((p.epsilon[m] - 1) / p.epsilon[m]) * ((share * Y) / x) ** (
-            1 / p.epsilon[m]
-        )
+        x_f = np.maximum(x, _FLOOR)
+        MPx = Z ** ((p.epsilon[m] - 1) / p.epsilon[m]) * (
+            (share * Y) / x_f
+        ) ** (1 / p.epsilon[m])
 
     return MPx
 
@@ -656,13 +674,17 @@ def solve_L(Y, K, K_g, p, method, m=-1):
     if epsilon == 1.0:
         L = (Y / (Z * K**gamma * K_g**gamma_g)) ** (1 / (1 - gamma - gamma_g))
     else:
+        Y_f = np.maximum(Y, _FLOOR)
+        K_f = np.maximum(K, _FLOOR)
+        K_g_f = np.maximum(K_g, _FLOOR)
+        numer = (
+            (Y_f / Z) ** ((epsilon - 1) / epsilon)
+            - gamma ** (1 / epsilon) * K_f ** ((epsilon - 1) / epsilon)
+            - gamma_g ** (1 / epsilon) * K_g_f ** ((epsilon - 1) / epsilon)
+        )
+        numer = np.maximum(numer, _FLOOR)
         L = (
-            (
-                (Y / Z) ** ((epsilon - 1) / epsilon)
-                - gamma ** (1 / epsilon) * K ** ((epsilon - 1) / epsilon)
-                - gamma_g ** (1 / epsilon) * K_g ** ((epsilon - 1) / epsilon)
-            )
-            / ((1 - gamma - gamma_g) ** (1 / epsilon))
+            numer / ((1 - gamma - gamma_g) ** (1 / epsilon))
         ) ** (epsilon / (epsilon - 1))
 
     return L