# crabbymetrics

> crabbymetrics is a compact Rust-backed Python econometrics package with a NumPy-facing, scikit-adjacent API. It provides linear, IV, panel causal, regularized, likelihood, moment/GMM, semiparametric, balancing-weight, transform, and optimizer utilities while keeping the runtime Python dependency footprint intentionally small: NumPy only.

This file is for coding agents and LLMs that need to use `crabbymetrics` correctly in applications. Prefer this file when you need a dense API map, data-shape conventions, estimator selection guidance, and links into the human documentation site.

Canonical project locations:

- Documentation site: https://apoorvalal.github.io/crabbymetrics/
- API reference: https://apoorvalal.github.io/crabbymetrics/api.html
- GitHub repository: https://github.com/apoorvalal/crabbymetrics
- PyPI package: https://pypi.org/project/crabbymetrics/

## Install and import

Install from PyPI:

```bash
uv pip install crabbymetrics
# or
pip install crabbymetrics
```

Typical import:

```python
import numpy as np
import crabbymetrics as cm
```

The package is built with PyO3/maturin and ships native wheels. Python package versioning is sourced from `Cargo.toml` in the repository. Runtime dependency: `numpy`.

## Global conventions

- Inputs are NumPy arrays or array-like objects that can be converted to NumPy arrays.
- Most feature matrices are 2D, shape `(n_observations, n_features)`.
- Most regression/classification targets are 1D, shape `(n_observations,)`.
- Estimators mutate in place: instantiate, call `fit(...)`, then call `predict(...)`, `summary(...)`, `bootstrap(...)`, or estimator-specific methods when available.
- `summary()` returns a plain Python `dict`, not a rich result object. Downstream code should index by keys.
- `bootstrap(n_bootstrap, seed=None)` returns NumPy arrays or dicts depending on estimator; if an estimator has an intercept, bootstrap draws place the intercept first.
- OLS, Ridge, ElasticNet, Logit, MultinomialLogit, Poisson, and TwoSLS fit an intercept internally.
- `FixedEffectsOLS`, `EPLM`, `AverageDerivative`, `PartiallyLinearDML`, `AIPW`, `GMM`, and `MEstimator` are estimation-first and intentionally do not expose `predict()`.
- Weighted fits exist for OLS, Ridge, FixedEffectsOLS, and TwoSLS.
- Robust covariance conventions vary by estimator; the main linear estimators share `summary(vcov="vanilla" | "hc1" | "newey_west" | "cluster", lags=None, clusters=None)`.
- Iterative likelihood and callback estimators expose `converged`, `iterations`, `termination_reason`, and `objective` in `summary()`. They reject budget-exhausted fits unless the estimator explicitly documents retention of a diagnostic best iterate, as `MatrixCompletion` does.
- The package does not require pandas, scikit-learn, scipy, statsmodels, or polars at runtime. If you start from a DataFrame, convert to NumPy explicitly before fitting.

## Implementation-detail references

Every estimator, transform, and optimizer class page under `docs/reference/` contains an implementation walkthrough in addition to the estimand, inference, and performance sections. These walkthroughs are source-audited descriptions of parameter layout, initialization, linear algebra, callback flow, stopping rules, prediction reconstruction, randomization, and retained fitted state.

Important ownership boundaries:

- `Logit`, `MultinomialLogit`, and `Poisson` own their likelihoods and analytic derivatives. Binary logit uses stable-softplus L-BFGS rather than IRLS; multinomial logit uses a full-class softmax parameterization; Poisson uses exact-Hessian Newton-CG.
- OLS, Ridge, TwoSLS, ABCOLS, BaggedPolynomialRegressor, the semiparametric estimators, GMM, balancing weights, survival models, and panel estimators own their estimator-specific numerical construction in Rust. Argmin or dense linear-algebra crates provide generic solver/factorization primitives where documented.
- `FixedEffectsOLS` delegates absorption to the Rust `within` crate but owns batch orchestration, the slope solve, absorbed degrees of freedom, covariance, and the Python contract.
- `ElasticNet` delegates coordinate descent to `linfa-elasticnet`; crabbymetrics owns input validation, the independently recomputed duality-gap convergence gate, diagnostics, and bootstrap policy.
- Exact `PCA` and exact training-Gram construction for `KernelBasis` delegate narrow primitives to Linfa. RandomizedPCA, NystromBasis, RandomFourierFeatures, out-of-sample cross-kernel construction, and all wrapper diagnostics are package-owned.
- There is no monolithic `KernelRidge` class. Exact or approximate kernel ridge is the explicit composition `KernelBasis` / `NystromBasis` / `RandomFourierFeatures` followed by `Ridge`; the two class pages document both halves of that pipeline.

Do not run `docs/_generate_class_reference.py --force` as a routine maintenance step. It creates first-pass scaffolds and would replace the audited hand-written sections. Without `--force`, it skips existing pages.

## Quick estimator selection

Use these defaults unless the application has a reason to do otherwise:

- Plain linear regression with robust standard errors: `cm.OLS()`.
- Anytime-valid confidence sequences for fitted OLS coefficients: `model.summary(vcov="vanilla", anytime_valid=True, g=cm.optimal_g(n, k, alpha))`.
- Linear regression with abundance-based constraints / weighted effect coding for categorical main effects and categorical modifiers: `cm.ABCOLS()`.
- Penalized linear regression / shrinkage / cross-validated ridge: `cm.Ridge(penalty=None_or_grid, cv=5)`.
- Nonlinear prediction with a bagged random-subspace polynomial ridge ensemble: `cm.BaggedPolynomialRegressor(...)`.
- Multi-way categorical fixed effects with slopes: `cm.FixedEffectsOLS()`.
- Instrumental variables / 2SLS: `cm.TwoSLS()`.
- Single treated time series with donor paths and simplex donor weights: `cm.SyntheticControl()`.
- Balanced panel causal ATT with absorbing treatment matrix: start with `cm.SyntheticDID()` or `cm.HorizontalPanelRidge()`; use `cm.MatrixCompletion()` when a low-rank completion model is appropriate.
- No-covariate interactive fixed-effects factor model on a balanced panel: `cm.InteractiveFixedEffects(rank=..., force=...)`.
- Binary treatment ATE with nuisance models: `cm.AIPW()`.
- Continuous treatment partially linear model: `cm.PartiallyLinearDML()` or `cm.EPLM()`.
- Continuous-treatment average derivative: `cm.AverageDerivative(method="dr")`.
- Covariate balancing/calibration weights: `cm.BalancingWeights()`.
- Binary classification: use `cm.Logit()` for a convergence-checked likelihood estimator.
- Multiclass classification: `cm.MultinomialLogit()`.
- Count regression: `cm.Poisson()`.
- User-supplied moment conditions: `cm.GMM(moment_fn, jacobian_fn=None, ...)`.
- User-supplied objective and scores: `cm.MEstimator(objective_fn, score_fn, ...)`.
- Feature transformations: `cm.PCA(...)`, `cm.RandomizedPCA(...)`, `cm.KernelBasis(...)`, `cm.NystromBasis(...)`, and `cm.RandomFourierFeatures(...)`.
- Regularized estimators expose point estimates for prediction but do not claim analytic standard errors unless the fitted objective and covariance match.
- Direct numerical optimization with Python callbacks: `cm.Optimizers.minimize_*`.

## Panel causal API: most important recent contract

The paved panel causal API is matrix-first:

```python
Y = np.asarray(Y, dtype=float)  # shape (n_units, n_periods)
W = np.asarray(W, dtype=float)  # same shape; binary 0/1 absorbing treatment indicator

model = cm.SyntheticDID()
model.fit(Y, W)
out = model.summary()
att = out["att"]
counterfactual = out["counterfactual"]
treatment_effect = out["treatment_effect"]
event_study = out["event_study"]
group_means = out["group_means"]
```

`Y` is the observed balanced outcome panel. `W` is a same-shaped binary absorbing treatment matrix where `W[i, t] == 1` means unit `i` is treated at time `t`. Never-treated donor units have all-zero rows in `W`. Do not pass long/dataframe inputs to the panel causal estimators.

First-class panel causal estimators using `fit(Y, W)`:

- `HorizontalPanelRidge(penalty=1.0)`: horizontal forecasting estimator. Uses never-treated donor paths as contemporaneous features, trains cohort-specific ridge forecasts on pre-treatment periods, and predicts treated-unit counterfactual paths. Works surprisingly well as a pure forecasting baseline.
- `SyntheticDID(zeta_omega=None, zeta_lambda=None, max_iterations=1000)`: synthetic difference-in-differences with cohort-specific donor and time weights under the common panel contract.
- `MatrixCompletion(lambda_l=None, lambda_fraction=0.25, fit_unit_effects=True, fit_time_effects=True, max_iterations=500, effect_iterations=2, tolerance=1e-6)`: matrix completion for panels where untreated cells (`W == 0`) are observed entries and treated cells are completed as counterfactuals.

Common panel summary keys:

- `att`: scalar average treatment effect on treated cells.
- `counterfactual`: matrix of counterfactual untreated outcomes for treated units/cells in the estimator's natural output shape.
- `treatment_effect`: observed-minus-counterfactual treatment effect matrix.
- `event_study`: event-time summaries, including unweighted and treated-count-weighted aggregates with confidence-interval-style columns where available.
- `group_means`: cohort/event-time rows plus nested unweighted and weighted event-time aggregates for plotting.

Panel caution:

- `HorizontalPanelRidge` and `SyntheticDID` currently assume at least one never-treated donor unit and use cohort-specific donor-based models for staggered timing.
- Treatment should be absorbing. If a unit is treated at time `t`, later periods for that unit should also be treated.
- The panel estimators infer treated units, first-treatment cohorts, never-treated donors, pre/post windows, event time, and aggregate summaries internally from `W`.

## Core estimator signatures

These signatures are the public Python surface exposed by the current package.

### Linear and IV estimators

```python
cm.OLS()
fit(x, y)
fit_weighted(x, y, sample_weight)
predict(x)
summary(vcov="hc1", lags=None, clusters=None)
bootstrap(n_bootstrap, seed=None)
```

`summary()` keys include `intercept`, `coef`, `intercept_se`, and `coef_se`. Use `vcov="hc1"` for heteroskedasticity-robust inference, `vcov="newey_west"` with `lags`, and `vcov="cluster"` with integer cluster labels.

```python
cm.ABCOLS()
fit(y, x, categories, cont_cat_interactions=None, cat_cat_interactions=None, center_continuous=True)
predict(x, categories)
summary()
fitted_values()
residuals()
design_matrix()
constraint_matrix()
column_names()
```

`ABCOLS` is an OLS reparameterization for categorical modifiers using abundance-based constraints / weighted effect coding. `x` is a 2D float matrix of continuous covariates; `categories` is a 2D `uint32`-like matrix of zero-based contiguous observed category codes. Categorical main effects use one coefficient per level with empirical-abundance-weighted mean zero. `cont_cat_interactions=[(j, k)]` adds and constrains `x_j:C_k` so the continuous main slope is the sample-abundance-weighted average group-specific slope. `cat_cat_interactions=[(a, b)]` adds all `C_a:C_b` cells with weighted zero-margin constraints. `summary()` keys include `coef`, `se`, `column_names`, `constraint_names`, `sigma2`, `df_resid`, `rank`, `max_constraint_violation`, `continuous_means`, and `n_levels`. Current inference is homoskedastic OLS in the constrained parameterization.

```python
cm.Ridge(penalty=None, cv=5)
fit(x, y)
fit_weighted(x, y, sample_weight)
predict(x)
summary(vcov="hc1", lags=None, clusters=None)
bootstrap(n_bootstrap, seed=None)
```

`penalty` may be a scalar or a grid. With a grid, Ridge cross-validates, stores coefficient paths, and exposes `penalty`, `penalties`, `best_penalty_index`, `cv_mse`, `intercept_path`, and `coef_path` in addition to coefficient summaries.

```python
cm.FixedEffectsOLS()
fit(x, fe, y)
fit_weighted(x, fe, y, sample_weight)
summary(vcov="hc1", lags=None, clusters=None)
bootstrap(n_bootstrap, seed=None)
```

`fe` is a 2D `uint32`-like matrix of zero-based category codes, shape `(n_observations, n_fixed_effect_dimensions)`. The estimator partials out fixed effects and estimates slopes without an intercept. `summary()` reports absorbed and residual degrees of freedom; one- and two-way absorbed rank is exact, while three-or-more-way rank is explicitly labeled conservative.

```python
cm.TwoSLS()
fit(x_endog, x_exog, z, y)
fit_weighted(x_endog, x_exog, z, y, sample_weight)
predict(x)
summary(vcov="hc1", lags=None, clusters=None)
bootstrap(n_bootstrap, seed=None)
```

`x_endog` is the endogenous regressor matrix, `x_exog` is the included exogenous regressor matrix, `z` is the excluded instrument matrix, and `y` is the outcome. Matrix-valued endogenous regressors and instruments are supported.

### Regularized and likelihood estimators

```python
cm.ElasticNet(penalty=1.0, l1_ratio=0.5, tolerance=1e-4, max_iterations=1000)
fit(x, y)
predict(x)
summary()
bootstrap(n_bootstrap, seed=None)
```

`l1_ratio=0.0` is ridge-style shrinkage; `l1_ratio=1.0` is lasso-style shrinkage. The implementation solves against `y - mean(y)` without centering `X`, then predicts `mean(y) + X @ coef`; this differs from jointly optimizing an unpenalized intercept unless `X` is column-centered. `summary()` exposes the final duality gap and fit diagnostics and marks analytic inference unavailable. A nonconverged fit raises.

```python
cm.Logit(alpha=0.0, max_iterations=100, gradient_tolerance=1e-4)
fit(x, y)
predict(x)
summary()
bootstrap(n_bootstrap, seed=None)
```

Binary labels must be coded 0/1. The native stable-softplus L-BFGS fit minimizes unnormalized Bernoulli negative log-likelihood plus L2 slope penalty and rejects nonconverged results. `predict()` returns probabilities for label 1 and `predict_label()` returns labels. Fisher inference is available only for the default unpenalized fit; `alpha>0` is prediction-only.

```python
cm.MultinomialLogit(alpha=0.0, max_iterations=100, gradient_tolerance=1e-4)
fit(x, y)
predict(x)
summary()
bootstrap(n_bootstrap, seed=None)
```

Multiclass labels should be integer-coded. The native stable-log-sum-exp L-BFGS fit minimizes full-class softmax negative log-likelihood plus L2 slope penalty and rejects nonconverged results. `predict()` returns class probabilities, `predict_lin()` returns logits, and `predict_label()` returns labels. The default unpenalized summary reports identified class-versus-last-class contrasts; `alpha>0` is prediction-only.

```python
cm.Poisson(alpha=0.0, max_iterations=100, tolerance=1e-4)
fit(x, y)
predict(x)
summary(vcov="vanilla")
bootstrap(n_bootstrap, seed=None)
```

Use Poisson for nonnegative count-like outcomes. Newton-CG minimizes unnormalized Poisson negative log-likelihood plus L2 slope penalty and rejects nonconverged results. Model-based and QMLE/sandwich inference are available only at `alpha=0`; penalized fits are prediction-only.

### Survival and recurrent-event estimators

```python
cm.ExponentialPH()
cm.WeibullPH()
fit(x, time, event)
```

The parametric PH classes maximize their right-censored likelihoods with analytic Newton steps. They expose absolute hazard, cumulative hazard, and survival predictions. `summary()` includes inverse-observed-information covariance and common fit diagnostics.

```python
cm.CoxPH()
fit(x, time, event)

cm.AndersenGill()
fit(x, start, stop, event)
```

The semiparametric classes maximize the Breslow-ties partial likelihood and predict relative risk; they do not estimate a baseline hazard. Andersen-Gill supports interval risk sets but has no subject identifier or subject-clustered robust covariance. All four survival estimators raise on iteration-budget exhaustion.

### Panel estimators and panel helpers

```python
cm.HorizontalPanelRidge(penalty=1.0)
fit(Y, W)
predict()
treatment_effect()
summary()
```

Important summary keys: `att`, `counterfactual`, `treatment_effect`, `event_study`, `group_means`, `cohort_intercepts`, `cohort_coef`, `pre_rmse`, `penalty`, `control_units`, `treated_units`, `cohorts`.

```python
cm.SyntheticControl(max_iterations=500)
fit(donors, treated)
predict(donors)
summary()
bootstrap(n_bootstrap, seed=None)
```

Low-level simplex donor-weight API for one treated path. `donors` is usually shaped `(n_periods, n_donors)` and `treated` is the treated pre-period vector. Summary keys include `weights` and `pre_rmse`.

```python
cm.SyntheticDID(zeta_omega=None, zeta_lambda=None, max_iterations=1000)
fit(Y, W)
predict()
treatment_effect()
summary()
```

Important summary keys: `att`, `unit_weights`, `time_weights`, `counterfactual`, `synthetic_outcome`, `treatment_effect`, `event_study`, `group_means`, `pre_rmse`, `unit_intercept`, `time_intercept`, `zeta_omega`, `zeta_lambda`, `control_units`, `treated_units`, `cohorts`.

```python
cm.MatrixCompletion(lambda_l=None, lambda_fraction=0.25, fit_unit_effects=True, fit_time_effects=True, max_iterations=500, effect_iterations=2, tolerance=1e-6)
fit(Y, W)
predict()
summary()
```

Important summary keys: `att`, `completed`, `counterfactual`, `treatment_effect`, `event_study`, `group_means`, `low_rank`, `unit_effects`, `time_effects`, `singular_values`, `lambda_l`, `converged`, `iterations`, `termination_reason`, `objective`, `history_objective`, `history_rmse`. Budget exhaustion retains the final iterate with `converged=False`.

```python
cm.InteractiveFixedEffects(rank=0, force=3)
fit(y)
predict()
summary()
```

No-covariate fect-style principal-components factor model helper for balanced panels. Lower-level helpers `cm.panel_factor(...)` and `cm.panel_fe(...)` are also exported.

### Balancing and semiparametric estimators

```python
cm.BalancingWeights(objective="quadratic", solver="auto", autoscale=False, min_weight=0.0, max_weight=1.0, l2_norm=0.0, max_iterations=200, tolerance=1e-8)
fit(covariates, target_covariates, baseline_weights=None, target_weights=None)
get_weights()
summary()
```

Low-level calibration/reweighting API. `solver_converged` and the generic convergence fields describe the scaled calibration system; `success` additionally checks normalized weight feasibility. `scaled_residual_norm` and `solver_objective` are numerical diagnostics, while `original_balance_l2` and `original_balance_max_abs` remain in user-supplied units.

```python
cm.EPLM(fd_eps=1e-6)
fit(y, d, w)
summary(vcov=None, lags=None, clusters=None)
```

Robins-Newey partially linear E-estimator for scalar continuous treatment `d` and controls `w`.

```python
cm.AverageDerivative(method="dr", fd_eps=1e-6)
fit(y, d, w)
summary(vcov=None, lags=None, clusters=None)
```

Average-derivative estimator with `method="ob"`, `"ipw"`, or `"dr"` under a scalar continuous-treatment working model.

```python
cm.PartiallyLinearDML(penalty=None, cv=5, n_folds=5, seed=42)
fit(y, d, x)
summary(vcov=None, lags=None, clusters=None)
```

Cross-fit ridge nuisance implementation of the partially linear Double ML score for a scalar continuous treatment. Folds use a deterministic seeded shuffle.

```python
cm.AIPW(penalty=None, cv=5, n_folds=5, propensity_clip=0.02, seed=42)
fit(y, d, x)
summary(vcov=None, lags=None, clusters=None)
```

Binary-treatment augmented inverse probability weighting estimator for the ATE. Treatment `d` should be coded `0/1`; fold assignment is deterministically shuffled within treatment strata.

### Bagged polynomial regression

```python
cm.BaggedPolynomialRegressor(n_estimators=100, degree=2, max_features=None, max_samples=None, bootstrap=True, penalty=1.0, seed=None)
fit(x, y)
predict(x)
summary()
```

Prediction-only ensemble with guarded polynomial expansion, internal column scaling, row bagging, random feature subspaces, and out-of-bag diagnostics when bootstrap sampling is enabled.

### Moment and callback-driven estimators

```python
cm.GMM(moment_fn, jacobian_fn=None, max_iterations=100, tolerance=1e-6, ridge=1e-8, fd_eps=1e-6)
fit(data, theta0, weighting="auto")
summary(vcov="sandwich", omega="iid", lags=None, clusters=None)
```

`moment_fn(theta, data)` should return a 2D array of per-observation moments. `jacobian_fn` is optional; if omitted, finite differences are used. Use this for custom just-identified or overidentified moment systems.

```python
cm.MEstimator(objective_fn, score_fn, max_iterations=100, tolerance=1e-6, derivative_step=1e-6)
fit(data, theta0)
compute_vcov()
summary()
bootstrap(n_bootstrap, seed=None)
```

Lowest-level callback M-estimation interface. Python supplies the objective and per-observation scores; Rust checks optimization convergence, rejects budget exhaustion, and estimates $A^{-1} B A^{-T}/n$ from a numerical mean-score Jacobian and empirical score meat.

### Transforms

```python
cm.PCA(n_components, whiten=False)
fit(x)
transform(x)
fit_transform(x)
inverse_transform(scores)
summary()
```

PCA transformer for design-matrix preprocessing. Explained variance uses $s_j^2/(n-1)$ and ratios use the full centered sample variance. Whitening is undone during inverse transformation.

```python
cm.RandomizedPCA(n_components, oversamples=10, power_iter=1, seed=None)
fit(x)
transform(x)
fit_transform(x)
inverse_transform(scores)
summary()
```

Randomized SVD PCA transformer for lower-rank approximate decomposition.

```python
cm.KernelBasis(kernel="gaussian", bandwidth=0.5, coef0=1.0, degree=2.0)
fit(x)
transform(x)
fit_transform(x)
summary()
```

Exact kernel feature basis transformer. Its Gaussian kernel is `exp(-||x-z||^2 / bandwidth)`.

```python
cm.NystromBasis(n_components, kernel="gaussian", bandwidth=0.5, coef0=1.0, degree=2.0, ridge=1e-10, seed=None)
cm.RandomFourierFeatures(n_components, bandwidth=0.5, seed=None)
```

Both expose `fit`, `transform`, `fit_transform`, and `summary`. Nystrom uses sampled landmark rows and a regularized kernel inverse square root. RandomFourierFeatures approximates the same Gaussian kernel with frequencies drawn from $N(0,2I/\text{bandwidth})$.

### Optimizers

`cm.Optimizers` exposes scipy-style static methods returning dictionaries with fields such as `x`, `fun`, `nit`, `success`, `message`, and `method`.

```python
cm.Optimizers.minimize_lbfgs(fun, x0, grad, max_iterations=100, tolerance=1e-6)
cm.Optimizers.minimize_bfgs(fun, x0, grad, max_iterations=100, tolerance=1e-6)
cm.Optimizers.minimize_nonlinear_cg(fun, x0, grad, max_iterations=100, restart_iters=10, restart_orthogonality=0.1, tolerance=1e-6)
cm.Optimizers.minimize_gauss_newton_ls(residual_fn, x0, jacobian_fn, max_iterations=100, tolerance=1e-6)
cm.Optimizers.minimize_simulated_annealing(fun, x0, lower=None, upper=None, temp=15.0, step_size=0.1, max_iterations=5000, seed=None)
```

## Minimal working examples

### OLS with robust standard errors

```python
import numpy as np
import crabbymetrics as cm

rng = np.random.default_rng(1)
x = rng.normal(size=(500, 3))
y = 0.5 + x @ np.array([1.0, -2.0, 0.25]) + rng.normal(size=500)

m = cm.OLS()
m.fit(x, y)
s = m.summary(vcov="hc1")
print(s["intercept"], s["coef"], s["coef_se"])
```

### Fixed effects OLS

```python
n = 1000
x = rng.normal(size=(n, 2))
unit = rng.integers(0, 50, size=n, dtype=np.uint32)
time = rng.integers(0, 10, size=n, dtype=np.uint32)
fe = np.column_stack([unit, time]).astype(np.uint32)
y = x @ np.array([1.0, -0.5]) + rng.normal(size=n)

m = cm.FixedEffectsOLS()
m.fit(x, fe, y)
print(m.summary(vcov="cluster", clusters=unit)["coef"])
```

### Synthetic DID panel ATT

```python
n, t = 30, 12
Y = rng.normal(size=(n, t))
W = np.zeros((n, t), dtype=float)
W[20:, 8:] = 1.0
Y[W == 1.0] += 1.0

m = cm.SyntheticDID()
m.fit(Y, W)
out = m.summary()
print(out["att"])
print(out["event_study"].keys())
```

### Matrix completion panel ATT

```python
m = cm.MatrixCompletion(lambda_fraction=0.05, max_iterations=500)
m.fit(Y, W)
out = m.summary()
print(out["att"])
cf = out["counterfactual"]
```

### Balancing weights

```python
source_x = rng.normal(size=(300, 4))
target_x = rng.normal(loc=0.25, size=(100, 4))

bw = cm.BalancingWeights(objective="quadratic", autoscale=True)
bw.fit(source_x, target_x)
weights = bw.get_weights()
print(bw.summary()["effective_sample_size"])
```

### Custom GMM

```python
def moments(theta, data):
    y = data[:, 0]
    x = data[:, 1:]
    resid = y - x @ theta
    return x * resid[:, None]

data = np.column_stack([y, np.column_stack([np.ones(len(y)), x])])
theta0 = np.zeros(data.shape[1] - 1)

g = cm.GMM(moments)
g.fit(data, theta0)
print(g.summary())
```

## Source tree map for agents

If you need to inspect or modify source, start here:

- `src/lib.rs`: PyO3 module registration; exported Python classes and functions.
- `src/estimators/linear.rs`: OLS, FixedEffectsOLS, anytime-valid linear inference, and shared linear covariance helpers.
- `src/estimators/iv.rs`: TwoSLS.
- `src/estimators/panel.rs`: InteractiveFixedEffects, MatrixCompletion, HorizontalPanelRidge, factor/fixed-effect panel helpers, and shared panel summaries.
- `src/estimators/synthetic.rs`: SyntheticControl and SyntheticDID.
- `src/estimators/regularized.rs`: Ridge, ElasticNet, BaggedPolynomialRegressor.
- `src/estimators/mle.rs`: Logit, MultinomialLogit, Poisson, MEstimator.
- `src/estimators/survival.rs`: ExponentialPH, WeibullPH, CoxPH, AndersenGill.
- `src/estimators/semiparametric.rs`: EPLM, AverageDerivative, PartiallyLinearDML, AIPW.
- `src/estimators/balancing.rs`: BalancingWeights.
- `src/estimators/gmm.rs`: GMM.
- `src/estimators/transforms.rs`: PCA, RandomizedPCA, KernelBasis, NystromBasis, and RandomFourierFeatures.
- `src/fit.rs`: shared convergence status and fit diagnostics.
- `src/validation.rs`: shared finite-value, binary-label, positivity, and sample-weight validation.
- `src/optimizers.rs`: Optimizers namespace.
- `src/utils.rs`: shared array conversion, linear algebra, covariance, bootstrap, and helper routines.
- `tests/`: Python tests for public behavior; use these as executable API examples when docs are thin.
- `docs/examples/*.qmd`: rendered example notebooks.
- `docs/ablations/*.qmd`: cached numerical experiments and panel case-study workflows.
- `docs/api.qmd`: generated/verified public-surface reference.

Development commands:

```bash
uv sync
uv run maturin develop
uv run pytest tests -q
uv sync --extra docs
uv run quarto render docs
```

## Documentation links

- [Home](https://apoorvalal.github.io/crabbymetrics/): overview, examples, ablations, and support pages.
- [API reference](https://apoorvalal.github.io/crabbymetrics/api.html): verified public surface, method signatures, summary schemas, optimizer catalog, and smoke checks.
- Estimator pages under the API reference document the implemented objective or estimating equations, covariance method, inferential limitations, and computational scaling.
- [Binding Crash Course: OLS](https://apoorvalal.github.io/crabbymetrics/mechanics-ols.html): shortest end-to-end Rust-to-Python wrapper walkthrough.
- [Binding Internals: Poisson](https://apoorvalal.github.io/crabbymetrics/mechanics-poisson.html): deeper built-in estimator walkthrough.
- [Binding Internals: MEstimator](https://apoorvalal.github.io/crabbymetrics/mechanics-mestimator.html): callback-heavy bridge where Rust owns optimization and Python supplies objective/scores.

## Estimator example pages

- [OLS](https://apoorvalal.github.io/crabbymetrics/examples/ols.html): baseline linear regression and covariance choices.
- [Ridge](https://apoorvalal.github.io/crabbymetrics/examples/ridge.html): L2 regularization and cross-validated penalty grids.
- [Bagged Polynomial Regression](https://apoorvalal.github.io/crabbymetrics/examples/bagged-polynomial-regression.html): leakage-free repeated prediction comparison and OOB diagnostics.
- [Fixed Effects OLS](https://apoorvalal.github.io/crabbymetrics/examples/fixed-effects-ols.html): one-way and multi-way fixed effects.
- [ElasticNet](https://apoorvalal.github.io/crabbymetrics/examples/elastic-net.html): mixed L1/L2 regularized linear regression.
- [Synthetic Control](https://apoorvalal.github.io/crabbymetrics/examples/synthetic-control.html): low-level simplex donor weights for one treated path.
- [Synthetic DID](https://apoorvalal.github.io/crabbymetrics/examples/synthetic-did.html): panel `fit(Y, W)` synthetic difference-in-differences workflow.
- [Logit](https://apoorvalal.github.io/crabbymetrics/examples/logit.html): binary logistic regression.
- [Multinomial Logit](https://apoorvalal.github.io/crabbymetrics/examples/multinomial-logit.html): multiclass classification.
- [Poisson](https://apoorvalal.github.io/crabbymetrics/examples/poisson.html): count regression and QMLE/sandwich inference.
- [Survival / Recurrent Events](https://apoorvalal.github.io/crabbymetrics/examples/survival-time-to-event-and-recurrent-events.html): exponential PH, Weibull PH, Cox PH, and Andersen--Gill counting-process models.
- [TwoSLS](https://apoorvalal.github.io/crabbymetrics/examples/twosls.html): instrumental variables regression.
- [GMM](https://apoorvalal.github.io/crabbymetrics/examples/gmm.html): custom moment estimators.
- [MEstimator Poisson](https://apoorvalal.github.io/crabbymetrics/examples/mestimator-poisson.html): callback M-estimation matched against built-in Poisson.
- [Balancing Weights](https://apoorvalal.github.io/crabbymetrics/examples/balancing-weights.html): entropy/quadratic calibration weights.
- [EPLM](https://apoorvalal.github.io/crabbymetrics/examples/eplm.html): partially linear E-estimator.
- [Average Derivative](https://apoorvalal.github.io/crabbymetrics/examples/average-derivative.html): Oaxaca-Blinder, IPW, and doubly robust average derivative.
- [Double ML And AIPW](https://apoorvalal.github.io/crabbymetrics/examples/double-ml.html): cross-fit nuisance models for DML and binary AIPW.
- [Richer Regression With Transformer Pipelines](https://apoorvalal.github.io/crabbymetrics/examples/richer-regression.html): KernelBasis/PCA features in downstream regressions.
- [PCA And Kernel Basis](https://apoorvalal.github.io/crabbymetrics/examples/pca-and-kernel-basis.html): transformer examples.
- [Optimizers](https://apoorvalal.github.io/crabbymetrics/examples/optimizers.html): direct optimizer usage.
- [GMM With Optimizers](https://apoorvalal.github.io/crabbymetrics/examples/gmm-with-optimizers.html): lower-level optimizer view of GMM-style estimation.

## Ablation and case-study links

- [Variance Estimators](https://apoorvalal.github.io/crabbymetrics/ablations/variance-estimators.html): Monte Carlo coverage for OLS, Poisson, and GMM variance estimators.
- [Semiparametric Estimator Comparisons](https://apoorvalal.github.io/crabbymetrics/ablations/semiparametric-estimators.html): EPLM versus DML and regression/balancing/AIPW comparisons.
- [Bridging Finite And Superpopulation](https://apoorvalal.github.io/crabbymetrics/ablations/bridging-finite-and-superpopulation.html): SATE/PATE coverage comparisons.
- [Panel Estimator DGP Comparisons](https://apoorvalal.github.io/crabbymetrics/ablations/panel-dgp-comparisons.html): TWFE, synthetic control, synthetic DID, and matrix completion across panel DGPs.
- [Same Root Panel Case Studies](https://apoorvalal.github.io/crabbymetrics/ablations/same-root-panel-cases.html): Basque Country and California Proposition 99 panel workflows with horizontal/vertical estimators, synthetic DID, matrix completion, HAC intervals, and donor-placebo inference.

## Optional background links

- [First Course Ding overview](https://apoorvalal.github.io/crabbymetrics/first-course-ding.html): translation map for Ding causal inference chapters.
- [Ding chapter pages](https://apoorvalal.github.io/crabbymetrics/first-course-ding.html): use the site navigation for chapter-specific causal inference notes and translated examples.

## Common mistakes to avoid

- Do not pass pandas DataFrames directly if you need predictable behavior; convert to NumPy arrays first.
- Do not use long panel data for `HorizontalPanelRidge`, `SyntheticDID`, or `MatrixCompletion`; use balanced `(n_units, n_periods)` matrices `Y` and `W`.
- Do not treat `SyntheticControl` and `SyntheticDID` as the same API. `SyntheticControl.fit(donors, treated)` is a low-level one-treated-path donor-weight API; `SyntheticDID.fit(Y, W)` is a panel causal estimator.
- Do not expect uniform `summary()` keys across estimators. Check the relevant estimator page or the API reference.
- Do not expect scipy result objects or sklearn estimators. The API is scikit-adjacent but intentionally lightweight.
- Do not assume every estimator has `predict()` or `bootstrap()`.
- For fixed effects, `fe` must contain zero-based integer category codes, not strings.
- For clustered standard errors, pass a 1D cluster-label array aligned with rows.
- For Newey-West covariance, pass `lags`.
- For binary treatment AIPW, code treatment as `0/1`.
- For classification labels, use integer-coded labels.