xjtulyc/factor-models

Factor Models

A comprehensive skill for estimating and interpreting asset pricing factor models in academic and professional research. Covers downloading Fama-French factor data from the Ken French Data Library, computing excess returns, OLS time-series regressions, the Carhart (1997) 4-factor momentum model, the q-factor model, rolling factor loadings, and the Gibbons-Ross-Shanken (1989) test for pricing errors.

---

Core Functions

1. Downloading Fama-French Factor Data

import os
import numpy as np
import pandas as pd
import pandas_datareader.data as web
import statsmodels.api as sm
from scipy import linalg
import matplotlib.pyplot as plt
import matplotlib.dates as mdates


# ── Fama-French dataset names ─────────────────────────────────────────────────
_FF_DATASETS = {
    "3factor": "F-F_Research_Data_Factors",
    "5factor": "F-F_Research_Data_5_Factors_2x3",
    "momentum": "F-F_Momentum_Factor",
    "qfactor": "q5_factors_daily_2022",   # from AQR or q-factor website
}


def download_ff_factors(
    model: str = "5factor",
    start: str = "2000-01-01",
    end: str = "2024-12-31",
    freq: str = "monthly",
) -> pd.DataFrame:
    """
    Download Fama-French factor returns from the Ken French Data Library.

    Parameters
    ----------
    model : str
        One of '3factor', '5factor', 'momentum'. '4factor' downloads 3factor
        + momentum and merges them automatically.
    start : str
        Start date in 'YYYY-MM-DD' format.
    end : str
        End date in 'YYYY-MM-DD' format.
    freq : str
        'monthly' or 'daily'. Ken French provides both.

    Returns
    -------
    pd.DataFrame
        DataFrame with factor columns in percent terms (divide by 100 for decimals).
        Always includes 'RF' (risk-free rate) and 'Mkt-RF' (market excess return).
    """
    if model == "4factor":
        ff3 = download_ff_factors("3factor", start, end, freq)
        mom = download_ff_factors("momentum", start, end, freq)
        combined = ff3.join(mom, how="inner")
        return combined

    dataset_key = _FF_DATASETS.get(model)
    if dataset_key is None:
        raise ValueError(f"Unsupported model: {model!r}. Choose from {list(_FF_DATASETS)}")

    # pandas_datareader fetches directly from Ken French's website
    factors_raw = web.DataReader(
        dataset_key,
        "famafrench",
        start=start,
        end=end,
    )
    # DataReader returns a dict; index 0 is monthly, index 1 is annual
    if freq == "monthly":
        df = factors_raw[0]
    elif freq == "daily":
        df = factors_raw[0]  # daily datasets only have one table
    else:
        raise ValueError(f"freq must be 'monthly' or 'daily', got {freq!r}")

    df.index = pd.to_datetime(df.index.astype(str), format="%Y-%m")
    df.index.name = "date"
    return df / 100.0  # convert from percent to decimal


def compute_excess_returns(
    prices: pd.DataFrame,
    rf: pd.Series,
    freq: str = "monthly",
) -> pd.DataFrame:
    """
    Compute simple excess returns for a portfolio or individual assets.

    Parameters
    ----------
    prices : pd.DataFrame
        Adjusted closing prices. Index = DatetimeIndex, columns = ticker symbols.
    rf : pd.Series
        Risk-free rate series aligned to the same period frequency (decimal form).
    freq : str
        'monthly' or 'daily'. Used to determine the resample frequency.

    Returns
    -------
    pd.DataFrame
        Excess return columns: R_i - RF for each asset.
    """
    if freq == "monthly":
        # Resample to month-end and compute simple returns
        monthly_prices = prices.resample("ME").last()
        returns = monthly_prices.pct_change().dropna()
    elif freq == "daily":
        returns = prices.pct_change().dropna()
    else:
        raise ValueError(f"freq must be 'monthly' or 'daily', got {freq!r}")

    # Align rf index to returns index
    rf_aligned = rf.reindex(returns.index).ffill()
    excess = returns.subtract(rf_aligned, axis=0)
    return excess.dropna()

2. OLS Factor Regression

def run_factor_regression(
    excess_ret: pd.Series,
    factors: pd.DataFrame,
    add_intercept: bool = True,
) -> dict:
    """
    Run an OLS time-series factor regression for a single asset or portfolio.

    Parameters
    ----------
    excess_ret : pd.Series
        Excess returns of the test asset (decimal).
    factors : pd.DataFrame
        Factor return data aligned to the same dates (decimal). Should NOT
        include the risk-free rate column.
    add_intercept : bool
        Whether to add a constant (alpha) to the regression.

    Returns
    -------
    dict with keys:
        alpha        : float  — annualised Jensen's alpha
        alpha_tstat  : float  — t-statistic for alpha
        alpha_pval   : float  — p-value for alpha
        betas        : dict   — factor loadings keyed by factor name
        r_squared    : float  — R² of the regression
        adj_r2       : float  — adjusted R²
        residuals    : pd.Series
        model        : RegressionResultsWrapper (full statsmodels result)
    """
    # Align on common dates
    aligned = pd.concat([excess_ret, factors], axis=1, join="inner").dropna()
    y = aligned.iloc[:, 0]
    X = aligned.iloc[:, 1:]
    if add_intercept:
        X = sm.add_constant(X)

    model = sm.OLS(y, X).fit(cov_type="HC3")

    alpha_monthly = model.params.get("const", 0.0)
    alpha_annual = alpha_monthly * 12  # annualise for monthly data

    factor_names = [c for c in model.params.index if c != "const"]
    betas = {name: model.params[name] for name in factor_names}

    return {
        "alpha": alpha_annual,
        "alpha_monthly": alpha_monthly,
        "alpha_tstat": model.tvalues.get("const", np.nan),
        "alpha_pval": model.pvalues.get("const", np.nan),
        "betas": betas,
        "r_squared": model.rsquared,
        "adj_r2": model.rsquared_adj,
        "residuals": model.resid,
        "model": model,
    }


def rolling_factor_loadings(
    excess_ret: pd.Series,
    factors: pd.DataFrame,
    window: int = 36,
) -> pd.DataFrame:
    """
    Estimate rolling OLS factor loadings using a fixed rolling window.

    Parameters
    ----------
    excess_ret : pd.Series
        Excess returns (decimal, monthly).
    factors : pd.DataFrame
        Factor returns aligned to same index.
    window : int
        Rolling window size in periods (default 36 months = 3 years).

    Returns
    -------
    pd.DataFrame
        Rolling coefficients. Columns: 'alpha', plus one per factor.
        First (window - 1) rows are NaN.
    """
    aligned = pd.concat([excess_ret, factors], axis=1, join="inner").dropna()
    y = aligned.iloc[:, 0]
    factor_cols = aligned.iloc[:, 1:]

    result_records = []
    for end in range(window, len(aligned) + 1):
        y_w = y.iloc[end - window : end]
        X_w = sm.add_constant(factor_cols.iloc[end - window : end])
        try:
            res = sm.OLS(y_w, X_w).fit()
            record = {"date": y.index[end - 1]}
            record["alpha"] = res.params.get("const", np.nan) * 12
            for col in factor_cols.columns:
                record[col] = res.params.get(col, np.nan)
        except Exception:
            record = {"date": y.index[end - 1], "alpha": np.nan}
            for col in factor_cols.columns:
                record[col] = np.nan
        result_records.append(record)

    df = pd.DataFrame(result_records).set_index("date")
    return df

3. GRS Test for Pricing Errors

def grs_test(
    excess_rets: pd.DataFrame,
    factors: pd.DataFrame,
) -> dict:
    """
    Gibbons, Ross, and Shanken (1989) test for joint pricing errors (alphas = 0).

    The GRS statistic tests the null hypothesis H0: alpha_i = 0 for all i
    simultaneously under the assumption of multivariate normality.

    Parameters
    ----------
    excess_rets : pd.DataFrame
        Excess returns of N test assets (T × N).
    factors : pd.DataFrame
        Factor returns (T × K). Must be aligned to same dates as excess_rets.

    Returns
    -------
    dict with keys:
        grs_stat   : float  — GRS F-statistic
        p_value    : float  — p-value under F(N, T-N-K) distribution
        alphas     : np.ndarray  — estimated alphas (N,)
        t_stats    : np.ndarray  — individual alpha t-statistics
    """
    from scipy.stats import f as f_dist

    aligned = pd.concat([excess_rets, factors], axis=1, join="inner").dropna()
    N = excess_rets.shape[1]
    K = factors.shape[1]
    T = len(aligned)

    alpha_hat = np.zeros(N)
    residuals = np.zeros((T, N))
    t_stats = np.zeros(N)

    factor_X = sm.add_constant(aligned[factors.columns])

    for j, col in enumerate(excess_rets.columns):
        y = aligned[col]
        res = sm.OLS(y, factor_X).fit()
        alpha_hat[j] = res.params["const"]
        residuals[:, j] = res.resid
        t_stats[j] = res.tvalues["const"]

    # Residual covariance matrix
    Sigma_hat = (residuals.T @ residuals) / (T - K - 1)

    # Mean factor returns and factor covariance
    mu_f = aligned[factors.columns].mean().values  # (K,)
    Omega_f = aligned[factors.columns].cov().values  # (K, K)

    # Sharpe ratio adjustment term
    sh_squared = float(mu_f @ linalg.inv(Omega_f) @ mu_f)

    # GRS F-statistic
    grs_stat = (
        (T / N)
        * ((T - N - K) / (T - K - 1))
        * (alpha_hat @ linalg.inv(Sigma_hat) @ alpha_hat)
        / (1 + sh_squared)
    )

    df1, df2 = N, T - N - K
    p_value = 1.0 - f_dist.cdf(grs_stat, df1, df2)

    return {
        "grs_stat": grs_stat,
        "p_value": p_value,
        "alphas": alpha_hat,
        "t_stats": t_stats,
        "df1": df1,
        "df2": df2,
    }


def factor_exposure_attribution(
    excess_ret: pd.Series,
    factors: pd.DataFrame,
) -> pd.DataFrame:
    """
    Decompose a portfolio's expected return into factor contributions.

    Returns
    -------
    pd.DataFrame
        Columns: 'factor_mean_return', 'loading', 'contribution', 'pct_explained'.
    """
    result = run_factor_regression(excess_ret, factors)
    factor_means = factors.mean() * 12  # annualise

    rows = []
    for factor_name, beta in result["betas"].items():
        mean_ret = factor_means.get(factor_name, np.nan)
        contribution = beta * mean_ret
        rows.append({
            "factor": factor_name,
            "factor_mean_return": mean_ret,
            "loading": beta,
            "contribution": contribution,
        })

    df = pd.DataFrame(rows).set_index("factor")
    total_explained = df["contribution"].sum()
    df["pct_explained"] = df["contribution"] / total_explained * 100
    return df

---

Example 1: Portfolio Alpha vs Fama-French 5-Factor Model

Estimate the risk-adjusted alpha (Jensen's alpha) of a custom equal-weight portfolio against the Fama-French 5-factor model.

import yfinance as yf
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

# ── Download factor data ──────────────────────────────────────────────────────
START = "2015-01-01"
END   = "2024-12-31"

ff5 = download_ff_factors(model="5factor", start=START, end=END, freq="monthly")
print("FF5 factors loaded:", ff5.columns.tolist())
print(ff5.tail(3))

# ── Download stock price data for a hypothetical portfolio ────────────────────
TICKERS = ["AAPL", "MSFT", "GOOGL", "AMZN", "META"]
raw_prices = yf.download(TICKERS, start=START, end=END, auto_adjust=True, progress=False)["Close"]

# ── Compute portfolio excess returns ─────────────────────────────────────────
rf = ff5["RF"]
# Equal-weight portfolio returns
asset_excess = compute_excess_returns(raw_prices, rf, freq="monthly")
portfolio_excess = asset_excess.mean(axis=1)
portfolio_excess.name = "EW_Portfolio"

# ── Extract Fama-French 5 factors (exclude RF) ────────────────────────────────
factor_cols = [c for c in ff5.columns if c != "RF"]
factors_df = ff5[factor_cols]

# ── Run 5-factor regression ───────────────────────────────────────────────────
result = run_factor_regression(portfolio_excess, factors_df)

print("\n=== Fama-French 5-Factor Regression ===")
print(f"Alpha (annualised):   {result['alpha']:.4f}  ({result['alpha']*100:.2f}%)")
print(f"Alpha t-stat:         {result['alpha_tstat']:.3f}")
print(f"Alpha p-value:        {result['alpha_pval']:.4f}")
print(f"R-squared:            {result['r_squared']:.4f}")
print("\nFactor loadings:")
for factor, beta in result["betas"].items():
    print(f"  {factor:12s}: {beta:.4f}")

# ── Factor attribution ────────────────────────────────────────────────────────
attr = factor_exposure_attribution(portfolio_excess, factors_df)
print("\nReturn Attribution:")
print(attr.round(4))

# ── Plot rolling 3-year alpha ─────────────────────────────────────────────────
rolling = rolling_factor_loadings(portfolio_excess, factors_df, window=36)

fig, axes = plt.subplots(2, 1, figsize=(12, 8), sharex=True)

axes[0].plot(rolling.index, rolling["alpha"] * 100, label="Rolling Alpha (annualised %)")
axes[0].axhline(0, color="black", linewidth=0.8, linestyle="--")
axes[0].set_ylabel("Alpha (%)")
axes[0].set_title("Rolling 36-Month Alpha — EW Tech Portfolio vs FF5")
axes[0].legend()
axes[0].grid(True)

axes[1].plot(rolling.index, rolling["Mkt-RF"], label="Market Beta", color="C1")
axes[1].set_ylabel("Beta")
axes[1].set_xlabel("Date")
axes[1].set_title("Rolling Market Beta")
axes[1].legend()
axes[1].grid(True)

plt.tight_layout()
plt.savefig("ff5_factor_regression.png", dpi=150, bbox_inches="tight")
plt.show()

---

Example 2: Compare Factor Loadings of Growth vs Value ETFs

Use the Carhart 4-factor model to compare factor exposures of a growth ETF (IWF) and a value ETF (IWD), and run the GRS test to check if both alphas are jointly zero.

import yfinance as yf
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

START = "2010-01-01"
END   = "2024-12-31"

# ── Download factors ──────────────────────────────────────────────────────────
ff4 = download_ff_factors(model="4factor", start=START, end=END, freq="monthly")
factor_cols_4 = [c for c in ff4.columns if c != "RF"]
factors_4 = ff4[factor_cols_4]

# ── Download ETF prices ───────────────────────────────────────────────────────
etfs = ["IWF", "IWD"]   # iShares Russell 1000 Growth / Value
prices = yf.download(etfs, start=START, end=END, auto_adjust=True, progress=False)["Close"]

# ── Compute excess returns ────────────────────────────────────────────────────
rf = ff4["RF"]
excess_etfs = compute_excess_returns(prices, rf, freq="monthly")

# ── Individual regressions ────────────────────────────────────────────────────
results_summary = {}
for etf in etfs:
    r = run_factor_regression(excess_etfs[etf], factors_4)
    results_summary[etf] = r
    print(f"\n=== {etf} — Carhart 4-Factor Model ===")
    print(f"  Alpha (ann.):  {r['alpha']*100:.2f}%  t={r['alpha_tstat']:.2f}  p={r['alpha_pval']:.4f}")
    print(f"  R²:            {r['r_squared']:.4f}")
    for factor, beta in r["betas"].items():
        print(f"  {factor:12s}: {beta:.4f}")

# ── GRS test: are both alphas jointly zero? ───────────────────────────────────
grs = grs_test(excess_etfs[etfs], factors_4)
print(f"\n=== GRS Test ===")
print(f"  GRS F-statistic: {grs['grs_stat']:.4f}")
print(f"  p-value:         {grs['p_value']:.4f}")
print(f"  Individual alpha t-stats: IWF={grs['t_stats'][0]:.3f}, IWD={grs['t_stats'][1]:.3f}")
if grs["p_value"] < 0.05:
    print("  -> Reject H0: at least one alpha is statistically significant.")
else:
    print("  -> Fail to reject H0: alphas are jointly indistinguishable from zero.")

# ── Rolling loadings comparison ───────────────────────────────────────────────
roll_iwf = rolling_factor_loadings(excess_etfs["IWF"], factors_4, window=36)
roll_iwd = rolling_factor_loadings(excess_etfs["IWD"], factors_4, window=36)

fig, axes = plt.subplots(3, 1, figsize=(13, 10), sharex=True)

factor_plot_pairs = [
    ("SMB", "Size (SMB) Exposure"),
    ("HML", "Value (HML) Exposure"),
    ("Mom", "Momentum (MOM) Exposure"),
]

for ax, (factor, title) in zip(axes, factor_plot_pairs):
    if factor in roll_iwf.columns:
        ax.plot(roll_iwf.index, roll_iwf[factor], label="IWF (Growth)", color="C0")
    if factor in roll_iwd.columns:
        ax.plot(roll_iwd.index, roll_iwd[factor], label="IWD (Value)", color="C3")
    ax.axhline(0, color="black", linewidth=0.6, linestyle="--")
    ax.set_title(title)
    ax.set_ylabel("Beta")
    ax.legend()
    ax.grid(True)

axes[-1].set_xlabel("Date")
plt.suptitle("Rolling 36-Month Factor Loadings: Growth vs Value ETFs (Carhart 4-Factor)", y=1.01)
plt.tight_layout()
plt.savefig("growth_vs_value_factor_loadings.png", dpi=150, bbox_inches="tight")
plt.show()

# ── Bar chart: factor exposure comparison at full-sample ─────────────────────
bar_factors = list(results_summary["IWF"]["betas"].keys())
iwf_betas = [results_summary["IWF"]["betas"][f] for f in bar_factors]
iwd_betas = [results_summary["IWD"]["betas"][f] for f in bar_factors]

x = np.arange(len(bar_factors))
width = 0.35
fig2, ax2 = plt.subplots(figsize=(9, 5))
ax2.bar(x - width / 2, iwf_betas, width, label="IWF (Growth)", color="steelblue")
ax2.bar(x + width / 2, iwd_betas, width, label="IWD (Value)", color="firebrick")
ax2.set_xticks(x)
ax2.set_xticklabels(bar_factors)
ax2.set_ylabel("Factor Loading (Beta)")
ax2.set_title("Full-Sample Factor Loadings: Growth vs Value (Carhart 4-Factor)")
ax2.legend()
ax2.axhline(0, color="black", linewidth=0.6)
ax2.grid(True, axis="y")
plt.tight_layout()
plt.savefig("factor_loadings_bar.png", dpi=150, bbox_inches="tight")
plt.show()

---

Notes and Best Practices

• Newey-West standard errors: For time-series regressions, heteroskedasticity and autocorrelation-consistent (HAC) standard errors improve inference. Pass cov_type='HAC' and cov_kwds={'maxlags': 6} to sm.OLS.fit().
• Frequency consistency: Ensure prices, factors, and risk-free rates all use the same frequency. Do not mix daily prices with monthly factor data.
• Stale RF alignment: When resampling prices to monthly frequency, use month-end (ME) and align the RF series by index before subtracting.
• GRS assumptions: The GRS test assumes i.i.d. multivariate normal residuals and requires T > N + K. Verify this constraint before applying the test.
• Look-ahead bias: When running rolling regressions for backtesting, ensure you are using only data available at each estimation date.
• Factor data license: Fama-French factors from Kenneth French's website are freely available for academic research. Commercial use may require separate licensing agreements.

---

Dependencies Installation

pip install numpy>=1.24.0 pandas>=2.0.0 pandas-datareader>=0.10.0 \
            statsmodels>=0.14.0 scipy>=1.11.0 matplotlib>=3.7.0 yfinance>=0.2.0