现在有2个DataFrame，第一个名为tesla_monthly，其中Adj Close列存储了特斯拉月度收盘价，第二个名为f5_df，是月度五因子数据，包括的列为Mkt-RF、SMB、HML、RMW、CMA、RF，需要用这两个DataFrame演示Fama-French五因子模型，请给出示例代码。

import pandas as pd
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
import seaborn as sns
from matplotlib.ticker import PercentFormatter

# 1. 数据准备与预处理
# 假设已有三个DataFrame:
#   tesla_monthly: 包含日期索引和'Adj Close'列
#   factors: 包含日期索引和列['Mkt-RF','SMB','HML','RF']
#   mom_df: 包含日期索引和动量因子数据

# 1.1 计算特斯拉月度收益率import pandas as pd
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
import seaborn as sns
from statsmodels.stats.stattools import jarque_bera
from statsmodels.stats.diagnostic import het_white

# 1. 数据准备与预处理
# 假设tesla_monthly包含日期索引和'Adj Close'列
# 假设f5_df包含日期索引和列['Mkt-RF','SMB','HML','RMW','CMA','RF']

# 计算特斯拉月度收益率
tesla_returns = tesla_monthly['Adj Close'].pct_change().dropna()
tesla_returns.name = 'TSLA_Return'

# 合并因子数据（确保日期对齐）
ff5_data = pd.concat([
    tesla_returns,
    f5_df
], axis=1).dropna()

# 计算特斯拉超额收益
ff5_data['TSLA_Excess'] = ff5_data['TSLA_Return'] - ff5_data['RF']

# 确认列名
print("数据包含的列:", ff5_data.columns.tolist())
print("\n前5行数据示例:")
print(ff5_data.head())

# 2. Fama-French五因子模型回归
def run_ff5_regression(data):
    """执行Fama-French五因子回归"""
    # 提取五因子：市场因子Mkt-RF，规模因子SMB，价值因子HML，盈利能力因子RMW，投资因子CMA
    X = data[['Mkt-RF', 'SMB', 'HML', 'RMW', 'CMA']]
    X = sm.add_constant(X)  # 添加截距项
    y = data['TSLA_Excess']
  
    model = sm.OLS(y, X).fit(cov_type='HAC', cov_kwds={'maxlags': 3})  # 使用异方差自相关稳健标准误
    return model

# 全样本回归
full_model = run_ff5_regression(ff5_data)

# 3. 滚动估计（60个月窗口）
WINDOW_SIZE = 60  # 5年滚动窗口
dates = []
alphas = []
betas_mkt = []
betas_smb = []
betas_hml = []
betas_rmw = []
betas_cma = []
r_squared = []

print("\n开始滚动五因子模型估计...")
for i in range(WINDOW_SIZE, len(ff5_data)):
    window_data = ff5_data.iloc[i-WINDOW_SIZE:i]
    model = run_ff5_regression(window_data)
  
    # 存储结果
    dates.append(ff5_data.index[i])
    alphas.append(model.params['const'])
    betas_mkt.append(model.params['Mkt-RF'])
    betas_smb.append(model.params['SMB'])
    betas_hml.append(model.params['HML'])
    betas_rmw.append(model.params['RMW'])
    betas_cma.append(model.params['CMA'])
    r_squared.append(model.rsquared)

# 创建结果DataFrame
results = pd.DataFrame({
    'Date': dates,
    'Alpha': alphas,
    'Beta_Mkt': betas_mkt,
    'Beta_SMB': betas_smb,
    'Beta_HML': betas_hml,
    'Beta_RMW': betas_rmw,
    'Beta_CMA': betas_cma,
    'R_squared': r_squared
}).set_index('Date')

# 4. 模型诊断与可视化
plt.figure(figsize=(18, 20))

# 4.1 Alpha系数变化
plt.subplot(4, 2, 1)
plt.plot(results.index, results['Alpha'], 'b-', lw=1.5)
plt.axhline(results['Alpha'].mean(), color='r', linestyle='--', label=f'平均Alpha={results["Alpha"].mean():.2f}%')
plt.axhline(0, color='k', linestyle='--', lw=0.8)
plt.title('特斯拉Alpha系数变化 (滚动60个月)', fontsize=16)
plt.ylabel('Alpha (%)', fontsize=14)
plt.legend()
plt.grid(alpha=0.3)

# 4.2 五因子暴露变化
plt.subplot(4, 2, 2)
plt.plot(results.index, results['Beta_Mkt'], 'b-', label='市场因子(Mkt-RF)')
plt.plot(results.index, results['Beta_SMB'], 'g-', label='规模因子(SMB)')
plt.plot(results.index, results['Beta_HML'], 'r-', label='价值因子(HML)')
plt.plot(results.index, results['Beta_RMW'], 'c-', label='盈利能力因子(RMW)')
plt.plot(results.index, results['Beta_CMA'], 'm-', label='投资因子(CMA)')
plt.axhline(0, color='k', linestyle='--', lw=0.8)
plt.title('五因子暴露变化', fontsize=16)
plt.ylabel('Beta系数', fontsize=14)
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.grid(alpha=0.3)

# 4.3 模型解释力变化
plt.subplot(4, 2, 3)
plt.plot(results.index, results['R_squared'], 'purple')
plt.title('模型解释力变化 (R²)', fontsize=16)
plt.ylabel('R²', fontsize=14)
plt.ylim(0, 1)
plt.grid(alpha=0.3)

# 4.4 残差诊断 - 正态性检验
from statsmodels.stats.stattools import jarque_bera

# 正确解包所有返回值
jb_stat, jb_pvalue, skewness, kurtosis = jarque_bera(full_model.resid)

plt.subplot(4, 2, 4)
sns.histplot(full_model.resid, kde=True, color='skyblue')
plt.title(f'残差分布 (偏度={skewness:.2f}, 峰度={kurtosis:.2f}, p值={jb_pvalue:.4f})', fontsize=12)  # 添加更多统计量
plt.xlabel('残差', fontsize=10)
plt.grid(alpha=0.3)

# 4.5 残差诊断 - 异方差检验
from statsmodels.stats.diagnostic import het_white

# 获取White检验结果（返回4个值）
white_test = het_white(full_model.resid, full_model.model.exog)
white_stat, white_pvalue, _, _ = white_test  # 正确解包

plt.subplot(4, 2, 5)
plt.scatter(full_model.fittedvalues, full_model.resid, alpha=0.5)
plt.axhline(0, color='red', linestyle='--')
plt.title(f'残差 vs 拟合值 (White检验 p值={white_pvalue:.4f})', fontsize=12)
plt.xlabel('拟合值', fontsize=10)
plt.ylabel('残差', fontsize=10)
plt.grid(alpha=0.3)

# 4.6 因子暴露相关性热图
plt.subplot(4, 2, 6)
factor_betas = results[['Beta_Mkt', 'Beta_SMB', 'Beta_HML', 'Beta_RMW', 'Beta_CMA']].corr()
sns.heatmap(factor_betas, annot=True, cmap='coolwarm', fmt=".2f", vmin=-1, vmax=1)
plt.title('因子暴露相关性热图', fontsize=16)

# 4.7 市场暴露 vs Alpha
plt.subplot(4, 2, 7)
sns.scatterplot(x=results['Beta_Mkt'], y=results['Alpha'], 
                size=np.abs(results['R_squared']), sizes=(20, 200),
                hue=results['Beta_CMA'], palette='viridis', alpha=0.7)
plt.axhline(0, color='grey', linestyle='--')
plt.axvline(1, color='grey', linestyle='--')
plt.title('市场风险暴露 vs Alpha', fontsize=16)
plt.xlabel('市场风险暴露(Beta_Mkt)', fontsize=14)
plt.ylabel('Alpha (%)', fontsize=14)
plt.grid(alpha=0.3)

# 4.8 模型预测 vs 实际值
plt.subplot(4, 2, 8)
plt.scatter(full_model.fittedvalues, ff5_data['TSLA_Excess'], alpha=0.5)
plt.plot([min(full_model.fittedvalues), max(full_model.fittedvalues)], 
         [min(full_model.fittedvalues), max(full_model.fittedvalues)], 
         color='red')
plt.title(f'模型预测 vs 实际值 (R²={full_model.rsquared:.2f})', fontsize=16)
plt.xlabel('预测超额收益(%)', fontsize=14)
plt.ylabel('实际超额收益(%)', fontsize=14)
plt.grid(alpha=0.3)

plt.tight_layout()
plt.savefig('tesla_fama_french_5_factor.png', dpi=300)
plt.show()

# 5. 输出全样本回归结果
print("\nFama-French五因子模型回归结果:")
print(full_model.summary())

# 6. 特斯拉特征分析
print("\n特斯拉股票特征分析:")
print(f"平均Alpha: {results['Alpha'].mean():.2f}%")
print(f"市场风险暴露: {results['Beta_Mkt'].mean():.2f}")
print(f"规模因子暴露: {results['Beta_SMB'].mean():.2f}")
print(f"价值因子暴露: {results['Beta_HML'].mean():.2f}")
print(f"盈利能力暴露: {results['Beta_RMW'].mean():.2f}")
print(f"投资模式暴露: {results['Beta_CMA'].mean():.2f}")
print(f"平均模型解释力(R²): {results['R_squared'].mean():.2%}")

# 7. 因子暴露的时间趋势分析
trend_df = results.rolling(24).mean().dropna()  # 2年移动平均

plt.figure(figsize=(14, 10))
for factor in ['Beta_Mkt', 'Beta_SMB', 'Beta_HML', 'Beta_RMW', 'Beta_CMA']:
    plt.plot(trend_df.index, trend_df[factor], label=factor[5:], lw=2)

plt.title('特斯拉因子暴露趋势 (2年移动平均)', fontsize=16)
plt.ylabel('Beta系数', fontsize=14)
plt.legend()
plt.grid(alpha=0.3)

# 添加重大事件标记
events = {
    '2018-07': 'Model 3量产',
    '2020-03': 'COVID冲击',
    '2021-10': '市值破万亿美元',
    '2023-01': '大幅降价'
}

for date, event in events.items():
    if date in trend_df.index:
        plt.axvline(pd.Timestamp(date), color='gray', linestyle='--', alpha=0.7)
        plt.text(pd.Timestamp(date), plt.ylim()[0]+0.05, event, rotation=45, fontsize=9)

plt.tight_layout()
plt.show()

# 8. 模型诊断结论
print("\n模型诊断结论:")
if jb_pvalue < 0.05:
    print("→ 残差显著偏离正态分布 (p={:.4f})，模型可能有遗漏变量".format(jb_pvalue))
else:
    print("→ 残差满足正态分布假设 (p={:.4f})".format(jb_pvalue))

if white_pvalue < 0.05:
    print("→ 存在显著的异方差问题 (p={:.4f})，结果需谨慎解读".format(white_pvalue))
else:
    print("→ 无异方差问题 (p={:.4f})".format(white_pvalue))

# 9. 投资含义解读
if results['Beta_RMW'].mean() > 0:
    print("\n→ 特斯拉对盈利能力因子有正向暴露，表明其具有高盈利增长特征")
else:
    print("\n→ 特斯拉对盈利能力因子有负向暴露，与传统高盈利能力公司特征不符")

if results['Beta_CMA'].mean() > 0:
    print("→ 特斯拉对投资因子有正向暴露，表明其保守的投资行为")
else:
    print("→ 特斯拉对投资因子有负向暴露，表明其积极的投资扩张策略")

if results['Alpha'].mean() > 0.5:
    print("→ 特斯拉持续产生显著的正Alpha，表明技术优势创造了超额收益")
elif results['Alpha'].mean() < -0.5:
    print("→ 特斯拉产生负Alpha，估值溢价未被基本面支撑")
else:
    print("→ 特斯拉未产生显著的超额收益")import pandas as pd
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
import seaborn as sns
from statsmodels.stats.stattools import jarque_bera
from statsmodels.stats.diagnostic import het_white

# 1. 数据准备与预处理
# 假设tesla_monthly包含日期索引和'Adj Close'列
# 假设f5_df包含日期索引和列['Mkt-RF','SMB','HML','RMW','CMA','RF']

# 计算特斯拉月度收益率
tesla_returns = tesla_monthly['Adj Close'].pct_change().dropna()
tesla_returns.name = 'TSLA_Return'

# 合并因子数据（确保日期对齐）
ff5_data = pd.concat([
    tesla_returns,
    f5_df
], axis=1).dropna()

# 计算特斯拉超额收益
ff5_data['TSLA_Excess'] = ff5_data['TSLA_Return'] - ff5_data['RF']

# 确认列名
print("数据包含的列:", ff5_data.columns.tolist())
print("\n前5行数据示例:")
print(ff5_data.head())

# 2. Fama-French五因子模型回归
def run_ff5_regression(data):
    """执行Fama-French五因子回归"""
    # 提取五因子：市场因子Mkt-RF，规模因子SMB，价值因子HML，盈利能力因子RMW，投资因子CMA
    X = data[['Mkt-RF', 'SMB', 'HML', 'RMW', 'CMA']]
    X = sm.add_constant(X)  # 添加截距项
    y = data['TSLA_Excess']
  
    model = sm.OLS(y, X).fit(cov_type='HAC', cov_kwds={'maxlags': 3})  # 使用异方差自相关稳健标准误
    return model

# 全样本回归
full_model = run_ff5_regression(ff5_data)

# 3. 滚动估计（60个月窗口）
WINDOW_SIZE = 60  # 5年滚动窗口
dates = []
alphas = []
betas_mkt = []
betas_smb = []
betas_hml = []
betas_rmw = []
betas_cma = []
r_squared = []

print("\n开始滚动五因子模型估计...")
for i in range(WINDOW_SIZE, len(ff5_data)):
    window_data = ff5_data.iloc[i-WINDOW_SIZE:i]
    model = run_ff5_regression(window_data)
  
    # 存储结果
    dates.append(ff5_data.index[i])
    alphas.append(model.params['const'])
    betas_mkt.append(model.params['Mkt-RF'])
    betas_smb.append(model.params['SMB'])
    betas_hml.append(model.params['HML'])
    betas_rmw.append(model.params['RMW'])
    betas_cma.append(model.params['CMA'])
    r_squared.append(model.rsquared)

# 创建结果DataFrame
results = pd.DataFrame({
    'Date': dates,
    'Alpha': alphas,
    'Beta_Mkt': betas_mkt,
    'Beta_SMB': betas_smb,
    'Beta_HML': betas_hml,
    'Beta_RMW': betas_rmw,
    'Beta_CMA': betas_cma,
    'R_squared': r_squared
}).set_index('Date')

# 4. 模型诊断与可视化
plt.figure(figsize=(18, 20))

# 4.1 Alpha系数变化
plt.subplot(4, 2, 1)
plt.plot(results.index, results['Alpha'], 'b-', lw=1.5)
plt.axhline(results['Alpha'].mean(), color='r', linestyle='--', label=f'平均Alpha={results["Alpha"].mean():.2f}%')
plt.axhline(0, color='k', linestyle='--', lw=0.8)
plt.title('特斯拉Alpha系数变化 (滚动60个月)', fontsize=16)
plt.ylabel('Alpha (%)', fontsize=14)
plt.legend()
plt.grid(alpha=0.3)

# 4.2 五因子暴露变化
plt.subplot(4, 2, 2)
plt.plot(results.index, results['Beta_Mkt'], 'b-', label='市场因子(Mkt-RF)')
plt.plot(results.index, results['Beta_SMB'], 'g-', label='规模因子(SMB)')
plt.plot(results.index, results['Beta_HML'], 'r-', label='价值因子(HML)')
plt.plot(results.index, results['Beta_RMW'], 'c-', label='盈利能力因子(RMW)')
plt.plot(results.index, results['Beta_CMA'], 'm-', label='投资因子(CMA)')
plt.axhline(0, color='k', linestyle='--', lw=0.8)
plt.title('五因子暴露变化', fontsize=16)
plt.ylabel('Beta系数', fontsize=14)
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.grid(alpha=0.3)

# 4.3 模型解释力变化
plt.subplot(4, 2, 3)
plt.plot(results.index, results['R_squared'], 'purple')
plt.title('模型解释力变化 (R²)', fontsize=16)
plt.ylabel('R²', fontsize=14)
plt.ylim(0, 1)
plt.grid(alpha=0.3)

# 4.4 残差诊断 - 正态性检验
from statsmodels.stats.stattools import jarque_bera

# 正确解包所有返回值
jb_stat, jb_pvalue, skewness, kurtosis = jarque_bera(full_model.resid)

plt.subplot(4, 2, 4)
sns.histplot(full_model.resid, kde=True, color='skyblue')
plt.title(f'残差分布 (偏度={skewness:.2f}, 峰度={kurtosis:.2f}, p值={jb_pvalue:.4f})', fontsize=12)  # 添加更多统计量
plt.xlabel('残差', fontsize=10)
plt.grid(alpha=0.3)

# 4.5 残差诊断 - 异方差检验
from statsmodels.stats.diagnostic import het_white

# 获取White检验结果（返回4个值）
white_test = het_white(full_model.resid, full_model.model.exog)
white_stat, white_pvalue, _, _ = white_test  # 正确解包

plt.subplot(4, 2, 5)
plt.scatter(full_model.fittedvalues, full_model.resid, alpha=0.5)
plt.axhline(0, color='red', linestyle='--')
plt.title(f'残差 vs 拟合值 (White检验 p值={white_pvalue:.4f})', fontsize=12)
plt.xlabel('拟合值', fontsize=10)
plt.ylabel('残差', fontsize=10)
plt.grid(alpha=0.3)

# 4.6 因子暴露相关性热图
plt.subplot(4, 2, 6)
factor_betas = results[['Beta_Mkt', 'Beta_SMB', 'Beta_HML', 'Beta_RMW', 'Beta_CMA']].corr()
sns.heatmap(factor_betas, annot=True, cmap='coolwarm', fmt=".2f", vmin=-1, vmax=1)
plt.title('因子暴露相关性热图', fontsize=16)

# 4.7 市场暴露 vs Alpha
plt.subplot(4, 2, 7)
sns.scatterplot(x=results['Beta_Mkt'], y=results['Alpha'], 
                size=np.abs(results['R_squared']), sizes=(20, 200),
                hue=results['Beta_CMA'], palette='viridis', alpha=0.7)
plt.axhline(0, color='grey', linestyle='--')
plt.axvline(1, color='grey', linestyle='--')
plt.title('市场风险暴露 vs Alpha', fontsize=16)
plt.xlabel('市场风险暴露(Beta_Mkt)', fontsize=14)
plt.ylabel('Alpha (%)', fontsize=14)
plt.grid(alpha=0.3)

# 4.8 模型预测 vs 实际值
plt.subplot(4, 2, 8)
plt.scatter(full_model.fittedvalues, ff5_data['TSLA_Excess'], alpha=0.5)
plt.plot([min(full_model.fittedvalues), max(full_model.fittedvalues)], 
         [min(full_model.fittedvalues), max(full_model.fittedvalues)], 
         color='red')
plt.title(f'模型预测 vs 实际值 (R²={full_model.rsquared:.2f})', fontsize=16)
plt.xlabel('预测超额收益(%)', fontsize=14)
plt.ylabel('实际超额收益(%)', fontsize=14)
plt.grid(alpha=0.3)

plt.tight_layout()
plt.savefig('tesla_fama_french_5_factor.png', dpi=300)
plt.show()

# 5. 输出全样本回归结果
print("\nFama-French五因子模型回归结果:")
print(full_model.summary())

# 6. 特斯拉特征分析
print("\n特斯拉股票特征分析:")
print(f"平均Alpha: {results['Alpha'].mean():.2f}%")
print(f"市场风险暴露: {results['Beta_Mkt'].mean():.2f}")
print(f"规模因子暴露: {results['Beta_SMB'].mean():.2f}")
print(f"价值因子暴露: {results['Beta_HML'].mean():.2f}")
print(f"盈利能力暴露: {results['Beta_RMW'].mean():.2f}")
print(f"投资模式暴露: {results['Beta_CMA'].mean():.2f}")
print(f"平均模型解释力(R²): {results['R_squared'].mean():.2%}")

# 7. 因子暴露的时间趋势分析
trend_df = results.rolling(24).mean().dropna()  # 2年移动平均

plt.figure(figsize=(14, 10))
for factor in ['Beta_Mkt', 'Beta_SMB', 'Beta_HML', 'Beta_RMW', 'Beta_CMA']:
    plt.plot(trend_df.index, trend_df[factor], label=factor[5:], lw=2)

plt.title('特斯拉因子暴露趋势 (2年移动平均)', fontsize=16)
plt.ylabel('Beta系数', fontsize=14)
plt.legend()
plt.grid(alpha=0.3)

# 添加重大事件标记
events = {
    '2018-07': 'Model 3量产',
    '2020-03': 'COVID冲击',
    '2021-10': '市值破万亿美元',
    '2023-01': '大幅降价'
}

for date, event in events.items():
    if date in trend_df.index:
        plt.axvline(pd.Timestamp(date), color='gray', linestyle='--', alpha=0.7)
        plt.text(pd.Timestamp(date), plt.ylim()[0]+0.05, event, rotation=45, fontsize=9)

plt.tight_layout()
plt.show()

# 8. 模型诊断结论
print("\n模型诊断结论:")
if jb_pvalue < 0.05:
    print("→ 残差显著偏离正态分布 (p={:.4f})，模型可能有遗漏变量".format(jb_pvalue))
else:
    print("→ 残差满足正态分布假设 (p={:.4f})".format(jb_pvalue))

if white_pvalue < 0.05:
    print("→ 存在显著的异方差问题 (p={:.4f})，结果需谨慎解读".format(white_pvalue))
else:
    print("→ 无异方差问题 (p={:.4f})".format(white_pvalue))

# 9. 投资含义解读
if results['Beta_RMW'].mean() > 0:
    print("\n→ 特斯拉对盈利能力因子有正向暴露，表明其具有高盈利增长特征")
else:
    print("\n→ 特斯拉对盈利能力因子有负向暴露，与传统高盈利能力公司特征不符")

if results['Beta_CMA'].mean() > 0:
    print("→ 特斯拉对投资因子有正向暴露，表明其保守的投资行为")
else:
    print("→ 特斯拉对投资因子有负向暴露，表明其积极的投资扩张策略")

if results['Alpha'].mean() > 0.5:
    print("→ 特斯拉持续产生显著的正Alpha，表明技术优势创造了超额收益")
elif results['Alpha'].mean() < -0.5:
    print("→ 特斯拉产生负Alpha，估值溢价未被基本面支撑")
else:
    print("→ 特斯拉未产生显著的超额收益")
tesla_returns = tesla_monthly['Adj Close'].pct_change().dropna()
tesla_returns.name = 'TSLA_Return'

# 1.2 合并所有因子数据
# 确保所有数据频率相同（月度）
carhart_data = pd.concat([
    tesla_returns, 
    factors,
    mom_df
], axis=1).dropna()

# 1.3 计算特斯拉超额收益
carhart_data['TSLA_Excess'] = carhart_data['TSLA_Return'] - carhart_data['RF']

# 1.4 确认因子列名
# 假设动量因子列名为'MOM'
print("数据包含的列:", carhart_data.columns.tolist())
print("\n前5行数据示例:")
print(carhart_data.head())

# 2. Carhart四因子模型回归
def run_carhart_regression(data):
    """执行Carhart四因子回归"""
    X = data[['Mkt-RF', 'SMB', 'HML', 'Mom']]
    X = sm.add_constant(X)  # 添加截距项
    y = data['TSLA_Excess']
  
    model = sm.OLS(y, X).fit(cov_type='HAC', cov_kwds={'maxlags': 3})
    return model

# 全样本回归
full_model = run_carhart_regression(carhart_data)

# 3. 滚动估计（36个月窗口）
WINDOW_SIZE = 36  # 3年滚动窗口
dates = []
alphas = []
betas_mkt = []
betas_smb = []
betas_hml = []
betas_mom = []
r_squared = []
pvalues_alpha = []

print("\n开始滚动四因子模型估计...")
for i in range(WINDOW_SIZE, len(carhart_data)):
    window_data = carhart_data.iloc[i-WINDOW_SIZE:i]
    model = run_carhart_regression(window_data)
  
    # 存储结果
    dates.append(carhart_data.index[i])
    alphas.append(model.params['const'])
    betas_mkt.append(model.params['Mkt-RF'])
    betas_smb.append(model.params['SMB'])
    betas_hml.append(model.params['HML'])
    betas_mom.append(model.params['Mom'])
    r_squared.append(model.rsquared)
    pvalues_alpha.append(model.pvalues['const'])

# 创建结果DataFrame
results = pd.DataFrame({
    'Date': dates,
    'Alpha': alphas,
    'Beta_Mkt': betas_mkt,
    'Beta_SMB': betas_smb,
    'Beta_HML': betas_hml,
    'Beta_MOM': betas_mom,
    'R_squared': r_squared,
    'Alpha_Pvalue': pvalues_alpha
}).set_index('Date')

# 4. 可视化分析
plt.figure(figsize=(18, 15))

# 4.1 Alpha系数变化 (标记统计显著性)
plt.subplot(4, 2, 1)
plt.plot(results.index, results['Alpha'], 'b-', lw=1.5, label='Alpha')

# 标记显著的Alpha点（p<0.05）
significant_dates = results[results['Alpha_Pvalue'] < 0.05].index
significant_alphas = results.loc[significant_dates, 'Alpha']
plt.scatter(significant_dates, significant_alphas, 
            color='red', s=30, label='显著Alpha (p<0.05)')

plt.axhline(0, color='black', linestyle='--', lw=0.8)
plt.title('特斯拉Alpha系数变化 (滚动36个月)', fontsize=14)
plt.ylabel('Alpha (%)', fontsize=12)
plt.legend()
plt.grid(alpha=0.3)

# 4.2 因子暴露变化
plt.subplot(4, 2, 2)
plt.plot(results.index, results['Beta_Mkt'], 'b-', label='市场因子(Beta_Mkt)')
plt.plot(results.index, results['Beta_SMB'], 'g-', label='规模因子(Beta_SMB)')
plt.plot(results.index, results['Beta_HML'], 'r-', label='价值因子(Beta_HML)')
plt.plot(results.index, results['Beta_MOM'], 'm-', label='动量因子(Beta_MOM)')
plt.axhline(0, color='black', linestyle='--', lw=0.8)
plt.title('因子暴露变化', fontsize=14)
plt.ylabel('Beta系数', fontsize=12)
plt.legend(loc='best')
plt.grid(alpha=0.3)

# 4.3 R²变化
plt.subplot(4, 2, 3)
plt.plot(results.index, results['R_squared'], 'purple')
plt.title('模型解释力变化 (R²)', fontsize=14)
plt.ylabel('R²', fontsize=12)
plt.ylim(0, 1)
plt.grid(alpha=0.3)

# 4.4 Beta_Mkt与Beta_MOM关系
plt.subplot(4, 2, 4)
sns.regplot(x=results['Beta_Mkt'], y=results['Beta_MOM'], 
            scatter_kws={'alpha':0.5}, line_kws={'color':'red'})
plt.title('市场风险暴露 vs 动量因子暴露', fontsize=14)
plt.xlabel('市场风险暴露(Beta_Mkt)')
plt.ylabel('动量因子暴露(Beta_MOM)')
plt.grid(alpha=0.3)

# 4.5 Alpha分布直方图
plt.subplot(4, 2, 5)
plt.hist(results['Alpha'], bins=20, color='skyblue', alpha=0.8)
plt.axvline(results['Alpha'].mean(), color='r', linestyle='--', 
            label=f'平均值={results["Alpha"].mean():.2f}%')
plt.title('Alpha分布', fontsize=14)
plt.xlabel('Alpha (%)')
plt.legend()
plt.grid(alpha=0.3)

# 4.6 动量因子暴露分布
plt.subplot(4, 2, 6)
plt.hist(results['Beta_MOM'], bins=20, color='purple', alpha=0.8)
plt.axvline(results['Beta_MOM'].mean(), color='r', linestyle='--',
            label=f'平均值={results["Beta_MOM"].mean():.2f}')
plt.title('动量因子暴露分布', fontsize=14)
plt.xlabel('Beta_MOM')
plt.legend()
plt.grid(alpha=0.3)

# 4.7 市场风险暴露 vs Alpha
plt.subplot(4, 2, 7)
sns.scatterplot(x=results['Beta_Mkt'], y=results['Alpha'], 
                hue=results['R_squared'], palette='viridis')
plt.axhline(0, color='grey', linestyle='--')
plt.axvline(1, color='grey', linestyle='--')
plt.title('市场风险暴露 vs Alpha', fontsize=14)
plt.xlabel('市场风险暴露(Beta_Mkt)')
plt.ylabel('Alpha (%)')
plt.grid(alpha=0.3)

# 4.8 最终模型预测vs实际值
plt.subplot(4, 2, 8)
plt.scatter(full_model.fittedvalues, carhart_data['TSLA_Excess'], alpha=0.5)
plt.plot([min(full_model.fittedvalues), max(full_model.fittedvalues)], 
         [min(full_model.fittedvalues), max(full_model.fittedvalues)], 
         color='red')
plt.title(f'模型预测 vs 实际值 (R²={full_model.rsquared:.2f})', fontsize=14)
plt.xlabel('预测超额收益(%)')
plt.ylabel('实际超额收益(%)')
plt.grid(alpha=0.3)

plt.tight_layout()
plt.savefig('tesla_carhart_four_factor.png', dpi=300)
plt.show()

# 5. 输出结果分析
print("\n全样本四因子模型结果:")
print(full_model.summary())

# 6. 投资策略分析
avg_alpha = results['Alpha'].mean()
avg_mom_beta = results['Beta_MOM'].mean()

print("\n特斯拉股票特征:")
print(f"平均Alpha: {avg_alpha:.2f}%")
print(f"动量因子暴露: {avg_mom_beta:.2f}")
print(f"市场风险暴露: {results['Beta_Mkt'].mean():.2f}")
print(f"规模因子暴露: {results['Beta_SMB'].mean():.2f}")
print(f"价值因子暴露: {results['Beta_HML'].mean():.2f}")

# 7. 因子暴露相关性分析
beta_corr = results[['Beta_Mkt', 'Beta_SMB', 'Beta_HML', 'Beta_MOM']].corr()
print("\n因子暴露相关性矩阵:")
print(beta_corr)

# 可视化相关性矩阵
plt.figure(figsize=(10, 8))
sns.heatmap(beta_corr, annot=True, cmap='coolwarm', vmin=-1, vmax=1)
plt.title('因子暴露相关性热图', fontsize=14)
plt.tight_layout()
plt.show()

# Python代码
# 合并两个数据框架，确保日期对齐
merged_data = pd.merge(tesla_monthly, 
                       f5_df, 
                       left_index=True, 
                       right_index=True)

# 计算超额收益
merged_data['Excess Return'] = (
    merged_data['Adj Close'].pct_change() - merged_data['RF']
)
merged_data = merged_data.dropna()

# 准备回归所需的数据
X = sm.add_constant(merged_data[['Mkt-RF', 'SMB', 'HML', 'RMW', 'CMA']])
y = merged_data['Excess Return']

# 使用OLS进行回归
model = sm.OLS(y, X)
results = model.fit()

# 打印回归结果
print(results.params)

const     0.039222
Mkt-RF    2.397336
SMB      -1.184509
HML      -0.724754
RMW      -1.441905
CMA      -0.119482
dtype: float64