perf(numerics): replace inv(A)@b with solve(A,b) in IV/SUR/timeseries hot paths

12 sites across regression/{iv,sur,advanced_iv}, inference/jackknife,
timeseries/{cointegration,var,structural_break}: np.linalg.inv(A) @ b ->
np.linalg.solve(A, b) (and inv(Z'Z)@z' -> solve(Z'Z, Z')). Numerically
equivalent (backward-stable, drops an explicit matrix inverse); 924 affected
tests pass unchanged. Parity-safe.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

Author: brycewang-stanford

src/statspai/inference/jackknife.py

@@ -114,7 +114,7 @@ def jackknife_se(
         X_g = X[mask]
         Y_g = Y[mask]
         try:
-            beta_loo[g_idx] = np.linalg.inv(X_g.T @ X_g) @ X_g.T @ Y_g
+            beta_loo[g_idx] = np.linalg.solve(X_g.T @ X_g, X_g.T @ Y_g)
         except np.linalg.LinAlgError:
             beta_loo[g_idx] = np.linalg.lstsq(X_g, Y_g, rcond=None)[0]
 

src/statspai/regression/advanced_iv.py

@@ -115,11 +115,11 @@ def liml(
     k_endog = X_endog.shape[1]
 
     # Residual maker for X_exog
-    Px_exog = X_exog @ np.linalg.inv(X_exog.T @ X_exog) @ X_exog.T
+    Px_exog = X_exog @ np.linalg.solve(X_exog.T @ X_exog, X_exog.T)
     Mx = np.eye(n) - Px_exog
 
     # Projection onto all instruments
-    Pz = Z_all @ np.linalg.inv(Z_all.T @ Z_all) @ Z_all.T
+    Pz = Z_all @ np.linalg.solve(Z_all.T @ Z_all, Z_all.T)
     Mz = np.eye(n) - Pz
 
     # Compute LIML κ via the Anderson (1951) generalized symmetric
@@ -306,7 +306,7 @@ def jive(
     Z_all = np.column_stack([X_exog, Z_excl])
 
     # Leave-one-out projection for each endogenous variable
-    Pz = Z_all @ np.linalg.inv(Z_all.T @ Z_all) @ Z_all.T
+    Pz = Z_all @ np.linalg.solve(Z_all.T @ Z_all, Z_all.T)
     h_ii = np.diag(Pz)
 
     X_endog_hat = np.zeros_like(X_endog)
@@ -454,7 +454,7 @@ def lasso_iv(
     X_exog_mat = np.column_stack([np.ones(n)] + [df[v].values for v in x_exog]) if x_exog else np.ones((n, 1))
 
     # Partial out exogenous regressors from instruments and endogenous vars
-    Px = X_exog_mat @ np.linalg.inv(X_exog_mat.T @ X_exog_mat) @ X_exog_mat.T
+    Px = X_exog_mat @ np.linalg.solve(X_exog_mat.T @ X_exog_mat, X_exog_mat.T)
     Mx = np.eye(n) - Px
 
     Z_tilde = Mx @ Z_candidates  # residualized instruments

src/statspai/regression/iv.py

@@ -195,8 +195,8 @@ def _liml_kappa(
     W_full = np.column_stack([X_exog, Z])  # all instruments
 
     # Projection matrices
-    P_exog = X_exog @ np.linalg.inv(X_exog.T @ X_exog) @ X_exog.T
-    P_full = W_full @ np.linalg.inv(W_full.T @ W_full) @ W_full.T
+    P_exog = X_exog @ np.linalg.solve(X_exog.T @ X_exog, X_exog.T)
+    P_full = W_full @ np.linalg.solve(W_full.T @ W_full, W_full.T)
 
     M_exog = np.eye(n) - P_exog
     M_full = np.eye(n) - P_full

src/statspai/regression/sur.py

@@ -325,7 +325,7 @@ def three_sls(
         instruments = list(set().union(*[set(exog) for _, (_, exog, _) in equations.items()]))
 
     Z = np.column_stack([np.ones(n), df[instruments].values.astype(float)])
-    Pz = Z @ np.linalg.inv(Z.T @ Z) @ Z.T  # projection onto instruments
+    Pz = Z @ np.linalg.solve(Z.T @ Z, Z.T)  # projection onto instruments
 
     # Stage 1: 2SLS for each equation
     Y_list, X_list, k_list, var_names_list = [], [], [], []

src/statspai/timeseries/cointegration.py

@@ -257,7 +257,7 @@ def johansen(
 
     # Concentrate out Z from dY and Y_lag
     if Z.shape[1] > 0:
-        Pz = Z @ np.linalg.inv(Z.T @ Z) @ Z.T
+        Pz = Z @ np.linalg.solve(Z.T @ Z, Z.T)
         Mz = np.eye(T_eff) - Pz
         R0 = Mz @ dY_trim
         R1 = Mz @ Y_lag_trim
@@ -275,7 +275,7 @@ def johansen(
     # |λ S11 - S10 S00^{-1} S01| = 0
     try:
         S00_inv = np.linalg.inv(S00)
-        M = np.linalg.inv(S11) @ S10 @ S00_inv @ S01
+        M = np.linalg.solve(S11, S10 @ S00_inv @ S01)
         eigenvalues, eigenvectors = np.linalg.eig(M)
     except np.linalg.LinAlgError:
         eigenvalues = np.zeros(k)

src/statspai/timeseries/structural_break.py

@@ -293,7 +293,7 @@ def cusum_test(
         bt = np.linalg.lstsq(Xt, yt, rcond=None)[0]
         pred = X_data[t] @ bt
         resid = y_data[t] - pred
-        ft = 1 + X_data[t] @ np.linalg.inv(Xt.T @ Xt) @ X_data[t]
+        ft = 1 + X_data[t] @ np.linalg.solve(Xt.T @ Xt, X_data[t])
         rec_resid.append(resid / np.sqrt(max(ft, 1e-10)))
 
     rec_resid = np.array(rec_resid)

src/statspai/timeseries/var.py

@@ -300,7 +300,7 @@ def granger_causality(
     V = sigma2 * np.eye(len(coefs_all))  # simplified
     mid = R @ V @ R.T
     try:
-        F_stat = (r @ np.linalg.inv(mid) @ r) / len(restrict_indices)
+        F_stat = (r @ np.linalg.solve(mid, r)) / len(restrict_indices)
     except np.linalg.LinAlgError:
         F_stat = np.nan