用Python和Pandas实战Fama-French三因子模型：一步步计算股票特质波动率
在量化投资领域，特质波动率（Idiosyncratic Volatility）作为衡量个股特有风险的重要指标，近年来受到学术界和业界的广泛关注。不同于系统性风险，特质波动率捕捉的是无法通过市场因子解释的个股特异性波动，这种"非系统性风险"的度量对于构建多因子模型、优化投资组合具有独特价值。本文将手把手教你如何用Python的Pandas和statsmodels库，从原始数据开始完整实现Fama-French三因子模型下的特质波动率计算流程。
1. 环境准备与数据获取
特质波动率计算的第一步是准备Python环境和获取必要数据。推荐使用Anaconda创建独立的Python 3.8+环境，确保各库版本兼容：
conda create -n ff3 python=3.8 conda activate ff3 pip install pandas numpy statsmodels openpyxl金融数据来源有多种选择，国内常用Wind、CSMAR或RESSET，国际常用CRSP或Kenneth French官网。假设我们从CSMAR获取了以下数据：
- 个股日收益率：包含股票代码(Stkcd)、交易日期(Date)、考虑现金红利再投资的日个股回报率(Dretnd)
- 无风险收益率：通常使用国债日收益率(Nrrdaydt)
- 三因子数据：市场风险溢价(RiskPremium)、规模因子(SMB)、价值因子(HML)
提示：实际应用中需注意数据频率一致性，所有数据应为日频且日期对齐。缺失值处理建议采用向前填充或删除整条记录。
将下载的Excel文件按sheet分开保存，结构如下：
Data/ ├── FF3_Factors.xlsx # 三因子数据 ├── Stock_Returns.xlsx # 个股收益率 └── RiskFree_Rate.xlsx # 无风险利率2. 数据清洗与合并
数据质量决定模型效果，金融数据处理需特别注意日期格式和缺失值：
import pandas as pd from datetime import datetime # 读取原始数据 factors = pd.read_excel('Data/FF3_Factors.xlsx', parse_dates=['Date']) returns = pd.read_excel('Data/Stock_Returns.xlsx', parse_dates=['Date']) rf = pd.read_excel('Data/RiskFree_Rate.xlsx', parse_dates=['Date']) # 统一日期格式 def format_date(df): df['Date'] = pd.to_datetime(df['Date']).dt.normalize() return df factors = format_date(factors) returns = format_date(returns) rf = format_date(rf)使用Pandas的merge函数进行表连接时，注意处理多对多关系：
# 三表合并 merged_data = pd.merge( pd.merge(factors, returns, on='Date', how='inner'), rf, on='Date', how='inner' ) # 检查合并后数据 print(f"合并后数据量：{len(merged_data)}") print(f"日期范围：{merged_data['Date'].min()} 至 {merged_data['Date'].max()}") print(f"包含股票数量：{merged_data['Stkcd'].nunique()}")常见数据问题及处理方法：
| 问题类型 | 检测方法 | 处理方案 |
|---|---|---|
| 日期不连续 | pd.date_range对比 | 填充或删除缺失日期 |
| 收益率异常值 | 3σ原则或分位数检查 | Winsorize处理或设为NaN |
| 因子缺失 | isnull().sum() | 删除或插值补全 |
3. Fama-French三因子模型实现
特质波动率的核心是通过三因子模型分离系统性风险和特异性风险。我们使用statsmodels的OLS进行回归：
import statsmodels.formula.api as smf def ff3_regression(df): """执行Fama-French三因子回归""" model = smf.ols('Dretnd - Nrrdaydt ~ RiskPremium + SMB + HML', data=df) results = model.fit() return results # 示例：单只股票单月回归 sample = merged_data[ (merged_data['Stkcd'] == '600000') & (merged_data['Date'].between('2020-01-01', '2020-01-31')) ].dropna() results = ff3_regression(sample) print(results.summary())回归结果解析要点：
- Adj. R-squared：模型解释力度，通常0.6-0.8为佳
- 系数显著性：p-value应小于0.05
- 残差正态性：Jarque-Bera检验p-value需大于0.05
计算残差和特质波动率的完整流程：
def calculate_residuals(model, df): """计算日频残差序列""" predicted = model.predict(df) actual = df['Dretnd'] - df['Nrrdaydt'] return actual - predicted def idiosyncratic_volatility(residuals, ddof=1): """计算特质波动率""" return residuals.std(ddof=ddof) * (len(residuals)**0.5)4. 批量计算与结果分析
实际研究中需要计算全市场股票的特质波动率，这里展示批量处理方法：
# 添加年月字段便于分组 merged_data['YearMonth'] = merged_data['Date'].dt.to_period('M') # 定义计算函数 def batch_compute_iv(group): try: model = ff3_regression(group) residuals = calculate_residuals(model, group) iv = idiosyncratic_volatility(residuals) return pd.Series({ 'IV': iv, 'Obs': len(group), 'R2': model.rsquared }) except: return pd.Series({ 'IV': None, 'Obs': len(group), 'R2': None }) # 分组计算 iv_results = merged_data.groupby(['Stkcd', 'YearMonth']).apply(batch_compute_iv) iv_results = iv_results.reset_index()结果分析时需关注以下质量指标：
- 观测值数量：单月交易天数应大于15天（假设22个交易日）
- R平方：过低可能预示模型设定问题
- IV极端值：检查是否数据错误或特殊事件
保存结果供后续研究：
# 保存到Excel iv_results.to_excel('Output/Idiosyncratic_Volatility.xlsx', index=False) # 也可保存到数据库 from sqlalchemy import create_engine engine = create_engine('sqlite:///Output/finance.db') iv_results.to_sql('IV_Results', engine, if_exists='replace')5. 高级应用与注意事项
实际应用中，特质波动率计算有几个进阶技巧：
滚动窗口计算：使用60个月滚动窗口提高稳定性
from tqdm import tqdm def rolling_iv(stock_code, window=60): stock_data = merged_data[merged_data['Stkcd'] == stock_code] dates = sorted(stock_data['Date'].unique()) results = [] for i in tqdm(range(window, len(dates))): window_data = stock_data[ stock_data['Date'].between(dates[i-window], dates[i-1]) ] if len(window_data) < window * 0.8: # 缺失值超过20%则跳过 continue model = ff3_regression(window_data) residuals = calculate_residuals(model, window_data) iv = idiosyncratic_volatility(residuals) results.append({ 'Stkcd': stock_code, 'Date': dates[i], 'IV': iv, 'R2': model.rsquared }) return pd.DataFrame(results)行业调整：先按行业分类再计算相对特质波动率
# 假设有行业分类数据 industry_data = pd.read_excel('Data/Industry_Classification.xlsx') iv_industry = pd.merge(iv_results, industry_data, on='Stkcd') # 计算行业中性化IV iv_industry['IV_Adj'] = iv_industry.groupby(['Industry', 'YearMonth'])['IV'].transform( lambda x: x - x.mean() )常见问题排查指南：
回归不收敛：
- 检查因子间多重共线性（VIF > 10需处理）
- 确保收益率数值合理（通常绝对值<0.2）
特质波动率异常高：
- 验证是否使用了对数收益率
- 检查是否包含停牌日数据
结果不稳定：
- 增加滚动窗口长度
- 考虑使用Robust Regression减少异常值影响
在实盘应用中，我们发现特质波动率因子通常需要与其他因子结合使用。例如，将IV与动量因子结合构建多空组合，或者作为风险模型中的控制变量。存储计算结果时，建议同时保存中间指标（如R平方、观测值数量）以便后续质量检验。
