Statsmodels: Statistical Modeling and Econometrics

Linear Regression (OLS)

import statsmodels.api as sm
import numpy as np
import pandas as pd

# Prepare data - ALWAYS add constant for intercept
X = sm.add_constant(X_data)

# Fit OLS model
model = sm.OLS(y, X)
results = model.fit()

# View comprehensive results
print(results.summary())

# Key results
print(f"R-squared: {results.rsquared:.4f}")
print(f"Coefficients:\\n{results.params}")
print(f"P-values:\\n{results.pvalues}")

# Predictions with confidence intervals
predictions = results.get_prediction(X_new)
pred_summary = predictions.summary_frame()
print(pred_summary)  # includes mean, CI, prediction intervals

# Diagnostics
from statsmodels.stats.diagnostic import het_breuschpagan
bp_test = het_breuschpagan(results.resid, X)
print(f"Breusch-Pagan p-value: {bp_test[1]:.4f}")

# Visualize residuals
import matplotlib.pyplot as plt
plt.scatter(results.fittedvalues, results.resid)
plt.axhline(y=0, color='r', linestyle='--')
plt.xlabel('Fitted values')
plt.ylabel('Residuals')
plt.show()

Logistic Regression (Binary Outcomes)

from statsmodels.discrete.discrete_model import Logit

# Add constant
X = sm.add_constant(X_data)

# Fit logit model
model = Logit(y_binary, X)
results = model.fit()

print(results.summary())

# Odds ratios
odds_ratios = np.exp(results.params)
print("Odds ratios:\\n", odds_ratios)

# Predicted probabilities
probs = results.predict(X)

# Binary predictions (0.5 threshold)
predictions = (probs > 0.5).astype(int)

# Model evaluation
from sklearn.metrics import classification_report, roc_auc_score

print(classification_report(y_binary, predictions))
print(f"AUC: {roc_auc_score(y_binary, probs):.4f}")

# Marginal effects
marginal = results.get_margeff()
print(marginal.summary())

Time Series (ARIMA)

from statsmodels.tsa.arima.model import ARIMA
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

# Check stationarity
from statsmodels.tsa.stattools import adfuller

adf_result = adfuller(y_series)
print(f"ADF p-value: {adf_result[1]:.4f}")

if adf_result[1] > 0.05:
    # Series is non-stationary, difference it
    y_diff = y_series.diff().dropna()

# Plot ACF/PACF to identify p, q
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8))
plot_acf(y_diff, lags=40, ax=ax1)
plot_pacf(y_diff, lags=40, ax=ax2)
plt.show()

# Fit ARIMA(p,d,q)
model = ARIMA(y_series, order=(1, 1, 1))
results = model.fit()

print(results.summary())

# Forecast
forecast = results.forecast(steps=10)
forecast_obj = results.get_forecast(steps=10)
forecast_df = forecast_obj.summary_frame()

print(forecast_df)  # includes mean and confidence intervals

# Residual diagnostics
results.plot_diagnostics(figsize=(12, 8))
plt.show()

Generalized Linear Models (GLM)

import statsmodels.api as sm

# Poisson regression for count data
X = sm.add_constant(X_data)
model = sm.GLM(y_counts, X, family=sm.families.Poisson())
results = model.fit()

print(results.summary())

# Rate ratios (for Poisson with log link)
rate_ratios = np.exp(results.params)
print("Rate ratios:\\n", rate_ratios)

# Check overdispersion
overdispersion = results.pearson_chi2 / results.df_resid
print(f"Overdispersion: {overdispersion:.2f}")

if overdispersion > 1.5:
    # Use Negative Binomial instead
    from statsmodels.discrete.count_model import NegativeBinomial
    nb_model = NegativeBinomial(y_counts, X)
    nb_results = nb_model.fit()
    print(nb_results.summary())

Core Statistical Modeling Capabilities

1. Linear Regression Models

Comprehensive suite of linear models for continuous outcomes with various error structures.

Available models:

• OLS: Standard linear regression with i.i.d. errors
• WLS: Weighted least squares for heteroskedastic errors
• GLS: Generalized least squares for arbitrary covariance structure
• GLSAR: GLS with autoregressive errors for time series
• Quantile Regression: Conditional quantiles (robust to outliers)
• Mixed Effects: Hierarchical/multilevel models with random effects
• Recursive/Rolling: Time-varying parameter estimation

Key features:

• Comprehensive diagnostic tests
• Robust standard errors (HC, HAC, cluster-robust)
• Influence statistics (Cook's distance, leverage, DFFITS)
• Hypothesis testing (F-tests, Wald tests)
• Model comparison (AIC, BIC, likelihood ratio tests)
• Prediction with confidence and prediction intervals

When to use: Continuous outcome variable, want inference on coefficients, need diagnostics

Reference: See references/linear_models.md for detailed guidance on model selection, diagnostics, and best practices.

2. Generalized Linear Models (GLM)

Flexible framework extending linear models to non-normal distributions.

Distribution families:

• Binomial: Binary outcomes or proportions (logistic regression)
• Poisson: Count data
• Negative Binomial: Overdispersed counts
• Gamma: Positive continuous, right-skewed data
• Inverse Gaussian: Positive continuous with specific variance structure
• Gaussian: Equivalent to OLS
• Tweedie: Flexible family for semi-continuous data

Link functions:

• Logit, Probit, Log, Identity, Inverse, Sqrt, CLogLog, Power
• Choose based on interpretation needs and model fit

Key features:

• Maximum likelihood estimation via IRLS
• Deviance and Pearson residuals
• Goodness-of-fit statistics
• Pseudo R-squared measures
• Robust standard errors

When to use: Non-normal outcomes, need flexible variance and link specifications

Reference: See references/glm.md for family selection, link functions, interpretation, and diagnostics.

3. Discrete Choice Models

Models for categorical and count outcomes.

Binary models:

• Logit: Logistic regression (odds ratios)
• Probit: Probit regression (normal distribution)

Multinomial models:

• MNLogit: Unordered categories (3+ levels)
• Conditional Logit: Choice models with alternative-specific variables
• Ordered Model: Ordinal outcomes (ordered categories)

Count models:

• Poisson: Standard count model
• Negative Binomial: Overdispersed counts
• Zero-Inflated: Excess zeros (ZIP, ZINB)
• Hurdle Models: Two-stage models for zero-heavy data

Key features:

• Maximum likelihood estimation
• Marginal effects at means or average marginal effects
• Model comparison via AIC/BIC
• Predicted probabilities and classification
• Goodness-of-fit tests

When to use: Binary, categorical, or count outcomes

Reference: See references/discrete_choice.md for model selection, interpretation, and evaluation.

4. Time Series Analysis

Comprehensive time series modeling and forecasting capabilities.

Univariate models:

• AutoReg (AR): Autoregressive models
• ARIMA: Autoregressive integrated moving average
• SARIMAX: Seasonal ARIMA with exogenous variables
• Exponential Smoothing: Simple, Holt, Holt-Winters
• ETS: Innovations state space models

Multivariate models:

• VAR: Vector autoregression
• VARMAX: VAR with MA and exogenous variables
• Dynamic Factor Models: Extract common factors
• VECM: Vector error correction models (cointegration)

Advanced models:

• State Space: Kalman filtering, custom specifications
• Regime Switching: Markov switching models
• ARDL: Autoregressive distributed lag

Key features:

• ACF/PACF analysis for model identification
• Stationarity tests (ADF, KPSS)
• Forecasting with prediction intervals
• Residual diagnostics (Ljung-Box, heteroskedasticity)
• Granger causality testing
• Impulse response functions (IRF)
• Forecast error variance decomposition (FEVD)

When to use: Time-ordered data, forecasting, understanding temporal dynamics

Reference: See references/time_series.md for model selection, diagnostics, and forecasting methods.

5. Statistical Tests and Diagnostics

Extensive testing and diagnostic capabilities for model validation.

Residual diagnostics:

• Autocorrelation tests (Ljung-Box, Durbin-Watson, Breusch-Godfrey)
• Heteroskedasticity tests (Breusch-Pagan, White, ARCH)
• Normality tests (Jarque-Bera, Omnibus, Anderson-Darling, Lilliefors)
• Specification tests (RESET, Harvey-Collier)

Influence and outliers:

• Leverage (hat values)
• Cook's distance
• DFFITS and DFBETAs
• Studentized residuals
• Influence plots

Hypothesis testing:

• t-tests (one-sample, two-sample, paired)
• Proportion tests
• Chi-square tests
• Non-parametric tests (Mann-Whitney, Wilcoxon, Kruskal-Wallis)
• ANOVA (one-way, two-way, repeated measures)

Multiple comparisons:

• Tukey's HSD
• Bonferroni correction
• False Discovery Rate (FDR)

Effect sizes and power:

• Cohen's d, eta-squared
• Power analysis for t-tests, proportions
• Sample size calculations

Robust inference:

• Heteroskedasticity-consistent SEs (HC0-HC3)
• HAC standard errors (Newey-West)
• Cluster-robust standard errors

When to use: Validating assumptions, detecting problems, ensuring robust inference

Reference: See references/stats_diagnostics.md for comprehensive testing and diagnostic procedures.

Formula API (R-style)

Statsmodels supports R-style formulas for intuitive model specification:

import statsmodels.formula.api as smf

# OLS with formula
results = smf.ols('y ~ x1 + x2 + x1:x2', data=df).fit()

# Categorical variables (automatic dummy coding)
results = smf.ols('y ~ x1 + C(category)', data=df).fit()

# Interactions
results = smf.ols('y ~ x1 * x2', data=df).fit()  # x1 + x2 + x1:x2

# Polynomial terms
results = smf.ols('y ~ x + I(x**2)', data=df).fit()

# Logit
results = smf.logit('y ~ x1 + x2 + C(group)', data=df).fit()

# Poisson
results = smf.poisson('count ~ x1 + x2', data=df).fit()

# ARIMA (not available via formula, use regular API)

Model Selection and Comparison

Information Criteria

# Compare models using AIC/BIC
models = {
    'Model 1': model1_results,
    'Model 2': model2_results,
    'Model 3': model3_results
}

comparison = pd.DataFrame({
    'AIC': {name: res.aic for name, res in models.items()},
    'BIC': {name: res.bic for name, res in models.items()},
    'Log-Likelihood': {name: res.llf for name, res in models.items()}
})

print(comparison.sort_values('AIC'))
# Lower AIC/BIC indicates better model

Likelihood Ratio Test (Nested Models)

# For nested models (one is subset of the other)
from scipy import stats

lr_stat = 2 * (full_model.llf - reduced_model.llf)
df = full_model.df_model - reduced_model.df_model
p_value = 1 - stats.chi2.cdf(lr_stat, df)

print(f"LR statistic: {lr_stat:.4f}")
print(f"p-value: {p_value:.4f}")

if p_value < 0.05:
    print("Full model significantly better")