Linearmodels

linearmodels Skill

linearmodels: panel data, IV/GMM, system regression, and asset pricing models in Python. Covers PanelOLS (FE/RE), BetweenOLS, FirstDifferenceOLS, Fama-MacBeth, IV2SLS/LIML/GMM, SUR, IV3SLS, and Driscoll-Kraay SEs. Use for random effects estimation, between or first-difference panel models, system estimation (SUR, 3SLS), LIML/GMM instrumental variables, Fama-MacBeth regressions, or Driscoll-Kraay standard errors. Complements pyfixest (high-dimensional FE + DiD) and statsmodels (GLM + time series).

Comprehensive skill for panel data estimation, instrumental variables, system regression, and asset pricing with linearmodels (Kevin Sheppard). Use decision trees below to find the right guidance, then load detailed references.

What is linearmodels?

linearmodels extends statsmodels with specialized model classes for structured data:

• Panel data: PanelOLS (fixed effects), RandomEffects, BetweenOLS, FirstDifferenceOLS, PooledOLS, FamaMacBeth
• Instrumental variables: IV2SLS, IVLIML (k-class), IVGMM, IVGMMCUE (continuously updating), AbsorbingLS
• System estimation: SUR (Seemingly Unrelated Regression), IV3SLS, IVSystemGMM
• Asset pricing: LinearFactorModel, LinearFactorModelGMM, TradedFactorModel
• Rich inference: Driscoll-Kraay, clustered (1- and 2-way), HAC kernels (Bartlett, Parzen, Quadratic Spectral)
• Dual API: Formula-based (via formulaic) and array-based interfaces

How to Use This Skill

Reference File Structure

Reading Order

1. New to linearmodels? Start with quickstart.md then panel-models.md 2. Need IV/GMM? Read quickstart.md then iv-models.md 3. System estimation (SUR/3SLS)? Read quickstart.md then system-models.md 4. Asset pricing? Read quickstart.md then asset-pricing.md 5. Choosing SEs? Read covariance-inference.md 6. Coming from pyfixest? Read quickstart.md then gotchas.md

Related Skills

Quick Decision Trees

"I need a panel model"

What panel estimation method?
├─ Fixed effects (within estimator)
│   ├─ 1-2 way FE, no IV → linearmodels PanelOLS or pyfixest feols
│   ├─ 3+ way FE → pyfixest (linearmodels max 2-way in PanelOLS)
│   ├─ FE + IV combined → pyfixest (linearmodels has no Panel IV)
│   └─ FE + DiD → pyfixest (linearmodels has no DiD)
├─ Random effects (GLS) → linearmodels RandomEffects
│   └─ → ./references/panel-models.md
├─ FE vs RE comparison → linearmodels (run both, compare)
│   └─ → ./references/panel-models.md
├─ Between estimator → linearmodels BetweenOLS
│   └─ → ./references/panel-models.md
├─ First difference → linearmodels FirstDifferenceOLS
│   └─ → ./references/panel-models.md
├─ Pooled OLS (panel-aware SEs) → linearmodels PooledOLS
│   └─ → ./references/panel-models.md
└─ Fama-MacBeth → linearmodels FamaMacBeth
    └─ → ./references/panel-models.md

"I need IV / GMM estimation"

What IV method?
├─ 2SLS (standard IV)
│   ├─ With fixed effects → pyfixest (linearmodels has no Panel IV)
│   └─ Without FE → linearmodels IV2SLS or pyfixest
│       └─ → ./references/iv-models.md
├─ LIML / k-class (better finite-sample) → linearmodels IVLIML
│   └─ → ./references/iv-models.md
├─ GMM-IV (efficient, overidentified) → linearmodels IVGMM
│   └─ → ./references/iv-models.md
├─ Continuously updating GMM → linearmodels IVGMMCUE
│   └─ → ./references/iv-models.md
└─ High-dimensional absorbed FE (OLS) → linearmodels AbsorbingLS
    └─ → ./references/iv-models.md

"I need system estimation"

System of equations?
├─ Multiple equations, correlated errors → SUR
│   └─ → ./references/system-models.md
├─ Multiple equations + endogenous variables → IV3SLS
│   └─ → ./references/system-models.md
├─ System GMM → IVSystemGMM
│   └─ → ./references/system-models.md
├─ Cross-equation parameter restrictions → LinearConstraint
│   └─ → ./references/system-models.md
└─ Not sure which → Start with SUR
    └─ → ./references/system-models.md

"Something isn't working"

Having issues?
├─ TypeError about DataFrame index → ./references/gotchas.md
├─ Need FE + IV in one model → ./references/gotchas.md
├─ Need 3+ way fixed effects → ./references/gotchas.md
├─ Constant term confusion → ./references/gotchas.md
├─ Formula parsing errors → ./references/gotchas.md
├─ Want to compare with pyfixest → ./references/gotchas.md
└─ SUR performance issues → ./references/gotchas.md

File-First Execution in Research Workflows

Important: In data research pipelines (see CLAUDE.md), linearmodels estimation is executed through script files, not interactively. This ensures auditability and reproducibility.

The pattern: 1. Write model code to scripts/stage8_analysis/{step}_{task-name}.py 2. Execute via Bash with automatic output capture wrapper script 3. Validation results get automatically embedded in scripts as comments 4. If failed, create versioned copy for fixes

Closely read agent_reference/SCRIPT_EXECUTION_REFERENCE.md for the mandatory file-first execution protocol covering complete code file writing, output capture, and file versioning rules.

The examples below show linearmodels syntax. In research workflows, wrap them in scripts following the file-first pattern.

---

Quick Reference

Essential Imports

from linearmodels.panel import PanelOLS, RandomEffects, BetweenOLS
from linearmodels.panel import FirstDifferenceOLS, PooledOLS, FamaMacBeth
from linearmodels.iv import IV2SLS, IVLIML, IVGMM, IVGMMCUE, AbsorbingLS
from linearmodels.system import SUR, IV3SLS, IVSystemGMM
from linearmodels.panel import compare  # Panel model comparison tables

Data Setup (Critical — MultiIndex Required for Panel Models)

import pandas as pd

# Panel data MUST have a MultiIndex with (entity, time)
df = df.set_index(["entity_id", "year"])

# Verify the index
print(f"Index names: {df.index.names}")
print(f"Index levels: {df.index.nlevels}")

Core Operations

Formula Syntax

# Panel FE keywords (appear in formula, not after |)
"y ~ x1 + x2 + EntityEffects"               # Entity FE
"y ~ x1 + x2 + EntityEffects + TimeEffects"  # Two-way FE
"y ~ x1 + x2 + TimeEffects"                  # Time FE only

# IV bracket notation
"y ~ 1 + exog + [endog ~ instrument1 + instrument2]"

# Suppress intercept
"y ~ x1 + x2 - 1"

Topic Index

Citation

When this library is used as a primary analytical tool, include in the report's Software & Tools references:

Sheppard, K. linearmodels: Econometric models for panel data, IV/GMM, and system regression [Computer software]. https://bashtage.github.io/linearmodels/

Cite when: linearmodels is used for panel estimation (RE, between), IV/GMM, Fama-MacBeth, or system regression (SUR, 3SLS). Do not cite when: Only imported but no estimation performed.

For method-specific citations (e.g., individual estimators or techniques), consult the reference files in this skill and agent_reference/CITATION_REFERENCE.md.

Asset Pricing Models

Reference for asset pricing estimators in linearmodels.asset_pricing. These models test whether a set of risk factors explains the cross-section of expected returns. For Fama-MacBeth regressions (the most commonly used asset pricing estimation method), see ./panel-models.md -- FamaMacBeth lives in linearmodels.panel. This file covers the specialized factor model classes that test asset pricing theories.

Contents

• Overview
• LinearFactorModel (Non-Traded Factors)
• LinearFactorModelGMM
• TradedFactorModel
• Interpreting Output
• Covariance Options
• Common Use Cases
• Data Format
• References and Further Reading

Overview

Three model classes for testing linear factor pricing models:

• LinearFactorModel -- two-step estimation for non-traded factors
• LinearFactorModelGMM -- GMM estimation for non-traded factors (more efficient)
• TradedFactorModel -- SUR-based estimation for traded factor portfolios

All test the same pricing restriction: E[R_i] = lambda_0 + lambda_1 beta_i1 + ... + lambda_K beta_iK. Test portfolios (rows) are assets or portfolios; factors explain cross-sectional return differences. The J-statistic tests whether the model correctly prices all test portfolios (null: all pricing errors are jointly zero).

LinearFactorModel (Non-Traded Factors)

Two-step estimation: (1) time-series regressions of each portfolio's excess returns on factors to get betas, (2) cross-sectional regression of mean excess returns on estimated betas to get risk premia. For factors that are NOT excess returns themselves (e.g., macro factors, consumption growth, labor income growth).

from linearmodels.asset_pricing import LinearFactorModel

# portfolios: T x N DataFrame of excess returns (time periods x test portfolios)
# factors: T x K DataFrame of factor values (not necessarily excess returns)
mod = LinearFactorModel(portfolios=excess_returns, factors=factors)
res = mod.fit(cov_type="kernel", kernel="bartlett", bandwidth=6)
print(res.summary)

# J-statistic: test of whether all pricing errors (alphas) are jointly zero
print(f"J-statistic: {res.j_statistic.stat:.3f} (p={res.j_statistic.pval:.3f})")

Risk-Free Rate Handling

# LinearFactorModel estimates a free intercept (lambda_0) in the cross-sectional
# regression. If the model is correct, lambda_0 should equal the risk-free rate
# (or zero if portfolios are already excess returns).
print(f"Estimated risk-free rate (lambda_0): {res.risk_premia.iloc[0]:.4f}")

LinearFactorModelGMM

GMM estimation of the same non-traded factor model. Jointly estimates factor betas and risk premia in a single system rather than sequentially. More efficient than two-step when the model is correctly specified. Same interface as LinearFactorModel.

from linearmodels.asset_pricing import LinearFactorModelGMM

mod = LinearFactorModelGMM(portfolios=excess_returns, factors=factors)
res = mod.fit(cov_type="kernel", kernel="bartlett", bandwidth=6)
print(res.summary)

# Same output structure as LinearFactorModel
print(f"J-statistic: {res.j_statistic.stat:.3f} (p={res.j_statistic.pval:.3f})")
print(res.risk_premia)

TradedFactorModel

For factors that ARE traded portfolios or excess returns (e.g., Fama-French HML, SMB, Mkt-RF). Uses SUR-based estimation: runs time-series regressions of each test portfolio on the factors, then tests whether the intercepts (alphas) are jointly zero. No cross-sectional regression needed because factor risk premia are directly observed as factor mean returns.

from linearmodels.asset_pricing import TradedFactorModel

# factors must be excess returns (e.g., Mkt-RF, SMB, HML from Ken French's library)
mod = TradedFactorModel(portfolios=excess_returns, factors=factor_returns)
res = mod.fit(cov_type="kernel", kernel="bartlett", bandwidth=6)
print(res.summary)

# GRS test (Gibbons, Ross, Shanken 1989) — tests if all alphas are jointly zero
print(f"J-statistic (GRS): {res.j_statistic.stat:.3f}")
print(f"p-value: {res.j_statistic.pval:.3f}")

# Individual alphas: pricing errors per portfolio
print(res.alphas)

Interpreting Alphas

# Large positive alpha = portfolio earns more than factors predict (underpriced)
# Large negative alpha = portfolio earns less than factors predict (overpriced)
# If model is correct, all alphas should be indistinguishable from zero
alphas = res.alphas
print("Alphas with largest absolute pricing errors:")
print(alphas.reindex(alphas.abs().sort_values(ascending=False).index).head(5))

Interpreting Output

All three model classes produce results with these key elements:

The J-statistic is the central test: rejection means the factor model is incomplete and fails to explain the cross-section of returns. Individual alphas identify which portfolios the model misprices.

res = mod.fit(cov_type="kernel", kernel="bartlett", bandwidth=6)

# Risk premia
print("Estimated risk premia (annualized if returns are monthly, multiply by 12):")
print(res.risk_premia)

# Factor betas (N x K matrix)
print(f"Betas shape: {res.betas.shape}")
print(res.betas.head())

# Full summary with standard errors, t-stats, p-values
print(res.summary)

Covariance Options

All asset pricing models support HAC covariance estimation for time-series dependence in returns.

# Newey-West with 6 lags (common for monthly data)
res = mod.fit(cov_type="kernel", kernel="bartlett", bandwidth=6)

# Parzen kernel
res = mod.fit(cov_type="kernel", kernel="parzen", bandwidth=6)

# Homoskedastic (no correction — rarely appropriate for asset returns)
res = mod.fit(cov_type="unadjusted")

Common Use Cases

CAPM test: Single market factor, test if alpha = 0 for portfolios sorted on size, value, momentum, etc.

mod = TradedFactorModel(portfolios=test_portfolios, factors=mkt_excess[["Mkt-RF"]])
res = mod.fit(cov_type="kernel", kernel="bartlett", bandwidth=6)

Fama-French 3-factor: Mkt-RF, SMB, HML -- all traded factors.

ff3 = factors[["Mkt-RF", "SMB", "HML"]]
mod = TradedFactorModel(portfolios=test_portfolios, factors=ff3)
res = mod.fit(cov_type="kernel", kernel="bartlett", bandwidth=6)

Fama-French 5-factor + momentum: Extend factor set.

ff6 = factors[["Mkt-RF", "SMB", "HML", "RMW", "CMA", "Mom"]]
mod = TradedFactorModel(portfolios=test_portfolios, factors=ff6)
res = mod.fit(cov_type="kernel", kernel="bartlett", bandwidth=6)

Macro factor models (consumption CAPM): Non-traded factors require LinearFactorModel.

macro_factors = factors[["consumption_growth", "labor_income_growth"]]
mod = LinearFactorModel(portfolios=test_portfolios, factors=macro_factors)
res = mod.fit(cov_type="kernel", kernel="bartlett", bandwidth=6)

Data Format

• portfolios: T x N pandas DataFrame. Rows = time periods. Columns = test portfolios (e.g., 25 size-value sorted portfolios). Values are excess returns (subtract risk-free rate before estimation).
• factors: T x K pandas DataFrame. Rows = time periods (must align with portfolios index). Columns = risk factors.
• Both DataFrames must share the same time index.
• No MultiIndex required -- these are NOT panel models.

import pandas as pd

# Typical setup with Fama-French data
portfolios = pd.read_parquet("data/raw/ff25_excess_returns.parquet")
factors = pd.read_parquet("data/raw/ff_factors.parquet")

# Verify alignment
assert portfolios.index.equals(factors.index), "Time indices must match"
print(f"Time periods: {len(portfolios)}")
print(f"Test portfolios: {portfolios.shape[1]}")
print(f"Factors: {factors.shape[1]}")

References and Further Reading

• Sheppard, K. linearmodels Asset Pricing documentation. https://bashtage.github.io/linearmodels/asset-pricing/introduction.html
• Fama, E.F. and French, K.R. (1993). "Common Risk Factors in the Returns on Stocks and Bonds." Journal of Financial Economics, 33(1), 3-56.
• Cochrane, J.H. (2005). Asset Pricing. Revised ed. Princeton University Press.
• Gibbons, M.R., Ross, S.A., and Shanken, J. (1989). "A Test of the Efficiency of a Given Portfolio." Econometrica, 57(5), 1121-1152.

Covariance Estimation and Inference

Reference for all standard error and covariance types in linearmodels across panel, IV, and system model families. For methodology guidance on when to use which SE type, see the data-scientist skill's statistical-modeling.md reference.

Contents

• Overview
• Panel Model Covariance Types
• IV Model Covariance Types
• System Model Covariance Types
• GMM Weight Matrices
• The debiased Parameter
• SE Type Decision Guide
• Comparison with pyfixest and statsmodels
• Key Differences from pyfixest
• References and Further Reading

Overview

linearmodels provides rich covariance estimation across all model types. Key points:

• Panel-aware clustering: entity and time dimensions are built into the API

(cluster_entity, cluster_time) -- no need to pass cluster variables manually

• Driscoll-Kraay SEs: kernel-based covariance robust to cross-sectional dependence,

available for panel models -- not offered by statsmodels

• Small-sample corrections: all covariance estimators accept a debiased parameter
• SE type is specified at `.fit()` time, not post-estimation -- unlike pyfixest's

.vcov() method, you must re-fit to change the covariance type

Panel Model Covariance Types

All panel estimators (PanelOLS, RandomEffects, BetweenOLS, FirstDifferenceOLS, PooledOLS, FamaMacBeth) accept cov_type in .fit().

Homoskedastic (Unadjusted)

res = mod.fit(cov_type="unadjusted")  # or "homoskedastic"

Assumes homoskedastic, uncorrelated errors. Rarely appropriate for panel data -- use only as a baseline or when you have strong theoretical justification.

Heteroskedasticity-Robust

res = mod.fit(cov_type="robust")  # or "heteroskedastic"

White's heteroskedasticity-consistent estimator. Corrects for non-constant error variance but does NOT account for within-entity serial correlation or cross-sectional dependence. Typically insufficient for panel data.

Clustered (One-Way and Two-Way)

# Cluster by entity (most common for panel data)
res = mod.fit(cov_type="clustered", cluster_entity=True)

# Cluster by time
res = mod.fit(cov_type="clustered", cluster_time=True)

# Two-way clustering (entity AND time)
res = mod.fit(cov_type="clustered", cluster_entity=True, cluster_time=True)

# Cluster by a custom variable (must align with the data index)
res = mod.fit(cov_type="clustered", clusters=df["state_id"])

• Entity clustering accounts for arbitrary within-entity serial correlation
• Time clustering accounts for cross-sectional correlation within time periods
• Two-way clustering accounts for both simultaneously (Cameron, Gelbach, Miller 2011)
• When clustering by the same variable used for FE, the df correction properly

accounts for absorbed effects

• Custom clusters variable must have the same index as the dependent variable

Driscoll-Kraay (Kernel)

res = mod.fit(cov_type="kernel", kernel="bartlett", bandwidth=5)

Robust to heteroskedasticity, serial correlation, AND cross-sectional dependence. This is the key covariance type that linearmodels offers beyond what most Python packages provide.

• Appropriate when T is large relative to N (unlike clustered SEs which assume

large N, fixed T)

• Kernel options: "bartlett" (Newey-West), "parzen", "qs" (Quadratic Spectral)
• bandwidth: number of lags included (larger = more serial correlation accounted for)
• Bandwidth rule of thumb for Bartlett: floor(4 * (T/100)^(2/9))

import math

# Automatic bandwidth selection (Bartlett kernel)
T = df.index.get_level_values(1).nunique()
bw = math.floor(4 * (T / 100) ** (2 / 9))
print(f"Suggested Bartlett bandwidth for T={T}: {bw}")

res = mod.fit(cov_type="kernel", kernel="bartlett", bandwidth=bw)

Fama-MacBeth Covariance

Only available for the FamaMacBeth estimator. Adjusts inference for time-series dependence in the averaged cross-sectional regression coefficients.

from linearmodels.panel import FamaMacBeth

mod = FamaMacBeth.from_formula("ret ~ 1 + beta + size + bm", data=df)

# Standard FM SEs (assumes independence across time periods)
res = mod.fit(cov_type="unadjusted")

# HAC-adjusted FM SEs (accounts for serial dependence in coefficient averages)
res = mod.fit(cov_type="kernel", kernel="bartlett", bandwidth=5)

IV Model Covariance Types

All IV estimators (IV2SLS, IVLIML, IVGMM, IVGMMCUE) accept cov_type in .fit(). IV models work with flat DataFrames (no MultiIndex).

Homoskedastic

res = mod.fit(cov_type="unadjusted")

Classical (homoskedastic) standard errors. Textbook settings only.

Heteroskedasticity-Robust

res = mod.fit(cov_type="robust")

Eicker-Huber-White heteroskedasticity-consistent SEs. The default choice for cross-sectional IV applications.

Clustered

res = mod.fit(cov_type="clustered", clusters=df["group"])

Cluster-robust SEs. Pass the grouping variable directly -- IV models do not have cluster_entity/cluster_time shortcuts (those are panel-only).

HAC / Kernel

res = mod.fit(cov_type="kernel", kernel="bartlett", bandwidth=5)

HAC (Newey-West style) covariance for time-series or spatial correlation in cross-sectional IV data. Same kernel options as panel models: "bartlett", "parzen", "qs".

System Model Covariance Types

System models (SUR, IV3SLS, IVSystemGMM) support covariance estimation on top of their cross-equation error structure. The system sigma matrix (cross-equation covariance) is always estimated; cov_type controls within-equation inference.

from linearmodels.system import SUR

mod = SUR(equations)

# Homoskedastic (GLS with estimated sigma)
res = mod.fit(cov_type="unadjusted")

# Heteroskedasticity-robust
res = mod.fit(cov_type="robust")

# HAC kernel
res = mod.fit(cov_type="kernel", kernel="bartlett", bandwidth=5)

# Clustered
res = mod.fit(cov_type="clustered", clusters=df["group"])

For IVSystemGMM, the weight_type parameter controls the GMM weight matrix independently of cov_type, following the same pattern as single-equation IVGMM.

GMM Weight Matrices

For IVGMM and IVSystemGMM, the weight matrix determines estimator efficiency while cov_type determines inference. These should generally match.

from linearmodels.iv import IVGMM

mod = IVGMM.from_formula(
    "np.log(wage) ~ 1 + exper + [educ ~ motheduc + fatheduc]", data=df
)

# Homoskedastic weight matrix
res = mod.fit(weight_type="unadjusted", cov_type="unadjusted")

# Robust weight matrix (heteroskedastic-efficient)
res = mod.fit(weight_type="robust", cov_type="robust")

# HAC weight matrix (serial-correlation-efficient)
res = mod.fit(weight_type="kernel", kernel="bartlett", bandwidth=5,
              cov_type="kernel")

# Clustered weight matrix
res = mod.fit(weight_type="clustered", clusters=df["group"],
              cov_type="clustered")

Mismatching weight_type and cov_type is legal but produces a less efficient estimator. The main case for intentional mismatch: use weight_type="robust" for point estimates but cov_type="clustered" for conservative inference.

The debiased Parameter

# Apply small-sample correction
res = mod.fit(cov_type="clustered", cluster_entity=True, debiased=True)

# No correction (default for some models)
res = mod.fit(cov_type="clustered", cluster_entity=True, debiased=False)

• When True: applies (N-1)/(N-K) type adjustments to the covariance matrix
• For clustered SEs with FE: the correction accounts for degrees of freedom

consumed by absorbed entity/time effects

• Default value varies by model type -- always specify explicitly for

reproducibility

# Explicit debiased across different SE types
res_robust = mod.fit(cov_type="robust", debiased=True)
res_cluster = mod.fit(cov_type="clustered", cluster_entity=True, debiased=True)
res_kernel = mod.fit(cov_type="kernel", kernel="bartlett", bandwidth=5, debiased=True)

The impact of debiased shrinks as sample size grows. For small panels (N < 50 entities or T < 20 periods), the correction can meaningfully affect inference.

SE Type Decision Guide

Default recommendation for most panel applications: entity-clustered SEs. They are robust to arbitrary within-entity correlation patterns and are the most widely reported in applied econometrics.

Comparison with pyfixest and statsmodels

Key Differences from pyfixest

• Fit-time vs post-estimation: linearmodels requires SE choice at .fit() time;

pyfixest allows post-estimation switching via .vcov(). To compare SE types in linearmodels, you must call .fit() multiple times.

• Kernel options: linearmodels offers three kernels for Driscoll-Kraay (Bartlett,

Parzen, Quadratic Spectral) with explicit bandwidth control; pyfixest has basic DK support with vcov="DK".

• Bootstrap: pyfixest offers CRV3 and wild bootstrap inference that linearmodels

does not provide. Use pyfixest when bootstrap-based inference is needed.

• Two-way clustering: both support it; linearmodels uses `cluster_entity=True,

cluster_time=True; pyfixest uses vcov={"CRV1": "entity+time"}`.

• Panel-aware shortcuts: linearmodels provides cluster_entity and cluster_time

booleans that automatically use the MultiIndex levels; pyfixest references variable names from the formula or data.

References and Further Reading

• Driscoll, J.C. and Kraay, A.C. (1998). "Consistent Covariance Matrix Estimation with Spatially Dependent Panel Data." Review of Economics and Statistics, 80(4), 549-560.
• Cameron, A.C. and Miller, D.L. (2015). "A Practitioner's Guide to Cluster-Robust Inference." Journal of Human Resources, 50(2), 317-372.
• Newey, W.K. and West, K.D. (1987). "A Simple, Positive Semi-definite, Heteroskedasticity and Autocorrelation Consistent Covariance Matrix." Econometrica, 55(3), 703-708.
• Petersen, M.A. (2009). "Estimating Standard Errors in Finance Panel Data Sets: Comparing Approaches." Review of Financial Studies, 22(1), 435-480.
• Sheppard, K. linearmodels Covariance documentation. https://bashtage.github.io/linearmodels/panel/reference.html

Common Gotchas and Troubleshooting

A reference for known sharp edges, silent failures, and library-specific behaviors in linearmodels. Each section describes a problem, why it happens, and how to fix it. No methodology content -- syntax and library guidance only.

---

Contents

1. MultiIndex Requirement for Panel Models 2. Maximum Two-Way Fixed Effects in PanelOLS 3. No Panel IV (FE + IV Combined) 4. No DiD Estimators 5. No Automatic Singleton Removal 6. No Post-Estimation SE Switching 7. Constant Term Handling 8. Formula Syntax Differences from pyfixest 9. SUR Performance with Large Datasets 10. Three-Way Package Boundary

---

1. MultiIndex Requirement for Panel Models

This is the most common linearmodels error. All panel model estimators (PanelOLS, RandomEffects, BetweenOLS, FirstDifferenceOLS, PooledOLS, FamaMacBeth) require a pandas DataFrame with a MultiIndex where level 0 is the entity identifier and level 1 is the time identifier. Passing a flat DataFrame raises an error immediately.

import pandas as pd
from linearmodels.panel import PanelOLS

# WRONG -- flat DataFrame with entity and time as regular columns
df = pd.read_csv("panel_data.csv")
mod = PanelOLS.from_formula("y ~ x + EntityEffects", data=df)
# ValueError: The index on the data must be a MultiIndex with 2 levels

# CORRECT -- set MultiIndex before estimation
df = df.set_index(["entity_id", "year"])
mod = PanelOLS.from_formula("y ~ x + EntityEffects", data=df)
res = mod.fit()

From Polars

df_pandas = df_polars.to_pandas()
df_pandas = df_pandas.set_index(["entity_id", "year"])

Checking Your Index

print(f"Index names: {df.index.names}")    # Should be ['entity_id', 'year']
print(f"Index levels: {df.index.nlevels}") # Should be 2
assert df.index.nlevels == 2, "MultiIndex required for panel models"

Important: IV models (IV2SLS, IVLIML, IVGMM, IVGMMCUE) do NOT require a MultiIndex -- they work with ordinary DataFrames. System models (SUR, 3SLS, SystemGMM) also do not require a MultiIndex.

---

2. Maximum Two-Way Fixed Effects in PanelOLS

PanelOLS supports at most two sets of absorbed effects: entity effects, time effects, or one of each plus an other_effects term. Attempting three or more dimensions raises an error.

from linearmodels.panel import PanelOLS

# Two-way -- OK
mod = PanelOLS(y, x, entity_effects=True, time_effects=True)

# Three-way -- ERROR: too many absorbed effects
mod = PanelOLS(y, x, entity_effects=True, time_effects=True, other_effects=z)

Workarounds

AbsorbingLS from linearmodels.iv supports unlimited absorbed dimensions via the pyhdfe backend:

from linearmodels.iv import AbsorbingLS

mod = AbsorbingLS(y, x, absorb=df[["entity", "year", "industry"]])
res = mod.fit()

pyfixest handles multi-way FE natively and is generally faster:

import pyfixest as pf

fit = pf.feols("y ~ x | entity + year + industry", data=df)

---

3. No Panel IV (FE + IV Combined)

linearmodels provides excellent panel estimators AND excellent IV estimators, but there is no estimator that absorbs fixed effects while simultaneously instrumenting endogenous variables. This is listed as planned but not yet implemented in the linearmodels roadmap.

# There is no way to do this in linearmodels:
# PanelOLS with entity_effects=True AND an instrument for an endogenous regressor

Workaround

Use pyfixest's three-part formula syntax:

import pyfixest as pf

# FE + IV in one call
fit = pf.feols("y ~ exog | entity + year | endog ~ instrument", data=df)

---

4. No DiD Estimators

linearmodels does not provide any difference-in-differences estimators. There is no TWFE convenience wrapper, no did2s, no event study, no local projections DiD, and no Sun-Abraham estimator.

Use pyfixest for all DiD work:

import pyfixest as pf

# TWFE
fit = pf.feols("y ~ treat | entity + year", data=df)

# did2s (Gardner 2022)
fit = pf.did2s(data=df, yname="y", first_stage="~ 0 | entity + year",
               second_stage="~ i(rel_year, ref=-1)", treatment="treat",
               cluster="entity")

---

5. No Automatic Singleton Removal

Unlike pyfixest (which drops singleton fixed-effect groups by default since v0.40 via fixef_rm="singleton"), linearmodels historically did not detect or remove singletons. As of recent versions, PanelOLS now has a singletons parameter (default True) that controls whether singleton groups are kept or dropped. Singleton groups -- entity-time combinations with exactly one observation -- can inflate degrees of freedom and produce misleading inference.

Check for Singletons Manually

# Check for entities with only one time period
group_sizes = df.groupby(level=0).size()
singletons = group_sizes[group_sizes == 1]
print(f"Singleton entities: {len(singletons)}")

# Remove if present
if len(singletons) > 0:
    df = df[~df.index.get_level_values(0).isin(singletons.index)]
    print(f"Removed {len(singletons)} singleton entities")

---

6. No Post-Estimation SE Switching

Unlike pyfixest where fit.vcov("hetero") switches standard errors without re-estimating, linearmodels requires re-fitting the model with a different cov_type argument:

from linearmodels.panel import PanelOLS

mod = PanelOLS.from_formula("y ~ x + EntityEffects", data=df)

# Must call .fit() separately for each SE type
res_robust = mod.fit(cov_type="robust")
res_clustered = mod.fit(cov_type="clustered", cluster_entity=True)
res_kernel = mod.fit(cov_type="kernel")

# pyfixest -- switch SEs post-estimation (no re-fit)
import pyfixest as pf
fit = pf.feols("y ~ x | entity", data=df)
fit.vcov("hetero")
fit.vcov({"CRV1": "entity"})

The computational cost of re-fitting in linearmodels is typically small since the estimation step itself is fast. The inconvenience is mainly syntactic.

---

7. Constant Term Handling

Array API

Like statsmodels, the linearmodels array API does NOT add a constant automatically. Omitting the constant silently fits a model through the origin with no warning.

import statsmodels.api as sm
from linearmodels.iv import IV2SLS

# WRONG -- no intercept (model through origin)
mod = IV2SLS(y, X, endog, instruments)

# CORRECT -- add constant explicitly
mod = IV2SLS(y, sm.add_constant(X), endog, instruments)

Formula API

The formula API adds a constant automatically (consistent with statsmodels formula behavior). To suppress the intercept: "y ~ x - 1".

Panel Models Exception

PanelOLS with entity_effects=True absorbs the intercept into the entity dummies -- you do not need to add or worry about a constant. However, RandomEffects DOES include an explicit intercept, so ensure the constant is present (the "1" in "y ~ 1 + x" via the formula API, or sm.add_constant(X) via the array API).

---

8. Formula Syntax Differences from pyfixest

linearmodels and pyfixest use different formula conventions. Mixing them up produces confusing errors.

Common Mistakes

# WRONG -- pyfixest pipe syntax in linearmodels
mod = PanelOLS.from_formula("y ~ x | entity", data=df)
# This does NOT absorb entity FE; it's parsed as a bitwise OR

# CORRECT -- use EntityEffects keyword
mod = PanelOLS.from_formula("y ~ x + EntityEffects", data=df)

# WRONG -- linearmodels IV bracket syntax in pyfixest
fit = pf.feols("y ~ 1 + exog + [endog ~ inst]", data=df)
# ParseError

# CORRECT -- pyfixest three-part formula
fit = pf.feols("y ~ exog | 0 | endog ~ inst", data=df)

---

9. SUR Performance with Large Datasets

The SUR (Seemingly Unrelated Regression) estimator uses dense matrix operations and does not scale well to large datasets. The cross-equation covariance matrix grows with the number of equations, and the system-level GLS step operates on the full stacked dataset.

Practical Thresholds

Alternative for Large Data

If cross-equation correlation is low, equation-by-equation OLS with clustered standard errors provides similar efficiency gains with much better scaling:

import pyfixest as pf

# Instead of SUR, estimate each equation separately
fit1 = pf.feols("y1 ~ x1 + x2", data=df, vcov={"CRV1": "cluster"})
fit2 = pf.feols("y2 ~ x1 + x3", data=df, vcov={"CRV1": "cluster"})

---

10. Three-Way Package Boundary

linearmodels, pyfixest, and statsmodels have complementary but overlapping capabilities. Choosing the wrong package wastes time and produces suboptimal results.

When to Use linearmodels

• Random effects, between effects, first difference, Fama-MacBeth
• IV without fixed effects (especially LIML, GMM, CUE-GMM)
• System estimation: SUR, 3SLS, system GMM
• Asset pricing factor models (linear factor models, risk premia)
• Driscoll-Kraay standard errors with fine-grained kernel control

When to Use pyfixest Instead

• Any regression with absorbed fixed effects (faster, more convenient)
• FE + IV combined (linearmodels cannot do this)
• Difference-in-differences / event study (linearmodels has none)
• Poisson / GLM with fixed effects (fepois)
• Wild bootstrap, CRV3, randomization inference
• Multiple estimation via sw() / csw()
• Publication tables (etable) and coefficient plots (coefplot)

When to Use statsmodels Instead

• GLM without fixed effects (logit, probit, negative binomial, zero-inflated)
• Time series (ARIMA, VAR, state space, GARCH via arch)
• Diagnostic tests (heteroskedasticity, normality, VIF, influence measures)
• Discrete choice (multinomial logit, ordered probit)
• Mixed / multilevel models (MixedLM)
• Robust regression (M-estimators, quantile regression)
• Prediction intervals and marginal effects

---

Quick Diagnostic Table

---

References and Further Reading

• Sheppard, K. linearmodels documentation: https://bashtage.github.io/linearmodels/
• Sheppard, K. linearmodels GitHub repository: https://github.com/bashtage/linearmodels
• Sheppard, K. linearmodels GitHub issues: https://github.com/bashtage/linearmodels/issues
• Cameron, A.C. and Trivedi, P.K. (2005). Microeconometrics: Methods and Applications. Cambridge University Press. (Panel data and IV theory)
• Wooldridge, J.M. (2010). Econometric Analysis of Cross Section and Panel Data. 2nd ed. MIT Press. (FE, RE, FD, Hausman test)

Instrumental Variable and GMM Models

Reference for all IV estimators in linearmodels.iv. For methodology guidance on when IV is appropriate, instrument validity, and identification, see the data-scientist skill's causal-inference.md reference.

Contents

• Data Requirements
• IV2SLS (Two-Stage Least Squares)
• IVLIML (Limited Information Maximum Likelihood)
• IVGMM (Generalized Method of Moments)
• IVGMMCUE (Continuously Updating Estimator)
• AbsorbingLS
• Diagnostics and Tests
• Covariance Options
• Estimator Selection Guide
• Common Pitfalls
• References and Further Reading

Data Requirements

• IV models do not require a MultiIndex (unlike panel models)
• Accept regular pandas DataFrames
• Terms: dependent, exog (exogenous regressors), endog (endogenous regressors),

instruments (excluded instruments)

• The model automatically distinguishes between included exogenous variables and

excluded instruments

• The order of identification requires at least as many excluded instruments as

endogenous regressors

IV2SLS (Two-Stage Least Squares)

The workhorse IV estimator. Nests OLS when there are no endogenous variables (endog=None). Consistent under heteroskedasticity with robust covariance.

Formula API

from linearmodels.iv import IV2SLS

# Basic IV: educ is endogenous, instrumented by motheduc and fatheduc
mod = IV2SLS.from_formula(
    "np.log(wage) ~ 1 + exper + exper_sq + [educ ~ motheduc + fatheduc]",
    data=df
)
res = mod.fit(cov_type="robust")
print(res.summary)

Array API

from linearmodels.iv import IV2SLS
import statsmodels.api as sm

dependent = df["log_wage"]
exog = sm.add_constant(df[["exper", "exper_sq"]])  # Included exogenous (with constant)
endog = df[["educ"]]                                 # Endogenous variable(s)
instruments = df[["motheduc", "fatheduc"]]           # Excluded instruments

mod = IV2SLS(dependent, exog, endog, instruments)
res = mod.fit(cov_type="robust")

Formula Bracket Syntax

Brackets enclose the endogenous ~ instruments mapping:

# Single endogenous, two instruments
"y ~ 1 + exog1 + exog2 + [endog1 ~ inst1 + inst2]"

# Multiple endogenous variables
"y ~ 1 + exog1 + [endog1 + endog2 ~ inst1 + inst2 + inst3]"

# Constants: include "1" in the exogenous part, NOT inside brackets
"y ~ 1 + exog1 + [endog1 ~ inst1]"        # Correct
"y ~ exog1 + [endog1 ~ 1 + inst1]"        # Wrong — constant goes outside brackets

# Suppress constant
"y ~ exog1 + [endog1 ~ inst1] - 1"

Compare with pyfixest pipe-separated syntax:

# pyfixest: "y ~ exog | 0 | endog ~ inst"
# linearmodels: "y ~ 1 + exog + [endog ~ inst]"

Interpreting IV2SLS Output

res.summary prints multiple sections:

• Parameter estimates: Coefficients, standard errors, t-stats, p-values, CIs
• First-stage diagnostics (per endogenous variable):
• Partial R-squared of excluded instruments
• Partial F-statistic — rule of thumb: F > 10 suggests instruments are

not weak (Staiger-Stock)

• Sargan/Hansen J-test: Test of overidentifying restrictions (only when

# instruments > # endogenous). H0: instruments are valid. Rejection suggests at least one instrument is endogenous.

• Wooldridge regression test for endogeneity: Durbin-Wu-Hausman equivalent.

H0: the endogenous variable is actually exogenous. Failure to reject means OLS may be consistent and IV is unnecessary.

IVLIML (Limited Information Maximum Likelihood)

Alternative to 2SLS with better finite-sample properties. Less biased than 2SLS when instruments are weak. LIML is a k-class estimator: 2SLS uses k=1; LIML uses k=lambda_min (smallest eigenvalue of a particular matrix).

from linearmodels.iv import IVLIML

# Standard LIML
mod = IVLIML.from_formula(
    "np.log(wage) ~ 1 + exper + [educ ~ motheduc + fatheduc]",
    data=df
)
res = mod.fit(cov_type="robust")

Fuller's Modified LIML

Subtracts alpha / (n - L) from k for finite-sample bias correction, where L is the number of instruments. The fuller parameter controls alpha.

# Fuller's modified LIML (alpha=1 is the standard choice)
mod = IVLIML.from_formula(
    "np.log(wage) ~ 1 + exper + [educ ~ motheduc + fatheduc]",
    data=df,
    fuller=1
)
res = mod.fit(cov_type="robust")

When to Use LIML vs 2SLS

• Many instruments: LIML is less biased (2SLS bias is proportional to

# instruments / n)

• Weak instruments: LIML is more robust (but still problematic with very weak

instruments)

• Fuller's LIML with alpha=1: approximately unbiased with very small efficiency

loss — a good default when instrument strength is uncertain

• Single endogenous variable + single instrument: 2SLS and LIML are numerically

identical (k=1 in both cases)

IVGMM (Generalized Method of Moments)

Most efficient IV estimator when instruments are valid and the model is overidentified. Estimates an optimal weight matrix from a first-step estimation, then re-estimates.

from linearmodels.iv import IVGMM

# 2-step efficient GMM with robust weight matrix
mod = IVGMM.from_formula(
    "np.log(wage) ~ 1 + exper + [educ ~ motheduc + fatheduc]",
    data=df
)
res = mod.fit(cov_type="robust", weight_type="robust")

Weight Matrix Options

Steps Parameter

# 1-step GMM (uses initial weight matrix only)
res = mod.fit(steps=1, weight_type="robust")

# 2-step GMM (default — efficient under correct specification)
res = mod.fit(steps=2, weight_type="robust")

# Iterative GMM (repeats until weight matrix converges)
res = mod.fit(iter_limit=100, weight_type="robust")

When to Use GMM vs 2SLS

• Just-identified (# instruments = # endogenous): 2SLS and efficient GMM produce

identical estimates. Use 2SLS for simplicity.

• Overidentified with heteroskedasticity: GMM with weight_type="robust" is more

efficient than 2SLS.

• Overidentified with homoskedasticity: 2SLS is already efficient; GMM gains nothing.

IVGMMCUE (Continuously Updating Estimator)

Jointly optimizes parameters and weight matrix in a single objective function (non-linear optimization). Can be more robust to weak instruments than 2-step GMM, but is more computationally expensive and may not converge.

from linearmodels.iv import IVGMMCUE

mod = IVGMMCUE.from_formula(
    "np.log(wage) ~ 1 + exper + [educ ~ motheduc + fatheduc]",
    data=df
)
res = mod.fit(cov_type="robust")

Use CUE when:

• Concerned about finite-sample bias from sequential 2-step estimation
• Moderate number of moment conditions (CUE becomes unstable with many)
• Willing to accept higher computation time and possible convergence issues

AbsorbingLS

OLS/WLS with high-dimensional absorbed fixed effects. Uses pyhdfe for FE absorption (similar to Stata's reghdfe or pyfixest). Lives in the IV module but is fundamentally an OLS estimator — does not support endogenous variables.

from linearmodels.iv import AbsorbingLS

# Absorb entity and year fixed effects
mod = AbsorbingLS(
    dependent=df["y"],
    exog=df[["x1", "x2"]],
    absorb=df[["entity_id", "year"]]
)
res = mod.fit()

Interacted Fixed Effects

from linearmodels.iv.absorbing import Interaction

# Entity-by-industry interacted FE
interact = Interaction(df["entity_id"], df["industry"])
mod = AbsorbingLS(
    dependent=df["y"],
    exog=df[["x1"]],
    absorb=interact
)
res = mod.fit()

AbsorbingLS vs PanelOLS vs pyfixest

Recommendation: For most absorbed-FE work, prefer pyfixest for speed and formula convenience. Use AbsorbingLS when you need to stay within the linearmodels ecosystem or need the Interaction class for complex FE structures that are not easily expressed in pyfixest's ^ syntax.

Diagnostics and Tests

First-Stage Statistics

res = mod.fit(cov_type="robust")

# Per-endogenous-variable first-stage results
for name, fs in res.first_stage.individual.items():
    print(f"--- First stage for {name} ---")
    print(f"  Partial R²: {fs.rsquared:.4f}")
    print(f"  Partial F:  {fs.f_statistic.stat:.2f} (p={fs.f_statistic.pval:.4f})")

Overidentification Test (Sargan / Hansen J-Test)

Available only when # instruments > # endogenous:

# Sargan test (assumes homoskedastic errors)
print(f"Sargan J-stat:  {res.sargan.stat:.3f}")
print(f"Sargan p-value: {res.sargan.pval:.3f}")

# Wooldridge overidentification test (robust to heteroskedasticity)
print(f"Wooldridge overid stat:  {res.wooldridge_overid.stat:.3f}")
print(f"Wooldridge overid p-val: {res.wooldridge_overid.pval:.3f}")

Endogeneity Test (Wu-Hausman)

Tests H0: the suspected endogenous variable is actually exogenous.

wh = res.wu_hausman()
print(f"Wu-Hausman F-stat: {wh.stat:.3f}")
print(f"Wu-Hausman p-val:  {wh.pval:.3f}")

Wooldridge Regression-Based Endogeneity Test

endo_test = res.wooldridge_regression
print(f"Wooldridge stat: {endo_test.stat:.3f}")
print(f"Wooldridge pval: {endo_test.pval:.3f}")

Weak Instrument Warning Signs

• First-stage partial F < 10 (Staiger-Stock rule of thumb)
• Large difference between 2SLS and LIML estimates (suggests weak instrument bias

in 2SLS — LIML is median-unbiased, so divergence from 2SLS indicates 2SLS bias)

• Wide confidence intervals relative to OLS
• Anderson-Rubin confidence sets (not built in — construct manually or use

alternative packages if needed)

Covariance Options

All IV estimators accept cov_type in .fit():

# Robust (heteroskedasticity-consistent)
res = mod.fit(cov_type="robust")

# Clustered by state
res = mod.fit(cov_type="clustered", clusters=df["state"])

# HAC with Bartlett kernel, bandwidth 5
res = mod.fit(cov_type="kernel", kernel="bartlett", bandwidth=5)

For GMM models, weight_type and cov_type are separate arguments. The weight matrix determines efficiency; the covariance type determines inference. They should generally match (e.g., weight_type="robust" with cov_type="robust").

Estimator Selection Guide

How to choose an IV estimator?
├─ Just-identified (# inst = # endog)
│   ├─ All estimators give identical point estimates
│   └─ Use IV2SLS for simplicity
├─ Overidentified (# inst > # endog)
│   ├─ Instruments strong (F > 10)?
│   │   ├─ Yes, homoskedastic → IV2SLS
│   │   ├─ Yes, heteroskedastic → IVGMM (2-step, weight_type="robust")
│   │   └─ Yes, clustered → IVGMM (2-step, weight_type="clustered")
│   └─ Instruments possibly weak (F ~ 10)?
│       ├─ Few instruments → IVLIML or Fuller (fuller=1)
│       └─ Many instruments → Fuller (fuller=1) preferred
├─ Concerned about finite-sample bias?
│   ├─ IVLIML (less biased than 2SLS)
│   ├─ Fuller LIML (approximately unbiased)
│   └─ IVGMMCUE (if willing to accept computation cost)
└─ Just want a robust default?
    └─ IV2SLS with cov_type="robust" (widely understood, easy to report)

Common Pitfalls

Constant term placement: The constant (1) goes in the exogenous part, never inside the brackets.

# Correct
"y ~ 1 + exper + [educ ~ motheduc + fatheduc]"

# Wrong — do not put constant inside brackets
"y ~ exper + [educ ~ 1 + motheduc + fatheduc]"

Forgetting to check first-stage strength: Always inspect res.first_stage before interpreting second-stage results. Weak first stages invalidate standard inference even with large samples.

Using Sargan test with heteroskedastic data: The Sargan test assumes homoskedasticity. Use res.wooldridge_overid for heteroskedasticity-robust overidentification testing.

Confusing weight_type and cov_type in GMM: weight_type determines the GMM weight matrix (affects point estimates in overidentified models). cov_type determines standard errors (affects inference only). They should generally be set to the same value.

Mixing up with panel IV: linearmodels does not have a combined panel-IV estimator. If you need IV with entity/time fixed effects, use pyfixest:

# pyfixest: IV with entity and time FE
import pyfixest as pf
fit = pf.feols("log_wage ~ 1 + exper | entity + year | educ ~ motheduc + fatheduc", data=df)

References and Further Reading

• Sheppard, K. linearmodels IV documentation. https://bashtage.github.io/linearmodels/iv/introduction.html
• Angrist, J.D. and Pischke, J.-S. (2009). Mostly Harmless Econometrics. Princeton University Press.
• Baum, C.F., Schaffer, M.E., and Stillman, S. (2007). "Enhanced Routines for Instrumental Variables/Generalized Method of Moments Estimation and Testing." Stata Journal, 7(4), 465-506.
• Stock, J.H. and Yogo, M. (2005). "Testing for Weak Instruments in Linear IV Regression." In Andrews, D.W.K. and Stock, J.H. (eds.), Identification and Inference for Econometric Models. Cambridge University Press.
• Hansen, L.P. (1982). "Large Sample Properties of Generalized Method of Moments Estimators." Econometrica, 50(4), 1029-1054.

Panel Data Models

Reference for all panel estimators in linearmodels.panel. For methodology guidance on when to use FE vs RE, assumption checking, and model selection, see the data-scientist skill's statistical-modeling.md reference.

Contents

• Data Requirements
• PanelOLS (Fixed Effects)
• RandomEffects
• BetweenOLS
• FirstDifferenceOLS
• PooledOLS
• FamaMacBeth
• Weighted Estimation
• Model Comparison with compare()
• Common Patterns
• References and Further Reading

Data Requirements

All panel models require a pandas DataFrame with a two-level MultiIndex of (entity, time). Balanced panels are not required -- all estimators handle unbalanced panels. Observations with NaN in the dependent or any independent variable are dropped automatically.

Standard Setup Pattern

import pandas as pd
from linearmodels.panel import PanelOLS

# Load data (Grunfeld investment dataset as example)
df = pd.read_parquet("data/raw/grunfeld.parquet")

# Set MultiIndex — entity first, time second
df = df.set_index(["firm_id", "year"])

# Verify structure before estimation
print(f"Index names: {df.index.names}")      # ['firm_id', 'year']
print(f"Index levels: {df.index.nlevels}")    # 2
print(f"Entities: {df.index.get_level_values(0).nunique()}")
print(f"Periods: {df.index.get_level_values(1).nunique()}")
print(f"Total obs: {len(df)}")

# Check balance
obs_per_entity = df.groupby(level=0).size()
print(f"Balanced: {obs_per_entity.nunique() == 1}")
print(f"Min periods per entity: {obs_per_entity.min()}")
print(f"Max periods per entity: {obs_per_entity.max()}")

PanelOLS (Fixed Effects)

Within estimator that absorbs entity and/or time fixed effects by demeaning. This is the workhorse panel model for controlling unobserved heterogeneity.

Effect Types

PanelOLS supports a maximum of 2 sets of effects. For 3+ way FE, use pyfixest.

Formula API

from linearmodels.panel import PanelOLS

# Entity fixed effects
mod = PanelOLS.from_formula("invest ~ value + capital + EntityEffects", data=df)
res = mod.fit()
print(res.summary)

# Time fixed effects only
mod = PanelOLS.from_formula("invest ~ value + capital + TimeEffects", data=df)

# Two-way fixed effects (entity + time)
mod = PanelOLS.from_formula(
    "invest ~ value + capital + EntityEffects + TimeEffects", data=df
)

Array API

mod = PanelOLS(
    dependent=df["invest"],
    exog=df[["value", "capital"]],
    entity_effects=True,
    time_effects=True,
)
res = mod.fit()

Other Effects (Arbitrary Groupings)

# Absorb a third grouping variable (array API only)
mod = PanelOLS(
    dependent=df["invest"],
    exog=df[["value", "capital"]],
    entity_effects=True,
    other_effects=df["industry_code"],
)
res = mod.fit()

check_rank Parameter

# Disable rank check when you know the model is identified
# Useful with large numbers of effects where rank checking is slow
mod = PanelOLS.from_formula(
    "invest ~ value + capital + EntityEffects", data=df, check_rank=False
)

Interpreting PanelOLS Output

The results object provides three R-squared measures:

Key output elements:

• F-test for poolability: Tests the null that all entity effects are jointly zero. Rejection means entity effects matter and pooled OLS is inappropriate.
• Estimated parameters do NOT include entity/time dummies -- they are absorbed via demeaning.
• Degrees of freedom account for absorbed effects.

res = mod.fit()

# Access fit statistics
print(f"Within R-sq: {res.rsquared_within:.4f}")
print(f"Between R-sq: {res.rsquared_between:.4f}")
print(f"Overall R-sq: {res.rsquared_overall:.4f}")

# F-test for poolability (entity effects jointly zero)
print(f"F-stat (poolable): {res.f_pooled.stat:.4f}")
print(f"F p-value: {res.f_pooled.pval:.4f}")

Extracting Fixed Effects

res = mod.fit()

# Estimated entity/time effects
effects = res.estimated_effects
print(effects.head(10))

# Merge effects back for inspection
# Effects are indexed by the same MultiIndex as the original data
print(f"Shape: {effects.shape}")
print(f"Mean effect (should be ~0 for entity FE): {effects.mean().values}")

RandomEffects

GLS estimator with quasi-demeaning. Assumes entity effects are uncorrelated with regressors -- more efficient than FE when this assumption holds, inconsistent when it does not.

Usage

from linearmodels.panel import RandomEffects

# RE requires an intercept — include 1 in the formula
mod = RandomEffects.from_formula("invest ~ 1 + value + capital", data=df)
res = mod.fit()
print(res.summary)

Array API

from linearmodels.panel import RandomEffects
import numpy as np

# Must add constant column manually in array API
exog = df[["value", "capital"]].copy()
exog.insert(0, "const", 1.0)

mod = RandomEffects(dependent=df["invest"], exog=exog)
res = mod.fit()

Variance Decomposition

# RE decomposes total variance into entity and idiosyncratic components
print(f"Sigma^2_entity: {res.variance_decomposition.Effects:.4f}")
print(f"Sigma^2_idiosyncratic: {res.variance_decomposition.Residual:.4f}")

# Theta: quasi-demeaning parameter (0 = pooled OLS, 1 = full within)
# Closer to 1 means more within-entity variation, RE approaches FE
print(f"Theta: {res.theta.iloc[0]:.4f}")

FE vs RE Comparison

linearmodels does NOT have a built-in Hausman test. Compare FE and RE results using the compare() function and inspect whether coefficients diverge substantively. For the methodology behind the FE vs RE decision, see the data-scientist skill's statistical-modeling.md.

from linearmodels.panel import PanelOLS, RandomEffects, compare

fe_mod = PanelOLS.from_formula(
    "invest ~ value + capital + EntityEffects", data=df
)
re_mod = RandomEffects.from_formula("invest ~ 1 + value + capital", data=df)

fe_res = fe_mod.fit()
re_res = re_mod.fit()

# Side-by-side comparison table
comp = compare({"Fixed Effects": fe_res, "Random Effects": re_res})
print(comp.summary)

# Manual coefficient comparison for Hausman-style inspection
fe_params = fe_res.params
re_params = re_res.params[fe_params.index]  # Align on common regressors
diff = fe_params - re_params
print("Coefficient differences (FE - RE):")
print(diff)
print("Large differences suggest RE assumption is violated → use FE")

BetweenOLS

Regresses entity time-averages on entity-averaged regressors. Uses only between-entity (cross-sectional) variation. Useful as a diagnostic or as part of variance decomposition analysis.

from linearmodels.panel import BetweenOLS

# Between estimator requires an intercept
mod = BetweenOLS.from_formula("invest ~ 1 + value + capital", data=df)
res = mod.fit()
print(res.summary)

# Number of effective observations = number of entities
print(f"Entities used: {res.nobs}")

FirstDifferenceOLS

Eliminates entity effects by first-differencing: delta_y_it = delta_x_it * beta + delta_epsilon_it. Alternative to within-estimation that can handle certain forms of serial correlation better. Loses one time period per entity.

from linearmodels.panel import FirstDifferenceOLS

# No intercept or effects keywords — differencing removes both
mod = FirstDifferenceOLS.from_formula("invest ~ value + capital", data=df)
res = mod.fit()
print(res.summary)

# First-difference requires at least 2 time periods per entity
# Entities with only 1 period are dropped automatically
print(f"Observations used: {res.nobs}")

PooledOLS

Standard OLS ignoring panel structure, but with panel-aware covariance estimation. Entity and time identifiers are used for clustering standard errors, not for absorbing effects. Useful as a baseline model.

from linearmodels.panel import PooledOLS

# Include intercept explicitly
mod = PooledOLS.from_formula("invest ~ 1 + value + capital", data=df)

# Without clustered SEs — equivalent to plain OLS
res_naive = mod.fit()

# With entity-clustered SEs — accounts for within-entity correlation
res_clustered = mod.fit(cov_type="clustered", cluster_entity=True)
print(res_clustered.summary)

FamaMacBeth

Two-step procedure: (1) run a cross-sectional regression for each time period, (2) average coefficients across periods. Standard approach in empirical asset pricing for testing factor models. See ./asset-pricing.md for factor model application context.

from linearmodels.panel import FamaMacBeth

# Typical asset pricing specification
mod = FamaMacBeth.from_formula("ret ~ 1 + beta + size + bm", data=df)

# Basic Fama-MacBeth SEs (assumes independence across periods)
res = mod.fit()
print(res.summary)

# HAC-adjusted SEs for time-series dependence in averaged coefficients
res_hac = mod.fit(cov_type="kernel", kernel="bartlett", bandwidth=5)
print(res_hac.summary)

# Access time-period-by-period coefficient estimates
# Not directly available in results — re-estimate manually if needed

Weighted Estimation

All panel models accept a weights parameter for WLS-style estimation. Weights must have the same shape as the dependent variable (one weight per observation).

# Weighted fixed effects
mod = PanelOLS.from_formula(
    "invest ~ value + capital + EntityEffects",
    data=df,
    weights=df["weight"],
)
res = mod.fit()

# Weighted random effects
mod = RandomEffects.from_formula(
    "invest ~ 1 + value + capital",
    data=df,
    weights=df["weight"],
)
res = mod.fit()

Common weight sources:

• Population weights: weight by group size for representative estimates
• Precision weights: inverse of known variance for heteroskedastic data
• Frequency weights: when data is pre-aggregated with counts

Model Comparison with compare()

The compare() function produces a side-by-side table showing coefficients, standard errors, and fit statistics for all models. It is the linearmodels equivalent of pyfixest's etable().

from linearmodels.panel import (
    PanelOLS, RandomEffects, BetweenOLS,
    FirstDifferenceOLS, PooledOLS, compare,
)

# Estimate all models on the same specification
pooled_res = PooledOLS.from_formula(
    "invest ~ 1 + value + capital", data=df
).fit(cov_type="clustered", cluster_entity=True)

fe_res = PanelOLS.from_formula(
    "invest ~ value + capital + EntityEffects", data=df
).fit()

re_res = RandomEffects.from_formula(
    "invest ~ 1 + value + capital", data=df
).fit()

between_res = BetweenOLS.from_formula(
    "invest ~ 1 + value + capital", data=df
).fit()

fd_res = FirstDifferenceOLS.from_formula(
    "invest ~ value + capital", data=df
).fit()

# Compare all
comp = compare({
    "Pooled": pooled_res,
    "FE": fe_res,
    "RE": re_res,
    "Between": between_res,
    "FD": fd_res,
})
print(comp.summary)

Common Patterns

Typical Panel Analysis Workflow

1. Estimate PooledOLS as baseline 2. Estimate PanelOLS (FE) -- check if entity effects matter via F-test for poolability 3. Estimate RandomEffects -- compare coefficients with FE 4. If FE and RE coefficients diverge substantially, use FE (RE assumption likely violated) 5. Consider FirstDifferenceOLS as a robustness check 6. Report preferred specification with appropriate clustered SEs

Adding Time Trends Instead of Time FE

# Linear time trend as a regressor (not as absorbed time effects)
# Useful when you want to model the trend parametrically
df["trend"] = (
    df.index.get_level_values("year")
    - df.index.get_level_values("year").min()
)

mod = PanelOLS.from_formula(
    "invest ~ value + capital + trend + EntityEffects", data=df
)
res = mod.fit()

Clustered Standard Errors with Panel Models

# Entity-clustered (default clustering dimension for panel data)
res = mod.fit(cov_type="clustered", cluster_entity=True)

# Time-clustered
res = mod.fit(cov_type="clustered", cluster_time=True)

# Two-way clustering (entity and time)
res = mod.fit(cov_type="clustered", cluster_entity=True, cluster_time=True)

For full covariance estimation options including Driscoll-Kraay and HAC kernels, see ./covariance-inference.md.

Accessing Results Programmatically

res = mod.fit()

# Coefficients and standard errors
print(res.params)           # pd.Series of point estimates
print(res.std_errors)       # pd.Series of standard errors
print(res.pvalues)          # pd.Series of p-values
print(res.conf_int())       # pd.DataFrame with 95% CI (default)
print(res.conf_int(0.90))   # 90% CI

# Residuals and fitted values
print(res.resids.head())    # Residuals
print(res.fitted_values.head())

# Number of observations and entities
print(f"N obs: {res.nobs}")
print(f"N entities: {res.entity_info['total']}")
print(f"N time periods: {res.time_info['total']}")

References and Further Reading

• Sheppard, K. linearmodels Panel Models documentation. https://bashtage.github.io/linearmodels/panel/introduction.html
• Wooldridge, J.M. (2010). Econometric Analysis of Cross Section and Panel Data. 2nd ed. MIT Press.
• Baltagi, B.H. (2021). Econometric Analysis of Panel Data. 6th ed. Springer.
• Mundlak, Y. (1978). "On the Pooling of Time Series and Cross Section Data." Econometrica, 46(1), 69-85.

linearmodels Quickstart

A technical reference for using linearmodels v7.0 in Python. Covers installation, the MultiIndex requirement, formula vs array APIs, first panel and IV models, reading output, and syntax comparison with pyfixest and statsmodels. No methodology content -- syntax and library guidance only.

Contents

• Installation
• The MultiIndex Requirement (CRITICAL)
• Formula API vs Array API
• Your First Panel Model
• Reading the Summary Output
• Key Results Attributes
• Your First IV Model
• Standard Errors
• Model Comparison Tables
• Quick Syntax Comparison: linearmodels vs pyfixest vs statsmodels
• References and Further Reading

---

Installation

pip install linearmodels

Key dependencies (installed automatically):

Optional dependencies:

Verify installation:

import linearmodels
print(linearmodels.__version__)  # Should print 7.x

Requires Python 3.10+.

---

The MultiIndex Requirement (CRITICAL)

This is the single most important thing to know about linearmodels. All panel models (PanelOLS, RandomEffects, BetweenOLS, FirstDifferenceOLS, PooledOLS, FamaMacBeth) require a pandas DataFrame with a two-level MultiIndex where:

• Level 0 = entity (firm, person, state, school, etc.)
• Level 1 = time (year, quarter, date, etc.)

Converting from a Flat DataFrame

import pandas as pd

# Flat data with entity and time as regular columns
df = pd.read_parquet("panel_data.parquet")
print(df.columns.tolist())
# ['firm_id', 'year', 'invest', 'value', 'capital']

# Set the MultiIndex — entity first, time second
df = df.set_index(["firm_id", "year"])

# Verify
print(f"Index names: {df.index.names}")    # ['firm_id', 'year']
print(f"Index levels: {df.index.nlevels}")  # 2

Converting from Polars

import polars as pl

df_polars = pl.read_parquet("panel_data.parquet")
df = df_polars.to_pandas().set_index(["firm_id", "year"])

What Happens If You Forget

Passing a DataFrame without a MultiIndex raises an error:

ValueError: Panel models require a MultiIndex with 2 levels that corresponds
to entities and time periods

If you see this error, set the index before passing data to any panel model.

Important: IV Models Do NOT Require MultiIndex

The IV2SLS, IVLIML, IVGMM, and IVGMMCUE classes work with regular (non-panel) DataFrames. Only panel models require the MultiIndex.

---

Formula API vs Array API

linearmodels exposes two parallel interfaces for every model class.

Formula API

Uses the .from_formula() class method. Formulas are parsed by formulaic (not patsy — formulaic is the modern replacement used by both linearmodels and pyfixest).

from linearmodels.panel import PanelOLS

mod = PanelOLS.from_formula("invest ~ value + capital + EntityEffects", data=df)
res = mod.fit()

Array API

Passes dependent variable, exogenous regressors, and optional parameters directly to the constructor.

from linearmodels.panel import PanelOLS

mod = PanelOLS(df["invest"], df[["value", "capital"]], entity_effects=True)
res = mod.fit()

Same Model, Both Ways

import pandas as pd
from linearmodels.panel import PanelOLS

df = pd.read_parquet("grunfeld.parquet")
df = df.set_index(["firm_id", "year"])

# --- Formula API ---
res_f = PanelOLS.from_formula(
    "invest ~ value + capital + EntityEffects", data=df
).fit(cov_type="clustered", cluster_entity=True)

# --- Array API ---
res_a = PanelOLS(
    df["invest"], df[["value", "capital"]], entity_effects=True
).fit(cov_type="clustered", cluster_entity=True)

# Both produce identical coefficient estimates
print(res_f.params)
print(res_a.params)

When to Prefer Each

Formula Syntax Unique to linearmodels

linearmodels formulas differ from statsmodels and pyfixest in two key ways:

1. Fixed effects as keywords in the formula (not after `|`)

# linearmodels — EntityEffects and TimeEffects are keywords IN the formula
"invest ~ value + capital + EntityEffects"
"invest ~ value + capital + EntityEffects + TimeEffects"

# pyfixest — fixed effects go AFTER the pipe
"invest ~ value + capital | firm_id"
"invest ~ value + capital | firm_id + year"

2. IV uses bracket notation (not pipe notation)

# linearmodels — endogenous variable and instruments in brackets
"np.log(wage) ~ 1 + exper + exper_sq + [educ ~ motheduc + fatheduc]"

# pyfixest — endogenous and instruments after second pipe
"np.log(wage) ~ 1 + exper + exper_sq | 0 | educ ~ motheduc + fatheduc"

---

Your First Panel Model

PanelOLS with Entity Effects (Formula)

import pandas as pd
from linearmodels.panel import PanelOLS

# Load and set MultiIndex
df = pd.read_parquet("grunfeld.parquet")
df = df.set_index(["firm_id", "year"])

# Estimate with entity fixed effects
mod = PanelOLS.from_formula(
    "invest ~ value + capital + EntityEffects", data=df
)
res = mod.fit(cov_type="clustered", cluster_entity=True)
print(res.summary)

PanelOLS with Entity Effects (Array)

mod = PanelOLS(
    dependent=df["invest"],
    exog=df[["value", "capital"]],
    entity_effects=True,
)
res = mod.fit(cov_type="clustered", cluster_entity=True)
print(res.summary)

Two-Way Fixed Effects

# Formula: add both keywords
mod = PanelOLS.from_formula(
    "invest ~ value + capital + EntityEffects + TimeEffects", data=df
)
res = mod.fit(cov_type="clustered", cluster_entity=True)

# Array: set both flags
mod = PanelOLS(
    df["invest"], df[["value", "capital"]],
    entity_effects=True, time_effects=True,
)
res = mod.fit(cov_type="clustered", cluster_entity=True)

---

Reading the Summary Output

res.summary (a property, not a method — no parentheses) prints a multi-panel text table.

Panel 1: Model Information

                          PanelOLS Estimation Summary
================================================================================
Dep. Variable:              invest   R-squared:                        0.7668
Estimator:                PanelOLS   R-squared (Between):              0.8194
No. Observations:              200   R-squared (Within):               0.7668
Date:                ...            R-squared (Overall):              0.8132
Time:                ...            Log-likelihood                   -1070.2
Cov. Estimator:          Clustered
                                     F-statistic:                      324.86
Entities:                       10   P-value                           0.0000
Avg Obs:                     20.00   Distribution:                  F(2,188)
Min Obs:                     20.00
Max Obs:                     20.00   F-statistic (robust):             160.84
                                     P-value                           0.0000
Time periods:                   20   Distribution:                  F(2,188)
Avg Obs:                     10.00
Min Obs:                     10.00
Max Obs:                     10.00

Key fields:

Panel 2: Coefficient Table

                Parameter Estimates
============================================
            Parameter  Std. Err.  T-stat  P-value  Lower CI  Upper CI
--------------------------------------------
value          0.1101    0.0155   7.1073   0.0000    0.0796    0.1406
capital        0.3100    0.0530   5.8488   0.0000    0.2055    0.4146

Entity effects are absorbed — they do not appear as coefficients. The F-test for effects tests whether entity effects are jointly significant.

Panel 3: Additional Statistics

F-test for Poolability:      49.177
P-value:                     0.0000
Distribution:                F(9,188)

The poolability F-test (also called the F-test for fixed effects) tests the null hypothesis that all entity intercepts are equal. A small p-value supports using fixed effects over pooled OLS.

---

Key Results Attributes

All attributes are available on the results object returned by .fit().

# Direct attribute access
print(res.params)
print(res.pvalues)
print(res.rsquared)
print(res.conf_int(level=0.99))  # 99% CI

---

Your First IV Model

IV models (IV2SLS, IVLIML, IVGMM) do not require a MultiIndex — they work with regular DataFrames on cross-sectional or any flat data.

IV2SLS (Formula)

from linearmodels.iv import IV2SLS

# Bracket notation: [endogenous ~ instruments]
mod = IV2SLS.from_formula(
    "np.log(wage) ~ 1 + exper + exper_sq + [educ ~ motheduc + fatheduc]",
    data=df,
)
res = mod.fit(cov_type="robust")
print(res.summary)

IV2SLS (Array)

from linearmodels.iv import IV2SLS
import numpy as np

mod = IV2SLS(
    dependent=np.log(df["wage"]),
    exog=df[["const", "exper", "exper_sq"]],  # exogenous regressors (include constant)
    endog=df[["educ"]],                        # endogenous regressor
    instruments=df[["motheduc", "fatheduc"]],  # excluded instruments
)
res = mod.fit(cov_type="robust")

IV Formula Notes

• The 1 in the formula adds an intercept (constant term)
• Brackets [endog ~ instruments] contain the endogenous variable on the left

and excluded instruments on the right

• Multiple endogenous variables: [endog1 + endog2 ~ inst1 + inst2 + inst3]
• Exogenous variables outside brackets are included in both stages automatically

---

Standard Errors

Standard errors are specified at .fit() time via cov_type and related keyword arguments.

Panel Models

from linearmodels.panel import PanelOLS

mod = PanelOLS.from_formula("invest ~ value + capital + EntityEffects", data=df)

# Unadjusted (homoskedastic) — default
res = mod.fit(cov_type="unadjusted")

# Heteroskedasticity-robust
res = mod.fit(cov_type="robust")

# Clustered by entity
res = mod.fit(cov_type="clustered", cluster_entity=True)

# Clustered by time
res = mod.fit(cov_type="clustered", cluster_time=True)

# Two-way clustered (entity and time)
res = mod.fit(cov_type="clustered", cluster_entity=True, cluster_time=True)

# Driscoll-Kraay (kernel-based, robust to cross-sectional dependence)
res = mod.fit(cov_type="kernel", kernel="bartlett", bandwidth=5)

IV Models

from linearmodels.iv import IV2SLS

mod = IV2SLS.from_formula(
    "np.log(wage) ~ 1 + exper + [educ ~ motheduc + fatheduc]", data=df
)

# Unadjusted — default
res = mod.fit(cov_type="unadjusted")

# Heteroskedasticity-robust
res = mod.fit(cov_type="robust")

# Kernel (HAC)
res = mod.fit(cov_type="kernel", kernel="bartlett", bandwidth=4)

# Clustered
res = mod.fit(cov_type="clustered", clusters=df["state"])

Quick SE Reference

---

Model Comparison Tables

linearmodels provides a compare() function for side-by-side model comparison.

from linearmodels.panel import PanelOLS, RandomEffects, BetweenOLS, PooledOLS, compare

df = df.set_index(["firm_id", "year"])

# Estimate multiple specifications
res_pooled = PooledOLS.from_formula("invest ~ 1 + value + capital", data=df).fit()
res_fe = PanelOLS.from_formula(
    "invest ~ value + capital + EntityEffects", data=df
).fit(cov_type="clustered", cluster_entity=True)
res_re = RandomEffects.from_formula(
    "invest ~ 1 + value + capital", data=df
).fit()
res_be = BetweenOLS.from_formula(
    "invest ~ 1 + value + capital", data=df
).fit()

# Compare in a single table
comp = compare({
    "Pooled": res_pooled,
    "FE": res_fe,
    "RE": res_re,
    "Between": res_be,
})
print(comp.summary)

The comparison table shows coefficients, standard errors, and fit statistics across models in aligned columns — useful for assessing sensitivity to estimation approach.

---

Quick Syntax Comparison: linearmodels vs pyfixest vs statsmodels

When to Use Which

Use linearmodels when:

• You need random effects, between estimation, or first difference estimators
• Running Fama-MacBeth cross-sectional regressions
• Estimating LIML, GMM-IV, or continuously updating GMM
• Building system models (SUR, 3SLS, system GMM)
• You need Driscoll-Kraay standard errors for panel data
• Asset pricing factor model tests

Use pyfixest when:

• High-dimensional fixed effects (3+ way FE, interacted FE)
• IV combined with fixed effects (linearmodels cannot do this)
• Difference-in-differences or event study designs
• Wild bootstrap or CRV3 clustered standard errors
• Publication-quality tables and coefficient plots

Use statsmodels when:

• GLMs (logit, probit, Poisson, negative binomial)
• Time series models (ARIMA, VAR, state space)
• Diagnostic tests (heteroskedasticity, normality, autocorrelation)
• Discrete choice models
• No panel structure or fixed effects needed

---

References and Further Reading

• Sheppard, K. linearmodels documentation: https://bashtage.github.io/linearmodels/
• Sheppard, K. linearmodels GitHub: https://github.com/bashtage/linearmodels
• formulaic documentation: https://matthewwardrop.github.io/formulaic/
• pyhdfe documentation: https://github.com/jeffgortmaker/pyhdfe

System Regression Models

Reference for system estimators in linearmodels.system: SUR, IV3SLS, IVSystemGMM, and cross-equation constraints. For methodology guidance on simultaneous equations and identification, see the data-scientist skill's statistical-modeling.md reference.

Contents

• Overview
• Data Setup
• SUR (Seemingly Unrelated Regression)
• IV3SLS (Three-Stage Least Squares)
• IVSystemGMM
• Cross-Equation Constraints (LinearConstraint)
• Interpreting System Output
• Performance Notes
• References and Further Reading

Overview

System models estimate multiple equations jointly, exploiting cross-equation error correlation for efficiency gains over equation-by-equation estimation.

Three estimators handle different settings:

Key differences from panel and IV models:

• Input format: dict or OrderedDict of equations, not a single formula
• No MultiIndex required: data is a flat DataFrame (or arrays per equation)
• No `.from_formula()` on the constructor: use a dict of formula strings instead
• Cross-equation covariance (sigma matrix) is estimated and used for GLS

Data Setup

Each equation is defined as a dict with keys for its components. Equations are collected into an OrderedDict (or regular dict on Python 3.7+, though OrderedDict makes equation ordering explicit).

Array Input

from collections import OrderedDict
from linearmodels.system import SUR
import statsmodels.api as sm

# Build equation specifications
equations = OrderedDict()
equations["earnings"] = {
    "dependent": data["hrearn"],
    "exog": data[["const", "exper", "tenure"]],
}
equations["benefits"] = {
    "dependent": data["hrbens"],
    "exog": data[["const", "exper", "union"]],
}

mod = SUR(equations)
res = mod.fit(cov_type="robust")
print(res.summary)

Each equation dict accepts these keys:

Formula Input

Pass a dict of formula strings to .from_formula():

from linearmodels.system import SUR

formulas = {
    "earnings": "hrearn ~ 1 + exper + tenure",
    "benefits": "hrbens ~ 1 + exper + union",
}
mod = SUR.from_formula(formulas, data=data)
res = mod.fit()
print(res.summary)

For IV formulas, use the same bracket syntax as single-equation IV models:

from linearmodels.system import IV3SLS

formulas = {
    "demand": "quantity ~ 1 + income + [price ~ cost + weather]",
    "supply": "quantity ~ 1 + cost + [price ~ income + weather]",
}
mod = IV3SLS.from_formula(formulas, data=data)
res = mod.fit(cov_type="robust")

SUR (Seemingly Unrelated Regression)

GLS estimator that exploits cross-equation error correlation for efficiency. All regressors must be exogenous. More efficient than equation-by-equation OLS when (a) errors are correlated across equations AND (b) regressors differ across equations. If all equations share identical regressors, SUR reduces to OLS (Zellner's invariance result).

from collections import OrderedDict
from linearmodels.system import SUR

equations = OrderedDict()
equations["earnings"] = {
    "dependent": data["hrearn"],
    "exog": data[["const", "exper", "tenure"]],
}
equations["benefits"] = {
    "dependent": data["hrbens"],
    "exog": data[["const", "exper", "union"]],
}

mod = SUR(equations)
res = mod.fit(cov_type="robust")
print(res.summary)

# System R-squared
print(f"System R-squared: {res.rsquared:.4f}")

# Cross-equation covariance (sigma) matrix
print("Sigma matrix:")
print(res.sigma)

Iterative SUR (FGLS)

Iterative SUR re-estimates the sigma matrix and GLS coefficients until convergence. Produces maximum likelihood estimates under normality.

# Iterative SUR — set iter_limit > 1
res_isur = mod.fit(iter_limit=100, cov_type="robust")
print(f"Iterations: {res_isur.iterations}")

Covariance Options

# Unadjusted (homoskedastic)
res = mod.fit(cov_type="unadjusted")

# Heteroskedasticity-robust
res = mod.fit(cov_type="robust")

# HAC (kernel-based)
res = mod.fit(cov_type="kernel", kernel="bartlett", bandwidth=5)

IV3SLS (Three-Stage Least Squares)

SUR extended to allow endogenous variables. Each equation can have its own set of endogenous regressors and excluded instruments. The three stages are: (1) 2SLS per equation to get consistent residuals, (2) estimate cross-equation covariance from those residuals, (3) GLS using the estimated covariance.

from collections import OrderedDict
from linearmodels.system import IV3SLS

equations = OrderedDict()
equations["demand"] = {
    "dependent": data["quantity"],
    "exog": data[["const", "income"]],
    "endog": data[["price"]],
    "instruments": data[["cost", "weather"]],
}
equations["supply"] = {
    "dependent": data["quantity"],
    "exog": data[["const", "cost"]],
    "endog": data[["price"]],
    "instruments": data[["income", "weather"]],
}

mod = IV3SLS(equations)
res = mod.fit(cov_type="robust")
print(res.summary)

Mixed Equations (Some IV, Some Not)

Not every equation needs endogenous variables. Equations without "endog" and "instruments" keys are treated as exogenous:

equations = OrderedDict()
equations["demand"] = {
    "dependent": data["quantity"],
    "exog": data[["const", "income"]],
    "endog": data[["price"]],
    "instruments": data[["cost", "weather"]],
}
equations["wage"] = {
    "dependent": data["wage"],
    "exog": data[["const", "educ", "exper"]],
}

mod = IV3SLS(equations)
res = mod.fit(cov_type="robust")

IVSystemGMM

System GMM estimation for multiple IV equations. More efficient than 3SLS under heteroskedasticity because it uses an optimal weight matrix rather than the parametric GLS covariance.

from linearmodels.system import IVSystemGMM

mod = IVSystemGMM(equations)
res = mod.fit(weight_type="robust")
print(res.summary)

Weight Matrix Options

When to Use IVSystemGMM vs IV3SLS

• Homoskedastic errors: 3SLS is efficient; GMM gains nothing
• Heteroskedastic errors: System GMM with weight_type="robust" dominates 3SLS
• Small samples: 3SLS may be more stable (GMM weight matrix estimation is noisy)
• Report both as a robustness check when feasible

Cross-Equation Constraints (LinearConstraint)

Test or impose linear restrictions across equation parameters using the R * beta = q format. Common use: constraining a coefficient to be equal across equations.

Parameter Ordering

System models stack parameters as a single vector: all parameters from equation 1, then all from equation 2, etc. Use res.params to inspect the ordering:

res = mod.fit()
print("Parameter names and indices:")
for i, name in enumerate(res.params.index):
    print(f"  [{i}] {name}")

Imposing an Equality Constraint

from linearmodels.system import SUR, LinearConstraint
import numpy as np

# Suppose res.params.index is:
#   [0] earnings_const
#   [1] earnings_exper
#   [2] earnings_tenure
#   [3] benefits_const
#   [4] benefits_exper
#   [5] benefits_union

# Constraint: exper coefficient equal across equations
# earnings_exper - benefits_exper = 0
total_params = len(res.params)
r = np.zeros((1, total_params))
r[0, 1] = 1.0    # earnings_exper
r[0, 4] = -1.0   # benefits_exper
q = np.array([0.0])

constraint = LinearConstraint(r, q)
res_constrained = mod.fit(constraints=constraint, cov_type="robust")
print(res_constrained.summary)

Multiple Simultaneous Constraints

# Constrain both exper AND constant to be equal across equations
r = np.zeros((2, total_params))
r[0, 0] = 1.0    # earnings_const
r[0, 3] = -1.0   # benefits_const
r[1, 1] = 1.0    # earnings_exper
r[1, 4] = -1.0   # benefits_exper
q = np.array([0.0, 0.0])

constraint = LinearConstraint(r, q)
res_constrained = mod.fit(constraints=constraint, cov_type="robust")

Interpreting System Output

Per-Equation Results

res.summary reports coefficients, standard errors, t-stats, and p-values grouped by equation. Each equation has its own R-squared.

res = mod.fit(cov_type="robust")

# Full summary with all equations
print(res.summary)

# Access equation-level results
for eq_name in res.equations:
    eq = res.equations[eq_name]
    print(f"\n--- {eq_name} ---")
    print(f"  R-squared: {eq.rsquared:.4f}")
    print(f"  Params:\n{eq.params}")
    print(f"  Std Errors:\n{eq.std_errors}")

System-Level Statistics

# System R-squared (weighted average across equations)
print(f"System R-squared: {res.rsquared:.4f}")

# Cross-equation covariance matrix (sigma)
# Off-diagonal elements measure error correlation across equations
print("Sigma (cross-equation covariance):")
print(res.sigma)

# Total number of observations and parameters
print(f"Total observations: {res.nobs}")
print(f"Total parameters: {len(res.params)}")

Joint Hypothesis Tests

Use the Wald test via constraints to test joint hypotheses across equations:

# Test H0: exper coefficient = 0 in BOTH equations simultaneously
r = np.zeros((2, total_params))
r[0, 1] = 1.0   # earnings_exper = 0
r[1, 4] = 1.0   # benefits_exper = 0
q = np.array([0.0, 0.0])

# Compare constrained vs unconstrained fit
constraint = LinearConstraint(r, q)
res_restricted = mod.fit(constraints=constraint)

Performance Notes

• SUR, 3SLS, and System GMM use dense matrix operations on the stacked system.

Memory and computation scale with (sum of N_i * K_i)^2 where N_i is the observations and K_i is the parameters per equation.

• For large N with few equations, consider equation-by-equation estimation with

cluster-robust standard errors as a scalable alternative (sacrificing cross-equation efficiency).

• Cross-equation constraints add minimal overhead -- the constraint matrix is small

relative to the system.

• Iterative SUR (iter_limit > 1) adds one matrix factorization per iteration but

typically converges in fewer than 20 iterations.

References and Further Reading

• Sheppard, K. linearmodels System documentation. https://bashtage.github.io/linearmodels/system/index.html
• Zellner, A. (1962). "An Efficient Method of Estimating Seemingly Unrelated Regressions and Tests for Aggregation Bias." JASA, 57(298), 348-368.
• Zellner, A. and Theil, H. (1962). "Three-Stage Least Squares: Simultaneous Estimation of Simultaneous Equations." Econometrica, 30(1), 54-78.
• Greene, W.H. (2018). Econometric Analysis. 8th ed. Pearson. Ch. 10 (Systems of Equations).
• Wooldridge, J.M. (2010). Econometric Analysis of Cross Section and Panel Data. 2nd ed. MIT Press. Ch. 7 (Simultaneous Equations Models).