Skip to content

Statistics API

asrquant.statistics

Regression, time-series tests, bootstrap inference, and factor analysis.

RegressionResult dataclass

Aligned regression outputs with compact plotting helpers.

Source code in src/asrquant/statistics.py 
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
@dataclass
class RegressionResult:
    """Aligned regression outputs with compact plotting helpers."""

    model: object
    coefficients: pd.Series
    confidence_intervals: pd.DataFrame
    fitted: pd.Series
    residuals: pd.Series
    diagnostics: pd.Series

    @property
    def summary(self) -> pd.Series:
        values = {f"coefficient:{name}": value for name, value in self.coefficients.items()}
        values.update({f"diagnostic:{name}": value for name, value in self.diagnostics.items()})
        return pd.Series(values)

    def to_frame(self) -> pd.DataFrame:
        return pd.concat({"fitted": self.fitted, "residual": self.residuals}, axis=1)

    def to_dict(self) -> dict[str, object]:
        return {
            "result_type": "regression",
            "summary": self.summary.to_dict(),
            "coefficients": self.coefficients.to_dict(),
            "diagnostics": self.diagnostics.to_dict(),
        }

    def plot(self, kind: str = "residuals"):
        from .viz import regression as viz

        if kind in {"residuals", "diagnostics"}:
            return viz.residual_diagnostics(self)
        if kind in {"coefficients", "intervals"}:
            return viz.coefficient_intervals(self)
        if kind in {"fitted", "actual_vs_fitted"}:
            return viz.actual_vs_fitted(self)
        raise ValueError("kind must be residuals, coefficients, or fitted")

ols 

ols(y: Series, x: Series | DataFrame, add_constant: bool = True, covariance: str = 'HAC', maxlags: int | None = None) -> RegressionResult

OLS with classical, HC, or Newey-West/HAC covariance.

Source code in src/asrquant/statistics.py 
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
def ols(
    y: pd.Series,
    x: pd.Series | pd.DataFrame,
    add_constant: bool = True,
    covariance: str = "HAC",
    maxlags: int | None = None,
) -> RegressionResult:
    """OLS with classical, HC, or Newey-West/HAC covariance."""
    y = pd.Series(y, dtype=float).rename("y")
    x = pd.DataFrame(x, dtype=float)
    data = pd.concat([y, x], axis=1).dropna()
    y_clean = data.iloc[:, 0]
    x_clean = data.iloc[:, 1:]
    if add_constant:
        x_clean = sm.add_constant(x_clean, has_constant="add")
    model = sm.OLS(y_clean, x_clean)
    cov_upper = covariance.upper()
    if cov_upper == "HAC":
        lag = maxlags if maxlags is not None else max(1, int(len(data) ** 0.25))
        fitted = model.fit(cov_type="HAC", cov_kwds={"maxlags": lag})
    elif cov_upper in {"HC0", "HC1", "HC2", "HC3"}:
        fitted = model.fit(cov_type=cov_upper)
    elif cov_upper in {"NONROBUST", "CLASSICAL"}:
        fitted = model.fit()
    else:
        raise ValueError("covariance must be HAC, HC0-HC3, or nonrobust")
    residuals = pd.Series(fitted.resid, index=data.index, name="residual")
    diagnostics = pd.Series(
        {
            "R2": fitted.rsquared,
            "Adjusted R2": fitted.rsquared_adj,
            "Durbin-Watson": sm.stats.stattools.durbin_watson(residuals),
            "Jarque-Bera p": sm.stats.stattools.jarque_bera(residuals)[1],
            "Ljung-Box p": float(acorr_ljungbox(residuals, lags=[min(10, max(1, len(residuals)//5))], return_df=True)["lb_pvalue"].iloc[0]),
            "Breusch-Pagan p": het_breuschpagan(residuals, fitted.model.exog)[1],
            "White p": het_white(residuals, fitted.model.exog)[1] if fitted.model.exog.shape[1] > 1 else np.nan,
        }
    )
    return RegressionResult(
        model=fitted,
        coefficients=pd.Series(fitted.params),
        confidence_intervals=fitted.conf_int().set_axis(["lower", "upper"], axis=1),
        fitted=pd.Series(fitted.fittedvalues, index=data.index, name="fitted"),
        residuals=residuals,
        diagnostics=diagnostics,
    )

rolling_regression 

rolling_regression(y: Series, x: Series | DataFrame, window: int = 63, add_constant: bool = True) -> pd.DataFrame

Rolling least-squares coefficients.

Source code in src/asrquant/statistics.py 
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
def rolling_regression(
    y: pd.Series,
    x: pd.Series | pd.DataFrame,
    window: int = 63,
    add_constant: bool = True,
) -> pd.DataFrame:
    """Rolling least-squares coefficients."""
    y = pd.Series(y, dtype=float).rename("y")
    x = pd.DataFrame(x, dtype=float)
    data = pd.concat([y, x], axis=1).dropna()
    x_clean = data.iloc[:, 1:]
    if add_constant:
        x_clean = sm.add_constant(x_clean, has_constant="add")
    model = RollingOLS(data.iloc[:, 0], x_clean, window=window)
    return model.fit().params

stationarity_tests 

stationarity_tests(series: Series) -> pd.Series

ADF and KPSS tests reported together to reduce one-test overinterpretation.

Source code in src/asrquant/statistics.py 
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
def stationarity_tests(series: pd.Series) -> pd.Series:
    """ADF and KPSS tests reported together to reduce one-test overinterpretation."""
    s = pd.Series(series, dtype=float).dropna()
    adf = adfuller(s, autolag="AIC")
    try:
        kpss_result = kpss(s, regression="c", nlags="auto")
        kpss_stat, kpss_p = kpss_result[0], kpss_result[1]
    except ValueError:
        kpss_stat, kpss_p = np.nan, np.nan
    return pd.Series(
        {
            "ADF statistic": adf[0],
            "ADF p-value": adf[1],
            "KPSS statistic": kpss_stat,
            "KPSS p-value": kpss_p,
        }
    )

block_bootstrap 

block_bootstrap(series: Series, statistic: Callable[[Series], float] = np.mean, n_boot: int = 1000, block_size: int | None = None, confidence: float = 0.95, random_state: int | None = 0) -> pd.Series

Moving-block bootstrap confidence interval for dependent data.

Source code in src/asrquant/statistics.py 
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
def block_bootstrap(
    series: pd.Series,
    statistic: Callable[[pd.Series], float] = np.mean,
    n_boot: int = 1_000,
    block_size: int | None = None,
    confidence: float = 0.95,
    random_state: int | None = 0,
) -> pd.Series:
    """Moving-block bootstrap confidence interval for dependent data."""
    s = pd.Series(series, dtype=float).dropna().reset_index(drop=True)
    n = len(s)
    if n < 2:
        raise ValueError("at least two observations are required")
    block_size = block_size or max(2, int(round(n ** (1 / 3))))
    if not 1 <= block_size <= n:
        raise ValueError("block_size must be between 1 and sample length")
    rng = np.random.default_rng(random_state)
    starts = np.arange(0, n - block_size + 1)
    estimates = np.empty(n_boot)
    for i in range(n_boot):
        chunks: list[np.ndarray] = []
        while sum(len(c) for c in chunks) < n:
            start = int(rng.choice(starts))
            chunks.append(s.iloc[start : start + block_size].to_numpy())
        sample = pd.Series(np.concatenate(chunks)[:n])
        estimates[i] = statistic(sample)
    alpha = 1 - confidence
    return pd.Series(
        {
            "estimate": float(statistic(s)),
            "lower": float(np.quantile(estimates, alpha / 2)),
            "upper": float(np.quantile(estimates, 1 - alpha / 2)),
            "bootstrap_std": float(np.std(estimates, ddof=1)),
        }
    )

permutation_test 

permutation_test(x: Series, y: Series, statistic: Callable[[ndarray, ndarray], float] | None = None, n_permutations: int = 5000, random_state: int | None = 0) -> pd.Series

Two-sided permutation test for a difference in means by default.

Source code in src/asrquant/statistics.py 
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
def permutation_test(
    x: pd.Series,
    y: pd.Series,
    statistic: Callable[[np.ndarray, np.ndarray], float] | None = None,
    n_permutations: int = 5_000,
    random_state: int | None = 0,
) -> pd.Series:
    """Two-sided permutation test for a difference in means by default."""
    x = pd.Series(x, dtype=float).dropna().to_numpy()
    y = pd.Series(y, dtype=float).dropna().to_numpy()
    statistic = statistic or (lambda a, b: float(np.mean(a) - np.mean(b)))
    observed = statistic(x, y)
    pooled = np.concatenate([x, y])
    rng = np.random.default_rng(random_state)
    count = 0
    for _ in range(n_permutations):
        shuffled = rng.permutation(pooled)
        value = statistic(shuffled[: len(x)], shuffled[len(x) :])
        count += abs(value) >= abs(observed)
    return pd.Series({"statistic": observed, "p_value": (count + 1) / (n_permutations + 1)})

benjamini_hochberg 

benjamini_hochberg(p_values: Iterable[float], alpha: float = 0.05) -> pd.DataFrame

Benjamini-Hochberg false-discovery-rate correction.

Source code in src/asrquant/statistics.py 
198
199
200
201
202
203
204
205
206
207
208
209
def benjamini_hochberg(p_values: Iterable[float], alpha: float = 0.05) -> pd.DataFrame:
    """Benjamini-Hochberg false-discovery-rate correction."""
    p = np.asarray(list(p_values), dtype=float)
    if np.any((p < 0) | (p > 1)):
        raise ValueError("p-values must lie in [0, 1]")
    n = len(p)
    order = np.argsort(p)
    ranked = p[order]
    adjusted_ranked = np.minimum.accumulate((ranked * n / np.arange(1, n + 1))[::-1])[::-1]
    adjusted = np.empty(n)
    adjusted[order] = np.clip(adjusted_ranked, 0, 1)
    return pd.DataFrame({"p_value": p, "adjusted_p": adjusted, "reject": adjusted <= alpha})

quantile_regression 

quantile_regression(y: Series, x: Series | DataFrame, quantile: float = 0.5, add_constant: bool = True) -> RegressionResult

Linear quantile regression for conditional-tail relationships.

Source code in src/asrquant/statistics.py 
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
def quantile_regression(
    y: pd.Series,
    x: pd.Series | pd.DataFrame,
    quantile: float = 0.5,
    add_constant: bool = True,
) -> RegressionResult:
    """Linear quantile regression for conditional-tail relationships."""
    if not 0 < quantile < 1:
        raise ValueError("quantile must lie in (0, 1)")
    y_s = pd.Series(y, dtype=float).rename("y"); x_f = pd.DataFrame(x, dtype=float)
    data = pd.concat([y_s, x_f], axis=1).dropna(); y_c = data.iloc[:, 0]; x_c = data.iloc[:, 1:]
    if add_constant:
        x_c = sm.add_constant(x_c, has_constant="add")
    fitted = sm.QuantReg(y_c, x_c).fit(q=quantile)
    residuals = pd.Series(y_c - fitted.fittedvalues, index=data.index, name="residual")
    return RegressionResult(
        model=fitted,
        coefficients=pd.Series(fitted.params),
        confidence_intervals=fitted.conf_int().set_axis(["lower", "upper"], axis=1),
        fitted=pd.Series(fitted.fittedvalues, index=data.index, name="fitted"),
        residuals=residuals,
        diagnostics=pd.Series({"Pseudo R2": getattr(fitted, "prsquared", np.nan), "quantile": quantile}),
    )

polynomial_regression 

polynomial_regression(y: Series, x: Series | DataFrame, degree: int = 2, covariance: str = 'HAC', maxlags: int | None = None) -> RegressionResult

Polynomial least-squares regression with named transformed features.

Source code in src/asrquant/statistics.py 
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
def polynomial_regression(
    y: pd.Series,
    x: pd.Series | pd.DataFrame,
    degree: int = 2,
    covariance: str = "HAC",
    maxlags: int | None = None,
) -> RegressionResult:
    """Polynomial least-squares regression with named transformed features."""
    if degree < 1:
        raise ValueError("degree must be at least 1")
    from sklearn.preprocessing import PolynomialFeatures
    x_f = pd.DataFrame(x, dtype=float)
    transformer = PolynomialFeatures(degree=degree, include_bias=False)
    values = transformer.fit_transform(x_f)
    transformed = pd.DataFrame(values, index=x_f.index, columns=transformer.get_feature_names_out(x_f.columns.astype(str)))
    result = ols(y, transformed, add_constant=True, covariance=covariance, maxlags=maxlags)
    result.model.asrquant_transformer = transformer
    return result

regularized_regression 

regularized_regression(y: Series, x: Series | DataFrame, method: str = 'ridge', alpha: float = 1.0, l1_ratio: float = 0.5, standardize: bool = True)

Fit Ridge, Lasso, or Elastic Net in a scikit-learn pipeline.

Source code in src/asrquant/statistics.py 
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
def regularized_regression(
    y: pd.Series,
    x: pd.Series | pd.DataFrame,
    method: str = "ridge",
    alpha: float = 1.0,
    l1_ratio: float = 0.5,
    standardize: bool = True,
):
    """Fit Ridge, Lasso, or Elastic Net in a scikit-learn pipeline."""
    from sklearn.linear_model import ElasticNet, Lasso, Ridge
    from sklearn.pipeline import make_pipeline
    from sklearn.preprocessing import StandardScaler
    x_f = pd.DataFrame(x, dtype=float); y_s = pd.Series(y, dtype=float).rename("y")
    data = pd.concat([y_s, x_f], axis=1).dropna(); y_c = data.iloc[:, 0]; x_c = data.iloc[:, 1:]
    key = method.lower().replace("-", "_")
    if key == "ridge":
        estimator = Ridge(alpha=alpha)
    elif key == "lasso":
        estimator = Lasso(alpha=alpha, max_iter=20_000)
    elif key in {"elastic_net", "elasticnet"}:
        estimator = ElasticNet(alpha=alpha, l1_ratio=l1_ratio, max_iter=20_000)
    else:
        raise ValueError("method must be ridge, lasso, or elastic_net")
    model = make_pipeline(StandardScaler(), estimator) if standardize else estimator
    model.fit(x_c, y_c)
    fitted = pd.Series(model.predict(x_c), index=data.index, name="fitted")
    residuals = (y_c - fitted).rename("residual")
    return {"model": model, "fitted": fitted, "residuals": residuals, "r2": float(model.score(x_c, y_c)), "features": list(x_c.columns)}

logistic_regression 

logistic_regression(y: Series, x: Series | DataFrame, add_constant: bool = True, covariance: str = 'HC1') -> RegressionResult

Binary logistic regression with robust covariance support.

Source code in src/asrquant/statistics.py 
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
def logistic_regression(
    y: pd.Series,
    x: pd.Series | pd.DataFrame,
    add_constant: bool = True,
    covariance: str = "HC1",
) -> RegressionResult:
    """Binary logistic regression with robust covariance support."""
    y_s = pd.Series(y, dtype=float).rename("y"); x_f = pd.DataFrame(x, dtype=float)
    data = pd.concat([y_s, x_f], axis=1).dropna(); y_c = data.iloc[:, 0]; x_c = data.iloc[:, 1:]
    if not set(y_c.unique()).issubset({0.0, 1.0}):
        raise ValueError("logistic target must contain only 0 and 1")
    if add_constant:
        x_c = sm.add_constant(x_c, has_constant="add")
    fitted = sm.Logit(y_c, x_c).fit(disp=False, cov_type=covariance)
    probabilities = pd.Series(fitted.predict(x_c), index=data.index, name="fitted")
    residuals = (y_c - probabilities).rename("residual")
    return RegressionResult(
        model=fitted,
        coefficients=pd.Series(fitted.params),
        confidence_intervals=fitted.conf_int().set_axis(["lower", "upper"], axis=1),
        fitted=probabilities,
        residuals=residuals,
        diagnostics=pd.Series({"Pseudo R2": fitted.prsquared, "AIC": fitted.aic, "BIC": fitted.bic}),
    )

factor_regression 

factor_regression(asset_returns: Series, factors: DataFrame, risk_free: float | Series = 0.0, covariance: str = 'HAC', maxlags: int | None = None) -> RegressionResult

Run a CAPM or multi-factor time-series regression on excess returns.

Source code in src/asrquant/statistics.py 
313
314
315
316
317
318
319
320
321
322
323
324
def factor_regression(
    asset_returns: pd.Series,
    factors: pd.DataFrame,
    risk_free: float | pd.Series = 0.0,
    covariance: str = "HAC",
    maxlags: int | None = None,
) -> RegressionResult:
    """Run a CAPM or multi-factor time-series regression on excess returns."""
    asset = pd.Series(asset_returns, dtype=float)
    rf = pd.Series(risk_free, index=asset.index, dtype=float) if np.isscalar(risk_free) else pd.Series(risk_free, dtype=float)
    excess = (asset - rf).rename("excess_return")
    return ols(excess, pd.DataFrame(factors, dtype=float), covariance=covariance, maxlags=maxlags)

cointegration_test 

cointegration_test(y: Series, x: Series, trend: str = 'c') -> pd.Series

Engle-Granger two-step cointegration test.

Source code in src/asrquant/statistics.py 
327
328
329
330
331
332
def cointegration_test(y: pd.Series, x: pd.Series, trend: str = "c") -> pd.Series:
    """Engle-Granger two-step cointegration test."""
    from statsmodels.tsa.stattools import coint
    data = pd.concat([pd.Series(y, dtype=float), pd.Series(x, dtype=float)], axis=1).dropna()
    statistic, p_value, critical = coint(data.iloc[:, 0], data.iloc[:, 1], trend=trend)
    return pd.Series({"statistic": statistic, "p_value": p_value, "critical_1%": critical[0], "critical_5%": critical[1], "critical_10%": critical[2]})

granger_causality 

granger_causality(x: Series, y: Series, maxlag: int = 5) -> pd.DataFrame

Report Granger-predictive F-test p-values; this is not structural causality.

Source code in src/asrquant/statistics.py 
335
336
337
338
339
340
341
342
343
344
def granger_causality(x: pd.Series, y: pd.Series, maxlag: int = 5) -> pd.DataFrame:
    """Report Granger-predictive F-test p-values; this is not structural causality."""
    from statsmodels.tsa.stattools import grangercausalitytests
    data = pd.concat([pd.Series(y, name="y"), pd.Series(x, name="x")], axis=1).dropna()
    results = grangercausalitytests(data, maxlag=maxlag, verbose=False)
    rows = []
    for lag, output in results.items():
        test = output[0]["ssr_ftest"]
        rows.append({"lag": lag, "F": test[0], "p_value": test[1], "df_denom": test[2], "df_num": test[3]})
    return pd.DataFrame(rows).set_index("lag")

arima_fit 

arima_fit(series: Series, order: tuple[int, int, int] = (1, 0, 0), trend: str | None = None)

Fit an ARIMA model and return the statsmodels result object.

Source code in src/asrquant/statistics.py 
347
348
349
350
351
def arima_fit(series: pd.Series, order: tuple[int, int, int] = (1, 0, 0), trend: str | None = None):
    """Fit an ARIMA model and return the statsmodels result object."""
    from statsmodels.tsa.arima.model import ARIMA
    s = pd.Series(series, dtype=float).dropna()
    return ARIMA(s, order=order, trend=trend).fit()

var_fit 

var_fit(data: DataFrame, lags: int = 1, trend: str = 'c')

Fit a vector autoregression to a multivariate stationary panel.

Source code in src/asrquant/statistics.py 
354
355
356
357
358
359
360
def var_fit(data: pd.DataFrame, lags: int = 1, trend: str = "c"):
    """Fit a vector autoregression to a multivariate stationary panel."""
    from statsmodels.tsa.api import VAR
    frame = pd.DataFrame(data, dtype=float).dropna()
    if frame.shape[1] < 2:
        raise ValueError("VAR requires at least two series")
    return VAR(frame).fit(lags, trend=trend)

autoregression_fit 

autoregression_fit(series: Series, lags: int | list[int] = 1, trend: str = 'c', *, old_names: bool = False)

Fit an explicit AR(p) model with statsmodels AutoReg.

Parameters:

Name Type Description Default
series Series

Univariate time series.

 required
lags int | list[int]

Maximum lag order or an explicit list of included lags.

 1
trend str

Deterministic terms: "n", "c", "t" or "ct".

 'c'
Source code in src/asrquant/statistics.py 
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
def autoregression_fit(
    series: pd.Series,
    lags: int | list[int] = 1,
    trend: str = "c",
    *,
    old_names: bool = False,
):
    """Fit an explicit AR(p) model with statsmodels AutoReg.

    Parameters
    ----------
    series:
        Univariate time series.
    lags:
        Maximum lag order or an explicit list of included lags.
    trend:
        Deterministic terms: ``"n"``, ``"c"``, ``"t"`` or ``"ct"``.
    """
    from statsmodels.tsa.ar_model import AutoReg

    values = pd.Series(series, dtype=float).dropna()
    if len(values) < 3:
        raise ValueError("at least three finite observations are required")
    return AutoReg(values, lags=lags, trend=trend, old_names=old_names).fit()