facmodCS: Cross-Sectional Factor Models in Python

A comprehensive Python implementation of cross-sectional fundamental factor models, providing tools for factor model estimation, risk decomposition, and portfolio analysis.

This package is a Python translation of the R facmodCS package, offering full feature parity with enhanced interactive visualizations using Plotly.

Features

Core Capabilities

Multiple Estimation Methods: OLS, WLS, Robust (MM-estimation), Weighted Robust
Risk Decomposition: Standard Deviation, Value-at-Risk (VaR), Expected Shortfall (ES)
Factor Analysis: Exposure analysis, factor returns, correlation matrices
Portfolio Analytics: Performance attribution, risk decomposition by factor
Interactive Visualizations: 10+ Plotly-based charts for model diagnostics
Comprehensive Reporting: Automated reports with statistics and insights

Key Features

✅ Cross-sectional regression at each time period ✅ Categorical variables (sectors, industries) ✅ Exposure standardization (cross-sectional, time-series) ✅ Variance estimation (EWMA, Robust EWMA, GARCH) ✅ Euler risk decomposition with marginal, component, and percentage contributions ✅ Production-ready with 86% test coverage and 223 passing tests

Installation

Prerequisites

Python 3.9 or higher
UV package manager (recommended) or pip

Install from Source

# Clone the repository
cd facmodcs_py

# Install with UV (recommended)
uv pip install -e ".[test]"

# Or install with pip
pip install -e ".[test]"

Dependencies

Core dependencies:

polars - Fast DataFrame operations
pandas - Data manipulation
numpy - Numerical computing
scipy - Scientific computing
statsmodels - Statistical models
plotly - Interactive visualizations
numba - JIT compilation for performance

Quick Start

Basic Usage

import polars as pl
from facmodcs import fit_ffm, fm_sd_decomp, port_sd_decomp
from facmodcs import plot_factor_returns, create_model_dashboard

# Load your panel data (balanced panel: Date, Asset, Return, Exposures)
data = pl.read_csv("panel_data.csv")

# Fit a cross-sectional factor model
model = fit_ffm(
    data=data,
    asset_var="Asset",
    ret_var="Return",
    date_var="Date",
    exposure_vars=["BP", "PM12M1M", "Beta60M"],  # Your factor exposures
    fit_method="LS",                              # LS, WLS, Rob, or W-Rob
    standardize="CrossSection",                   # Standardize exposures
    lag_exposures_flag=True,                      # Lag exposures by 1 period
)

# View model summary
print(model.summary())

# Access fitted components
beta_matrix = model.beta              # (N × K) factor exposures
factor_returns = model.factor_returns # (T × K) factor returns
residuals = model.residuals           # (N × T) residual matrix
r_squared = model.r_squared           # (T,) R-squared by period

# Risk decomposition (asset-level)
sd_decomp = fm_sd_decomp(model)
first_asset = model.asset_names[0]
print(sd_decomp[first_asset].summary())

# Portfolio risk decomposition
import numpy as np
weights = np.ones(model.n_assets) / model.n_assets  # Equal weights
port_sd = port_sd_decomp(model, weights)
print(f"Portfolio SD: {port_sd.total_risk:.4f}")
print(f"Factor contributions: {port_sd.percentage}")

Interactive Visualization

from facmodcs import (
    plot_factor_returns,
    plot_cumulative_factor_returns,
    plot_sd_decomposition,
    create_model_dashboard,
)

# Plot factor returns over time
fig = plot_factor_returns(model)
fig.show()

# Plot cumulative factor returns
fig = plot_cumulative_factor_returns(model)
fig.show()

# Plot SD decomposition for top 10 risky assets
fig = plot_sd_decomposition(model, top_n=10)
fig.show()

# Create comprehensive dashboard
fig = create_model_dashboard(model, weights)
fig.show()

# Save plots to HTML
fig.write_html("model_dashboard.html")

Main Entry Points

1. Model Fitting

`fit_ffm()` - Main fitting function

from facmodcs import fit_ffm

model = fit_ffm(
    data: pl.DataFrame,           # Input panel data
    asset_var: str,               # Asset identifier column
    ret_var: str,                 # Return column
    date_var: str,                # Date column
    exposure_vars: list[str],     # List of exposure variables
    fit_method: str = "LS",       # LS, WLS, Rob, W-Rob
    standardize: str = "CrossSection",  # CrossSection, TimeSeries, or none
    lag_exposures_flag: bool = True,    # Lag exposures by 1 period
    var_model: str = "EWMA",      # For WLS: StdDev, EWMA, RobustEWMA, GARCH
    **kwargs
) -> FactorModel

Returns: FactorModel object with:

beta: Factor exposures (N × K)
factor_returns: Factor returns (T × K)
residuals: Residuals (N × T)
r_squared: R² values (T,)
factor_cov: Factor covariance matrix (K × K)
resid_var: Residual variances (N,)

2. Risk Decomposition

Asset-Level Decomposition

from facmodcs import fm_sd_decomp, fm_var_decomp, fm_es_decomp

# Standard Deviation decomposition
sd_decomp = fm_sd_decomp(model)
# Returns: dict[str, DecompositionResult]
# Each asset has: total_risk, component, marginal, percentage

# Value-at-Risk decomposition
var_decomp = fm_var_decomp(model, p=0.05, type="np")
# type: "np" (non-parametric) or "normal" (parametric)

# Expected Shortfall decomposition
es_decomp = fm_es_decomp(model, p=0.05, type="np")

Portfolio-Level Decomposition

from facmodcs import port_sd_decomp, port_var_decomp, port_es_decomp

weights = np.ones(model.n_assets) / model.n_assets  # Portfolio weights

# Portfolio SD decomposition
port_sd = port_sd_decomp(model, weights)
print(f"Total Portfolio SD: {port_sd.total_risk:.4f}")
print(f"Factor contributions (%): {port_sd.percentage}")

# Portfolio VaR decomposition
port_var = port_var_decomp(model, weights, p=0.05, type="np")

# Portfolio ES decomposition
port_es = port_es_decomp(model, weights, p=0.05, type="np")

3. Post-Fitting Analysis

from facmodcs import fm_cov, fm_rsq, fm_tstats, vif, summary_statistics

# Covariance matrix
cov_matrix = fm_cov(model)  # (N × N) asset return covariance

# R-squared statistics
rsq_stats = fm_rsq(model, adj_rsq=True)
print(f"Mean R²: {rsq_stats['mean_r_squared']:.4f}")
print(f"Mean Adj. R²: {rsq_stats['mean_adj_r_squared']:.4f}")

# T-statistics for factor significance
tstats = fm_tstats(model, z_alpha=1.96)  # 95% confidence
print(f"Significant factors: {tstats['total_significant']}")

# Variance Inflation Factors (multicollinearity check)
vif_results = vif(model)
for factor, mean_vif in vif_results['mean_vif'].items():
    print(f"{factor}: VIF = {mean_vif:.2f}")

# Comprehensive summary
summary_df = summary_statistics(model)
print(summary_df)

4. Reporting Functions

from facmodcs import (
    rep_exposures,
    rep_return,
    rep_risk,
    model_summary_report,
    factor_correlation_report,
    performance_attribution_summary,
)

# Exposure summary (mean, std, min, max, VIF)
exposure_report = rep_exposures(model)
print(exposure_report)

# Return attribution (factor contributions to returns)
asset = model.asset_names[0]
return_attr = rep_return(model, asset=asset)
print(return_attr.head(10))

# Risk decomposition report (top N risky assets)
risk_reports = rep_risk(model, risk_measures=["SD", "VaR", "ES"], p=0.05, top_n=10)
print(risk_reports['SD'])  # Top 10 by SD
print(risk_reports['VaR']) # Top 10 by VaR

# Model summary (comprehensive text report)
summary_text = model_summary_report(model)
print(summary_text)

# Factor correlation matrix
corr_matrix = factor_correlation_report(model)
print(corr_matrix)

# Performance attribution (portfolio level)
perf_attr = performance_attribution_summary(model, weights)
print(perf_attr)

5. Plotting Functions

All plotting functions return Plotly Figure objects with interactive features (zoom, pan, hover).

from facmodcs import (
    plot_factor_returns,            # Time series of factor returns
    plot_cumulative_factor_returns, # Cumulative returns
    plot_r_squared,                 # R² over time
    plot_factor_distribution,       # Histogram with normal curve
    plot_sd_decomposition,          # Stacked bar chart (top N assets)
    plot_portfolio_risk_decomposition,  # Waterfall chart
    plot_risk_comparison,           # Grouped bar chart (SD/VaR/ES)
    plot_exposure_heatmap,          # Factor exposure heatmap
    plot_residual_diagnostics,      # 4-panel diagnostic plots
    create_model_dashboard,         # Comprehensive 6-panel dashboard
)

# Factor returns time series
fig = plot_factor_returns(model, factors=["BP", "PM12M1M"])
fig.show()

# Cumulative returns (compound returns)
fig = plot_cumulative_factor_returns(model)
fig.write_html("cumulative_returns.html")

# R-squared over time
fig = plot_r_squared(model)
fig.show()

# Factor distribution with normal overlay
fig = plot_factor_distribution(model, factor="BP")
fig.show()

# SD decomposition (top 10 assets)
fig = plot_sd_decomposition(model, top_n=10)
fig.show()

# Portfolio risk decomposition (waterfall chart)
fig = plot_portfolio_risk_decomposition(model, weights, risk_measure="SD")
fig.show()

# Risk comparison (SD vs VaR vs ES)
fig = plot_risk_comparison(model, weights, p=0.05)
fig.show()

# Exposure heatmap
fig = plot_exposure_heatmap(model, n_assets=20)
fig.show()

# Residual diagnostics (time series, histogram, Q-Q, ACF)
fig = plot_residual_diagnostics(model, asset=model.asset_names[0])
fig.show()

# Comprehensive dashboard (6 panels)
fig = create_model_dashboard(model, weights)
fig.write_html("dashboard.html")

Example Visualizations:

Data Requirements

Input data must be a balanced panel with the following structure:

Date	Asset	Return	BP	PM12M1M	Beta60M	Sector
2020-01-31	AAPL	0.0523	1.23	0.45	1.05	Tech
2020-01-31	MSFT	0.0412	1.45	0.38	0.98	Tech
...	...	...	...	...	...	...
2020-02-29	AAPL	0.0234	1.25	0.48	1.06	Tech

Requirements:

Balanced panel: All assets present in all time periods
Date column: Convertible to date type
Asset column: Unique identifiers
Return column: Numeric returns
Exposure columns: Numeric (continuous) or categorical (sectors)

Data Formats Supported:

Polars DataFrame (preferred for performance)
Pandas DataFrame (automatically converted)
CSV files (via pl.read_csv())

Estimation Methods

1. Ordinary Least Squares (LS)

model = fit_ffm(data, ..., fit_method="LS")

Standard OLS cross-sectional regression at each time period.

Pros: Fast, unbiased estimates Cons: Sensitive to outliers

2. Weighted Least Squares (WLS)

model = fit_ffm(data, ..., fit_method="WLS", var_model="EWMA")

Two-pass estimation with inverse variance weights.

Variance Models:

StdDev: Sample variance
EWMA: Exponentially weighted moving average (default λ=0.9)
RobustEWMA: Robust EWMA with outlier rejection
GARCH: GARCH(1,1) conditional variance

Pros: Accounts for heteroskedasticity Cons: Requires first-pass estimates

3. Robust Regression (Rob)

model = fit_ffm(data, ..., fit_method="Rob")

MM-estimation with Huber loss function.

Pros: Robust to outliers Cons: Slower convergence

4. Weighted Robust (W-Rob)

model = fit_ffm(data, ..., fit_method="W-Rob", var_model="RobustEWMA")

Combines WLS weighting with robust estimation.

Pros: Best of both worlds Cons: Most computationally intensive

Examples

See the examples/ directory for complete demonstrations:

examples/phase5_demo.py - Comprehensive demo of reporting & visualization
examples/end_to_end_analysis.ipynb - Jupyter notebook with full workflow
examples/plots/ - 10 interactive HTML plots

Run the Demo

python examples/phase5_demo.py

This generates:

6 console reports (model summary, exposures, correlations, etc.)
10 interactive HTML plots in examples/plots/

Testing

Run All Tests

pytest tests/ -v --cov=facmodcs --cov-report=html

Current Status:

✅ 223 tests passing
✅ 86% code coverage
✅ 0 failures

Test Categories

Unit tests: Individual function testing
Integration tests: End-to-end workflows
Regression tests: All method combinations
Edge cases: Missing data, outliers, numerical stability
Performance tests: Benchmarks for speed and memory
R comparison (Real Data): Validation against R package using stocks_factors.csv ⭐ NEW
R comparison (Synthetic Data): Validation against R package using generated data

Testing with Real Data

The test suite includes comprehensive R-Python comparison tests using a common real-world dataset (stocks_factors.csv). This ensures both implementations produce identical results on the same data:

Dataset Details:

294 US stocks from 2006-2010 (monthly)
Real-world patterns: Same data used in R package examples
No seed differences: Eliminates random number generation issues
Comprehensive: Validates fitting, decomposition, and portfolio analysis

Why Real Data?

✅ Both implementations handle real-world patterns identically
✅ No synthetic data generation differences (R vs Python seeds differ)
✅ Reproducible results across languages
✅ Validation against published R results

Run R Comparison Tests on Real Data:

# Requires: R installed, facmodCS R package, rpy2
pytest tests/test_r_comparison_real_data.py --run-r-comparison -v

# Install R package (if needed)
R -e "install.packages('facmodCS', repos='https://cloud.r-project.org')"

See: tests/README_TESTING.md for comprehensive testing guide

Run Specific Tests

# Real data R comparison (RECOMMENDED for validation)
pytest tests/test_r_comparison_real_data.py --run-r-comparison -v

# Synthetic data R comparison
pytest tests/test_r_comparison.py --run-r-comparison -v

# Seed behavior and reproducibility
pytest tests/test_seed_equivalence.py -v

# Plotting tests
pytest tests/test_plotting.py -v

# Reporting tests
pytest tests/test_reporting.py -v

# Performance benchmarks
pytest tests/test_performance.py -v -s

Random Seed Handling

IMPORTANT: R's set.seed() and Python's np.random.seed() produce different values even with the same seed number due to different RNG algorithms.

Best Practice:

✅ Use real data (stocks_factors.csv) for R-Python comparison
✅ Use synthetic data for Python-only edge case tests
❌ Don't compare R and Python synthetic data directly

See: docs/SEED_HANDLING.md for detailed explanation

Performance

Benchmarks (on medium dataset: 100 assets, 60 periods):

Operation	Time	Target
LS Fitting	1.2s	< 5s
SD Decomposition	0.3s	< 1s
Portfolio SD	0.02s	< 0.1s
Full Workflow	3.5s	< 10s

Memory Usage: < 500 MB for large models (500 assets, 120 periods)

Optimizations:

Polars for fast data operations
Numba JIT compilation for hot loops
Efficient matrix operations with NumPy
Lazy evaluation where possible

Project Structure

facmodcs_py/
├── facmodcs/              # Main package
│   ├── __init__.py        # Public API
│   ├── core.py            # Fitting functions (~380 lines)
│   ├── covariance.py      # Post-fitting analysis (~100 lines)
│   ├── decomposition.py   # Risk decomposition (~250 lines)
│   ├── plotting.py        # Visualizations (~180 lines)
│   ├── reporting.py       # Reports (~170 lines)
│   ├── models.py          # Data classes (~115 lines)
│   └── utils.py           # Utilities (~75 lines)
├── tests/                 # Test suite (~2,000 lines)
│   ├── conftest.py        # Test configuration
│   ├── test_core.py       # Core fitting tests
│   ├── test_decomposition.py
│   ├── test_plotting.py
│   ├── test_reporting.py
│   ├── test_regression_suite.py
│   ├── test_edge_cases.py
│   ├── test_performance.py
│   ├── test_r_comparison.py
│   └── fixtures/
│       └── synthetic_data_generator.py
├── examples/              # Examples and demos
│   ├── phase5_demo.py
│   ├── end_to_end_analysis.ipynb
│   └── plots/             # Generated plots
├── pyproject.toml         # Package configuration
├── requirements.txt       # Dependencies
├── README.md             # This file
└── IMPLEMENTATION_STATUS.md  # Detailed implementation notes

Advanced Usage

Custom Preprocessing

from facmodcs import lag_exposures, standardize_exposures

# Lag exposures manually
lagged_data = lag_exposures(data, asset_var="Asset", date_var="Date",
                             exposures=["BP", "PM12M1M"])

# Standardize with custom method
standardized_data = standardize_exposures(
    lagged_data,
    date_var="Date",
    exposures_num=["BP", "PM12M1M"],
    method="TimeSeries",  # or CrossSection
    lambda_ewma=0.9,      # EWMA decay
)

Model Customization

# Fit with specific options
model = fit_ffm(
    data=data,
    asset_var="Asset",
    ret_var="Return",
    date_var="Date",
    exposure_vars=["Sector", "BP", "PM12M1M"],  # Mix categorical + numeric
    fit_method="W-Rob",
    var_model="RobustEWMA",
    lambda_ewma=0.94,           # Custom EWMA decay
    robust_scale_threshold=2.5, # Outlier threshold
    max_iter=50,                # Robust regression iterations
    add_intercept=True,         # Include intercept
    standardize="CrossSection",
    lag_exposures_flag=True,
)

Batch Processing

# Fit multiple models
models = {}
for method in ["LS", "WLS", "Rob"]:
    models[method] = fit_ffm(data, ..., fit_method=method)

# Compare model performance
for method, model in models.items():
    rsq_stats = fm_rsq(model, adj_rsq=True)
    print(f"{method}: Mean R² = {rsq_stats['mean_r_squared']:.4f}")

Contributing

Contributions are welcome! Please follow these guidelines:

Fork the repository
Create a feature branch
Write tests for new functionality
Ensure all tests pass: pytest tests/ -v
Submit a pull request

Development Setup:

# Install in development mode with test dependencies
uv pip install -e ".[test]"

# Run tests with coverage
pytest tests/ -v --cov=facmodcs --cov-report=html

# View coverage report
open htmlcov/index.html

Citation

If you use this package in academic research please cite the original authors. The python implementation is below:

@software{facmodcs_python,
  title = {facmodCS: Cross-Sectional Factor Models in Python},
  author = {Hayes, Drew},
  year = {2025},
  url = {https://github.com/Druhayes/facmodcs}
}

Include the original authors of the R package as well

Original R package:

@manual{facmodCS_R,
  title = {facmodCS: Cross-Sectional Factor Models},
  author = {PCRA Team},
  year = {2023},
  note = {R package}
}

Acknowledgments

Original Authors

This codebase was written by Claude Sonnet 4.5 with Claude Code.

This package is based on the original R facmodCS package (version 1.0.3), developed by the PCRA Team.

Special thanks to:

Chicago Research on Security Prices, LLC for the cross-section of about 300 CRSP stocks data
S&P GLOBAL MARKET INTELLIGENCE for contributing 14 factor scores (alpha factors and factor exposures) fundamental data

Technology Stack

Built with Polars, NumPy, SciPy, Statsmodels, and Plotly
Testing framework using pytest and rpy2

Support

For questions, issues, or feature requests:

Issues: GitHub Issues
Documentation: See IMPLEMENTATION_STATUS.md for detailed implementation notes
Examples: See examples/ directory

Roadmap

Future enhancements:

Version: 0.1.0 Status: Production Ready Python: 3.9+ Last Updated: January 2025

Name		Name	Last commit message	Last commit date
Latest commit History 12 Commits
.github/workflows		.github/workflows
docs		docs
examples		examples
facmodcs		facmodcs
tests		tests
.gitignore		.gitignore
README.md		README.md
pyproject.toml		pyproject.toml
requirements.txt		requirements.txt
uv.lock		uv.lock

Folders and files

Latest commit

History

Repository files navigation

facmodCS: Cross-Sectional Factor Models in Python

Features

Core Capabilities

Key Features

Installation

Prerequisites

Install from Source

Dependencies

Quick Start

Basic Usage

Interactive Visualization

Main Entry Points

1. Model Fitting

fit_ffm() - Main fitting function

2. Risk Decomposition

Asset-Level Decomposition

Portfolio-Level Decomposition

3. Post-Fitting Analysis

4. Reporting Functions

5. Plotting Functions

Data Requirements

Estimation Methods

1. Ordinary Least Squares (LS)

2. Weighted Least Squares (WLS)

3. Robust Regression (Rob)

4. Weighted Robust (W-Rob)

Examples

Run the Demo

Testing

Run All Tests

Test Categories

Testing with Real Data

Run Specific Tests

Random Seed Handling

Performance

Project Structure

Advanced Usage

Custom Preprocessing

Model Customization

Batch Processing

Contributing

Citation

Acknowledgments

Original Authors

Technology Stack

Support

Roadmap

About

Resources

Uh oh!

Stars

Watchers

Forks

Releases

Packages 0

Uh oh!

Contributors

Uh oh!

Languages

`fit_ffm()` - Main fitting function

Packages