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
}
defdownload_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).
"""
model == :
ff3 = download_ff_factors(, start, end, freq)
mom = download_ff_factors(, start, end, freq)
combined = ff3.join(mom, how=)
combined
dataset_key = _FF_DATASETS.get(model)
dataset_key :
ValueError()
factors_raw = web.DataReader(
dataset_key,
,
start=start,
end=end,
)
freq == :
df = factors_raw[]
freq == :
df = factors_raw[]
:
ValueError()
df.index = pd.to_datetime(df.index.astype(), =)
df.index.name =
df /
() -> pd.DataFrame:
freq == :
monthly_prices = prices.resample().last()
returns = monthly_prices.pct_change().dropna()
freq == :
returns = prices.pct_change().dropna()
:
ValueError()
rf_aligned = rf.reindex(returns.index).ffill()
excess = returns.subtract(rf_aligned, axis=)
excess.dropna()
if
"4factor"
"3factor"
"momentum"
"inner"
return
if
is
None
raise
f"Unsupported model: {model!r}. Choose from {list(_FF_DATASETS)}"
# pandas_datareader fetches directly from Ken French's website
"famafrench"
# DataReader returns a dict; index 0 is monthly, index 1 is annual
if
"monthly"
0
elif
"daily"
0
# daily datasets only have one table
else
raise
f"freq must be 'monthly' or 'daily', got {freq!r}"
"""
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
"monthly"
# Resample to month-end and compute simple returns
"ME"
elif
"daily"
else
raise
f"freq must be 'monthly' or 'daily', got {freq!r}"
# Align rf index to returns index
0
return
2. OLS Factor Regression
defrun_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,
}
defrolling_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 inrange(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) * 12for 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
defgrs_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 inenumerate(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,
}
deffactor_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 * 100return 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.
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) inzip(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.
