07 超越线性关系

# 若 order_book_id 在索引中则重置为普通列
if 'order_book_id' in all_stock_data.index.names:  # 检查是否为多级索引
    all_stock_data = all_stock_data.reset_index()  # 重置索引为普通列
# 统一日期列名为 trade_date
if 'date' in all_stock_data.columns and 'trade_date' not in all_stock_data.columns:  # 兼容不同数据源
    all_stock_data = all_stock_data.rename(columns={'date': 'trade_date'})  # 重命名日期列

# 筛选海康威视并按日期排序
hikvision_data = all_stock_data[all_stock_data['order_book_id'] == '002415.XSHE'].sort_values('trade_date').copy()  # 提取海康威视数据

# 构造时间序号与收盘价
time_index = np.arange(len(hikvision_data))  # 生成0到N-1的时间序号作为自变量
close_price = hikvision_data['close'].values  # 提取收盘价数组作为因变量

# 绘制股价走势
plt.figure(figsize=(10, 4))  # 创建10x4英寸画布
plt.plot(time_index, close_price, linewidth=0.8, color='steelblue')  # 绑制收盘价折线
plt.xlabel('交易日序号')  # 设置横轴标签
plt.ylabel('收盘价（元）')  # 设置纵轴标签
plt.title(f'海康威视收盘价（共{len(hikvision_data)}个交易日）')  # 设置标题
plt.tight_layout()  # 自动调整布局
plt.show()  # 显示图表

from sklearn.pipeline import Pipeline  # 导入Pipeline用于构建建模流水线
from sklearn.preprocessing import PolynomialFeatures  # 导入多项式特征变换器
from sklearn.linear_model import LinearRegression  # 导入线性回归模型

degree_candidates = [1, 3, 5, 10, 15, 20]  # 候选多项式阶数列表
fig, axes = plt.subplots(2, 3, figsize=(12, 6))  # 创建2行3列子图网格
axes = axes.ravel()  # 将二维子图数组展平为一维便于迭代

time_features = time_index.reshape(-1, 1)  # 将时间序号变为列向量（sklearn要求二维输入）

for idx, degree in enumerate(degree_candidates):  # 遍历每个候选阶数
    # 构建多项式回归流水线：先做特征变换再拟合线性回归
    poly_pipeline = Pipeline([
        ('poly', PolynomialFeatures(degree=degree)),  # 生成1到d次多项式特征
        ('lr', LinearRegression())  # 对多项式特征做最小二乘拟合
    ])
    poly_pipeline.fit(time_features, close_price)  # 在全量数据上拟合模型
    predicted_price = poly_pipeline.predict(time_features)  # 预测拟合值

    axes[idx].scatter(time_index, close_price, s=1, alpha=0.3, color='gray')  # 绑制原始数据散点
    axes[idx].plot(time_index, predicted_price, color='red', linewidth=2)  # 绑制拟合曲线
    axes[idx].set_title(f'阶数 d={degree}')  # 设置子图标题

plt.tight_layout()  # 自动调整子图间距
plt.show()  # 显示图表

from sklearn.model_selection import TimeSeriesSplit, cross_val_score  # 导入时间序列CV工具

degree_list = list(range(1, 16))  # 候选阶数1到15
cv_mean_mse = []  # 存储各阶数的平均交叉验证MSE

tscv = TimeSeriesSplit(n_splits=5)  # 5折时间序列交叉验证（保持时间顺序）

for degree in degree_list:  # 遍历每个候选阶数
    poly_pipeline = Pipeline([
        ('poly', PolynomialFeatures(degree=degree)),  # 多项式特征变换
        ('lr', LinearRegression())  # 线性回归
    ])
    # scoring='neg_mean_squared_error'：返回负MSE（sklearn约定越大越好）
    cv_scores = cross_val_score(poly_pipeline, time_features, close_price, cv=tscv, scoring='neg_mean_squared_error')
    cv_mean_mse.append(-cv_scores.mean())  # 取负号还原为正MSE并求平均

optimal_degree = degree_list[np.argmin(cv_mean_mse)]  # 取MSE最小的阶数为最优

plt.figure(figsize=(8, 4))  # 创建画布
plt.plot(degree_list, cv_mean_mse, 'o-', color='steelblue')  # 绑制各阶数的CV MSE折线
plt.xlabel('多项式阶数')  # 横轴标签
plt.ylabel('交叉验证 MSE')  # 纵轴标签
plt.title(f'最优阶数 d = {optimal_degree}')  # 标题显示最优阶数
plt.tight_layout()  # 自动调整布局
plt.show()  # 显示图表

# 将时间序号等频分为8个区间（每个区间约含相同数量的观测值）
hikvision_data['Time'] = time_index  # 添加时间序号列
hikvision_data['Time_bin'] = pd.qcut(time_index, q=8, labels=False)  # 等频分箱为0-7共8组

# 为每个分箱生成哑变量并拟合OLS模型
step_dummies = pd.get_dummies(hikvision_data['Time_bin'], prefix='bin', drop_first=True, dtype=float)  # 生成虚拟变量
step_features = sm.add_constant(step_dummies)  # 添加截距项
step_model = sm.OLS(close_price, step_features).fit()  # 拟合阶梯函数OLS模型

plt.figure(figsize=(10, 4))  # 创建画布
plt.scatter(time_index, close_price, s=1, alpha=0.3, color='gray')  # 绑制原始数据
plt.plot(time_index, step_model.fittedvalues, color='red', linewidth=2)  # 绑制阶梯函数拟合值
plt.xlabel('交易日序号')  # 横轴标签
plt.ylabel('收盘价（元）')  # 纵轴标签
plt.title('阶梯函数拟合（8个等频分箱）')  # 图表标题
plt.tight_layout()  # 自动调整布局
plt.show()  # 显示图表

from patsy import dmatrix  # 导入patsy用于根据公式生成设计矩阵

df_candidates = list(range(3, 16))  # 候选自由度3到15
cv_mse_spline = []  # 存储各自由度的交叉验证MSE

for df_val in df_candidates:  # 遍历每个候选自由度
    fold_mse_list = []  # 当前自由度下各折的MSE
    for train_idx, test_idx in tscv.split(time_features):  # 时间序列5折交叉验证
        x_train_fold = time_index[train_idx]  # 当前折的训练集时间序号
        x_test_fold = time_index[test_idx]  # 当前折的测试集时间序号
        y_train_fold = close_price[train_idx]  # 训练集目标值
        y_test_fold = close_price[test_idx]  # 测试集目标值

        # 用patsy cr()生成自然样条基矩阵
        basis_train = dmatrix(f'cr(x, df={df_val})', {'x': x_train_fold}, return_type='dataframe')  # 训练集基函数矩阵
        basis_test = dmatrix(f'cr(x, df={df_val})', {'x': x_test_fold}, return_type='dataframe')  # 测试集基函数矩阵

        spline_fold_model = sm.OLS(y_train_fold, basis_train).fit()  # 拟合当前折的样条回归
        y_pred_fold = spline_fold_model.predict(basis_test)  # 在测试折上预测
        fold_mse_list.append(np.mean((y_test_fold - y_pred_fold) ** 2))  # 计算并记录MSE

    cv_mse_spline.append(np.mean(fold_mse_list))  # 求5折平均MSE

optimal_df = df_candidates[np.argmin(cv_mse_spline)]  # 选取MSE最小的自由度

plt.figure(figsize=(8, 4))  # 创建画布
plt.plot(df_candidates, cv_mse_spline, 'o-', color='steelblue')  # 绑制CV MSE随自由度变化的折线
plt.xlabel('自然样条自由度')  # 横轴标签
plt.ylabel('交叉验证 MSE')  # 纵轴标签
plt.title(f'最优自由度 df = {optimal_df}（CV MSE ≈ {min(cv_mse_spline):.2f}）')  # 显示最优结果
plt.tight_layout()  # 自动调整布局
plt.show()  # 显示图表

# 拟合15阶多项式
poly15 = Pipeline([
    ('poly', PolynomialFeatures(degree=15)),  # 15阶多项式特征
    ('lr', LinearRegression())  # 线性回归
])
poly15.fit(time_features, close_price)  # 拟合模型
pred_poly15 = poly15.predict(time_features)  # 预测

# 拟合df=15的自然样条
basis_full = dmatrix('cr(x, df=15)', {'x': time_index}, return_type='dataframe')  # 生成全量样条基
spline15_model = sm.OLS(close_price, basis_full).fit()  # 拟合样条模型
pred_spline15 = spline15_model.fittedvalues  # 获取拟合值

fig, axes = plt.subplots(1, 2, figsize=(12, 4))  # 创建1行2列子图

axes[0].scatter(time_index, close_price, s=1, alpha=0.3, color='gray')  # 左图：原始数据
axes[0].plot(time_index, pred_poly15, color='red', linewidth=2)  # 15阶多项式拟合线
axes[0].set_title('15阶多项式（边界剧烈振荡）')  # 子图标题
axes[0].set_ylim(-50, 100)  # 限制纵轴范围避免极端值

axes[1].scatter(time_index, close_price, s=1, alpha=0.3, color='gray')  # 右图：原始数据
axes[1].plot(time_index, pred_spline15, color='blue', linewidth=2)  # 自然样条拟合线
axes[1].set_title(f'自然样条 df=15（CV MSE ≈ {min(cv_mse_spline):.2f}）')  # 子图标题

plt.tight_layout()  # 自动调整布局
plt.show()  # 显示图表

from pygam import LinearGAM  # 导入pygam线性GAM模型
from pygam import s as s_gam  # 导入平滑项函数（避免与pygam.s冲突）

if not hasattr(np, 'int'):  # 兼容新版NumPy中移除np.int别名而旧版pygam仍会引用该名称的情况
    np.int = int  # 将np.int临时指向Python内置int以保证pygam内部样条基函数构造正常执行

lambda_values = [0.1, 10, 1000]  # 三个代表性的平滑参数值
fig, axes = plt.subplots(1, 3, figsize=(14, 4))  # 创建1行3列子图

for idx, lam in enumerate(lambda_values):  # 遍历三个λ值
    # s_gam(0, lam=lam)：对第0个特征施加平滑样条，指定λ
    gam_model = LinearGAM(s_gam(0, lam=lam)).fit(time_features, close_price)  # 拟合平滑样条GAM
    x_grid = np.linspace(0, len(time_index) - 1, 500).reshape(-1, 1)  # 生成500个均匀网格点用于绘图
    y_pred_grid = gam_model.predict(x_grid)  # 在网格上预测

    axes[idx].scatter(time_index, close_price, s=1, alpha=0.3, color='gray')  # 绑制原始数据散点
    axes[idx].plot(x_grid, y_pred_grid, color='red', linewidth=2)  # 绑制平滑样条拟合曲线
    axes[idx].set_title(f'λ = {lam}')  # 子图标题

plt.tight_layout()  # 自动调整布局
plt.show()  # 显示图表

# 使用pygam内置的网格搜索功能自动选择最优λ
optimal_gam = LinearGAM(s_gam(0)).gridsearch(  # 对第0个特征施加平滑样条
    time_features, close_price,  # 输入时间特征和目标股价
    lam=np.logspace(-3, 4, 30)  # 在[0.001, 10000]范围内搜索30个候选λ值
)

optimal_lambda = optimal_gam.lam[0][0]  # 提取最优λ值
effective_dof = optimal_gam.statistics_['edof']  # 提取有效自由度

x_grid = np.linspace(0, len(time_index) - 1, 500).reshape(-1, 1)  # 生成绘图用网格
y_pred_optimal = optimal_gam.predict(x_grid)  # 用最优模型在网格上预测

plt.figure(figsize=(10, 4))  # 创建画布
plt.scatter(time_index, close_price, s=1, alpha=0.3, color='gray')  # 原始数据散点
plt.plot(x_grid, y_pred_optimal, color='red', linewidth=2)  # 绑制最优平滑样条
plt.title(f'最优平滑样条（λ ≈ {optimal_lambda:.4f}，有效自由度 ≈ {effective_dof:.2f}）')  # 显示选中参数
plt.xlabel('交易日序号')  # 横轴标签
plt.ylabel('收盘价（元）')  # 纵轴标签
plt.tight_layout()  # 自动调整布局
plt.show()  # 显示图表

from statsmodels.nonparametric.smoothers_lowess import lowess  # 导入LOWESS平滑函数

frac_values = [0.01, 0.05, 0.2]  # 三个候选带宽（邻域比例）
fig, axes = plt.subplots(1, 3, figsize=(14, 4))  # 创建1行3列子图

for idx, frac_val in enumerate(frac_values):  # 遍历每个带宽
    # lowess返回(x排序, y平滑)的二维数组
    smoothed_result = lowess(close_price, time_index, frac=frac_val, return_sorted=True)  # LOWESS平滑

    axes[idx].scatter(time_index, close_price, s=1, alpha=0.3, color='gray')  # 原始数据
    axes[idx].plot(smoothed_result[:, 0], smoothed_result[:, 1], color='red', linewidth=2)  # 平滑曲线
    axes[idx].set_title(f'frac = {frac_val}')  # 设置子图标题

plt.tight_layout()  # 自动调整布局
plt.show()  # 显示图表

fig, axes = plt.subplots(1, 3, figsize=(14, 4))  # 创建1行3列子图

# 偏依赖展示：保持其他变量在训练集分布上平均化，展示单变量的边际效应
PartialDependenceDisplay.from_estimator(
    gam_pipeline,  # 训练好的GAM管道模型
    gam_feature_matrix,  # 训练数据
    features=['Momentum', 'Volatility', 'Month'],  # 要展示偏依赖的三个特征
    ax=axes,  # 绑制到预设的子图轴上
    subsample=50,  # 随机下采样50个样本以加速计算（注：生产环境建议≥200）
    grid_resolution=20  # 网格精度为20个点
)

plt.tight_layout()  # 自动调整布局
plt.show()  # 显示图表