多项式回归的基本思想

多项式回归通过添加高次项来扩展线性模型：

\[ y_i = \beta_0 + \beta_1 x_i + \beta_2 x_i^2 + \cdots + \beta_d x_i^d + \epsilon_i \]

• 本质上仍是线性模型——对参数 \(\beta\) 线性
• 阶数 \(d\) 控制灵活性：\(d\) 越大越灵活，但也越容易过拟合
• Runge 现象：高阶多项式在端点处可能剧烈振荡

多项式回归的数学本质

d阶多项式回归模型：

\[ \large{ y_i = \beta_0 + \beta_1 x_i + \beta_2 x_i^2 + \cdots + \beta_d x_i^d + \varepsilon_i } \tag{3}\]

等价表述：创建新变量 \(x_i^2, x_i^3, \ldots, x_i^d\)，然后做普通多元线性回归。

为什么d通常 ≤ 5？

• 高阶多项式在数据边界处剧烈振荡（龙格现象）
• 多项式系数随阶数增大变得极度不稳定
• 阶数过高 → 严重过拟合

补充说明：共线性问题

\(x, x^2, x^3\) 之间往往高度相关，导致系数估计不稳定。 解决方案：使用正交多项式（Orthogonal Polynomials）：

from numpy.polynomial import polynomial as P
# 正交多项式避免了共线性问题

多项式回归的数学本质

d阶多项式回归模型：

\[ \large{ y_i = \beta_0 + \beta_1 x_i + \beta_2 x_i^2 + \cdots + \beta_d x_i^d + \varepsilon_i } \tag{4}\]

等价表述：创建新变量 \(x_i^2, x_i^3, \ldots, x_i^d\)，然后做普通多元线性回归。

为什么d通常 ≤ 5？

• 高阶多项式在数据边界处剧烈振荡（龙格现象）
• 多项式系数随阶数增大变得极度不稳定
• 阶数过高 → 严重过拟合

补充说明：共线性问题

\(x, x^2, x^3\) 之间往往高度相关，导致系数估计不稳定。 解决方案：使用正交多项式（Orthogonal Polynomials）：

from numpy.polynomial import polynomial as P
# 正交多项式避免了共线性问题

不同阶数多项式拟合对比

Code

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()  # 显示图表

交叉验证选择最优阶数

Code

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()  # 显示图表

ANOVA逐步检验多项式阶数

除了交叉验证，可以用F检验逐步检验是否需要更高阶项：

import statsmodels.api as sm  # 导入statsmodels用于统计检验
from scipy.stats import f as f_dist  # 导入F分布用于计算p值

anova_results = []  # 存储各对嵌套模型的F检验结果
for degree in range(1, 6):  # 比较1阶vs2阶、2阶vs3阶...4阶vs5阶
    # 拟合当前阶数模型
    poly_current = PolynomialFeatures(degree=degree)  # 当前阶数特征变换
    features_current = poly_current.fit_transform(time_features)  # 生成当前阶特征矩阵
    model_current = sm.OLS(close_price, features_current).fit()  # 拟合OLS模型

    # 拟合下一阶数模型
    poly_next = PolynomialFeatures(degree=degree + 1)  # 下一阶特征变换
    features_next = poly_next.fit_transform(time_features)  # 生成下一阶特征矩阵
    model_next = sm.OLS(close_price, features_next).fit()  # 拟合OLS模型

    # 计算F统计量
    rss_current = model_current.ssr  # 当前模型的残差平方和
    rss_next = model_next.ssr  # 下一阶模型的残差平方和
    df_diff = model_next.df_model - model_current.df_model  # 自由度差异
    f_stat = ((rss_current - rss_next) / df_diff) / (rss_next / model_next.df_resid)  # F统计量
    p_val = 1 - f_dist.cdf(f_stat, df_diff, model_next.df_resid)  # 右尾p值

    anova_results.append({'比较': f'{degree}阶 vs {degree+1}阶', 'F统计量': f'{f_stat:.2f}', 'p值': f'{p_val:.4f}'})

pd.DataFrame(anova_results)  # 输出ANOVA比较结果表

ANOVA与F检验的直觉理解

ANOVA的核心问题：增加模型复杂度后，拟合改善是否显著？

\[ \large{ F = \frac{(\text{RSS}_1 - \text{RSS}_2) / (p_2 - p_1)}{\text{RSS}_2 / (n - p_2)} } \tag{5}\]

• \(\text{RSS}_1\)：简单模型（如1阶）的残差平方和
• \(\text{RSS}_2\)：复杂模型（如2阶）的残差平方和
• 分子：增加的复杂度带来多少拟合改善
• 分母：复杂模型的残差方差（基准噪声水平）

决策逻辑：

• \(F\) 值大 → 拟合改善显著 → 值得增加复杂度
• \(F\) 值小 → 改善可能只是偶然 → 保持简单模型
• 逐步检验：1阶 vs 2阶 → 2阶 vs 3阶 → … → 第一次p值不显著时停止

多项式回归的局限性

为什么需要更灵活的方法？

问题一：全局影响

• 修改一处数据 → 影响整条拟合曲线
• 高杠杆点可能导致曲线在远离该点处也发生扭曲

问题二：边界振荡

• 多项式在数据范围边界处容易剧烈振荡
• 龙格现象（Runge’s Phenomenon）：阶数越高，边界振荡越严重

问题三：不适合非光滑关系

• 如果真实关系有”尖角”或”阶梯”，多项式无法精确拟合
• 需要分段（Piecewise）的方法

过渡：这些局限引出了本章后续的方法——阶梯函数和样条，它们用局部拟合替代全局多项式。

ANOVA与F检验的直觉理解

ANOVA的核心问题：增加模型复杂度后，拟合改善是否显著？

\[ \large{ F = \frac{(\text{RSS}_1 - \text{RSS}_2) / (p_2 - p_1)}{\text{RSS}_2 / (n - p_2)} } \tag{6}\]

• \(\text{RSS}_1\)：简单模型（如1阶）的残差平方和
• \(\text{RSS}_2\)：复杂模型（如2阶）的残差平方和
• 分子：增加的复杂度带来多少拟合改善
• 分母：复杂模型的残差方差（基准噪声水平）

决策逻辑：

• \(F\) 值大 → 拟合改善显著 → 值得增加复杂度
• \(F\) 值小 → 改善可能只是偶然 → 保持简单模型
• 逐步检验：1阶 vs 2阶 → 2阶 vs 3阶 → … → 第一次p值不显著时停止

多项式回归的局限性

为什么需要更灵活的方法？

问题一：全局影响

• 修改一处数据 → 影响整条拟合曲线
• 高杠杆点可能导致曲线在远离该点处也发生扭曲

问题二：边界振荡

• 多项式在数据范围边界处容易剧烈振荡
• 龙格现象（Runge’s Phenomenon）：阶数越高，边界振荡越严重

问题三：不适合非光滑关系

• 如果真实关系有”尖角”或”阶梯”，多项式无法精确拟合
• 需要分段（Piecewise）的方法

过渡：这些局限引出了本章后续的方法——阶梯函数和样条，它们用局部拟合替代全局多项式。

阶梯函数：将连续变量离散化

阶梯函数将连续变量 \(X\) 分成 \(K\) 个区间，在每个区间内拟合常数：

\[ y_i = \beta_0 + \beta_1 \cdot \mathbf{1}(c_1 \le x_i < c_2) + \cdots + \beta_K \cdot \mathbf{1}(c_K \le x_i < c_{K+1}) + \epsilon_i \]

• 切点 \(c_1, \ldots, c_K\) 将 \(x\) 的取值范围分为 \(K+1\) 段
• 每段内的预测值是一个常数（该段的均值）
• 优点：简单直观，易于解释
• 缺点：对切点位置敏感，不连续的”阶梯”可能不自然

阶梯函数的数学表述

将连续变量 \(X\) 划分为 \(K+1\) 个区间：

\[ \large{ y_i = \beta_0 + \beta_1 C_1(x_i) + \beta_2 C_2(x_i) + \cdots + \beta_K C_K(x_i) + \varepsilon_i } \tag{7}\]

其中 \(C_k(x) = I(c_k \leq x < c_{k+1})\) 是指示函数。

系数解读：

• \(\beta_0\)：基准区间（\(x < c_1\)）的均值
• \(\beta_k\)：第 \(k\) 个区间相对于基准的平均差异
• 等价于：对分组变量做单因子ANOVA

金融应用：

• 将日成交量分为”低/中/高”三档，分析各档下的收益率特征
• 将市盈率按分位数分组，比较各组的未来表现
• 本质上就是”分组比较”——在量化投资中非常常见

阶梯函数的数学表述

将连续变量 \(X\) 划分为 \(K+1\) 个区间：

\[ \large{ y_i = \beta_0 + \beta_1 C_1(x_i) + \beta_2 C_2(x_i) + \cdots + \beta_K C_K(x_i) + \varepsilon_i } \tag{8}\]

其中 \(C_k(x) = I(c_k \leq x < c_{k+1})\) 是指示函数。

系数解读：

• \(\beta_0\)：基准区间（\(x < c_1\)）的均值
• \(\beta_k\)：第 \(k\) 个区间相对于基准的平均差异
• 等价于：对分组变量做单因子ANOVA

金融应用：

• 将日成交量分为”低/中/高”三档，分析各档下的收益率特征
• 将市盈率按分位数分组，比较各组的未来表现
• 本质上就是”分组比较”——在量化投资中非常常见

阶梯函数拟合与交叉验证

Code

# 将时间序号等频分为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()  # 显示图表

从分段多项式到样条

回归样条在每段区间拟合低阶多项式，并在节点处施加连续性约束：

• 分段多项式：在 \(K\) 个节点处分段，每段独立拟合多项式
• 连续性约束：要求函数值、一阶导数、二阶导数在节点处连续
• 三次样条：每添加一个节点增加 1 个自由度（共 \(K + 4\) 个参数）
• 截断幂基：\(h(x, \xi) = (x - \xi)^3_+\) 是节点 \(\xi\) 处的基函数

节点数量与位置的选择

节点选择策略：

经验法则：

• 通常 3-5 个节点已经足够
• 节点放在数据密度高的区域（分位数方法）
• 通过CV选择最优节点数量

过度节点的风险：

• 节点太多 → 过拟合 → CV误差上升
• 极端情况：\(K = n\)（每个数据点一个节点）→ 完全插值 → 严重过拟合

节点数量与位置的选择

节点选择策略：

经验法则：

• 通常 3-5 个节点已经足够
• 节点放在数据密度高的区域（分位数方法）
• 通过CV选择最优节点数量

过度节点的风险：

• 节点太多 → 过拟合 → CV误差上升
• 极端情况：\(K = n\)（每个数据点一个节点）→ 完全插值 → 严重过拟合

自然样条：边界行为更稳健

自然样条在边界区域（最外侧两个节点之外）强制为线性，有效减少端点振荡：

• 边界约束额外减少 4 个参数 → 自由度 = \(K\)
• 直觉：边界外数据稀少，线性外推比多项式外推更稳健
• Python 实现：patsy.cr() 生成自然样条基函数

自然样条的边界约束

三次样条的边界问题：数据范围之外的外推非常不可靠。

自然样条的额外约束：

在边界节点之外强制函数为线性（而非三次）。

• 等价于：\(f''(\xi_1) = f'''(\xi_1) = 0\) 且 \(f''(\xi_K) = f'''(\xi_K) = 0\)
• 每个边界多2个约束 → 自由度减少4个
• 自然样条自由度 = \(K\)（\(K\) 个节点）

对比：

金融应用的重要性：

• 预测未来股价时，边界外推不可避免
• 自然样条的线性外推比三次外推更保守、更安全

自然样条的边界约束

三次样条的边界问题：数据范围之外的外推非常不可靠。

自然样条的额外约束：

在边界节点之外强制函数为线性（而非三次）。

• 等价于：\(f''(\xi_1) = f'''(\xi_1) = 0\) 且 \(f''(\xi_K) = f'''(\xi_K) = 0\)
• 每个边界多2个约束 → 自由度减少4个
• 自然样条自由度 = \(K\)（\(K\) 个节点）

对比：

金融应用的重要性：

• 预测未来股价时，边界外推不可避免
• 自然样条的线性外推比三次外推更保守、更安全

自然样条拟合与CV选择自由度

Code

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()  # 显示图表

样条 vs 多项式：CV MSE 对比

Code

# 拟合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()  # 显示图表

样条方法的四大优势

平滑样条：拟合与光滑性的权衡

平滑样条通过最小化带惩罚的目标函数来拟合数据：

\[ \min_g \sum_{i=1}^{n} (y_i - g(x_i))^2 + \lambda \int g''(t)^2 \, dt \]

• 第一项：拟合误差（越小越贴合数据）
• 第二项：粗糙度惩罚（\(g''(t)^2\) 衡量曲率，越小越光滑）
• \(\lambda\) 控制两者的权衡：
  • \(\lambda = 0\)：完美插值（过拟合）
  • \(\lambda \to \infty\)：退化为直线（欠拟合）

平滑样条的优化问题

平滑样条求解的是一个带惩罚的最优化问题：

\[ \large{ \min_f \sum_{i=1}^n (y_i - f(x_i))^2 + \lambda \int [f''(t)]^2 dt } \tag{11}\]

两项的权衡：

• 第一项（拟合度）：残差平方和 → 希望尽量小 → 鼓励复杂模型
• 第二项（光滑度）：二阶导数的积分 → 惩罚曲率 → 鼓励简单模型
• \(\lambda \geq 0\)：调节旋钮，平衡两者

\(\lambda\) 的极端情况：

平滑样条的优化问题

平滑样条求解的是一个带惩罚的最优化问题：

\[ \large{ \min_f \sum_{i=1}^n (y_i - f(x_i))^2 + \lambda \int [f''(t)]^2 dt } \tag{12}\]

两项的权衡：

• 第一项（拟合度）：残差平方和 → 希望尽量小 → 鼓励复杂模型
• 第二项（光滑度）：二阶导数的积分 → 惩罚曲率 → 鼓励简单模型
• \(\lambda \geq 0\)：调节旋钮，平衡两者

\(\lambda\) 的极端情况：

有效自由度与LOOCV

• 平滑样条的解是在所有 \(n\) 个数据点处的自然样条，但通过惩罚实现收缩
• 有效自由度（Effective Degrees of Freedom）：\(\text{df}_\lambda = \text{tr}(S_\lambda)\)
  • \(S_\lambda\) 是平滑矩阵，\(\hat{y} = S_\lambda y\)
  • \(\text{df}_\lambda\) 随 \(\lambda\) 增大而减小
• LOOCV 快捷公式（无需反复拟合）：

\[ \text{RSS}_{cv}(\lambda) = \sum_{i=1}^{n} \left(\frac{y_i - \hat{g}_\lambda(x_i)}{1 - \{S_\lambda\}_{ii}}\right)^2 \]

pygam实现平滑样条

Code

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

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()  # 显示图表

网格搜索最优平滑参数

Code

# 使用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()  # 显示图表

平滑样条 vs 回归样条

LOWESS的优势与局限

优势：

• ✓ 完全非参数：不假设任何函数形式
• ✓ 局部自适应：不同区域可以有不同的趋势
• ✓ 对离群值稳健（使用稳健权重迭代时）
• ✓ 是探索性数据分析的利器

局限：

• ✗ 计算量大：每个预测点都需要做一次回归 → \(O(n^2)\)
• ✗ 维度诅咒：高维空间中”邻域”概念失效
• ✗ 无法给出参数估计或系数检验
• ✗ 不易推广到多变量场景

在金融中的典型用途：

• 探索性分析：绘制”非参数趋势线”
• 与参数模型对比：LOWESS曲线偏离线性 → 证据表明非线性
• 不适合正式建模，更多用于数据可视化

LOWESS的优势与局限

优势：

• ✓ 完全非参数：不假设任何函数形式
• ✓ 局部自适应：不同区域可以有不同的趋势
• ✓ 对离群值稳健（使用稳健权重迭代时）
• ✓ 是探索性数据分析的利器

局限：

• ✗ 计算量大：每个预测点都需要做一次回归 → \(O(n^2)\)
• ✗ 维度诅咒：高维空间中”邻域”概念失效
• ✗ 无法给出参数估计或系数检验
• ✗ 不易推广到多变量场景

在金融中的典型用途：

• 探索性分析：绘制”非参数趋势线”
• 与参数模型对比：LOWESS曲线偏离线性 → 证据表明非线性
• 不适合正式建模，更多用于数据可视化

LOWESS实现与带宽选择

Code

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()  # 显示图表

GAM的基本思想

广义加性模型（GAM）将线性模型中的每个线性项替换为灵活的非线性函数：

\[ y_i = \beta_0 + f_1(x_{i1}) + f_2(x_{i2}) + \cdots + f_p(x_{ip}) + \epsilon_i \]

• 每个 \(f_j\) 可以是样条、局部回归等非线性函数
• 可加性假设：各特征的效应独立相加（不含交互项）
• 这是GAM最大的优势也是限制：保持了可解释性，但可能遗漏交互效应

GAM的数学表述

广义可加模型（Generalized Additive Model）：

\[ \large{ y_i = \beta_0 + f_1(x_{i1}) + f_2(x_{i2}) + \cdots + f_p(x_{ip}) + \varepsilon_i } \tag{13}\]

每个 \(f_j\) 是单变量非参数函数（通常用平滑样条实现）。

与普通线性回归的对比：

关键假设——可加性：

• 各变量的效应独立相加，没有交互
• 这是一个强假设，但保证了可解释性
• 如需交互效应：可添加 \(f_{jk}(x_j, x_k)\) 项（二维平滑面）

GAM的数学表述

广义可加模型（Generalized Additive Model）：

\[ \large{ y_i = \beta_0 + f_1(x_{i1}) + f_2(x_{i2}) + \cdots + f_p(x_{ip}) + \varepsilon_i } \tag{14}\]

每个 \(f_j\) 是单变量非参数函数（通常用平滑样条实现）。

与普通线性回归的对比：

关键假设——可加性：

• 各变量的效应独立相加，没有交互
• 这是一个强假设，但保证了可解释性
• 如需交互效应：可添加 \(f_{jk}(x_j, x_k)\) 项（二维平滑面）

GAM的优势与局限

GAM在A股多因子选股中的应用

传统线性多因子模型：

\[ r_i = \alpha + \beta_1 \cdot \text{PE}_i + \beta_2 \cdot \text{ROE}_i + \beta_3 \cdot \text{市值}_i + \varepsilon_i \]

GAM增强版：

\[ r_i = \alpha + f_1(\text{PE}_i) + f_2(\text{ROE}_i) + f_3(\text{市值}_i) + \varepsilon_i \]

为什么GAM更适合？

• PE效应非线性：PE极低可能是价值陷阱，PE适中有价值溢价
• 规模效应非线性：小盘股和超大盘股的收益特征不同
• ROE效应非线性：ROE过高可能不可持续

偏依赖图的威力：

• 可视化每个因子的非线性效应
• 发现数据中的真实模式，而非强加线性假设
• 是因子研究的重要探索性工具

GAM在A股多因子选股中的应用

传统线性多因子模型：

\[ r_i = \alpha + \beta_1 \cdot \text{PE}_i + \beta_2 \cdot \text{ROE}_i + \beta_3 \cdot \text{市值}_i + \varepsilon_i \]

GAM增强版：

\[ r_i = \alpha + f_1(\text{PE}_i) + f_2(\text{ROE}_i) + f_3(\text{市值}_i) + \varepsilon_i \]

为什么GAM更适合？

• PE效应非线性：PE极低可能是价值陷阱，PE适中有价值溢价
• 规模效应非线性：小盘股和超大盘股的收益特征不同
• ROE效应非线性：ROE过高可能不可持续

偏依赖图的威力：

• 可视化每个因子的非线性效应
• 发现数据中的真实模式，而非强加线性假设
• 是因子研究的重要探索性工具

sklearn实现GAM：海康威视多因子模型

构建三因子模型预测海康威视未来5日收益率：

• 动量因子：过去20日累计收益率
• 波动率因子：过去20日收益率标准差
• 月份因子：捕捉A股日历效应

from sklearn.compose import ColumnTransformer  # 导入列变换器用于混合特征处理
from sklearn.preprocessing import SplineTransformer, OneHotEncoder  # 样条变换和独热编码
from sklearn.linear_model import Ridge  # 岭回归模型
from sklearn.pipeline import Pipeline  # 模型管道
from sklearn.inspection import PartialDependenceDisplay  # 偏依赖图工具

# 构造因子特征
hikvision_data['Momentum'] = hikvision_data['close'].pct_change(20)  # 20日动量（累计收益率）
hikvision_data['Volatility'] = hikvision_data['close'].pct_change().rolling(20).std()  # 20日波动率
hikvision_data['Month'] = hikvision_data['trade_date'].dt.month  # 交易月份
hikvision_data['Future_Return'] = hikvision_data['close'].pct_change(5).shift(-5)  # 未来5日收益率（目标变量）

# 删除缺失值（滚动窗口和前瞻产生的NA）
gam_analysis_data = hikvision_data[['Momentum', 'Volatility', 'Month', 'Future_Return']].dropna()  # 提取并清洗数据

GAM模型拟合与评估

# 准备特征和目标
gam_feature_matrix = gam_analysis_data[['Momentum', 'Volatility', 'Month']]  # 三因子特征矩阵
gam_target = gam_analysis_data['Future_Return']  # 目标变量

# 定义混合特征处理器：连续变量用样条，类别变量用独热编码
preprocessor = ColumnTransformer([
    ('spline', SplineTransformer(n_knots=5, degree=3), ['Momentum', 'Volatility']),  # 对动量和波动率做三次样条（5个节点）
    ('onehot', OneHotEncoder(drop='first', sparse_output=False), ['Month'])  # 月份独热编码（去掉第一个避免多重共线性）
])

# 构建GAM管道：预处理 → 岭回归
gam_pipeline = Pipeline([
    ('preprocess', preprocessor),  # 特征预处理步骤
    ('ridge', Ridge(alpha=1.0))  # 带L2正则化的线性回归
])

gam_pipeline.fit(gam_feature_matrix, gam_target)  # 在全量数据上拟合GAM模型
r2_score = gam_pipeline.score(gam_feature_matrix, gam_target)  # 计算R²拟合优度
print(f'GAM模型 R² ≈ {r2_score:.4f}')  # 打印模型评估指标

GAM模型 R² ≈ 0.0290

偏依赖图：可视化非线性效应

Code

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()  # 显示图表

逻辑回归GAM：预测涨跌方向

from pygam import LogisticGAM, s as s_gam, f  # 导入逻辑回归GAM、平滑项和因子项
from sklearn.metrics import accuracy_score  # 导入准确率评估函数

# 构造方向变量：未来5日收益率>0为上涨(1)，否则为下跌(0)
gam_analysis_data = gam_analysis_data.copy()  # 创建副本避免SettingWithCopyWarning
gam_analysis_data['Direction'] = (gam_analysis_data['Future_Return'] > 0).astype(int)  # 构造二分类标签

# 准备特征矩阵（pygam需要numpy数组）
pygam_features = gam_analysis_data[['Momentum', 'Volatility', 'Month']].values  # 提取为numpy数组
direction_target = gam_analysis_data['Direction'].values  # 目标标签数组

# 拟合逻辑回归GAM
# s_gam(0): 对动量施加非线性平滑项
# s_gam(1): 对波动率施加非线性平滑项
# f(2): 对月份施加因子项（离散类别变量）
logistic_gam = LogisticGAM(s_gam(0) + s_gam(1) + f(2))  # 定义模型结构
logistic_gam.gridsearch(pygam_features, direction_target)  # 网格搜索最优平滑参数

# 评估模型
predicted_direction = logistic_gam.predict(pygam_features)  # 预测涨跌方向
accuracy = accuracy_score(direction_target, predicted_direction)  # 计算预测准确率
print(f'逻辑回归GAM准确率 ≈ {accuracy:.4f}')  # 打印准确率指标
print(logistic_gam.summary())  # 输出模型摘要（包含Pseudo R²、各项显著性）

逻辑回归GAM准确率 ≈ 0.5760
LogisticGAM                                                                                               
=============================================== ==========================================================
Distribution:                      BinomialDist Effective DoF:                                     35.6438
Link Function:                        LogitLink Log Likelihood:                                 -2526.1802
Number of Samples:                         3764 AIC:                                              5123.648
                                                AICc:                                            5124.3883
                                                UBRE:                                               3.3688
                                                Scale:                                                 1.0
                                                Pseudo R-Squared:                                   0.0316
==========================================================================================================
Feature Function                  Lambda               Rank         EDoF         P > x        Sig. Code   
================================= ==================== ============ ============ ============ ============
s(0)                              [0.2512]             20           12.7         1.57e-04     ***         
s(1)                              [0.2512]             20           12.1         2.27e-08     ***         
f(2)                              [0.2512]             12           10.9         3.89e-07     ***         
intercept                                              1            0.0          6.59e-03     **          
==========================================================================================================
Significance codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

WARNING: Fitting splines and a linear function to a feature introduces a model identifiability problem
         which can cause p-values to appear significant when they are not.

WARNING: p-values calculated in this manner behave correctly for un-penalized models or models with
         known smoothing parameters, but when smoothing parameters have been estimated, the p-values
         are typically lower than they should be, meaning that the tests reject the null too readily.
None

ANOVA检验：月份变量的最佳形式

月份变量应以哪种形式入模？我们通过三个嵌套模型的偏差分析来判断：

from pygam import l  # 导入线性项函数

# 模型1: 月份作为因子项 f(2)
gam_factor = LogisticGAM(s_gam(0) + s_gam(1) + f(2)).fit(pygam_features, direction_target)  # 拟合因子项模型
# 模型2: 月份作为线性项 l(2)
gam_linear = LogisticGAM(s_gam(0) + s_gam(1) + l(2)).fit(pygam_features, direction_target)  # 拟合线性项模型
# 模型3: 月份作为平滑项 s_gam(2)
gam_smooth = LogisticGAM(s_gam(0) + s_gam(1) + s_gam(2)).fit(pygam_features, direction_target)  # 拟合平滑项模型

# 提取各模型的Deviance和有效自由度
dev_factor = float(np.ravel(gam_factor.statistics_['deviance'])[0])  # 因子项模型偏差
dev_linear = float(np.ravel(gam_linear.statistics_['deviance'])[0])  # 线性项模型偏差
dev_smooth = float(np.ravel(gam_smooth.statistics_['deviance'])[0])  # 平滑项模型偏差

comparison_table = pd.DataFrame({
    '模型': ['f(Month) - 因子项', 'l(Month) - 线性项', 's(Month) - 平滑项'],
    'Deviance': [f'{dev_factor:.2f}', f'{dev_linear:.2f}', f'{dev_smooth:.2f}'],
})
comparison_table  # 输出比较表

分类GAM中的链接函数

从回归GAM到分类GAM：

\[ \large{ \log\frac{P(Y=1|X)}{1-P(Y=1|X)} = \beta_0 + f_1(X_1) + f_2(X_2) + \cdots } \tag{15}\]

链接函数（Link Function）将线性预测值映射到概率：

分类GAM的优势：

• 比逻辑回归更灵活（允许非线性效应）
• 比黑箱模型更可解释（有偏依赖图）
• 适合中等复杂度的预测任务

分类GAM中的链接函数

从回归GAM到分类GAM：

\[ \large{ \log\frac{P(Y=1|X)}{1-P(Y=1|X)} = \beta_0 + f_1(X_1) + f_2(X_2) + \cdots } \tag{16}\]

链接函数（Link Function）将线性预测值映射到概率：

分类GAM的优势：

• 比逻辑回归更灵活（允许非线性效应）
• 比黑箱模型更可解释（有偏依赖图）
• 适合中等复杂度的预测任务

六种方法比较总结

实践建议

• 从简单开始：先尝试低阶多项式或基础样条
• 交叉验证：始终用 CV 选择模型复杂度（阶数、节点数、\(\lambda\)）
• 检查边界：特别关注拟合曲线在数据边界处的行为
• 多变量问题：优先考虑 GAMs，必要时加入交互项
• 与线性基线比较：非线性模型不一定优于线性——用 CV 验证
• 中国市场特殊性：A股日历效应和政策驱动使得非线性建模尤为重要

本章核心公式速查表

选择优先级：自然样条 > GAM > 平滑样条 > 多项式 > 阶梯 > LOWESS（用于建模）

本章核心公式速查表

选择优先级：自然样条 > GAM > 平滑样条 > 多项式 > 阶梯 > LOWESS（用于建模）