From 00ff6fcd78d9e378880673aa4250234f1005aff6 Mon Sep 17 00:00:00 2001 From: Vahid Ahmadi Date: Tue, 31 Mar 2026 13:34:04 +0100 Subject: [PATCH] Add numerical guards for CES production functions with heterogeneous epsilon MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit When epsilon differs from 1.0 across sectors, the CES production function raises K, L, and K_g to (epsilon-1)/epsilon, which is negative for epsilon < 1. If any input is zero or near-zero (common during initial solver iterations), this produces inf or NaN, causing the SS and TPI solvers to fail. Changes: - Add _FLOOR (1e-12) constant and apply np.maximum() guards before all CES power operations in get_Y, get_MPx, get_KLratio_KLonly, and solve_L - Fix bug in get_MPx: np.any(x) == 0 should be np.any(x == 0) - Guard denominator in get_KLratio_KLonly CES branch against negative values These guards only affect the general CES code paths (epsilon != 1); the Cobb-Douglas special case (epsilon == 1) is unchanged. Tested with OG-UK 8-sector model using heterogeneous epsilon values (0.40-1.30) — SS solver now converges successfully. Co-Authored-By: Claude Opus 4.6 --- ogcore/firm.py | 70 +++++++++++++++++++++++++++++++++----------------- 1 file changed, 46 insertions(+), 24 deletions(-) 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