facmodcs 0.1.1
Cross-Section Factor Models - High-performance Python implementation
License: MIT
Requires Python: >=3.12
Keywords: factor-models, finance, quantitative-finance, risk-analysis

# facmodCS: Cross-Sectional Factor Models in Python

[![Python 3.9+](https://img.shields.io/badge/python-3.9+-blue.svg)](https://www.python.org/downloads/)
[![Code Coverage](https://img.shields.io/badge/coverage-86%25-brightgreen.svg)](htmlcov/index.html)
[![Tests](https://img.shields.io/badge/tests-223%20passing-brightgreen.svg)](tests/)

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

```bash
# 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

```python
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

```python
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

```python
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

```python
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

```python
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

```python
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

```python
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).

```python
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**:

![Factor Returns](examples/plots/factor_returns.html)
![SD Decomposition](examples/plots/sd_decomposition.html)
![Model Dashboard](examples/plots/model_dashboard.html)

---

## 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)

```python
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)

```python
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)

```python
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)

```python
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

```bash
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

```bash
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**:
```bash
# 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

```bash
# 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

```python
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

```python
# 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

```python
# 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:

1. Fork the repository
2. Create a feature branch
3. Write tests for new functionality
4. Ensure all tests pass: `pytest tests/ -v`
5. Submit a pull request

**Development Setup**:

```bash
# 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:

```bibtex
@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:
```bibtex
@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**.

The original R `facmodCS` package (version 1.0.3) was authored by:
- **Mido Shammaa** (Author, Creator) - midoshammaa@yahoo.com
- **Doug Martin** (Contributor, Author) - martinrd3d@gmail.com
- **Kirk Li** (Author, Contributor) - cocokecoli@gmail.com
- **Avinash Acharya** (Contributor) - Jon.Spinney@vestcor.org
- **Lingjie Yi** (Contributor) - lingjy.uw@gmail.com

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](https://github.com/yourusername/facmodCS/issues)
- **Documentation**: See `IMPLEMENTATION_STATUS.md` for detailed implementation notes
- **Examples**: See `examples/` directory

---

## Roadmap

Future enhancements:

- [ ] Rolling window estimation
- [ ] Out-of-sample prediction
- [ ] Additional risk measures (CVaR, drawdown)
- [ ] Style analysis tools
- [ ] Performance dashboards
- [ ] API for real-time data feeds

---

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