用Python和Pandas实现Fama-French三因子模型与股票特质波动率计算

金融数据分析中，Fama-French三因子模型是衡量股票收益与风险关系的重要工具。本文将带你从零开始，用Python和Pandas实现这一经典模型，并计算股票特质波动率——这一指标在量化投资和学术研究中被广泛用于衡量个股的非系统性风险。

1. 环境准备与数据获取

在开始之前，确保你的Python环境已安装以下库：

pip install pandas numpy statsmodels matplotlib

金融数据通常来自专业数据库。假设你已从锐思数据库(RESSET)获取以下CSV文件：

• stock_returns.csv ：包含各股票日收益率数据
• ff3_factors.csv ：包含市场因子(MKT)、规模因子(SMB)、价值因子(HML)的日频数据

提示：实际应用中可能遇到数据缺失或格式不一致问题，建议先用Excel或文本编辑器快速检查原始文件结构。

2. 数据预处理关键步骤

2.1 读取并清洗股票收益率数据

import pandas as pd

# 读取股票收益率数据
stock_ret = pd.read_csv('stock_returns.csv', parse_dates=['date'])
stock_ret.set_index('date', inplace=True)

# 处理缺失值
print(f"原始数据缺失值占比: {stock_ret.isnull().mean().mean():.2%}")
stock_ret = stock_ret.dropna(axis=1, thresh=int(0.8*len(stock_ret)))  # 删除缺失超过20%的股票
stock_ret = stock_ret.fillna(0)  # 简单用0填充剩余缺失

2.2 处理三因子数据

# 读取Fama-French三因子数据
ff3 = pd.read_csv('ff3_factors.csv', parse_dates=['date'])
ff3.set_index('date', inplace=True)

# 确保日期范围匹配
start_date = max(stock_ret.index.min(), ff3.index.min())
end_date = min(stock_ret.index.max(), ff3.index.max())
stock_ret = stock_ret.loc[start_date:end_date]
ff3 = ff3.loc[start_date:end_date]

3. 实现Fama-French三因子模型

3.1 单只股票回归分析

我们先以单只股票为例演示回归过程：

import statsmodels.api as sm

# 选择示例股票
sample_stock = '600000.SH'  
excess_ret = stock_ret[sample_stock] - ff3['RF']  # 计算超额收益

# 准备解释变量
X = ff3[['MKT', 'SMB', 'HML']]
X = sm.add_constant(X)  # 添加常数项

# 执行OLS回归
model = sm.OLS(excess_ret, X, missing='drop')
results = model.fit()
print(results.summary())

3.2 批量处理所有股票

实际应用中我们需要计算所有股票的特质波动率：

residuals_dict = {}

for stock in stock_ret.columns:
    try:
        excess_ret = stock_ret[stock] - ff3['RF']
        X = sm.add_constant(ff3[['MKT', 'SMB', 'HML']])
        model = sm.OLS(excess_ret, X, missing='drop')
        results = model.fit()
        residuals_dict[stock] = results.resid
    except:
        print(f"处理股票 {stock} 时出错")
        continue

residuals_df = pd.DataFrame(residuals_dict)

4. 计算特质波动率

特质波动率是回归残差的标准差，反映无法被三因子解释的风险：

# 计算年化特质波动率
idiosyncratic_vol = residuals_df.std() * np.sqrt(252)  # 年化处理

# 结果可视化
import matplotlib.pyplot as plt

plt.figure(figsize=(10,6))
idiosyncratic_vol.sort_values().plot(kind='bar')
plt.title('各股票特质波动率比较')
plt.ylabel('年化特质波动率')
plt.show()

5. 实际应用中的注意事项

1. 数据对齐问题 ：

   • 确保股票收益率与因子数据的日期完全匹配
   • 处理停牌股票的特殊情况

2. 回归稳定性优化 ：

   • 考虑使用滚动窗口回归捕捉时变特征
   • 对极端值进行Winsorize处理

3. 性能优化技巧 ：

   • 对大规模股票池，考虑并行计算
   • 使用 joblib 或 multiprocessing 加速循环

# 并行计算示例
from joblib import Parallel, delayed

def calculate_ivol(stock):
    excess_ret = stock_ret[stock] - ff3['RF']
    X = sm.add_constant(ff3[['MKT', 'SMB', 'HML']])
    model = sm.OLS(excess_ret, X, missing='drop')
    return model.fit().resid

results = Parallel(n_jobs=4)(delayed(calculate_ivol)(stock) for stock in stock_ret.columns)
residuals_df = pd.concat(results, axis=1)

6. 结果分析与应用

计算得到的特质波动率可用于：

• 构建低波动率投资组合
• 作为风险因子加入多因子模型
• 研究"特质波动率之谜"现象

比较不同行业特质波动率的差异：

# 假设我们有行业分类数据
industry_data = pd.read_csv('industry_classification.csv', index_col=0)
idiosyncratic_vol.name = 'ivol'
result = pd.concat([idiosyncratic_vol, industry_data], axis=1)

# 按行业分组统计
industry_ivol = result.groupby('industry')['ivol'].mean().sort_values()

# 可视化展示
plt.figure(figsize=(12,6))
industry_ivol.plot(kind='bar')
plt.title('各行业平均特质波动率比较')
plt.ylabel('年化特质波动率')
plt.xticks(rotation=45)
plt.show()

在最近一个实际项目中，我们发现科技类股票的特质波动率普遍高于公用事业类股票，这与传统金融理论相符。但有趣的是，某些传统行业也表现出较高的特质波动，这可能反映了行业转型期的特殊风险特征。