|
| 1 | +""" |
| 2 | +Stochastic Frontier Analysis (SFA). |
| 3 | +
|
| 4 | +Estimates production/cost frontiers with composed error: |
| 5 | + y_i = x_i'β + v_i - u_i (production) |
| 6 | + y_i = x_i'β + v_i + u_i (cost) |
| 7 | +
|
| 8 | +where v_i ~ N(0, σ_v²) is noise and u_i ≥ 0 is inefficiency. |
| 9 | +
|
| 10 | +Equivalent to Stata's ``frontier`` and R's ``sfa::sfa()``. |
| 11 | +
|
| 12 | +References |
| 13 | +---------- |
| 14 | +Aigner, D., Lovell, C.A.K. & Schmidt, P. (1977). |
| 15 | +"Formulation and Estimation of Stochastic Frontier Production |
| 16 | +Function Models." *Journal of Econometrics*, 6(1), 21-37. |
| 17 | +
|
| 18 | +Battese, G.E. & Coelli, T.J. (1995). |
| 19 | +"A Model for Technical Inefficiency Effects in a Stochastic |
| 20 | +Frontier Production Function for Panel Data." |
| 21 | +*Empirical Economics*, 20(2), 325-332. |
| 22 | +""" |
| 23 | + |
| 24 | +from typing import Optional, List, Dict, Any |
| 25 | +import numpy as np |
| 26 | +import pandas as pd |
| 27 | +from scipy import stats |
| 28 | +from scipy.optimize import minimize |
| 29 | + |
| 30 | +from ..core.results import EconometricResults |
| 31 | + |
| 32 | + |
| 33 | +class FrontierResult(EconometricResults): |
| 34 | + """Extended results for stochastic frontier models.""" |
| 35 | + |
| 36 | + def efficiency(self) -> pd.Series: |
| 37 | + """Technical efficiency estimates E[exp(-u_i) | ε_i].""" |
| 38 | + if 'efficiency' in self.diagnostics: |
| 39 | + return pd.Series(self.diagnostics['efficiency'], name='efficiency') |
| 40 | + return pd.Series(dtype=float) |
| 41 | + |
| 42 | + |
| 43 | +def frontier( |
| 44 | + data: pd.DataFrame = None, |
| 45 | + y: str = None, |
| 46 | + x: List[str] = None, |
| 47 | + dist: str = "half-normal", |
| 48 | + cost: bool = False, |
| 49 | + maxiter: int = 200, |
| 50 | + tol: float = 1e-8, |
| 51 | + alpha: float = 0.05, |
| 52 | +) -> FrontierResult: |
| 53 | + """ |
| 54 | + Stochastic frontier model. |
| 55 | +
|
| 56 | + Estimates a production or cost frontier with composed error. |
| 57 | +
|
| 58 | + Equivalent to Stata's ``frontier y x, dist(hnormal)`` and |
| 59 | + R's ``sfa::sfa()``. |
| 60 | +
|
| 61 | + Parameters |
| 62 | + ---------- |
| 63 | + data : pd.DataFrame |
| 64 | + y : str |
| 65 | + Output (production) or cost variable. |
| 66 | + x : list of str |
| 67 | + Input/cost variables. |
| 68 | + dist : str, default 'half-normal' |
| 69 | + Inefficiency distribution: 'half-normal', 'exponential', 'truncated-normal'. |
| 70 | + cost : bool, default False |
| 71 | + If True, estimate cost frontier (u_i enters positively). |
| 72 | + maxiter : int, default 200 |
| 73 | + alpha : float, default 0.05 |
| 74 | +
|
| 75 | + Returns |
| 76 | + ------- |
| 77 | + FrontierResult |
| 78 | + With .efficiency() method for unit-level TE estimates. |
| 79 | +
|
| 80 | + Examples |
| 81 | + -------- |
| 82 | + >>> import statspai as sp |
| 83 | + >>> result = sp.frontier(df, y='log_output', x=['log_labor', 'log_capital']) |
| 84 | + >>> print(result.summary()) |
| 85 | + >>> eff = result.efficiency() |
| 86 | + """ |
| 87 | + df = data.dropna(subset=[y] + x) |
| 88 | + n = len(df) |
| 89 | + |
| 90 | + y_data = df[y].values.astype(float) |
| 91 | + X_data = np.column_stack([np.ones(n), df[x].values.astype(float)]) |
| 92 | + k = X_data.shape[1] |
| 93 | + var_names = ['_cons'] + list(x) |
| 94 | + |
| 95 | + sign = 1 if cost else -1 # sign of u in composed error |
| 96 | + |
| 97 | + if dist == 'half-normal': |
| 98 | + def neg_log_lik(theta): |
| 99 | + beta = theta[:k] |
| 100 | + ln_sigma_v = theta[k] |
| 101 | + ln_sigma_u = theta[k + 1] |
| 102 | + sigma_v = np.exp(ln_sigma_v) |
| 103 | + sigma_u = np.exp(ln_sigma_u) |
| 104 | + sigma = np.sqrt(sigma_v**2 + sigma_u**2) |
| 105 | + lam = sigma_u / sigma_v |
| 106 | + |
| 107 | + eps = y_data - X_data @ beta |
| 108 | + z = sign * eps * lam / sigma |
| 109 | + |
| 110 | + ll = np.sum( |
| 111 | + -0.5 * np.log(2 * np.pi) - np.log(sigma) |
| 112 | + - 0.5 * (eps / sigma)**2 |
| 113 | + + np.log(2 * stats.norm.cdf(z)) |
| 114 | + ) |
| 115 | + return -ll |
| 116 | + |
| 117 | + n_extra = 2 # ln_sigma_v, ln_sigma_u |
| 118 | + |
| 119 | + elif dist == 'exponential': |
| 120 | + def neg_log_lik(theta): |
| 121 | + beta = theta[:k] |
| 122 | + ln_sigma_v = theta[k] |
| 123 | + ln_sigma_u = theta[k + 1] |
| 124 | + sigma_v = np.exp(ln_sigma_v) |
| 125 | + sigma_u = np.exp(ln_sigma_u) |
| 126 | + |
| 127 | + eps = y_data - X_data @ beta |
| 128 | + mu_star = -sign * eps - sigma_v**2 / sigma_u |
| 129 | + sigma_star = sigma_v |
| 130 | + |
| 131 | + ll = np.sum( |
| 132 | + -np.log(sigma_u) + 0.5 * (sigma_v / sigma_u)**2 |
| 133 | + + sign * eps / sigma_u |
| 134 | + + np.log(np.clip(stats.norm.cdf(mu_star / sigma_star), 1e-20, None)) |
| 135 | + ) |
| 136 | + return -ll |
| 137 | + |
| 138 | + n_extra = 2 |
| 139 | + |
| 140 | + elif dist == 'truncated-normal': |
| 141 | + def neg_log_lik(theta): |
| 142 | + beta = theta[:k] |
| 143 | + ln_sigma_v = theta[k] |
| 144 | + ln_sigma_u = theta[k + 1] |
| 145 | + mu = theta[k + 2] # mean of truncated normal |
| 146 | + sigma_v = np.exp(ln_sigma_v) |
| 147 | + sigma_u = np.exp(ln_sigma_u) |
| 148 | + sigma = np.sqrt(sigma_v**2 + sigma_u**2) |
| 149 | + lam = sigma_u / sigma_v |
| 150 | + |
| 151 | + eps = y_data - X_data @ beta |
| 152 | + mu_star = (sign * eps * sigma_u**2 - mu * sigma_v**2) / sigma**2 * (-1) |
| 153 | + # Actually: mu_i* = (-sign*eps*sigma_u^2 + mu*sigma_v^2) / sigma^2 |
| 154 | + mu_star2 = (-sign * eps * sigma_u**2 + mu * sigma_v**2) / sigma**2 |
| 155 | + sigma_star = sigma_v * sigma_u / sigma |
| 156 | + |
| 157 | + ll = np.sum( |
| 158 | + -0.5 * np.log(2 * np.pi) - np.log(sigma) |
| 159 | + - 0.5 * ((eps + sign * mu) / sigma)**2 |
| 160 | + + np.log(np.clip(stats.norm.cdf(mu_star2 / sigma_star), 1e-20, None)) |
| 161 | + - np.log(np.clip(stats.norm.cdf(mu / sigma_u), 1e-20, None)) |
| 162 | + ) |
| 163 | + return -ll |
| 164 | + |
| 165 | + n_extra = 3 |
| 166 | + else: |
| 167 | + raise ValueError(f"Unknown distribution: {dist}") |
| 168 | + |
| 169 | + # Initialize with OLS |
| 170 | + beta_init = np.linalg.lstsq(X_data, y_data, rcond=None)[0] |
| 171 | + resid = y_data - X_data @ beta_init |
| 172 | + sigma_init = np.std(resid) |
| 173 | + theta0 = np.concatenate([ |
| 174 | + beta_init, |
| 175 | + [np.log(sigma_init * 0.7), np.log(sigma_init * 0.7)], |
| 176 | + ]) |
| 177 | + if dist == 'truncated-normal': |
| 178 | + theta0 = np.concatenate([theta0, [0.0]]) |
| 179 | + |
| 180 | + result = minimize(neg_log_lik, theta0, method='L-BFGS-B', |
| 181 | + options={'maxiter': maxiter, 'ftol': tol}) |
| 182 | + theta_hat = result.x |
| 183 | + |
| 184 | + beta_hat = theta_hat[:k] |
| 185 | + sigma_v = np.exp(theta_hat[k]) |
| 186 | + sigma_u = np.exp(theta_hat[k + 1]) |
| 187 | + sigma = np.sqrt(sigma_v**2 + sigma_u**2) |
| 188 | + lam = sigma_u / sigma_v |
| 189 | + |
| 190 | + # SE via numerical Hessian |
| 191 | + k_total = len(theta_hat) |
| 192 | + eps_h = 1e-5 |
| 193 | + H = np.zeros((k_total, k_total)) |
| 194 | + f0 = neg_log_lik(theta_hat) |
| 195 | + for i in range(k_total): |
| 196 | + ei = np.zeros(k_total) |
| 197 | + ei[i] = eps_h |
| 198 | + fp = neg_log_lik(theta_hat + ei) |
| 199 | + fm = neg_log_lik(theta_hat - ei) |
| 200 | + H[i, i] = (fp - 2 * f0 + fm) / eps_h**2 |
| 201 | + for j in range(i + 1, k_total): |
| 202 | + ej = np.zeros(k_total) |
| 203 | + ej[j] = eps_h |
| 204 | + fpp = neg_log_lik(theta_hat + ei + ej) |
| 205 | + fpm = neg_log_lik(theta_hat + ei - ej) |
| 206 | + fmp = neg_log_lik(theta_hat - ei + ej) |
| 207 | + fmm = neg_log_lik(theta_hat - ei - ej) |
| 208 | + H[i, j] = H[j, i] = (fpp - fpm - fmp + fmm) / (4 * eps_h**2) |
| 209 | + |
| 210 | + try: |
| 211 | + var_cov = np.linalg.inv(H) |
| 212 | + se = np.sqrt(np.abs(np.diag(var_cov))) |
| 213 | + except np.linalg.LinAlgError: |
| 214 | + se = np.full(k_total, np.nan) |
| 215 | + |
| 216 | + # Technical efficiency: E[exp(-u_i) | ε_i] (Jondrow et al. 1982) |
| 217 | + eps_hat = y_data - X_data @ beta_hat |
| 218 | + if dist == 'half-normal': |
| 219 | + mu_star = -sign * eps_hat * sigma_u**2 / sigma**2 |
| 220 | + sigma_star = sigma_v * sigma_u / sigma |
| 221 | + # E[u|ε] = μ* + σ* × φ(μ*/σ*) / Φ(μ*/σ*) |
| 222 | + ratio = mu_star / sigma_star |
| 223 | + E_u = mu_star + sigma_star * stats.norm.pdf(ratio) / np.clip(stats.norm.cdf(ratio), 1e-20, None) |
| 224 | + efficiency = np.exp(-E_u) |
| 225 | + else: |
| 226 | + efficiency = np.full(n, np.nan) |
| 227 | + |
| 228 | + efficiency = np.clip(efficiency, 0, 1) |
| 229 | + |
| 230 | + # Build results |
| 231 | + extra_names = ['ln_sigma_v', 'ln_sigma_u'] |
| 232 | + if dist == 'truncated-normal': |
| 233 | + extra_names.append('mu') |
| 234 | + |
| 235 | + all_names = var_names + extra_names |
| 236 | + params = pd.Series(theta_hat, index=all_names) |
| 237 | + std_errors = pd.Series(se, index=all_names) |
| 238 | + |
| 239 | + ll_val = -neg_log_lik(theta_hat) |
| 240 | + |
| 241 | + return FrontierResult( |
| 242 | + params=params, |
| 243 | + std_errors=std_errors, |
| 244 | + model_info={ |
| 245 | + 'model_type': f"Stochastic Frontier ({'Cost' if cost else 'Production'})", |
| 246 | + 'inefficiency_dist': dist, |
| 247 | + 'sigma_v': sigma_v, |
| 248 | + 'sigma_u': sigma_u, |
| 249 | + 'lambda': lam, |
| 250 | + 'mean_efficiency': float(np.nanmean(efficiency)), |
| 251 | + 'converged': result.success, |
| 252 | + }, |
| 253 | + data_info={ |
| 254 | + 'n_obs': n, |
| 255 | + 'dep_var': y, |
| 256 | + 'df_resid': n - k_total, |
| 257 | + }, |
| 258 | + diagnostics={ |
| 259 | + 'log_likelihood': ll_val, |
| 260 | + 'sigma_v': sigma_v, |
| 261 | + 'sigma_u': sigma_u, |
| 262 | + 'sigma': sigma, |
| 263 | + 'lambda': lam, |
| 264 | + 'aic': -2 * ll_val + 2 * k_total, |
| 265 | + 'bic': -2 * ll_val + np.log(n) * k_total, |
| 266 | + 'efficiency': efficiency, |
| 267 | + 'mean_efficiency': float(np.nanmean(efficiency)), |
| 268 | + }, |
| 269 | + ) |
0 commit comments