7  超越线性关系 (Moving Beyond Linearity)

7.1 本章学习契约

  • 先修:第 3、5、6 章的线性基线、前向验证、管道与正则化;能解释偏差—方差权衡。
  • 目标 O7.1:区分多项式、阶梯、回归样条、平滑样条与局部回归的基函数、支持域和外推行为。
  • 目标 O7.2:复用训练期样条设计信息,比较同一前向切分上的线性与非线性模型误差。
  • 目标 O7.3:正确区分 GCV/UBRE 准则与显式时间验证,并输出候选 \(\lambda\) 及日期边界。
  • 目标 O7.4:对冻结 20 日回撤标签设置等跨度 purge,确保训练标签窗口严格早于下一预测原点。
  • 目标 O7.5:把偏依赖限定为模型预测形状,结合支持密度、基线、AUC/Brier 判断证据边界。

入口检查(先作答再展开):结点数与平滑参数是否是同一概念?

展开入口检查答案与补修路径

不是;前者控制回归样条基的局部切分,后者通过粗糙度惩罚控制平滑样条有效自由度。答错者先复习 小节 7.2 和平滑样条目标函数,再画一张“基函数—惩罚—复杂度”关系图。

低风险检索:闭卷写出自然样条的边界行为,并解释为何 pyGAM gridsearch(objective='GCV') 不是时间交叉验证。

渐隐链:正文给完整的标签终点定义;练习 6 只保留二十日边界断言;项目题要求自行选 LogisticGAM 平滑度且只使用 M05 训练—验证键;练习 8 完全独立比较二维与可加结构。

陌生迁移:将活跃度—未来波动率迁移到长三角上市公司“杠杆率—下一期现金流风险”,重新定义支持域、时间窗口和线性基线。

到目前为止,本书主要关注线性模型。线性模型相对简单,易于描述和实现,在解释和推断方面具有优势。然而,标准线性回归在预测能力方面存在显著局限性。这是因为线性假设几乎总是一种近似,有时是一种较差的近似。在 章节 6 中,我们看到可以通过岭回归、Lasso、主成分回归等技术改进最小二乘法。在那个设置中,改进是通过降低线性模型的复杂性,从而降低估计的方差来实现的。但我们仍然使用线性模型,其改进空间是有限的!

在本章中,我们将放宽线性假设,同时尽可能保持可解释性。我们通过检查线性模型的非常简单的扩展来实现这一点,例如多项式回归和阶梯函数,以及更复杂的方法,如样条、局部回归和广义加性模型。

多项式回归(Polynomial Regression)通过添加额外的预测变量来扩展线性模型,这些预测变量是通过将每个原始预测变量提升到幂次获得的。例如,三次回归使用三个变量 \(X, X^2, X^3\) 作为预测变量。这种方法提供了对数据进行非线性拟合的简单方式。

阶梯函数(Step Functions)将变量的范围切割成 \(K\) 个不同的区域以产生定性变量。其效果是拟合分段常数函数。

回归样条(Regression Splines)比多项式和阶梯函数更灵活,实际上是两者的扩展。它们涉及将 \(X\) 的范围分成 \(K\) 个不同的区域。在每个区域内,对数据拟合多项式函数。然而,这些多项式被约束为在区域边界或节点(knots)处平滑连接。只要区间被分成足够的区域,这可以产生极其灵活的拟合。

平滑样条(Smoothing Splines)与回归样条类似,但产生于稍微不同的情况。平滑样条源于最小化残差平方和准则,并施加平滑性惩罚。

局部回归(Local Regression)与样条类似,但在一个重要方面有所不同。区域允许重叠,并且实际上以非常平滑的方式重叠。

广义加性模型(Generalized Additive Models, GAMs)允许我们将上述方法扩展以处理多个预测变量。

小节 7.2小节 7.7 中,我们提出了多种灵活建模响应变量 \(Y\) 与单个预测变量 \(X\) 之间关系的方法。在 小节 7.8 中,我们将展示这些方法可以无缝集成,以将 \(Y\) 建模为多个预测变量 \(X_1, \ldots, X_p\) 的函数。

7.2 多项式回归 (Polynomial Regression)

从历史上看,将线性回归扩展到预测变量与响应变量之间关系为非线性设置的标准方法是将标准线性模型

\[ y_i = \beta_0 + \beta_1 x_i + \epsilon_i \]

替换为多项式函数

\[ y_i = \beta_0 + \beta_1 x_i + \beta_2 x_i^2 + \beta_3 x_i^3 + \cdots + \beta_d x_i^d + \epsilon_i \tag{7.1}\]

其中 \(\epsilon_i\) 是误差项。这种方法称为多项式回归。事实上,我们在 章节 3 的多项式回归小节中已经看到了这种方法的例子。对于足够大的次数 \(d\),多项式回归允许我们产生高度非线性的曲线。请注意,式 7.1 中的系数可以很容易地使用最小二乘线性回归来估计,因为这只是一个标准的线性模型,预测变量为 \(x_i, x_i^2, x_i^3, \ldots, x_i^d\)。一般来说,使用大于3或4的 \(d\) 是不常见的,因为对于大的 \(d\) 值,多项式曲线会变得过于灵活,可能会呈现一些非常奇怪的形状。这在 \(X\) 变量的边界附近尤其如此。

Note: 多项式回归的数学原理

多项式回归本质上仍然是线性回归,因为它关于参数是线性的。我们将原始变量 \(X\) 转换为多个特征 \(X, X^2, X^3, \ldots, X^d\),然后对这些特征拟合线性模型。这意味着我们可以使用所有标准线性回归的工具和推断方法。

多项式回归的关键选择是次数 \(d\)

  • \(d=1\): 线性回归
  • \(d=2\): 二次回归,可以捕捉U形或倒U形关系
  • \(d=3\): 三次回归,可以捕捉更复杂的单峰关系
  • \(d \geq 4\): 高次多项式,可能过拟合,特别是在边界处

高次多项式在边界处的行为往往不稳定,这是龙格现象(Runge’s phenomenon)的表现之一。

7.2.1 中国案例:海康威视股价的长期趋势分析

非线性模型在金融时间序列分析中有着广泛的应用。虽然股票收益率通常被假设为难以预测(有效市场假说),但在长期视角下,股价往往表现出明显的非线性趋势周期性波动

让我们使用海康威视 (002415.XSHE) 的历史股价数据,探索如何使用超越线性的方法来提取其长期价格趋势。我们将尝试使用多项式回归、样条和GAM来拟合股价随时间的变化。

下面用海康威视(002415.XSHE)的前复权收盘价演示时间趋势拟合。前复权口径便于跨公司行为日期保持价格序列连续,但它不是未经调整的成交价格;复权规则、样本终点和价格非平稳性都会改变曲线。这里的时间索引只用于描述当前样本形状,不把趋势拟合解释为收益可预测性或市场机制。

列表 7.1: 导入必要的库与数据
# ── 导入科学计算与可视化核心库 ──
import numpy as np              # 数值计算(多项式特征矩阵、网格生成等)
import pandas as pd             # 表格数据处理(读取 h5、列筛选、缺失值处理)
import matplotlib.pyplot as plt # 绑定 Matplotlib 绘图接口
import statsmodels.api as sm    # 统计建模(OLS、GLM、置信区间等)
from sklearn.preprocessing import PolynomialFeatures  # 将 X 升维为 [1, X, X², ..., X^d]
from sklearn.linear_model import LinearRegression      # 普通最小二乘线性回归器
from sklearn.pipeline import Pipeline                  # 将"特征变换 → 模型拟合"打包成统一管道
from sklearn.model_selection import cross_val_score    # 自动化交叉验证并返回各折评分
import warnings  # 关闭冗余警告,保持输出整洁
warnings.filterwarnings('ignore')  # 关闭冗余警告,保持输出简洁
if not hasattr(np, 'int'):  # 旧版 pygam 仍调用 np.int,这里为兼容 NumPy 1.20+ 补回别名
    np.int = int  # 将已废弃的 np.int 映射到内置 int,避免 pygam 在构造样条基时崩溃

# 设置中文字体(思源黑体),使图表能正确显示中文标题和标签
plt.rcParams['font.family'] = ['Source Han Serif SC']  # 统一使用出版字体思源宋体
plt.rcParams['axes.unicode_minus'] = False  # 修复负号显示为方块的问题

接下来,我们加载并准备海康威视的历史行情数据。

# ── 1. 加载海康威视(002415.XSHE)历史日线行情 ──
import os  # 导入 os 用于跨平台路径拼接
import platform  # 导入 platform 用于严格按操作系统选择数据根目录
# 根据操作系统设置数据根目录路径
DATA_ROOT = os.path.abspath(os.path.expanduser(os.environ['BOOK_DATA_DIR']))  # 从必需环境变量解析数据根目录
assert os.path.isdir(DATA_ROOT), f'BOOK_DATA_DIR 不存在: {DATA_ROOT}'  # 缺少挂载时立即失败
path = os.path.join(DATA_ROOT, 'stock/stock_price_pre_adjusted.h5')  # 拼接前复权股价文件路径
market_stock_data = pd.read_hdf(path)  # 读取前复权日线数据(含全部 A 股)

# 若 h5 以 MultiIndex 存储,则先展平为普通列
if 'order_book_id' in market_stock_data.index.names:  # 获取索引
    market_stock_data = market_stock_data.reset_index()  # 将索引列还原为普通列
# 统一日期列名为 trade_date,兼容不同数据源的命名差异
# 获取列名列表
if 'date' in market_stock_data.columns and 'trade_date' not in market_stock_data.columns:
    market_stock_data = market_stock_data.rename(columns={'date': 'trade_date'})  # 重命名日期列
market_stock_data['trade_date'] = pd.to_datetime(market_stock_data['trade_date'])  # 将交易日期统一转换为 datetime 以确保日期轴正确显示

# 筛选海康威视,并按交易日升序排列
haikang_data = market_stock_data[market_stock_data['order_book_id'] == '002415.XSHE'].copy()  # 过滤海康威视数据
haikang_data = haikang_data.sort_values('trade_date').reset_index(drop=True)  # 按日期排序并重置索引

# 构造回归所需的两列:
#   Time — 等距整数索引 (0, 1, 2, ...),即"第几个交易日"
#   Price — 当日收盘价(前复权),作为被预测的响应变量 Y
haikang_data['Time'] = np.arange(len(haikang_data))  # 生成从 0 开始的等距时间索引
haikang_data['Price'] = haikang_data['close']  # 将收盘价列赋值给 Price 列

# 打印数据概况,帮助确认加载是否正常
print(f"数据范围: {haikang_data['trade_date'].min()}{haikang_data['trade_date'].max()}")  # 打印起止日期
print(f"样本数量: {len(haikang_data)}")  # 打印总交易日数
数据范围: 2010-05-28 00:00:00 至 2025-12-31 00:00:00
样本数量: 3789

运行时应从上表读取实际起止日和观测数,因为本地数据快照可能更新。该样本只代表可用交易区间,不能据跨度本身断言已覆盖所有市场状态。

在正式建模之前,我们先对数据的基本统计特征进行初步审视。下面的描述性统计摘要展示了时间索引(Time)与收盘价(Price)的均值、标准差及分位数分布,有助于我们把握数据的集中趋势和离散程度。

表 7.1: 海康威视股价数据摘要
# 打印 Time(交易日索引)和 Price(收盘价)的描述统计:
# 包含均值、标准差、四分位数等,帮助了解价格分布
print(haikang_data[['Time', 'Price']].describe())  # 生成描述性统计摘要
              Time        Price
count  3789.000000  3789.000000
mean   1894.000000    21.505153
std    1093.934413    13.992262
min       0.000000     2.771200
25%     947.000000     7.247600
50%    1894.000000    24.863300
75%    2841.000000    31.207300
max    3788.000000    61.037800

描述性统计量应以当次表格为准。均值、标准差、分位数与极差描述这份价格水平样本,但价格非平稳且受复权口径影响;均值高于中位数至多提示样本右尾较长,不能单凭两者确定总体分布形状。

现在让我们尝试拟合不同次数的多项式回归模型,看看它们通过时间预测股价的效果:

在完成了基础的时间序列构建后,这段代码正式开始了一场“用多项式堆砌复杂性”的视觉实验。我们利用 sklearn.pipeline 优雅地将特征升维器(PolynomialFeatures)和基础线性回归器(LinearRegression)打包成了一个组合管道。紧接着,代码写了一个 for 循环,强行让拟合多项式的最高次幂(\(d\))从非常朴素的 1 次(纯线性直线)、3 次、5 次,一路狂飙突进到惊悚的 10 次、15 次甚至是 20 次。在这张包含六个子图的矩阵中,你可以非常直观地感受到什么叫做“非线性灵活性带来的失控”:当 \(d=1\) 时,这只是一条呆板且无力的上升红线;当 \(d=5\) 时,曲线开始随着海康威视股价中期的巨大波动展现出了优美的 S 型跟随;但是当你把目光移向 \(d=15\) 甚至 \(d=20\) 时,虽然图表标题上的 MSE(训练误差)被压榨到了一个极低的值,但那条红色的预测线在时间序列的两端(尤其是最右侧的未来端)却如同疯魔般向天空或深渊剧烈发散、扭曲拉扯。这就是统计学上臭名昭著的龙格现象(Runge’s phenomenon),它极其生动且惨痛地向我们宣告了单纯堆砌全局多项式阶数在预测长尾趋势时的无力与危险。

列表 7.2: 拟合不同次数的多项式回归模型
# ── 准备回归输入 ──
time_features = haikang_data[['Time']].values  # X: 交易日索引,shape=(n,1)
price_target = haikang_data['Price'].values     # Y: 当天收盘价,shape=(n,)

# 创建等间距时间网格,用于绘制平滑的拟合曲线(共 500 个点)
time_grid = np.linspace(time_features.min(), time_features.max(), 500).reshape(-1, 1)  # 取最大值
trade_date_numeric = haikang_data['trade_date'].astype('int64').to_numpy()  # 将交易日期转换为整数时间戳,供插值生成平滑日期网格
date_grid = pd.to_datetime(np.interp(time_grid.flatten(), time_features.flatten(), trade_date_numeric))  # 将数值时间网格映射回真实交易日期,避免图形横轴退回 1970 年

# ── 尝试 6 种不同次数的多项式拟合 ──
polynomial_degrees = [1, 3, 5, 10, 15, 20]  # d=1 是纯线性;d=20 是极端高次
fitted_models = []  # 初始化列表,用于存储各次多项式的已训练模型

接下来,图 7.1 可视化不同次数多项式的拟合效果,直观对比模型复杂度的影响。

plt.figure(figsize=(15, 10))  # 创建 2×3 子图矩阵

# 遍历每种多项式次数,依次拟合并绘制子图
for i, degree in enumerate(polynomial_degrees, 1):  # 遍历每个多项式阶数
    # Pipeline 将 "升维 → 回归" 两步封装:
    #   PolynomialFeatures(d) 把 X 变成 [1, X, X², ..., X^d]
    #   LinearRegression() 对升维后的矩阵做最小二乘
    current_model = Pipeline([  # 构建current_model模型Pipeline
        ('poly', PolynomialFeatures(degree=degree)),  # 生成多项式特征
        ('linear', LinearRegression())  # 初始化线性回归模型
    ])  # 完成构建

    current_model.fit(time_features, price_target)  # 在全部训练数据上拟合
    fitted_models.append(current_model)  # 将当前模型存入列表,供后续 ANOVA 等分析复用

    # 计算训练集 MSE(均方误差)——仅反映拟合精度,不代表泛化能力
    predicted_price_train = current_model.predict(time_features)  # 使用模型进行预测
    current_mse = np.mean((price_target - predicted_price_train) ** 2)  # 计算当前模型在训练集上的均方误差

    # ── 绘制第 i 个子图 ──
    plt.subplot(2, 3, i)  # 设置当前子图位置
    # 黑色半透明线 = 原始股价走势
    # 绑制折线图
    plt.plot(haikang_data['trade_date'], price_target, 'k-', alpha=0.3, linewidth=1, label='实际股价')
    
    # 红色线 = 多项式在时间网格上的预测
    predicted_price_grid = current_model.predict(time_grid)  # 在连续网格上生成预测曲线
    plt.plot(date_grid, predicted_price_grid, 'r-', linewidth=2, label=f'd={degree}')  # 红色曲线:在真实交易日期网格上的多项式拟合线
    
    plt.title(f'Degree = {degree}\nMSE = {current_mse:.0f}', fontsize=12)  # 子图标题:显示次数与 MSE
    plt.legend()  # 添加图例
    plt.grid(True, alpha=0.3)  # 添加半透明网格线
    if i > 3:  # 仅底部一行子图显示 x 轴标签,减少视觉冗余
        plt.xlabel('交易日期')  # 仅底部行显示真实交易日期标签

plt.tight_layout()  # 自动调整子图间距,防止标签重叠
plt.show()  # 渲染并显示图表
六幅子图按一、三、五、十、十五和二十次展示海康威视股价与全样本多项式拟合曲线,标题报告训练MSE。
图 7.1: 同一股价时间样本上六种全局多项式的训练拟合

图 7.1 只报告同一样本上的训练拟合。若高次曲线在边界附近明显摆动,可把它作为数值与外推风险的诊断线索;是否过拟合仍须由 图 7.2 的前向验证损失确认,不能由训练 MSE 或曲线外观单独判定。

# ── 交叉验证选择最优多项式次数 ──
cross_validation_scores = []  # 初始化交叉验证分数存储列表
cv_polynomial_degrees = list(range(1, 11))  # 候选次数 1~10

# 金融时间序列不能随机打乱!必须用 TimeSeriesSplit:
# 前 N 天训练 → 第 N+1 段验证,逐步推进(模拟真实投资决策)
from sklearn.model_selection import TimeSeriesSplit  # 时间序列专用交叉验证(前N天训练→后续验证)

time_series_cv = TimeSeriesSplit(n_splits=5)  # 5 折时间序列切分

# 对每个候选多项式次数执行交叉验证
for degree in cv_polynomial_degrees:  # 遍历循环
    # 构建多项式回归 Pipeline(特征升维 + 线性回归)
    current_model = Pipeline([  # 构建current_model模型Pipeline
        ('poly', PolynomialFeatures(degree=degree)),  # 生成多项式特征
        ('linear', LinearRegression())  # 初始化线性回归模型
    ])  # 完成构建
    
    # cross_val_score 返回负 MSE(sklearn 约定"越大越好"),取负号还原为正 MSE
    # 执行交叉验证评估
    fold_mse_scores = -cross_val_score(current_model, time_features, price_target, cv=time_series_cv, scoring='neg_mean_squared_error')
    cross_validation_scores.append(fold_mse_scores.mean())  # 5 折平均 MSE

# ── 绘制 CV 误差曲线 ──
plt.figure(figsize=(10, 6))  # 创建画布
plt.plot(cv_polynomial_degrees, cross_validation_scores, 'bo-', linewidth=2)  # 蓝色圆点折线:各次数的 CV MSE
plt.xlabel('多项式次数', fontsize=14)  # x 轴标签
plt.ylabel('CV MSE (TS-Split)', fontsize=14)  # y 轴标签
plt.title('多项式次数选择 (时间序列交叉验证)', fontsize=16)  # 图表标题
plt.grid(True, alpha=0.3)  # 半透明网格线
plt.show()  # 渲染并显示图表

# 选出使 CV MSE 最小的那个次数
optimal_degree = cv_polynomial_degrees[np.argmin(cross_validation_scores)]  # 获取最小值的索引
print(f'最优多项式次数: {optimal_degree}')  # 输出最优多项式次数
折线展示一至十次多项式的前向验证均方误差,最低点对应当次验证选中的次数。
图 7.2: 不同多项式次数的交叉验证误差
最优多项式次数: 1

候选次数及其前向验证误差必须从当次输出读取;被选次数不是预先固定的模型排名。若低次数胜出,只能说明它在这一数据快照、候选集和验证窗口中的误差较小;若高次数胜出,也仍须检查边界振荡与锁定测试期。该比较不能单凭一次排名“证明”龙格现象或金融序列普遍偏好线性趋势。

Tip: 关于多项式拟合金融时间序列的龙格现象

观察上图中的高次多项式(如 \(d=15, 20\)),你会发现一个严重的问题:尽管曲线在样本中间部分可能拟合得不错,但在两端(早期和近期),曲线表现出剧烈的振荡

这就是著名的龙格现象(Runge’s phenomenon)。这告诉我们,试图通过增加多项式次数来捕捉长期的复杂趋势通常是一个坏主意。不仅解释性差,而且在边界处的预测极不稳定(这一点对于预测未来股价尤为致命)。

这正是我们需要引入样条(Splines)的动机。

7.2.2 方差分析

对于多项式回归,我们可以使用方差分析(ANOVA)来检验是否需要更高次的项。

列表 7.3: 多项式回归的方差分析 (Haikang)
from scipy.stats import f  # F 分布 CDF,用于计算 ANOVA 的 p 值

# ── 逐次拟合多项式(1 次 → 5 次),记录残差平方和 RSS 与自由度 ──
max_polynomial_degree = 5  # 设置多项式回归的阶数
residual_sum_squares = []       # 各次数模型的 RSS
degrees_of_freedom_values = []  # 各模型的残差自由度 (n - 参数个数)

# 从 1 次到最高次逐步拟合多项式模型
for degree in range(1, max_polynomial_degree + 1):  # 遍历循环
    # PolynomialFeatures 生成设计矩阵 [1, X, X², ..., X^d]
    # 拟合并转换数据
    standardized_time = (time_features - np.mean(time_features)) / np.std(time_features)
    polynomial_features = PolynomialFeatures(degree=degree).fit_transform(standardized_time)
    ols_model = sm.OLS(price_target, polynomial_features).fit()  # statsmodels OLS 拟合
    current_rss = np.sum(ols_model.resid ** 2)  # 残差平方和 = Σ(y - ŷ)²
    residual_sum_squares.append(current_rss)  # 将当前模型 RSS 存入列表
    degrees_of_freedom_values.append(ols_model.df_resid)  # 残差自由度 = n - (d+1)

# ── F 检验:逐级对比"低次 vs 高一次"是否显著改善拟合 ──
print('多项式次数的方差分析:')  # 打印表头标题
print('-' * 80)  # 打印分隔线
print(f'{"模型":<20} {"残差自由度":<15} {"RSS (x1e6)":<15} {"F统计量":<15} {"p值":<15}')  # 打印列名
print('-' * 80)  # 打印分隔线

# 逐级比较相邻两个模型的拟合改善程度
for i in range(1, len(residual_sum_squares)):  # 计算元素数量
    # 自由度之差 = 新增参数个数(高次模型比低次模型多了几个参数)
    dof_difference = degrees_of_freedom_values[i-1] - degrees_of_freedom_values[i]  # 计算dof_difference的有效自由度
    # RSS 之差 = 低次 RSS - 高次 RSS(高次模型拟合更好,RSS 更小)
    rss_difference = residual_sum_squares[i-1] - residual_sum_squares[i]  # 高次多项式的 RSS 应更小
    
    # F = (RSS 改善量 / 新增自由度) / (当前模型 RSS / 当前残差自由度)
    # 完成模型/Pipeline构建
    f_statistic = max(0.0, (rss_difference / dof_difference) / (residual_sum_squares[i] / degrees_of_freedom_values[i]))
    # p 值:在 F 分布下,观测到的 F 统计量以上的尾部面积
    # 完成模型/Pipeline构建
    anova_p_value = f.sf(f_statistic, dof_difference, degrees_of_freedom_values[i])

    print(f'{i}次 vs {i+1}{"":<12} {dof_difference:<15} {rss_difference/1e6:<15.2f} {f_statistic:<15.4f} {anova_p_value:<15.6f}')  # 输出每对模型的 F 检验结果
多项式次数的方差分析:
--------------------------------------------------------------------------------
模型                   残差自由度           RSS (x1e6)      F统计量            p值             
--------------------------------------------------------------------------------
1次 vs 2次             1.0             0.04            825.1820        0.000000       
2次 vs 3次             1.0             0.07            2344.0388       0.000000       
3次 vs 4次             1.0             0.00            6.5605          0.010465       
4次 vs 5次             1.0             0.01            461.9279        0.000000       

表中使用标准化时间和 f.sf 计算有限的尾部概率。每一行只检验在经典同方差、独立误差假设下新增一阶多项式项是否改善样本内拟合;应从当次输出读取统计量与 \(p\) 值。股价时间序列通常不满足这些经典误差假设,因此该表是嵌套模型诊断,不是未来预测证据;复杂度仍须用后文的前向验证选择。

# ── 用 5 次多项式拟合海康威视股价,并绘制 95% 置信区间 ──
degree = 5  # 选择的多项式次数

# 构建 sklearn Pipeline:先生成多项式特征,再做线性回归
poly_pipeline = Pipeline([  # 构建多项式回归Pipeline(标准化+多项式特征+线性回归)
    ('poly', PolynomialFeatures(degree=degree)),   # 生成 [1, X, X², ..., X⁵]
    ('linear', LinearRegression())                  # 在多项式特征上做 OLS
])  # 完成构建
poly_pipeline.fit(time_features, price_target)      # 训练模型

# ── 转用 statsmodels OLS 以获取统计推断(置信区间) ──
polynomial_features = PolynomialFeatures(degree=degree).fit_transform(time_features)  # 生成多项式特征矩阵
statsmodels_poly_model = sm.OLS(price_target, polynomial_features).fit()  # 用 statsmodels OLS 拟合以获取统计推断

# 在连续时间网格上做预测,使曲线平滑
grid_poly_features = PolynomialFeatures(degree=degree).fit_transform(time_grid)  # 拟合并转换数据
model_predictions = statsmodels_poly_model.get_prediction(grid_poly_features)  # 获取预测及统计信息
predicted_mean_price = model_predictions.predicted_mean                         # 预测均值 ŷ
predicted_confidence_interval = model_predictions.conf_int(alpha=0.05)          # 95% 置信区间

# ── 绘图:观测值 + 拟合曲线 + 置信带 ──
plt.figure(figsize=(12, 8))  # 创建新图形
plt.plot(haikang_data['trade_date'], price_target, 'k-', alpha=0.3, linewidth=1, label='观测股价')  # 黑色半透明线:原始股价走势
plt.plot(date_grid, predicted_mean_price, 'b-', linewidth=2, label=f'{degree}次多项式拟合')  # 蓝色线:在真实交易日期网格上的多项式拟合曲线
# fill_between 用半透明蓝色填充上下界之间的区域,形成置信带
# 填充区域
plt.fill_between(date_grid, predicted_confidence_interval[:, 0], predicted_confidence_interval[:, 1],
                 alpha=0.2, color='blue', label='95%置信区间')  # 定义alpha变量
plt.xlabel('时间', fontsize=14)  # x 轴标签
plt.ylabel('股价 (元)', fontsize=14)  # y 轴标签
plt.title('海康威视股价的多项式回归拟合', fontsize=16)  # 图表标题
plt.legend(fontsize=12)  # 添加图例
plt.grid(True, alpha=0.3)  # 半透明网格线
plt.show()  # 渲染并显示图表
海康威视价格散点、五次多项式拟合曲线及其置信带随交易日期变化。
图 7.3: 带置信区间的多项式回归拟合 (Degree=5)

图 7.3 的蓝色曲线展示5次多项式拟合,浅蓝色阴影是模型设定下的95%点态均值置信区间。边界处区间若变宽,提示该处估计不确定性较大;由于股价误差存在时序相关且模型使用全样本,该区间不应当作未来价格预测区间。

7.2.3 逻辑回归的多项式扩展

多项式回归不仅适用于连续响应变量,也可以扩展到分类问题。为了让分类阈值与海康威视的真实价格分布相匹配,下面我们将’高价股’定义为前复权收盘价达到样本75%分位数的交易日,并使用多项式逻辑回归来刻画该高价状态随时间变化的概率:

列表 7.4: 多项式逻辑回归
# ── 构造二元响应变量:股价达到样本 75% 分位数记为"高价"(1),否则为"低价"(0) ──
high_price_threshold = float(np.quantile(price_target, 0.75))  # 使用样本 75% 分位数定义高价阈值,避免标签完全失衡
is_high_price = (price_target >= high_price_threshold).astype(int)  # 将高于或等于阈值的交易日标记为高价状态

# ── 用 3 次多项式逻辑回归 (GLM-Binomial) 拟合高价概率 ──
degree = 3  # 设置多项式回归的阶数
polynomial_features = PolynomialFeatures(degree=degree).fit_transform(time_features)  # 生成 3 次多项式特征矩阵
# GLM + Binomial 族 = 逻辑回归,链接函数为 logit
logistic_regression_model = sm.GLM(is_high_price, polynomial_features,  # 拟合广义线性模型(Logistic回归)
                    family=sm.families.Binomial()).fit()  # 训练/拟合模型

# 在连续时间网格上预测高价概率
grid_poly_features = PolynomialFeatures(degree=degree).fit_transform(time_grid)  # 拟合并转换数据
predicted_probabilities = logistic_regression_model.predict(grid_poly_features)  # 输出在 [0,1] 之间

# 获取预测概率的 95% 置信区间
prediction_summary = logistic_regression_model.get_prediction(grid_poly_features).summary_frame()  # 在响应概率尺度上提取预测均值及置信区间
predicted_probabilities = prediction_summary['mean'].to_numpy()  # 读取预测概率均值序列
predicted_confidence_lower = prediction_summary['mean_ci_lower'].clip(lower=0).to_numpy()  # 获取并截断概率下界,确保不低于 0
predicted_confidence_upper = prediction_summary['mean_ci_upper'].clip(upper=1).to_numpy()  # 获取并截断概率上界,确保不高于 1

图 7.4 将二元状态事件带、概率曲线与点态区间放在同一坐标系中。

# ── 绘图:低价日 / 高价日标记 + 概率曲线 + 置信带 ──
plt.figure(figsize=(12, 8))  # 创建新图形

# 将低价日与高价日分别映射到真实交易日期,便于与概率曲线共用同一日期横轴
low_price_dates = haikang_data.loc[is_high_price == 0, 'trade_date']  # 提取低价状态对应的交易日期
high_price_dates = haikang_data.loc[is_high_price == 1, 'trade_date']  # 提取高价状态对应的交易日期
low_price_offsets = 0.01 + 0.012 * (np.arange(len(low_price_dates)) % 4)  # 通过分层错位将低价日分散到 y=0 附近的小带状区域
high_price_offsets = 0.99 - 0.012 * (np.arange(len(high_price_dates)) % 4)  # 通过分层错位将高价日分散到 y=1 附近的小带状区域
plt.scatter(low_price_dates, low_price_offsets, alpha=0.25, s=18, color='gray', marker='|', linewidths=0.8, label='低价交易日')  # 用分层错位后的短竖线标记低价日
plt.scatter(high_price_dates, high_price_offsets, alpha=0.25, s=18, color='gray', marker='|', linewidths=0.8, label='高价交易日')  # 用分层错位后的短竖线标记高价日

# 绘制预测概率曲线及其置信带
plt.plot(date_grid, predicted_probabilities, 'b-', linewidth=3, label='预测概率')  # 在真实交易日期网格上绘制预测概率曲线
# 填充区域
plt.fill_between(date_grid, predicted_confidence_lower, predicted_confidence_upper,
                 alpha=0.2, color='blue', label='95%置信区间')  # 在概率尺度上填充 95% 置信带

plt.xlabel('交易日期', fontsize=14)  # x 轴标签
plt.ylabel(f'P(股价 ≥ {high_price_threshold:.2f}元)', fontsize=14)  # y 轴标签
plt.title('时间对高价概率的影响(多项式逻辑回归)', fontsize=16)  # 图表标题
plt.ylim([0, 1])    # 概率范围 [0, 1]
plt.legend(fontsize=12)  # 添加图例
plt.grid(True, alpha=0.3)  # 半透明网格线
plt.show()  # 渲染并显示图表
低价与高价交易日事件带叠加预测概率曲线和百分之九十五置信区间。
图 7.4: 多项式逻辑回归拟合的高价状态概率及置信带

这里我们不再使用脱离样本分布的固定阈值,而是将高价状态定义为’进入样本价格上四分位区间’。这样构造后,逻辑回归能够基于真实样本识别出高价状态在时间上的集中阶段:蓝色曲线给出进入高价区间的估计概率,浅蓝色阴影表示该概率的95%置信区间。灰色短竖线通过分层错位后形成了上下两个清晰的事件带,更便于观察高价日与低价日在时间轴上的聚集分布。在分类建模中,阈值必须贴合样本分布,否则图形诊断和模型解释都会被失真。

7.3 阶梯函数 (Step Functions)

使用特征的多项式函数作为线性模型中的预测变量会对 \(X\) 的非线性函数施加全局结构。我们可以使用阶梯函数来避免施加这种全局结构。这里我们将 \(X\) 的范围分成箱子(bins),并在每个箱子中拟合不同的常数。这相当于将连续变量转换为有序分类变量。

更详细地说,我们在 \(X\) 的范围内创建切点 \(c_1, c_2, \ldots, c_K\),然后构造 \(K+1\) 个新变量:

\[ \begin{aligned} C_0(X) &= I(X < c_1), \\ C_1(X) &= I(c_1 \leq X < c_2), \\ C_2(X) &= I(c_2 \leq X < c_3), \\ &\vdots \\ C_{K-1}(X) &= I(c_{K-1} \leq X < c_K), \\ C_K(X) &= I(c_K \leq X), \end{aligned} \]

其中 \(I(\cdot)\)指示函数(indicator function),如果条件为真则返回1,否则返回0。这些有时称为虚拟变量(dummy variables)。请注意,对于 \(X\) 的任何值,\(C_0(X) + C_1(X) + \cdots + C_K(X) = 1\),因为 \(X\) 必须恰好位于 \(K+1\) 个区间之一。

然后我们使用最小二乘法拟合使用 \(C_1(X), C_2(X), \ldots, C_K(X)\) 作为预测变量的线性模型:

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

对于给定的 \(X\) 值,\(C_1, C_2, \ldots, C_K\) 中最多有一个可以非零。请注意,当 \(X < c_1\) 时,式 7.2 中的所有预测变量都为零,因此 \(\beta_0\) 可以解释为 \(X < c_1\)\(Y\) 的平均值。相比之下,式 7.2 预测对于 \(c_j \leq X < c_{j+1}\),响应为 \(\beta_0 + \beta_j\),因此 \(\beta_j\) 表示相对于 \(X < c_1\)\(X\)\(c_j \leq X < c_{j+1}\) 时的平均增加。

Caution: 阶梯函数的局限性

阶梯函数方法有一些重要限制:

  1. 离散化损失信息:将连续变量离散化会损失信息,特别是在箱子数量较少时
  2. 边界处的不连续:拟合函数在切点处是不连续的,这在许多应用中是不现实的
  3. 箱内常数假设:在每个箱子内假设响应为常数,这可能过于简单

何时使用阶梯函数

  • 预测变量与响应变量之间的关系确实是分段常数
  • 数据中存在自然的断点(如不同政策阶段的分界点)
  • 需要简单易解释的模型
  • 作为其他非线性方法的基线比较

何时避免使用阶梯函数

  • 变量与响应之间的关系平滑变化
  • 需要在边界处进行预测
  • 数据量较少且箱子数量较多时

7.3.1 阶梯函数的Python实现

让我们使用阶梯函数来分析股价随时间的变化(例如,按时间分段):

既然单纯依靠多项式的全局曲线很难驯服时间序列的桀骜不驯,下面的 Python 代码转向了人类最古老也最易懂的分段策略——“阶梯函数”(Step Functions)。在这段代码中,我们将海康威视这段连续的时间历史,通过 pandas.qcut 量化切割成了 8 个基于分位数的等距时间“箱子”(Bins)。这意味着我们假设在这 8 个固定阶段内,海康威视的股价被强行认定为是某种恒定不变的“常数”。接着代码通过 get_dummies 独热编码技巧,把这些时间箱子转化为了一组互相排斥的虚拟变量,并丢进了一个最基础的 OLS 线性模型中。

列表 7.5: 阶梯函数回归 (Haikang)
# ── 将连续时间轴按分位数切成 8 个等频"箱子" ──
num_time_bins = 8  # 设置时间分箱的数量
# pd.qcut: 基于分位数切分,每个箱子含大致相同数量的交易日
haikang_data['time_cut'] = pd.qcut(haikang_data['Time'], num_time_bins)  # 对连续变量进行等频分箱

# ── 将分箱结果转换为虚拟变量(独热编码) ──
# drop_first=True: 丢弃第 1 个箱子作为基准,避免完全多重共线性
# dtype=float: 显式指定输出为浮点型,避免布尔数组与 numpy/statsmodels 不兼容
# 生成时间分箱的虚拟变量
time_dummy_variables = pd.get_dummies(haikang_data['time_cut'], drop_first=True, dtype=float)

# ── 拟合 OLS 阶梯回归 ──
step_features = time_dummy_variables.values  # 提取为NumPy数组
# sm.add_constant 添加截距列(即基准箱的平均值)
step_features_sm = sm.add_constant(step_features)  # 添加常数项(截距项)
step_regression_model = sm.OLS(price_target, step_features_sm).fit()  # OLS 拟合阶梯回归模型

# 输出系数表:每个虚拟变量的系数代表该箱相对于基准箱的价格差异
print(step_regression_model.summary().tables[1])  # 输出模型摘要统计
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          3.5015      0.241     14.551      0.000       3.030       3.973
x1             2.8333      0.340      8.326      0.000       2.166       3.500
x2             6.3669      0.340     18.699      0.000       5.699       7.034
x3            17.0705      0.340     50.162      0.000      16.403      17.738
x4            23.8551      0.340     70.061      0.000      23.188      24.523
x5            38.8566      0.340    114.180      0.000      38.189      39.524
x6            29.1395      0.340     85.581      0.000      28.472      29.807
x7            25.9185      0.340     76.162      0.000      25.251      26.586
==============================================================================

系数表量化各时间箱相对基准箱的样本均价差;截距是基准箱样本均价。具体数值、最高箱及对应日期必须从当次箱边界和输出读取。这只是历史价格的分段描述,不是“价值”估计或未来状态预测。

图 7.5 把阶梯回归映射回交易日期:蓝线在每个训练期分箱内为常数,并在切点处跳跃。它描述各箱样本均价,不把历史分段解释成经济状态或未来制度断点。

# ── 绘制阶梯函数拟合效果 ──
plt.figure(figsize=(12, 8))  # 创建新图形
plt.plot(haikang_data['trade_date'], price_target, 'k-', alpha=0.3, linewidth=1, label='观测股价')  # 黑色半透明线:原始股价走势

# 预测所有原始数据点 → 生成阶梯状的拟合值
predicted_price_step = step_regression_model.predict(step_features_sm)  # 用阶梯模型预测股价
plt.plot(haikang_data['trade_date'], predicted_price_step, 'b-', linewidth=2, label=f'{num_time_bins}阶梯拟合')  # 蓝色线:阶梯拟合曲线

# ── 用红色虚线标记各分箱的分界点 ──
time_quantiles_indices = pd.qcut(haikang_data['Time'], num_time_bins).unique()  # 获取唯一值
# retbins=True 返回分位数的切点位置
time_bins_edges = pd.qcut(haikang_data['Time'], num_time_bins, retbins=True)[1]  # 对连续变量进行等频分箱
for bin_edge in time_bins_edges[1:-1]:       # 跳过首尾边界
    # 将 Time 索引映射回真实日期
    bin_date_value = haikang_data.iloc[int(bin_edge)]['trade_date']  # 根据分箱边界索引获取对应日期
    plt.axvline(x=bin_date_value, color='r', linestyle='--', alpha=0.3)  # 垂直虚线

plt.xlabel('时间', fontsize=14)  # x 轴标签
plt.ylabel('股价 (元)', fontsize=14)  # y 轴标签
plt.title('海康威视股价的阶梯函数拟合 (分段常数)', fontsize=16)  # 图表标题
plt.legend(fontsize=12)  # 添加图例
plt.grid(True, alpha=0.3)  # 半透明网格线
plt.show()  # 渲染并显示图表
海康威视价格散点叠加按时间箱分段的阶梯状预测线。
图 7.5: 阶梯函数回归拟合

那么,我们该如何科学地决定究竟要切出多少个时间箱子(Bins)呢?这段代码再一次请出了交叉验证(Cross-Validation)这个万能裁判。我们构建了一个包含 5、10、15、20 甚至 30 个箱子的备选清单。为了避免直接使用 qcut 在划分测试集外推时可能导致的逻辑报错,代码非常严谨地引入了专为机器学习流水线设计的 sklearn.preprocessing.KBinsDiscretizer 来处理分箱特征工程。通过外层的 5 折循环切割,我们记录下不同装箱粒度对应的验证集 MSE 误差曲线。在输出的蓝色 CV MSE 折线图中,你会寻找那个误差最低洼的底部节点,从而让数据自己发声,选择出一个在“拟合历史”与“外推未来”之间达成最完美平衡的阶梯数量切分方案。

列表 7.6: 使用交叉验证选择最优箱子数量
from sklearn.model_selection import TimeSeriesSplit  # 导入保持时间先后的前向验证器

# ── 候选分箱数量:从 5 到 30 ──
candidate_bin_counts = [5, 10, 15, 20, 30]  # 定义候选参数列表
step_cv_mse_scores = []  # 存放每种分箱数的平均 CV MSE

# 五折扩展窗口严格保证训练日早于验证日
kfold_iterator = TimeSeriesSplit(n_splits=5, gap=1)  # 留一日 gap 避免边界紧邻泄漏

# 对每种候选分箱数量执行交叉验证
for n_bins in candidate_bin_counts:  # 遍历循环
    current_fold_mse = []  # 当前分箱数下各折的 MSE
    
    time_index_features = haikang_data[['Time']].values  # 特征:时间索引

    # 对每个交叉验证折进行训练和验证
    # 遍历K折交叉验证的训练集和验证集索引
    for train_indices, validation_indices in kfold_iterator.split(time_index_features):
        # ── 分割训练集和验证集 ──
        # 设置features_train参数
        features_train, features_val = time_index_features[train_indices], time_index_features[validation_indices]
        target_train, target_val = price_target[train_indices], price_target[validation_indices]  # 分割目标变量

        # ── 使用 KBinsDiscretizer 进行分箱特征工程 ──
        # 相比 pd.qcut,它能在 sklearn Pipeline 中正确处理训练/测试边界
        from sklearn.preprocessing import KBinsDiscretizer  # 分箱离散化工具
        
        # strategy='quantile': 按分位数分箱;encode='onehot-dense': 输出独热编码矩阵
        # 拟合discretizer_model模型
        discretizer_model = KBinsDiscretizer(n_bins=n_bins, encode='onehot-dense', strategy='quantile')
        binned_features_train = discretizer_model.fit_transform(features_train)   # 在训练集上拟合分箱并转换
        binned_features_val = discretizer_model.transform(features_val)           # 在验证集上仅转换
        
        linear_regressor = LinearRegression().fit(binned_features_train, target_train)  # 在分箱特征上拟合线性回归
        predicted_prices_step = linear_regressor.predict(binned_features_val)  # 在验证集上预测
        current_fold_mse.append(np.mean((target_val - predicted_prices_step) ** 2))  # 验证集 MSE

    # 汇总当前分箱数的平均 CV MSE
    if current_fold_mse:  # 条件判断
        step_cv_mse_scores.append(np.mean(current_fold_mse))  # 5 折平均 MSE
    else:  # 无有效折时的安全处理
        step_cv_mse_scores.append(np.inf)  # 无有效折时设为无穷大

# 只从已经计算的前向验证误差中选择候选
optimal_bin_count = candidate_bin_counts[np.argmin(step_cv_mse_scores)]  # 获取最小验证误差对应的箱数
print(f'前向验证选中的箱子数量: {optimal_bin_count}')  # 输出运行时选择而不写死赢家
前向验证选中的箱子数量: 30

接下来,图 7.6 可视化 CV MSE 随分箱数量的变化;最低点只是当前候选集内的运行时选择。

# ── 绘制 CV MSE 随分箱数量的变化曲线 ──
plt.figure(figsize=(10, 6))  # 创建画布
plt.plot(candidate_bin_counts, step_cv_mse_scores, 'bo-', linewidth=2, markersize=8)  # 蓝色圆点折线
plt.xlabel('分箱数量 (Bins)', fontsize=14)  # x 轴标签
plt.ylabel('CV MSE', fontsize=14)  # y 轴标签
plt.title('阶梯函数:箱子数量选择', fontsize=16)  # 图表标题
plt.grid(True, alpha=0.3)  # 半透明网格线
plt.show()  # 渲染并显示图表
折线以箱数为横轴、前向验证MSE为纵轴,展示五个候选分箱复杂度的误差。
图 7.6: 阶梯函数候选箱数的前向验证均方误差

图 7.6 与当次前向验证输出读取所选箱数,并同时报告各折误差。该选择只是在候选集合中的预测损失最小者;箱宽不等于经济周期,也不能据此推断均值回复。

7.4 基函数 (Basis Functions)

多项式和分段常数回归模型实际上是基函数(basis function)方法的特例。其思想是拥有一族可以应用于变量 \(X\) 的函数或变换:\(b_1(X), b_2(X), \ldots, b_K(X)\)。我们不拟合 \(X\) 的线性模型,而是拟合模型:

\[ y_i = \beta_0 + \beta_1 b_1(x_i) + \beta_2 b_2(x_i) + \beta_3 b_3(x_i) + \cdots + \beta_K b_K(x_i) + \epsilon_i \tag{7.3}\]

请注意,基函数 \(b_1(\cdot), b_2(\cdot), \ldots, b_K(\cdot)\) 是固定且已知的(换句话说,我们提前选择这些函数)。对于多项式回归,基函数是 \(b_j(x_i) = x_i^j\),对于分段常数函数,它们是 \(b_j(x_i) = I(c_j \leq x_i < c_{j+1})\)。我们可以将 式 7.3 视为以预测变量 \(b_1(x_i), b_2(x_i), \ldots, b_K(x_i)\) 的标准线性模型。因此,我们可以使用最小二乘法来估计 式 7.3 中的未知回归系数。重要的是,这意味着 章节 3 中讨论的所有线性模型推断工具(如系数估计的标准误差和模型整体显著性的F统计量)在这种设置下都是可用的。

到目前为止,我们考虑了将多项式函数和分段常数函数用作基函数;然而,许多其他选择也是可能的。例如,我们可以使用小波或傅里叶级数来构造基函数。在下一节中,我们将研究一种非常常见的基函数选择:回归样条(regression splines)。

Note: 基函数方法的灵活性

基函数框架为我们提供了统一的视角来理解各种非线性方法:

  1. 统一表示:多项式回归、阶梯函数、样条等都可以表示为基函数的线性组合
  2. 推断工具:由于本质上是线性模型,所有线性模型的推断工具都适用
  3. 灵活扩展:可以轻松组合不同类型的基函数

常见基函数类型

  • 多项式基函数\(1, x, x^2, x^3, \ldots\)
  • 分段常数基函数:指示函数
  • B样条基函数:分段多项式,具有良好性质
  • 傅里叶基函数\(\sin(2\pi k x), \cos(2\pi k x)\),用于周期性数据
  • 小波基函数:适用于多尺度分析
  • 径向基函数\(K(x, x_i) = \exp(-\gamma \|x - x_i\|^2)\)

选择基函数时需要考虑:

  • 数据的性质(平滑性、周期性等)
  • 计算复杂度
  • 可解释性
  • 边界行为

7.5 回归样条 (Regression Splines)

现在我们讨论一类灵活的基函数,它扩展了我们刚刚看到的多项式回归和分段常数回归方法。

7.5.1 分段多项式 (Piecewise Polynomials)

不在整个 \(X\) 范围上拟合高次多项式,分段多项式回归(piecewise polynomial regression)涉及在 \(X\) 的不同区域拟合不同的低次多项式。例如,分段三次多项式通过拟合形式为

\[ y_i = \beta_0 + \beta_1 x_i + \beta_2 x_i^2 + \beta_3 x_i^3 + \epsilon_i \tag{7.4}\]

的三次回归模型来工作,其中系数 \(\beta_0, \beta_1, \beta_2, \beta_3\)\(X\) 范围的不同部分中不同。系数变化的点称为节点(knots)。例如,没有节点的分段三次多项式只是标准的三次多项式,如 式 7.1\(d=3\) 的情况。在点 \(c\) 处有单个节点的分段三次多项式采用以下形式:

\[ y_i = \begin{cases} \beta_{01} + \beta_{11} x_i + \beta_{21} x_i^2 + \beta_{31} x_i^3 + \epsilon_i, & \text{如果 } x_i < c \\ \beta_{02} + \beta_{12} x_i + \beta_{22} x_i^2 + \beta_{32} x_i^3 + \epsilon_i, & \text{如果 } x_i \geq c \end{cases} \]

换句话说,我们对数据拟合两个不同的多项式函数,一个在 \(x_i < c\) 的观测子集上,另一个在 \(x_i \geq c\) 的观测子集上。

使用更多节点会产生更灵活的分段多项式。一般来说,如果我们在 \(X\) 的范围内放置 \(K\) 个不同的节点,最终将拟合 \(K+1\) 个不同的三次多项式。

7.5.2 约束与样条

如果没有对分段多项式施加约束,拟合曲线可能看起来不合理。为了解决这个问题,我们可以在约束下拟合分段多项式,要求拟合曲线必须是连续的。换句话说,在节点处不能有跳跃。

我们还可以添加额外约束:分段多项式的一阶和二阶导数在节点处连续。换句话说,我们要求分段多项式不仅在节点处连续,而且非常平滑。每我们对分段三次多项式施加一个约束,实际上都会释放一个自由度,通过降低所得分段多项式拟合的复杂性。

在左下图中,我们施加了三个约束(连续性、一阶导数连续性和二阶导数连续性),因此剩下5个自由度。下图中的曲线称为三次样条(cubic spline)。一般来说,具有 \(K\) 个节点的三次样条总共使用 \(4+K\) 个自由度。

在图7.3中,右下图是线性样条(linear spline),在 age=50 处连续。\(d\) 次样条的一般定义是:它是分段 \(d\) 次多项式,在每个节点处直到 \(d-1\) 阶的导数连续。因此,线性样条通过在每个节点定义的预测变量空间区域中拟合线,要求在每个节点处连续来获得。

7.5.3 样条基表示

我们刚刚看到的回归样条可能看起来有些复杂:如何拟合分段 \(d\) 次多项式,并约束它(可能还有它的前 \(d-1\) 阶导数)连续?事实证明,我们可以使用基模型 式 7.3 来表示回归样条。具有 \(K\) 个节点的三次样条可以建模为:

\[ y_i = \beta_0 + \beta_1 b_1(x_i) + \beta_2 b_2(x_i) + \cdots + \beta_{K+3} b_{K+3}(x_i) + \epsilon_i \tag{7.5}\]

对于基函数 \(b_1, b_2, \ldots, b_{K+3}\) 的适当选择。然后可以使用最小二乘法拟合模型 式 7.5

就像有多种表示多项式的方法一样,使用不同的基函数选择在 式 7.5 中表示三次样条也有许多等效方法。使用 式 7.5 表示三次样条的最直接方法是从三次多项式的基开始——即 \(x, x^2, x^3\)——然后为每个节点添加一个截断幂基函数(truncated power basis function)。

截断幂基函数定义为:

\[ h(x, \xi) = (x-\xi)_+^3 = \begin{cases} (x-\xi)^3, & \text{如果 } x > \xi \\ 0, & \text{否则} \end{cases} \tag{7.6}\]

其中 \(\xi\) 是节点。可以证明,将形式为 \(\beta_4 h(x, \xi)\) 的项添加到三次多项式模型 式 7.4 中,将仅在 \(\xi\) 处的三阶导数中产生不连续性;该函数将在每个节点处保持连续,具有连续的一阶和二阶导数。

换句话说,为了对具有 \(K\) 个节点的数据集拟合三次样条,我们对具有截距和 \(3+K\) 个预测变量的形式为 \(X, X^2, X^3, h(X, \xi_1), h(X, \xi_2), \ldots, h(X, \xi_K)\) 的预测变量执行最小二乘回归,其中 \(\xi_1, \ldots, \xi_K\) 是节点。这相当于估计总共 \(K+4\) 个回归系数;因此,拟合具有 \(K\) 个节点的三次样条使用 \(K+4\) 个自由度。

不幸的是,样条在预测变量的外围范围可能具有高方差——即当 \(X\) 取非常小或非常大的值时。自然样条(natural spline)是具有额外边界约束的回归样条:函数被要求在边界处是线性的(在 \(X\) 小于最小节点或大于最大节点的区域中)。这个额外约束意味着自然样条通常在边界处产生更稳定的估计。

数学推导:自然样条的边界条件

为什么自然三次样条(Natural Cubic Splines)相比于普通三次样条在边界区域(即小于最小节点 \(\xi_1\) 和大于最大节点 \(\xi_K\) 的区域)具有更稳定的估计和更小的方差?这在数学上是如何约束的?

回忆三次回归样条可以使用截断幂基函数(Truncated Power Basis)表示: \[ f(x) = \beta_0 + \beta_1 x + \beta_2 x^2 + \beta_3 x^3 + \sum_{k=1}^K \theta_k (x - \xi_k)_+^3 \] 对于自然样条,我们要求其在边界 \([-\infty, \xi_1]\)\([\xi_K, \infty]\) 内严格保持线性。这意味着在这些外围区域内,多项式的二次和三次项的总系数必须化简为 0。

  1. 左边界条件 (\(x < \xi_1\)): 当 \(x < \xi_1\) 时,所有包含 \((x-\xi_k)_+^3\) 的节点截断项直接为 0。此时基函数退化为: \[ f(x) = \beta_0 + \beta_1 x + \beta_2 x^2 + \beta_3 x^3 \] 为了保证其呈严格线性状态,必须有: \[ \beta_2 = 0 \quad \text{和} \quad \beta_3 = 0 \]

  2. 右边界条件 (\(x > \xi_K\)): 当 \(x > \xi_K\) 时,所有截断幂项 \((x-\xi_k)_+^3\) 都被激活并展开为 \((x-\xi_k)^3\)\[ (x-\xi_k)^3 = x^3 - 3\xi_k x^2 + 3\xi_k^2 x - \xi_k^3 \] 我们将这些展开项放回原函数,并提取 \(x^3\)\(x^2\) 的同类项。为了让函数在右边界外也保持为线性,其二次和三次的总合成系数必须同时为 0(结合左边界得出的 \(\beta_2=0, \beta_3=0\)):

  • 三次项限制为 0\(\sum_{k=1}^K \theta_k = 0\)
  • 二次项限制为 0\(-3 \sum_{k=1}^K \theta_k \xi_k = 0 \Rightarrow \sum_{k=1}^K \theta_k \xi_k = 0\)

综上所述,左边界的 2 个约束和右边界的 2 个约束总共构成了 4 个额外的线性约束条件。这解释了为什么拥有 \(K\) 个内部节点的自然三次样条恰巧能节省出 4 个参数位,总共只需要 \(K\) 个自由度(而普通的一次函数带 K 个节点的三次样条需要 \(K+4\) 个自由度)。因为它们在数据稀疏的外围区间被“束缚”为相对平稳的线性延伸,这种更受限的刚性数学结构直接阻止了曲线对极端离群值的非理性振荡,从而大幅压缩了边界处的预测方差。

7.5.4 回归样条的Python实现

让我们使用 patsystatsmodels 来拟合回归样条,以平滑海康威视股价的长期走势:

列表 7.7: 导入样条相关库
from patsy import dmatrix                    # patsy:用 R 风格公式构建设计矩阵(样条基函数)
from sklearn.linear_model import LinearRegression  # 线性回归器
import statsmodels.api as sm                  # statsmodels:提供统计推断(p 值、置信区间)

回归三次样条是在节点之间使用分段三次多项式,但各段并非独立估计。patsy.bs 构造的 B-spline 基共享一组回归系数,函数、一阶导和二阶导在内部节点处的连续性由样条空间(等价地,由线性约束)内建。下面的 patsy.bs + OLS未惩罚回归样条:它没有平滑参数 \(\lambda\),也没有粗糙度惩罚。\(\lambda\int [g''(t)]^2dt\) 只在后面的平滑样条/GAM 语境中出现;因此不能把节点连续性称为“连续性惩罚”。

列表 7.8: 拟合三次样条 (Haikang)
# ── 在时间轴的 20%/40%/60%/80% 分位数处放置 4 个内部节点 ──
import numpy as np  # 数值计算库(分位数节点计算等)
spline_knots = np.quantile(haikang_data['Time'], [0.2, 0.4, 0.6, 0.8])  # 在 20%/40%/60%/80% 分位数处设置内部节点

# ── 用 patsy 的 bs() 构建三次 B-Spline 基函数矩阵 ──
# degree=3 → 三次样条;include_intercept=False → 截距交给 OLS
# 构建三次样条(Cubic Spline)基函数设计矩阵
cubic_spline_features = dmatrix('bs(Time, knots=spline_knots, degree=3, include_intercept=False)',
                   # 提取为NumPy数组
                   {'Time': haikang_data['Time'].values, 'spline_knots': spline_knots},
                   return_type='dataframe')  # 指定返回DataFrame格式

# 使用 statsmodels OLS 拟合(以便后续做统计推断)
cubic_spline_model = sm.OLS(price_target, cubic_spline_features).fit()  # 训练/拟合模型

模型已经拟合完成,接下来在 图 7.7 的连续时间网格上叠加三次样条拟合与原始股价,并明确标出训练样本确定的节点。

# ── 在连续时间网格上生成预测,使样条曲线平滑 ──
flattened_time_grid = time_grid.flatten()  # dmatrix 需要一维数组
# 用与训练时完全相同的 bs() 公式和节点,在网格上构建基函数
# 为预测网格构建三次样条设计矩阵
grid_cubic_spline_features = dmatrix('bs(time_grid, knots=spline_knots, degree=3, include_intercept=False)',
                          # 传入预测网格数据
                          {'time_grid': flattened_time_grid, 'spline_knots': spline_knots},
                          return_type='dataframe')  # 指定返回DataFrame格式

# 在网格上预测
predicted_mean_spline = cubic_spline_model.predict(grid_cubic_spline_features)  # 使用模型进行预测

# ── 绘图:观测值 + 三次样条拟合 + 节点位置 ──
plt.figure(figsize=(12, 8))  # 创建新图形
# 绑制折线图
plt.plot(haikang_data['trade_date'], price_target, 'k-', alpha=0.3, linewidth=1, label='观测股价')
plt.plot(date_grid, predicted_mean_spline, 'b-', linewidth=2, label='三次样条拟合')  # 在真实交易日期网格上绘制三次样条拟合曲线

# 在节点处画红色垂直虚线,直观展示"分段拼接"的位置
for current_knot in spline_knots:  # 遍历循环
    bin_date_value = haikang_data.iloc[int(current_knot)]['trade_date']  # 节点索引 → 日期
    plt.axvline(x=bin_date_value, color='r', linestyle='--', alpha=0.7, linewidth=1)  # 添加垂直参考线
    
plt.xlabel('时间', fontsize=14)  # x 轴标签
plt.ylabel('股价 (元)', fontsize=14)  # y 轴标签
plt.title('三次样条回归拟合 (B-Spline)', fontsize=16)  # 图表标题
plt.legend(fontsize=12)  # 添加图例
plt.grid(True, alpha=0.3)  # 半透明网格线
plt.show()  # 渲染并显示图表
海康威视价格散点叠加三次回归样条曲线,并标出内部结点位置。
图 7.7: 三次样条回归拟合

普通三次样条的边界形状可能不稳定。下面改用自然三次样条(patsy.cr):它在边界之外施加线性约束,从而降低高阶多项式式的边界振荡风险。这里的 df=10 只是候选复杂度,是否改善样本外误差仍须由同一时间验证和支持域证据判断。

列表 7.9: 拟合自然样条 (Natural Spline)
# ── 自然三次样条:cr() 函数,df=10 指定 10 个自由度 ──
# cr() 在两端施加"线性"约束,防止边界发散
# 提取为NumPy数组
natural_spline_features = dmatrix('cr(Time, df=10)', {'Time': haikang_data['Time'].values},
                    return_type='dataframe')  # 指定返回DataFrame格式

# OLS 拟合
natural_spline_model = sm.OLS(price_target, natural_spline_features).fit()  # 训练/拟合模型

自然样条模型拟合完毕后,我们同样在时间网格上生成预测值,并将其与观测股价叠加绘图。通过与前文三次样条的对比,读者可以清楚地观察到自然样条在序列两端的边界约束如何有效抑制了曲线的过度发散。

图 7.8 展示自然样条在训练支持域内的拟合形状;边界线性约束不等于样本外预测有效。

# ── 在网格上构建自然样条基函数并预测 ──
grid_natural_spline_features = dmatrix('cr(time_grid, df=10)',  # 为预测网格构建自然样条设计矩阵
                          {'time_grid': flattened_time_grid},  # 传入预测网格数据
                          return_type='dataframe')  # 指定返回DataFrame格式

predicted_mean_natural = natural_spline_model.predict(grid_natural_spline_features)  # 自然样条在网格上的预测值

# ── 绘图:自然样条拟合曲线(绿色) ──
plt.figure(figsize=(12, 8))  # 创建新图形
plt.plot(haikang_data['trade_date'], price_target, 'k-', alpha=0.3, linewidth=1, label='观测股价')  # 黑色半透明线:原始股价走势
plt.plot(date_grid, predicted_mean_natural, 'g-', linewidth=2, label='自然样条 (df=10)')  # 绿色线:在真实交易日期网格上的自然样条拟合曲线

plt.xlabel('时间', fontsize=14)  # x 轴标签
plt.ylabel('股价 (元)', fontsize=14)  # y 轴标签
plt.title('自然样条回归拟合', fontsize=16)  # 图表标题
plt.legend(fontsize=12)  # 添加图例
plt.grid(True, alpha=0.3)  # 半透明网格线
plt.show()  # 渲染并显示图表
海康威视价格散点叠加自然样条曲线,边界附近曲线趋于线性。
图 7.8: 自然样条回归拟合

图中绿色曲线是当前前复权样本上、df=10 候选的拟合。自然样条在边界外施加线性约束,通常比同阶全局高次多项式少一种端点振荡来源;但本次端点是否更稳定、哪些波段被跟踪以及 df=10 是否合适,都必须由现场图、前向验证和支持密度判断,不能在输出前写成固定结论。

7.5.5 选择节点数量和位置

当拟合样条时,我们应该在哪里放置节点?回归样条在包含大量节点的区域最灵活,因为在这些区域中多项式系数可以快速变化。因此,一种选择是在函数似乎变化最快的地方放置更多节点,而在它看起来更稳定的地方放置较少节点。

虽然这个选项可以很好地工作,但在实践中,通常以均匀的方式放置节点。一种方法是指定所需的自由度,然后让软件自动将相应数量的节点放置在数据的均匀分位数处。

下面在 3—15 个自由度候选中使用前向折比较验证 MSE。最低点只是在当前切分与候选集内的选择证据,不是全局最优或稳定机制;设计矩阵的节点和边界必须只由相应训练折确定。

列表 7.10: 使用交叉验证选择样条的自由度
from patsy import build_design_matrices  # 用训练期 design_info 转换验证数据
from sklearn.model_selection import TimeSeriesSplit  # 导入前向时间验证器

# ── 候选自由度列表:从 3(最简单)到 15(最灵活) ──
candidate_degrees_of_freedom = [3, 4, 5, 6, 7, 8, 10, 12, 15]  # 定义候选参数列表
spline_cv_mse_scores = []  # 存放每种 df 的平均 CV MSE

kfold_iterator = TimeSeriesSplit(n_splits=5, gap=1)  # 严格保持训练折早于验证折

# 对每个候选自由度执行交叉验证
for current_dof in candidate_degrees_of_freedom:  # 遍历循环
    current_fold_mse = []  # 当前 df 下各折 MSE

    # 对每个交叉验证折进行训练和验证
    for train_indices, validation_indices in kfold_iterator.split(haikang_data):  # 遍历K折交叉验证的训练集和验证集索引
        data_train, data_validation = haikang_data.iloc[train_indices], haikang_data.iloc[validation_indices]  # 分割训练集和验证集

        # 分别在训练集和验证集上构建自然样条基函数
        features_train = dmatrix(f'cr(Time, df={current_dof})',  # 为训练集构建自然样条设计矩阵
                         {'Time': data_train['Time'].values},  # 提取为NumPy数组
                         return_type='dataframe')  # 指定返回DataFrame格式
        features_val = build_design_matrices([features_train.design_info], {'Time': data_validation['Time'].values})[0]  # 复用训练期节点、边界与列顺序

        # 在训练集上拟合 OLS
        # 提取为NumPy数组
        current_spline_model = sm.OLS(data_train['Price'].values, features_train).fit()

        # 在验证集上预测并计算 MSE
        predicted_prices_spline = current_spline_model.predict(features_val)  # 使用模型进行预测
        current_mse = np.mean((data_validation['Price'].values - predicted_prices_spline) ** 2)  # 计算当前折的验证集 MSE
        current_fold_mse.append(current_mse)  # 将当前折 MSE 存入列表

    spline_cv_mse_scores.append(np.mean(current_fold_mse))  # 5 折平均 MSE

接下来,图 7.9 绘制 CV MSE 随自由度的变化,并从实际输出中选择当前候选集内的自由度。

# ── 绘制 CV MSE 随自由度变化的曲线 ──
plt.figure(figsize=(10, 6))  # 创建画布
plt.plot(candidate_degrees_of_freedom, spline_cv_mse_scores, 'bo-', linewidth=2, markersize=8)  # 蓝色圆点折线
plt.xlabel('自由度', fontsize=14)  # x 轴标签
plt.ylabel('交叉验证均方误差', fontsize=14)  # y 轴标签
plt.title('自然样条:自由度选择', fontsize=16)  # 图表标题
plt.grid(True, alpha=0.3)  # 半透明网格线

# 用红色虚线标出最优自由度的位置
# 获取最小值的索引
optimal_degrees_of_freedom = candidate_degrees_of_freedom[np.argmin(spline_cv_mse_scores)]
minimum_cv_mse = min(spline_cv_mse_scores)  # 取出最小的 CV MSE 值
plt.axvline(x=optimal_degrees_of_freedom, color='r', linestyle='--',  # 添加垂直参考线
            label=f'最优自由度 = {optimal_degrees_of_freedom}')  # 设置图例标签文本
plt.legend(fontsize=12)  # 添加图例
plt.show()  # 渲染并显示图表

print(f'最优自由度: {optimal_degrees_of_freedom}')  # 打印最优 df
print(f'对应的最小交叉验证MSE: {minimum_cv_mse:.4f}')  # 打印最小 CV MSE
折线以自然样条自由度为横轴、前向验证MSE为纵轴,并以竖虚线标出当次最低误差对应的自由度。
图 7.9: 自然样条候选自由度的前向验证均方误差
最优自由度: 4
对应的最小交叉验证MSE: 129.3479

从当次前向验证输出读取所选自由度和 MSE,并检查各折是否一致。自然样条自由度控制设计矩阵维数,不等同于节点数或“年度级信息维度”;若曲线在最优点附近平坦,应把多个复杂度视为近似等价并优先较简单者。

7.5.6 与多项式回归的比较

让我们比较一下自然样条和高次多项式的表现:

为了彻底为你呈现多项式回归和自然样条的实力鸿沟,这段收尾代码安排了一场令人窒息的同级别“自由度”生死决战。代码左手训练了一个极其复杂狂暴的 15 次多项式,右手则构造了一个同样消耗 15 个自由度顶额配置的自然样条。我们将这俩放进了同样的 5 折交叉验证炼丹炉里。 下图只展示当前训练样本中的两种拟合形状。若 15 次全局多项式在边界振荡而自然样条更平稳,这与自然样条的边界约束一致;但图形本身不证明它捕捉了真实信号,也不证明样条普遍胜出。模型选择必须读取同一前向折误差、离散度、支持域和基线;若样条未改善验证损失,应明确报告“无增量证据”。

列表 7.11: 比较自然样条与多项式回归
# ── 对手 A:拟合 15 次多项式 ──
polynomial_model_15deg = Pipeline([  # 构建多项式回归Pipeline(标准化+多项式特征+线性回归)
    ('poly', PolynomialFeatures(degree=15)),  # 生成 16 维特征 [1, X, ..., X¹⁵]
    ('linear', LinearRegression())  # 初始化线性回归模型
])  # 完成构建
polynomial_model_15deg.fit(time_features, price_target)  # 在全部数据上训练 15 次多项式模型

# ── 对手 B:拟合 15 自由度的自然样条 ──
natural_spline_features_15dof = dmatrix('cr(Time, df=15)',  # 构建自然样条(Natural Spline)设计矩阵
                       {'Time': haikang_data['Time'].values},  # 提取为NumPy数组
                       return_type='dataframe')  # 指定返回DataFrame格式
natural_spline_model_15dof = sm.OLS(haikang_data['Price'].values, natural_spline_features_15dof).fit()  # OLS 拟合 15 自由度自然样条

接下来,通过 5 折交叉验证对比两种方法的泛化能力。

# ── 5 折 CV 对决:多项式 vs 自然样条 ──
# sklearn cross_val_score 可直接用于 Pipeline
# 执行交叉验证评估
cv_mse_poly_15deg = -cross_val_score(polynomial_model_15deg, time_features, price_target, cv=kfold_iterator,
                               scoring='neg_mean_squared_error').mean()  # 计算均值

# 对样条需手动循环(statsmodels 不兼容 sklearn cross_val_score)
fold_mse_natural_15dof = []  # 初始化各折MSE存储列表
# 手动循环 5 折交叉验证,计算自然样条的 CV MSE
for train_indices, validation_indices in kfold_iterator.split(haikang_data):  # 遍历K折交叉验证的训练集和验证集索引
    data_train, data_validation = haikang_data.iloc[train_indices], haikang_data.iloc[validation_indices]  # 分割训练/验证集
    features_train = dmatrix('cr(Time, df=15)',  # 为训练集构建自然样条设计矩阵
                     {'Time': data_train['Time'].values},  # 提取为NumPy数组
                     return_type='dataframe')  # 训练集样条基函数
    features_val = build_design_matrices(  # 复用训练折的样条基定义
        [features_train.design_info],  # 训练折节点、边界和列顺序
        {'Time': data_validation['Time'].values}  # 验证折时间
    )[0]  # 取出验证设计矩阵
    current_spline_model = sm.OLS(data_train['Price'].values, features_train).fit()  # 在训练集上拟合 OLS
    predicted_prices_spline = current_spline_model.predict(features_val)  # 在验证集上预测
    fold_mse_natural_15dof.append(np.mean((data_validation['Price'].values - predicted_prices_spline) ** 2))  # 计算并存储当前折 MSE
cv_mse_natural_15dof = np.mean(fold_mse_natural_15dof)  # 5 折平均 MSE

print(f'15次多项式的交叉验证MSE: {cv_mse_poly_15deg:.4f}')  # 打印多项式 CV MSE
print(f'15自由度自然样条的交叉验证MSE: {cv_mse_natural_15dof:.4f}')  # 打印样条 CV MSE
15次多项式的交叉验证MSE: 888782456260.5305
15自由度自然样条的交叉验证MSE: 508.9858

比较表应报告当次前向验证中高次多项式与自然样条的逐折 MSE。若高次多项式在某些折发生边界发散,这说明该参数化在此数据与外推范围下数值不稳;它不是对所有数据集的普遍禁令。自然样条的边界线性约束通常更稳,但仍须在同一时间切分和预先冻结的候选集上评估。

图 7.10 只比较两个在全样本拟合后的支持域内形状,模型选择另读前向验证误差。

# ── 在连续网格上分别生成两种模型的预测 ──
grid_poly_features_15deg = PolynomialFeatures(degree=15).fit_transform(time_grid)  # 拟合并转换数据
predicted_price_poly_15deg = polynomial_model_15deg.predict(time_grid)         # 15 次多项式预测

grid_natural_features_15dof = dmatrix('cr(time_grid, df=15)',  # 为预测网格构建自然样条设计矩阵
                              {'time_grid': flattened_time_grid},  # 传入预测网格数据
                              return_type='dataframe')  # 网格上 15df 自然样条基函数
predicted_price_natural_15dof = natural_spline_model_15dof.predict(grid_natural_features_15dof)  # 自然样条预测

# ── 左右两幅子图并排对比 ──
fig, axes = plt.subplots(1, 2, figsize=(18, 7))  # 创建子图布局

# 左图:15 次多项式 — 注意边界处的发散
# 绘制散点图
axes[0].scatter(haikang_data['Time'], haikang_data['Price'], alpha=0.3, s=20, label='观测数据')
axes[0].plot(time_grid, predicted_price_poly_15deg, 'r-', linewidth=3, label='15次多项式')  # 红色线:多项式拟合
axes[0].set_xlabel('时间 (交易日)', fontsize=14)  # x 轴标签
axes[0].set_ylabel('股价 (元)', fontsize=14)  # y 轴标签
axes[0].set_title('15次多项式拟合', fontsize=16)  # 子图标题
axes[0].legend(fontsize=12)  # 添加图例
axes[0].grid(True, alpha=0.3)  # 半透明网格线

# 右图:15 自由度自然样条 — 边界稳定
# 绘制散点图
axes[1].scatter(haikang_data['Time'], haikang_data['Price'], alpha=0.3, s=20, label='观测数据')
axes[1].plot(time_grid, predicted_price_natural_15dof, 'b-', linewidth=3, label='15自由度自然样条')  # 蓝色线:自然样条拟合
axes[1].set_xlabel('时间 (交易日)', fontsize=14)  # x 轴标签
axes[1].set_ylabel('股价 (元)', fontsize=14)  # y 轴标签
axes[1].set_title('15自由度自然样条拟合', fontsize=16)  # 子图标题
axes[1].legend(fontsize=12)  # 添加图例
axes[1].grid(True, alpha=0.3)  # 半透明网格线

plt.tight_layout()  # 自动调整子图间距
plt.show()  # 渲染并显示图表
同一坐标系比较自然样条与多项式拟合在样本内部及边界处的曲线形状。
图 7.10: 自然样条与多项式回归的比较

回归样条通常比多项式回归提供更优的结果。这是因为:

  1. 局部灵活性:样条通过增加节点数量但保持次数固定来引入灵活性,而多项式必须使用高次来产生灵活拟合

  2. 稳定性:样条通常产生更稳定的估计,因为它们避免了高次多项式在边界处的振荡行为

  3. 适应性:样条允许我们在函数 \(f\) 似乎变化快速的区域放置更多节点(因此更灵活),而在 \(f\) 看起来更稳定的区域放置较少节点

  4. 边界行为:自然样条在边界处强制线性约束,避免了多项式的边界振荡问题

实践建议

  • 对于大多数应用,优先选择样条而非高次多项式
  • 使用自然样条来避免边界问题
  • 通过交叉验证选择自由度(通常4-10个自由度足够)
  • 如果数据量充足,可以考虑放置更多节点以捕捉局部特征

7.6 平滑样条 (Smoothing Splines)

在上一节中,我们讨论了回归样条,它是通过指定一组节点、产生一系列基函数,然后使用最小二乘法来估计样条系数而创建的。我们现在介绍一种略有不同的方法,它也产生样条。

7.6.1 平滑样条概述

在对一组数据拟合平滑曲线时,我们真正想要做的是找到某个函数,比如 \(g(x)\),它能很好地拟合观测数据:也就是说,我们希望

\[ \text{RSS} = \sum_{i=1}^n (y_i - g(x_i))^2 \]

很小。然而,这种方法有一个问题。如果不对 \(g(x_i)\) 施加任何约束,我们总是可以通过选择 \(g\) 使得它插值所有 \(y_i\) 来使RSS为零。这样的函数会严重过拟合数据——它太灵活了。我们真正想要的是一个使RSS小,但也平滑的函数 \(g\)

如何确保 \(g\) 是平滑的?有多种方法可以做到这一点。一种自然的方法是找到最小化以下函数的 \(g\)

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

其中 \(\lambda\) 是非负调整参数(tuning parameter)。最小化 式 7.7 的函数 \(g\) 称为平滑样条(smoothing spline)。

式 7.7 是什么意思?式 7.7 采用”损失+惩罚”的形式,类似于 章节 6 中岭回归和Lasso的背景。项 \(\sum_{i=1}^n (y_i - g(x_i))^2\) 是鼓励 \(g\) 很好拟合数据的损失函数,项 \(\lambda \int g''(t)^2 dt\) 是惩罚 \(g\) 变异性的惩罚项。符号 \(g''(t)\) 表示函数 \(g\) 的二阶导数。一阶导数 \(g'(t)\) 测量函数在 \(t\) 处的斜率,二阶导数对应斜率的变化量。因此,广义地说,函数的二阶导数是其粗糙度(roughness)的度量:如果 \(g(t)\)\(t\) 附近非常波动,它的绝对值就大;否则接近零。(直线的二阶导数为零;注意直线是完全平滑的。)

积分符号 \(\int\) 是一个积分,我们可以将其视为在整个 \(t\) 范围内的总和。换句话说,\(\int g''(t)^2 dt\) 只是函数 \(g'(t)\) 在其整个范围内的总变化的度量。如果 \(g\) 非常平滑,那么 \(g'(t)\) 将接近常数,\(\int g''(t)^2 dt\) 将取小值。相反,如果 \(g\) 跳跃且可变,那么 \(g'(t)\) 将显著变化,\(\int g''(t)^2 dt\) 将取大值。因此,在 式 7.7 中,\(\lambda \int g''(t)^2 dt\) 鼓励 \(g\) 平滑。\(\lambda\) 的值越大,\(g\) 越平滑。

\(\lambda = 0\) 时,式 7.7 中的惩罚项没有影响,因此函数 \(g\) 将非常跳跃,并且将精确插值训练观测值。当 \(\lambda \to \infty\) 时,\(g\) 将完全平滑——它只是一条尽可能接近训练点的直线。事实上,在这种情况下,\(g\) 将是线性最小二乘线,因为 式 7.7 中的损失函数相当于最小化残差平方和。对于 \(\lambda\) 的中间值,\(g\) 将近似训练观测值但会有所平滑。我们看到 \(\lambda\) 控制平滑样条的偏差-方差权衡。

可以证明,最小化 式 7.7 的函数 \(g(x)\) 具有一些特殊性质:它是一个分段三次多项式,在唯一值 \(x_1, \ldots, x_n\) 处有节点,并且在每个节点处具有连续的一阶和二阶导数。此外,它在极端节点之外的区域中是线性的。换句话说,最小化 式 7.7 的函数 \(g(x)\) 是在 \(x_1, \ldots, x_n\) 处有节点的自然三次样条!然而,它不是如果应用上述样条基表示小节中描述的基函数方法在 \(x_1, \ldots, x_n\) 处有节点所获得的自然三次样条——相反,它是这样一个自然三次样条的收缩版本,其中 式 7.7 中调整参数 \(\lambda\) 的值控制收缩水平。

7.6.2 选择平滑参数 \(\lambda\)

我们已经看到,平滑样条只是在每个唯一值 \(x_i\) 处有节点的自然三次样条。平滑样条似乎有太多的自由度,因为在每个数据点都有一个节点允许很大的灵活性。但调整参数 \(\lambda\) 控制平滑样条的粗糙度,从而控制有效自由度(effective degrees of freedom)。可以证明,随着 \(\lambda\) 从0增加到 \(\infty\),有效自由度(我们写为 \(df_\lambda\))从 \(n\) 减少到2。

在平滑样条的背景下,为什么我们讨论有效自由度而不是自由度?通常,自由度指的是自由参数的数量,例如多项式或三次样条中拟合的系数数量。虽然平滑样条有 \(n\) 个参数因此有 \(n\) 个名义自由度,但这些 \(n\) 个参数受到严重约束或收缩。因此 \(df_\lambda\) 是平滑样条灵活性的度量——它越高,平滑样条越灵活(且低偏差但高方差)。有效自由度的定义有些技术性。我们可以写成

\[ \hat{g}_\lambda = S_\lambda y \tag{7.8}\]

其中 \(\hat{g}_\lambda\) 是对于特定 \(\lambda\) 选择对 式 7.7 的解——也就是说,它是一个 \(n\) 向量,包含平滑样条在训练点 \(x_1, \ldots, x_n\) 处的拟合值。式 7.8 表明,将平滑样条应用于数据时拟合值的向量可以写成 \(n \times n\) 矩阵 \(S_\lambda\)(对此有公式)乘以响应向量 \(y\)。然后有效自由度定义为

\[ df_\lambda = \sum_{i=1}^n \{S_\lambda\}_{ii} \tag{7.9}\]

即矩阵 \(S_\lambda\) 的对角元素之和。

在拟合平滑样条时,我们不需要选择节点的数量或位置——每个训练观测值 \(x_1, \ldots, x_n\) 都会有一个节点。相反,我们有另一个问题:我们需要选择 \(\lambda\) 的值。毫不奇怪,这个问题的可能解决方案之一是交叉验证。换句话说,我们可以找到使交叉验证RSS最小的 \(\lambda\) 值。事实证明,平滑样条的留一交叉验证误差(LOOCV)可以非常高效地计算,成本基本上与单次拟合相同,使用以下公式:

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

符号 \(\hat{g}_{\lambda}^{(-i)}(x_i)\) 表示此平滑样条在 \(x_i\) 处的拟合值,其中拟合使用除第 \(i\) 个观测值 \((x_i, y_i)\) 之外的所有训练观测值。相比之下,\(\hat{g}_\lambda(x_i)\) 表示拟合到所有训练观测值的平滑样条函数在 \(x_i\) 处的评估。

这个非凡的公式说,我们可以仅使用 \(\hat{g}_\lambda\)(对所有数据的原始拟合)来计算这些留一拟合!我们在 章节 5 中针对最小二乘线性回归有一个非常相似的公式 式 5.2。使用 式 5.2,我们可以非常快速地对之前在本章中讨论的回归样条以及使用任意基函数的最小二乘回归执行LOOCV。

7.6.3 平滑样条的Python实现

让我们使用 pygam 或简单的平滑方法。虽然Python的标准库对平滑样条的支持不如R的 smooth.spline 那么直接(pygam 更通用),但我们可以演示其效果。或者使用 scipy.interpolate.UnivariateSpline

这里我们继续使用 pygam 库中的线性GAM(只包含s(0))来演示平滑样条,因为平滑样条实际上是GAM的一种特例。

平滑样条用 \(\lambda\) 控制拟合误差与二阶导数粗糙度之间的权衡,不需要手工把少量节点作为唯一复杂度旋钮。图 7.11 比较三个候选 \(\lambda\);较小惩罚通常允许更高有效自由度,较大惩罚通常得到更平滑的曲线,但具体形状、过拟合与欠拟合只能从本次输出和时间验证判断。图例中的 EDF 表示当前平滑器的有效自由度。

from pygam import LinearGAM, s as s_gam  # pygam: GAM 专用库;s_gam 代表平滑样条项

# ── 准备特征矩阵和目标向量 ──
time_features = haikang_data[['Time']].values  # shape=(n,1)
price_target = haikang_data['Price'].values     # shape=(n,)

# ── 对比三种 λ(惩罚强度)的平滑样条效果 ──
# λ 越大 → 惩罚越重 → 曲线越平滑;λ 越小 → 曲线越扭曲追逐噪音
smoothing_lambdas = [0.1, 10, 1000]  # 三种不同惩罚强度的候选 λ 值

plt.figure(figsize=(14, 8))  # 创建 14×8 英寸画布
plt.plot(haikang_data['trade_date'], price_target, 'k-', alpha=0.3, linewidth=1, label='观测股价')  # 绘制原始股价散点

colors = ['r', 'g', 'b']  # 三种 λ 对应的线条颜色:红/绿/蓝
for i, lam in enumerate(smoothing_lambdas):  # 遍历三种惩罚强度
    # s_gam(0, lam=lam): 对第 0 个特征施加平滑样条,指定惩罚 λ
    smoothing_gam = LinearGAM(s_gam(0, lam=lam))  # 构建指定 λ 的平滑样条 GAM
    smoothing_gam.fit(time_features, price_target)  # 拟合模型到训练数据

    predicted_price_smooth = smoothing_gam.predict(time_grid)  # 在均匀时间网格上生成预测值

    # EDF (Effective Degrees of Freedom): 有效自由度,衡量模型实际复杂度
    effective_dof = smoothing_gam.statistics_['edof']  # 从拟合统计信息中获取有效自由度

    plt.plot(date_grid, predicted_price_smooth, color=colors[i], linewidth=2,  # 在真实交易日期网格上绘制当前 λ 的拟合曲线
             label=f'λ={lam}, EDF={effective_dof:.1f}')  # 图例标注 λ 和有效自由度

plt.xlabel('时间', fontsize=14)  # x 轴标签
plt.ylabel('股价 (元)', fontsize=14)  # y 轴标签
plt.title('不同平滑参数的平滑样条拟合', fontsize=16)  # 图表标题
plt.legend(fontsize=11)  # 添加图例
plt.grid(True, alpha=0.3)  # 半透明网格线
plt.show()  # 渲染并显示图表
海康威视价格灰线叠加三个不同lambda的平滑样条曲线,图例同时报告各自有效自由度。
图 7.11: 三个惩罚参数下的平滑样条训练拟合

下面用 pygam.gridsearch\(10^{-3}\)\(10^4\) 的 30 个候选 \(\lambda\) 上比较 GCV,并在 图 7.12 展示所选曲线。这条教学路线没有构造验证折,也不保持金融时间顺序;GCV 选中的曲线只是在该准则与候选网格下的拟合复杂度折衷,不能称为时间交叉验证结果或普遍最优模型。用于未来收益预测时,应采用本章后文显式的按时间候选遍历与标签窗口 purge。

# ── pygam 内置网格搜索:高斯未知尺度模型默认比较 GCV,不构造验证折 ──
smoothing_gam = LinearGAM(s_gam(0))  # 初始化仅含一个平滑项的 GAM
# logspace(-3,4,30): 在 10⁻³ 到 10⁴ 之间等比取 30 个候选 λ
smoothing_gam.gridsearch(time_features, price_target, lam=np.logspace(-3, 4, 30), objective='GCV')  # 明确用GCV比较候选λ

optimal_lambda = float(np.ravel(smoothing_gam.lam)[0])  # 提取最优lambda标量值(pygam返回嵌套list)
print(f'GCV准则选中的lambda值: {optimal_lambda:.4f}')  # 打印准则选中的惩罚参数
# 从拟合统计信息中获取最优模型的有效自由度
effective_dof = float(np.ravel(smoothing_gam.statistics_['edof'])[0])  # 提取有效自由度标量值
print(f'对应的有效自由度: {effective_dof:.2f}')  # 打印有效自由度

# ── 绘制 GCV 选中 λ 对应的拟合曲线 ──
plt.figure(figsize=(12, 8))  # 创建 12×8 英寸画布
plt.plot(haikang_data['trade_date'], price_target, 'k-', alpha=0.3, linewidth=1, label='观测股价')  # 绘制原始股价

optimal_predicted_price = smoothing_gam.predict(time_grid)  # 最优模型网格预测
plt.plot(date_grid, optimal_predicted_price, 'b-', linewidth=3,  # 在真实交易日期网格上绘制最优拟合蓝色曲线
         label=f'GCV选中拟合 (λ={optimal_lambda:.2f})')  # 图例标注准则选中的 λ 值

plt.xlabel('时间', fontsize=14)  # x 轴标签
plt.ylabel('股价', fontsize=14)  # y 轴标签
plt.title('通过GCV准则选择的平滑样条', fontsize=16)  # 图表标题
plt.legend(fontsize=12)  # 添加图例
plt.grid(True, alpha=0.3)  # 半透明网格线
plt.show()  # 渲染并显示图表
  0% (0 of 30) |                         | Elapsed Time: 0:00:00 ETA:  --:--:--

 13% (4 of 30) |###                      | Elapsed Time: 0:00:00 ETA:   0:00:00

 26% (8 of 30) |######                   | Elapsed Time: 0:00:00 ETA:   0:00:00

 40% (12 of 30) |#########               | Elapsed Time: 0:00:00 ETA:   0:00:00

 50% (15 of 30) |############            | Elapsed Time: 0:00:00 ETA:   0:00:00

 63% (19 of 30) |###############         | Elapsed Time: 0:00:00 ETA:   0:00:00

 70% (21 of 30) |################        | Elapsed Time: 0:00:00 ETA:   0:00:00

 83% (25 of 30) |####################    | Elapsed Time: 0:00:00 ETA:   0:00:00

 96% (29 of 30) |####################### | Elapsed Time: 0:00:00 ETA:   0:00:00

100% (30 of 30) |########################| Elapsed Time: 0:00:00 Time:  0:00:00

GCV准则选中的lambda值: 0.0092
对应的有效自由度: 19.76
海康威视价格灰线叠加GCV在候选网格内选出的平滑样条蓝线,图例报告运行时lambda。
图 7.12: GCV候选准则所选平滑参数对应的训练拟合

从当次网格搜索输出读取 \(\lambda\) 与有效自由度。较小的 \(\lambda\) 只表示在该候选网格和损失函数下允许更多弯曲;EDF 是平滑器复杂度摘要,不是数据“有效信息维度”,也不能与自然样条自由度作机械的一一对应。

Note: 平滑样条与回归样条的区别

虽然平滑样条和回归样条都产生类似的结果,但它们在几个重要方面有所不同:

  1. 节点选择
    • 回归样条:用户选择节点数量和位置
    • 平滑样条:在每个数据点处都有节点
  2. 控制复杂性的方式
    • 回归样条:通过节点数量控制
    • 平滑样条:通过平滑参数 \(\lambda\) 控制
  3. 计算复杂度
    • 回归样条:\(O(n)\)(通常节点数远小于 \(n\)
    • 平滑样条:\(O(n^2)\)\(O(n^3)\)
  4. 灵活性
    • 回归样条:更适合大数据集
    • 平滑样条:更自动,但对于大数据集计算成本高

实践建议

  • 对于大数据集(\(n > 1000\)),优先使用回归样条
  • 对于中小数据集,平滑样条更方便(不需要选择节点)
  • 平滑样条的LOOCV计算非常高效,适合调参
  • 两者在适当调整后通常产生相似结果

7.7 局部回归 (Local Regression)

局部回归(local regression)是拟合灵活非线性函数的另一种方法,它涉及仅使用附近的训练观测值在目标点 \(x_0\) 处计算拟合。

7.7.1 局部回归算法

局部回归的算法如下(算法7.1):

算法7.1 在 \(X=x_0\) 处的局部回归

  1. 收集最接近 \(x_0\) 的训练点分数 \(s = k/n\)
  2. 为这个邻域中的每个点分配权重 \(K_{i0} = K(x_i, x_0)\),使得离 \(x_0\) 最远的点权重为零,最近的点权重最高。除了这 \(k\) 个最近邻之外,所有点的权重都为零。
  3. 使用上述权重,通过找到最小化以下目标的 \(\hat{\beta}_0\)\(\hat{\beta}_1\),对 \(y_i\)\(x_i\) 进行加权最小二乘回归:

\[ \sum_{i=1}^n K_{i0} (y_i - \beta_0 - \beta_1 x_i)^2 \tag{7.10}\]

  1. \(x_0\) 处的拟合值为 \(\hat{f}(x_0) = \hat{\beta}_0 + \hat{\beta}_1 x_0\)

请注意,在算法7.1的步骤3中,权重 \(K_{i0}\) 对每个 \(x_0\) 值都会不同。换句话说,为了在新点处获得局部回归拟合,我们需要通过针对一组新权重最小化 式 7.10 来拟合新的加权最小二乘回归模型。局部回归有时被称为基于内存的过程(memory-based procedure),因为像最近邻一样,每次我们希望计算预测时都需要所有训练数据。

7.7.2 局部回归的Python实现

让我们使用 statsmodelslowess (Locally Weighted Scatterplot Smoothing) 函数来平滑股价。这是技术分析中这类平滑方法的统计学基础。

图 7.13statsmodels 的 LOWESS 在同一历史价格样本上比较 1%、5% 和 20% 三种窗口比例。较小窗口通常更贴近局部波动,较大窗口通常更平滑;哪一个具有预测价值仍需时间验证,LOWESS 也不等同于移动平均。

from statsmodels.nonparametric.smoothers_lowess import lowess  # LOWESS 局部加权回归

# ── 对比三种不同窗口比例(frac)的局部回归效果 ──
# frac: 每次局部拟合使用的数据占总量的比例
local_regression_fractions = [0.01, 0.05, 0.2]  # 三种局部窗口比例:1%、5%、20%

plt.figure(figsize=(14, 8))  # 创建 14×8 英寸画布
plt.plot(haikang_data['trade_date'], haikang_data['Price'], 'k-', alpha=0.2, linewidth=1, label='观测股价')  # 绘制原始股价走势

for frac in local_regression_fractions:  # 遍历三种窗口比例
    # lowess(y, x, frac): 对 (x, y) 做局部加权回归
    # 返回 shape=(n,2) 的数组,第 0 列为排序后的 x,第 1 列为拟合值
    fitted_local_regression = lowess(haikang_data['Price'].values, haikang_data['Time'].values, frac=frac)  # 执行 LOWESS 平滑
    fitted_local_regression_dates = pd.to_datetime(np.interp(fitted_local_regression[:, 0], time_features.flatten(), trade_date_numeric))  # 将 LOWESS 返回的数值时间索引映射回真实交易日期

    # 标签中附带实际使用的交易日窗口大小
    label_text = f'frac={frac} (Window ≈ {int(frac*len(haikang_data))} days)'  # 构造图例文字
    plt.plot(fitted_local_regression_dates, fitted_local_regression[:, 1], linewidth=2, label=label_text)  # 在真实交易日期横轴上绘制局部回归拟合线

plt.xlabel('交易日期', fontsize=14)  # x 轴标签
plt.ylabel('股价 (元)', fontsize=14)  # y 轴标签
plt.title('不同 Span 值的局部线性回归 (Lowess)', fontsize=16)  # 图表标题
plt.legend(fontsize=12)  # 添加图例
plt.grid(True, alpha=0.3)  # 半透明网格线
plt.show()  # 渲染并显示图表
海康威视观测价格灰线叠加frac为0.01、0.05和0.20的三条LOWESS平滑曲线。
图 7.13: 三个窗口比例下的LOWESS历史价格平滑
列表 7.12: 局部回归不同 frac 参数的残差平方和比较
from statsmodels.nonparametric.smoothers_lowess import lowess  # LOWESS 局部加权回归
import numpy as np  # 数值计算库

# 候选的 frac 参数列表(实际生产中可更细粒度搜索,如 np.arange(0.01, 0.30, 0.01))
candidate_fractions = [0.01, 0.02, 0.05, 0.08, 0.10, 0.15, 0.20]  # 7 个候选窗口比例
rss_results = []  # 存储每个 frac 对应的残差平方和

time_values = haikang_data['Time'].values  # 提取时间索引数组
price_values = haikang_data['Price'].values  # 提取股价数组

准备工作就绪后,下面我们遍历所有候选的窗口比例 frac,对每个值分别执行 LOWESS 拟合,并计算其在训练集上的残差平方和(RSS)。RSS 越小意味着拟合越贴近实际股价,但过小的 frac 也可能导致过拟合。

列表 7.13: 遍历候选 frac 计算残差平方和
for frac_candidate in candidate_fractions:  # 遍历候选窗口比例
    # 执行 LOWESS 拟合,返回 (x_sorted, y_fitted) 的 n×2 数组
    lowess_fit = lowess(price_values, time_values, frac=frac_candidate)  # 局部回归拟合
    fitted_prices = lowess_fit[:, 1]  # 提取拟合值列
    residual_sum_squares = np.sum((price_values - fitted_prices) ** 2)  # 计算残差平方和 RSS
    rss_results.append({'frac': frac_candidate, 'RSS': residual_sum_squares,  # 将当前模型的误差存入结果列表
                        'Window': int(frac_candidate * len(haikang_data))})  # 记录结果

rss_comparison_df = pd.DataFrame(rss_results)  # 将结果转为 DataFrame
print(rss_comparison_df.to_string(index=False))  # 打印比较表
 frac          RSS  Window
 0.01  3382.680941      37
 0.02  6918.735934      75
 0.05 14696.373851     189
 0.08 24599.718525     303
 0.10 32919.581348     378
 0.15 58052.363062     568
 0.20 81108.332640     757

从当次表格读取各 frac 的训练 RSS。较小窗口通常更贴近训练观测,因此训练 RSS 偏低是机械现象,不能据此选择部署窗口或宣称泛化更好。

图 7.14 将各候选窗口比例对应的训练 RSS绘制成折线,并标出训练 RSS 最低点;该点不是样本外最优窗口。

plt.figure(figsize=(10, 5))  # 创建画布
plt.plot(rss_comparison_df['frac'], rss_comparison_df['RSS'], 'o-', linewidth=2, markersize=8)  # 绘制 RSS 曲线
plt.xlabel('frac (窗口比例)', fontsize=14)  # x 轴标签
plt.ylabel('残差平方和 (RSS)', fontsize=14)  # y 轴标签
plt.title('局部回归 frac 参数选择', fontsize=16)  # 图标题
plt.grid(True, alpha=0.3)  # 添加网格线

# 标注最优 frac
best_frac_idx = rss_comparison_df['RSS'].idxmin()  # 找到 RSS 最小的索引
best_frac_value = rss_comparison_df.loc[best_frac_idx, 'frac']  # 对应的最优 frac
best_rss_value = rss_comparison_df.loc[best_frac_idx, 'RSS']  # 对应的最小 RSS
plt.annotate(f'训练RSS最低 frac={best_frac_value}', xy=(best_frac_value, best_rss_value),
             xytext=(best_frac_value + 0.03, best_rss_value * 1.3),  # 设置标注文字的偏移位置
             arrowprops=dict(arrowstyle='->', color='red'), fontsize=12, color='red')  # 箭头标注
plt.show()  # 渲染并显示图表
print(f'\n训练RSS最低的候选 frac = {best_frac_value};不得据此选择预测窗口')
折线比较多个LOESS带宽frac对应的残差平方和,标示误差最低候选。
图 7.14: 不同 frac 参数下的残差平方和 (RSS) 对比

训练RSS最低的候选 frac = 0.01;不得据此选择预测窗口

图中标出的只是训练 RSS 最小的候选值,不是推荐窗口。用于预测时必须把 frac 的选择放进训练期内部,并用后续时间块评估;最终外层测试块不得参与窗口选择。

Tip: 局部回归的特点

局部回归有一些独特的特点:

  1. 局部适应性:能够适应数据的局部特征,特别适合在不同区域有不同行为的数据

  2. 直观性:方法直观易懂——在每个点附近拟合局部模型

  3. 灵活性控制:通过span参数控制灵活性:

    • 小span:更局部、更波动的拟合
    • 大span:更全局、更平滑的拟合
  4. 计算成本:需要存储所有训练数据,预测时需要重新计算局部模型

  5. 高维问题:在多维情况下表现不佳(维度灾难)

适用场景

  • 探索性数据分析
  • 数据在不同区域有不同的行为模式
  • 需要直观理解数据结构
  • 低维问题(通常 \(p < 4\)

不适用场景

  • 高维数据
  • 需要快速预测的在线应用
  • 数据量非常大的情况

7.8 广义加性模型 (Generalized Additive Models)

小节 7.2小节 7.7 中,我们提出了几种基于单个预测变量 \(X\) 灵活预测响应变量 \(Y\) 的方法。这些方法可以看作是简单线性回归的扩展。在这里,我们探索基于几个预测变量 \(X_1, \ldots, X_p\) 灵活预测 \(Y\) 的问题。这相当于多重线性回归的扩展。

广义加性模型(GAMs)提供了一个扩展标准线性模型的一般框架,允许每个变量的非线性函数,同时保持可加性;这一经典框架由 Hastie 和 Tibshirani (1986年) 系统提出。就像线性模型一样,GAMs可以应用于定量和定性响应。我们首先在下文回归问题的GAMs小节检查定量响应的GAM,然后在分类问题的GAMs小节检查定性响应。

7.8.1 回归问题的GAMs

扩展多重线性回归模型

\[ y_i = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \cdots + \beta_p x_{ip} + \epsilon_i \]

以允许每个特征和响应变量之间的非线性关系的一种自然方法是将每个线性分量 \(\beta_j x_{ij}\) 替换为(平滑)非线性函数 \(f_j(x_{ij})\)。我们将模型写为:

\[ \begin{aligned} y_i &= \beta_0 + \sum_{j=1}^p f_j(x_{ij}) + \epsilon_i \\ &= \beta_0 + f_1(x_{i1}) + f_2(x_{i2}) + \cdots + f_p(x_{ip}) + \epsilon_i \end{aligned} \tag{7.11}\]

这是GAM的一个例子。它被称为加性模型(additive model),因为我们为每个 \(X_j\) 计算一个单独的 \(f_j\),然后将它们的所有贡献加在一起。

小节 7.2小节 7.7 中,我们讨论了许多拟合单变量函数的方法。GAMs的美妙之处在于我们可以使用这些方法作为拟合加性模型的构建块。事实上,对于我们在本章中看到的大多数方法,这可以相当容易地完成。

以自然样条为例,考虑使用自然样条对year和age、将education作为定性预测变量来拟合模型:

\[ \text{wage} = \beta_0 + f_1(\text{year}) + f_2(\text{age}) + f_3(\text{education}) + \epsilon \tag{7.12}\]

这里year和age是定量变量,而变量education是定性的,有五个水平:<HS, HS, <Coll, Coll, >Coll,指的是个人完成的高中或大学教育的数量。我们使用自然样条拟合前两个函数。我们通过 章节 3 中介绍的虚拟变量方法,为每个水平拟合一个单独的常数来拟合第三个函数。

使用最小二乘拟合模型 式 7.12 很容易,因为如 小节 7.5 所讨论,自然样条可以使用适当选择的基函数来构造。因此,整个模型只是对样条基变量和虚拟变量的大型回归,所有这些都打包在一个大型回归矩阵中。

7.8.2 GAMs的优缺点

在继续之前,让我们总结GAM的优点和局限性:

优点

  • GAMs允许我们对每个 \(X_j\) 拟合非线性 \(f_j\),因此我们可以自动建模标准线性回归会错过的非线性关系。这意味着我们不需要手动尝试对每个变量进行许多不同的变换
  • 非线性拟合可能对响应变量 \(Y\) 做出更准确的预测
  • 由于模型是可加的,我们可以在保持所有其他变量固定的同时检查每个 \(X_j\)\(Y\) 的单独影响
  • 变量 \(X_j\) 的函数 \(f_j\) 的平滑度可以通过自由度来总结

缺点

  • GAMs的主要限制是模型被限制为可加的。对于许多变量,可能会错过重要的交互作用。然而,像线性回归一样,我们可以通过包含形式为 \(X_j \times X_k\) 的额外预测变量手动将交互项添加到GAM模型中。此外,我们可以添加形式为 \(f_{jk}(X_j, X_k)\) 的低维交互函数到模型中;此类项可以使用二维平滑器(如局部回归)或二维样条来拟合(此处不涵盖)

对于完全一般的模型,我们必须寻找更灵活的方法,如 章节 8 中描述的随机森林和提升。GAMs在线性和完全非参数模型之间提供了有用的折衷。

7.8.3 GAMs的Python实现

虽然Python中有专门的 pygam 库,但在 scikit-learn 中,我们可以通过 SplineTransformer 与线性模型的组合,非常优雅地构建广义加性模型。这种方法利用了线性模型的”可加性”和样条的”非线性基函数”,本质上就是GAM。

我们将构建如下金融因子模型: \[ \text{Future Return} = \beta_0 + f_1(\text{Momentum}) + f_2(\text{Volatility}) + f_3(\text{Month}) + \epsilon \]

下面以过去 20 日动量、过去 20 日波动率和月份为预测变量,以未来 5 日收益率为目标,演示加性样条管道。SplineTransformer 处理两个连续变量,OneHotEncoder 处理月份,Ridge 约束基函数系数。月份在这里是候选日历特征,不预设存在日历效应;预测价值由带隔离带的未来测试期决定。

列表 7.14: 使用Scikit-Learn构建GAM (Haikang)
# ── 导入 sklearn 模块 ──
from sklearn.pipeline import make_pipeline          # 串联"预处理+模型"的流水线工具
from sklearn.compose import ColumnTransformer        # 按列名分别施加不同的变换
from sklearn.preprocessing import SplineTransformer, OneHotEncoder, StandardScaler  # 样条变换、独热编码、标准化
from sklearn.linear_model import Ridge               # 岭回归——自带 L2 正则化防止过拟合

# ── 1. 构造量化因子 ──
# 动量 (Momentum): 过去 20 个交易日的累计涨跌幅,反映中短期趋势惯性
haikang_data['Momentum'] = haikang_data['Price'].pct_change(20)  # 计算20日动量指标
# 波动率 (Volatility): 过去 20 个交易日"日收益率"的标准差,衡量价格剧烈程度
haikang_data['Volatility'] = haikang_data['Price'].pct_change().rolling(20).std()  # 计算百分比变化(收益率)
# 月份 (Month): 1~12 的整数,作为不预设效应存在的候选类别特征
haikang_data['Month'] = haikang_data['trade_date'].dt.month  # 提取交易日期中的月份信息
# 目标变量 (Future Return): 未来 5 个交易日的收益率(shift(-5) 表示向前看 5 天)
haikang_data['Future_Return'] = haikang_data['Price'].pct_change(5).shift(-5)  # 计算未来5日收益率作为预测目标
haikang_data['Target_End_Date'] = haikang_data['trade_date'].shift(-5)  # 保存标签窗口终点用于切分审计

# 删除因滚动窗口 / 前瞻产生的缺失值
complete_analysis_data = haikang_data.dropna().copy()  # 创建数据副本避免修改原始数据

# 准备特征矩阵 X 和目标向量 y,并按时间先后切分
gam_features = complete_analysis_data[['Momentum', 'Volatility', 'Month']]  # 提取GAM模型的输入特征
gam_target_returns = complete_analysis_data['Future_Return']  # 提取目标变量
gam_split_point = int(len(complete_analysis_data) * 0.8)  # 用前80%训练、后20%测试
gam_purge_days = 5  # 未来5日标签要求切分边界隔离5个交易日
gam_train_end = gam_split_point - gam_purge_days  # 排除标签窗口跨越测试起点的训练行
gam_features_train = gam_features.iloc[:gam_train_end]  # 仅使用边界前且标签不重叠的样本训练
gam_features_test = gam_features.iloc[gam_split_point:]  # 保留较晚样本测试
gam_target_train = gam_target_returns.iloc[:gam_train_end]  # 训练期连续目标
gam_target_test = gam_target_returns.iloc[gam_split_point:]  # 测试期连续目标
assert complete_analysis_data.iloc[gam_train_end - 1]['Target_End_Date'] < complete_analysis_data.iloc[gam_split_point]['trade_date']  # 审计标签窗口不跨界
# ── 2. 构建预处理管道(GAM 的核心思想)──
# 对连续变量 → SplineTransformer: 自动生成 B 样条基函数,引入非线性
# 对分类变量 → OneHotEncoder: 月份转为 11 个哑变量(drop='first' 避免共线性)
gam_preprocessor = ColumnTransformer([  # 构建特征预处理器:对数值和分类特征分别变换
    ('splines', SplineTransformer(n_knots=5, degree=3), ['Momentum', 'Volatility']),  # 对数值特征应用样条基函数变换
    ('categorical', OneHotEncoder(drop='first', handle_unknown='ignore'), ['Month'])  # 测试期新类别按全零编码处理
])  # 完成构建

# ── 3. 组合成完整的 GAM 流水线 ──
# Ridge(alpha=1.0): 岭回归的正则化类似于平滑惩罚,防止样条基过拟合
ridge_gam_pipeline = make_pipeline(gam_preprocessor, Ridge(alpha=1.0))  # 初始化岭回归模型

# 拟合模型并只在未来测试期评估
from sklearn.metrics import mean_squared_error, r2_score  # 回归测试指标
ridge_gam_pipeline.fit(gam_features_train, gam_target_train)  # 仅拟合较早样本
gam_test_prediction = ridge_gam_pipeline.predict(gam_features_test)  # 预测未来测试期
gam_baseline_prediction = np.repeat(gam_target_train.mean(), len(gam_target_test))  # 训练均值基线
print(f'GAM 测试期 R2: {r2_score(gam_target_test, gam_test_prediction):.4f}')  # 输出样本外R2
print(f'GAM 测试期 MSE: {mean_squared_error(gam_target_test, gam_test_prediction):.6f}')  # 输出样本外MSE
print(f'均值基线 MSE: {mean_squared_error(gam_target_test, gam_baseline_prediction):.6f}')  # 输出基线MSE
GAM 测试期 R2: -0.0647
GAM 测试期 MSE: 0.001964
均值基线 MSE: 0.001863

这里的结论必须来自按时间留出的未来测试期,而不是训练集拟合优度。测试期 \(R^2\)、MSE 与训练期均值基线共同回答“模型是否改善未来预测”;即使统计误差略有改善,也不能自动推出可交易价值,还需计入换手、冲击成本和信号稳定性。

图 7.15 用于审计已拟合 GAM 的预测形状。三幅图分别对应动量、波动率和月份;曲线方向、弯曲程度与月份差异都必须从当次图形读取。它们不剥离混杂,也不证明作用机制、日历异象或可交易收益。

from sklearn.inspection import PartialDependenceDisplay  # 偏依赖图工具

# ── 绘制偏依赖图——揭示每个因子对未来收益的"独立边际效应" ──
fig, ax = plt.subplots(figsize=(14, 5))  # 创建 14×5 英寸画布

# 偏依赖图的计算参数
partial_dependence_params = {  # 定义参数字典
    'subsample': 50,        # 随机抽样 50 条训练样本加速计算(ICE 曲线数量)
    'n_jobs': 2,            # 并行 2 核
    'grid_resolution': 20,  # 特征取值离散化为 20 个网格点
    'random_state': 0,      # 固定随机种子保证可复现
}  # 执行数据处理操作

# 关注的三个特征名称列表
selected_gam_features = ['Momentum', 'Volatility', 'Month']  # 提取GAM模型的输入特征

# 一次性生成 3 张偏依赖图(average = 平均 ICE 曲线)
partial_dependence_display = PartialDependenceDisplay.from_estimator(  # 绘制偏依赖图展示各特征的边际效应
    ridge_gam_pipeline,      # 已拟合的 GAM 管道
    gam_features_train,      # 只用训练期经验分布确定网格并做边际平均
    selected_gam_features,   # 要展示的特征
    kind='average',          # 只画均值线,不画每条 ICE 曲线
    ax=ax,  # 定义ax变量
    **partial_dependence_params,  # 展开参数字典传入函数
)  # 完成构建

# 添加总标题,y=1.05 将标题上移避免被子图遮挡
# 执行数据处理操作
partial_dependence_display.figure_.suptitle('因子对未来收益的偏效应 (Partial Dependence)', fontsize=16, y=1.05)
plt.subplots_adjust(top=0.9)  # 为总标题留出空间
plt.show()  # 渲染并显示图表
三个面板分别显示动量、波动率和月份在训练支持域内的平均预测偏效应。
图 7.15: GAM模型的偏依赖图 (因子效应)

偏依赖曲线把某一特征设为网格值,并对训练期其余特征的联合经验分布上的预测取平均。它描述的是已拟合模型的边际预测形状,不是因果效应;若特征高度相关,网格还可能包含现实中稀少的组合,因此应结合样本密度和未来测试误差谨慎解读。

这种 SplineTransformer + LinearModel 的方法非常强大。它不仅保持了模型的可解释性(每个特征对预测的贡献是可加的),而且利用了 scikit-learn 强大的生态系统(如交叉验证、管道、评估指标)。

7.8.4 分类问题的GAMs

GAMs也可以用于 \(Y\) 是定性的情况。例如,我们预测明天股价是涨还是跌

\[ \log\frac{p(X)}{1-p(X)} = \beta_0 + f_1(\text{Momentum}) + f_2(\text{Volatility}) + f_3(\text{Month}) \]

这可以通过将线性回归器替换为逻辑回归器(LogisticRegression)来实现。

列表 7.15: 逻辑回归GAM (预测股价涨跌)
from sklearn.linear_model import LogisticRegression  # 逻辑回归分类器

# ── 构造二分类目标:未来 5 日收益率 > 0 → 上涨(1);否则 → 下跌(0) ──
# 转换数据类型
complete_analysis_data['Direction'] = (complete_analysis_data['Future_Return'] > 0).astype(int)
target_direction = complete_analysis_data['Direction']  # 提取目标变量
direction_train = target_direction.iloc[:gam_train_end]  # 使用同一5日隔离带后的较早训练期标签
direction_test = target_direction.iloc[gam_split_point:]  # 较晚测试期标签

# ── 构建逻辑回归 GAM 管道 ──
# 复用与回归 GAM 完全相同的预处理器(SplineTransformer + OneHotEncoder)
# 仅将末端的 Ridge 替换为 LogisticRegression
# penalty='l2', C=1.0: L2 正则化强度的倒数——C 越小正则越强
# 初始化逻辑回归模型
logistic_gam_pipeline = make_pipeline(gam_preprocessor, LogisticRegression(penalty='l2', C=1.0))

# 拟合较早样本,并在未来测试期与多数类基线比较
from sklearn.metrics import accuracy_score, brier_score_loss, roc_auc_score  # 分类、排序与概率校准指标
logistic_gam_pipeline.fit(gam_features_train, direction_train)  # 只拟合训练期
direction_probability = logistic_gam_pipeline.predict_proba(gam_features_test)[:, 1]  # 测试期上涨概率
direction_prediction = (direction_probability >= 0.5).astype(int)  # 按0.5阈值分类
majority_class = int(direction_train.mean() >= 0.5)  # 训练期多数类
majority_prediction = np.repeat(majority_class, len(direction_test))  # 多数类测试基线
print(f'逻辑GAM测试准确率: {accuracy_score(direction_test, direction_prediction):.4f}')  # 测试准确率
print(f'多数类基线准确率: {accuracy_score(direction_test, majority_prediction):.4f}')  # 基线准确率
print(f'逻辑GAM测试Brier分数: {brier_score_loss(direction_test, direction_probability):.4f}')  # 概率误差
逻辑GAM测试准确率: 0.4861
多数类基线准确率: 0.4847
逻辑GAM测试Brier分数: 0.2491

分类结论以未来测试期为准,并同时报告训练期多数类基线与 Brier 分数。准确率超过50%并不必然超过多数类基线,更不等于可交易价值;后者还要求稳定的概率校准、阈值选择和交易成本检验。

模型训练完成后,图 7.16 展示各特征与模型预测上涨概率之间的边际关系。计算时把目标特征设为网格值,并在训练样本中其余特征的联合经验分布上取平均,而不是简单地把其他变量固定在均值。

# ── 逻辑回归 GAM 的偏依赖图(概率尺度)──
fig, ax = plt.subplots(figsize=(14, 5))  # 创建 14×5 英寸画布

partial_dependence_display = PartialDependenceDisplay.from_estimator(  # 绘制偏依赖图展示各特征的边际效应
    logistic_gam_pipeline,      # 已拟合的逻辑回归 GAM 管道
    gam_features_train,         # 训练数据(用于计算网格和边际平均)
    selected_gam_features,      # 要展示的三个特征
    kind='average',             # 只画均值线,不画每条 ICE 曲线
    response_method='predict_proba',   # 输出预测"上涨概率"而非 log-odds
    ax=ax,                      # 绑定到当前画布
    **partial_dependence_params, # 复用偏依赖图参数
)  # 完成构建

# 总标题——展示各因子如何影响上涨概率
# 执行数据处理操作
partial_dependence_display.figure_.suptitle('因子对上涨概率的偏效应 (P(Up))', fontsize=16, y=1.05)
plt.subplots_adjust(top=0.9)  # 为总标题留出空间

# 在每个子图中添加 0.5 参考线——涨跌分界
for axis in partial_dependence_display.axes_[0]:  # 遍历所有子图坐标轴
    axis.axhline(0.5, color='gray', linestyle='--', alpha=0.5)  # 50% 分界线

plt.show()  # 渲染并显示图表
三个面板显示动量、波动率和月份与模型预测上涨概率的平均偏依赖及0.5参考线。
图 7.16: 逻辑回归GAM的偏依赖图 (涨跌概率)

逻辑回归 GAM 的偏依赖图展示各特征与模型预测上涨概率之间的边际预测关联,灰色虚线是 0.5 决策参考。倒 U 形、单调抑制或月份差异都只能在图中确实出现时描述,并须结合测试期 Brier、基线和样本密度判断稳定性。

7.8.5 选读:pygam 实现与完整预测器偏依赖

GAMs也可以用于 \(Y\) 是定性的情况。为简单起见,这里我们假设 \(Y\) 取值为0或1,令 \(p(X) = \Pr(Y=1|X)\) 为响应等于1的条件概率(给定预测变量)。回想逻辑回归模型@eq-multiple-logistic:

\[ \log\frac{p(X)}{1-p(X)} = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \cdots + \beta_p X_p \tag{7.13}\]

左边是 \(P(Y=1|X)\)\(P(Y=0|X)\) 的比率的对数,式 7.13 将其表示为预测变量的线性函数。扩展 式 7.13 以允许非线性关系的一种自然方法是使用模型:

\[ \log\frac{p(X)}{1-p(X)} = \beta_0 + f_1(X_1) + f_2(X_2) + \cdots + f_p(X_p) \tag{7.14}\]

式 7.14 是一个逻辑回归GAM。它具有前几节中讨论的定量响应的所有相同的优缺点。

例如,在股票分析中,我们可能只关心明天股价是涨还是跌,而不是具体的收益率数值。

\(p(X) = \Pr(Y=1|X)\) 为股价上涨 (\(Y=1\)) 的概率。我们可以构建逻辑回归GAM:

\[ \log\frac{p(X)}{1-p(X)} = \beta_0 + f_1(\text{Momentum}) + f_2(\text{Volatility}) + f_3(\text{Month}) \]

列表 7.16: 逻辑回归GAM (预测股价涨跌)
from pygam import LogisticGAM, s, f  # pygam 的逻辑回归 GAM;s: 平滑项, f: 因子项

# pygam 要求 NumPy 数组;沿用相同的时间切分
gam_features_train_array = gam_features_train.to_numpy()  # 较早训练期特征
gam_features_test_array = gam_features_test.to_numpy()  # 较晚测试期特征
direction_train_array = direction_train.to_numpy()  # 训练期标签
direction_test_array = direction_test.to_numpy()  # 测试期标签

# ── 在较早训练期内另留按时间靠后的验证段,并按5日标签跨度purge ──
# s(0): 对第 0 列 Momentum 施加平滑样条
# s(1): 对第 1 列 Volatility 施加平滑样条
# f(2): 对第 2 列 Month 施加因子项(类别变量)
gam_validation_start = int(len(gam_features_train_array) * 0.8)  # 较后20%训练期作为候选λ验证段
gam_tuning_end = gam_validation_start - gam_purge_days  # 五日标签窗口不得伸入验证段
assert gam_tuning_end > 0 and gam_validation_start < len(gam_features_train_array)  # 保证切分非空
assert complete_analysis_data.iloc[gam_tuning_end - 1]['Target_End_Date'] < complete_analysis_data.iloc[gam_validation_start]['trade_date']  # 审计标签窗口
candidate_lambdas = np.logspace(-2, 3, 8)  # 预先声明有限候选平滑参数
validation_auc_by_lambda = []  # 保存每个候选在较后验证段的AUC
for candidate_lambda in candidate_lambdas:  # 显式遍历,不调用GCV/UBRE gridsearch
    candidate_gam = LogisticGAM(s(0, lam=candidate_lambda) + s(1, lam=candidate_lambda) + f(2, lam=candidate_lambda))  # 固定候选λ
    candidate_gam.fit(gam_features_train_array[:gam_tuning_end], direction_train_array[:gam_tuning_end])  # 仅拟合较早调参段
    validation_probability = candidate_gam.predict_proba(gam_features_train_array[gam_validation_start:])  # 预测较后验证段
    validation_auc_by_lambda.append(roc_auc_score(direction_train_array[gam_validation_start:], validation_probability))  # 记录排序误差
selected_lambda = candidate_lambdas[int(np.argmax(validation_auc_by_lambda))]  # 只按验证段选择一次
print('候选lambda及时间验证AUC:', dict(zip(candidate_lambdas, validation_auc_by_lambda)))  # 显示真实准则与结果
print('调参段/验证段日期:', complete_analysis_data.iloc[0]['trade_date'], complete_analysis_data.iloc[gam_tuning_end - 1]['trade_date'], complete_analysis_data.iloc[gam_validation_start]['trade_date'], complete_analysis_data.iloc[gam_train_end - 1]['trade_date'])  # 输出折边界
linear_gam_model = LogisticGAM(s(0, lam=selected_lambda) + s(1, lam=selected_lambda) + f(2, lam=selected_lambda))  # 锁定验证选出的λ
linear_gam_model.fit(gam_features_train_array, direction_train_array)  # 用完整训练期重拟合锁定模型
pygam_test_probability = linear_gam_model.predict_proba(gam_features_test_array)  # 未来测试概率
print(f'pygam 测试 AUC: {roc_auc_score(direction_test_array, pygam_test_probability):.4f}')  # 测试AUC
print(f'pygam 测试 Brier: {brier_score_loss(direction_test_array, pygam_test_probability):.4f}')  # 测试Brier
候选lambda及时间验证AUC: {0.01: 0.5514296775620817, 0.05179474679231213: 0.5495381892594435, 0.2682695795279726: 0.5413306786129086, 1.3894954943731375: 0.5293180686908909, 7.196856730011521: 0.5240639345168961, 37.27593720314942: 0.5144626956473647, 193.06977288832496: 0.4959460206846966, 1000.0: 0.48112383164647976}
调参段/验证段日期: 2010-06-30 00:00:00 2020-05-13 00:00:00 2020-05-21 00:00:00 2022-11-10 00:00:00
pygam 测试 AUC: 0.4963
pygam 测试 Brier: 0.2712

这一实现只是对前一节的补充:模型仍须按时间外推评估。下面直接用完整线性预测器计算概率偏依赖;不能把单个项的 log-odds 贡献单独做 sigmoid,因为截距和其他项缺失时得到的并不是上涨概率。

# ── 在训练期联合经验分布上计算完整预测器的概率偏依赖 ──
# 三个特征的标签与中文标题
labels = ['Momentum (20d)', 'Volatility (20d)', 'Month']  # 横轴标签
titles = ['动量对上涨概率的影响',      # 第 1 张子图标题
          '波动率对上涨概率的影响',    # 第 2 张子图标题
          '季节性对上涨概率的影响']    # 第 3 张子图标题

partial_data = []  # 存储每个特征的网格与平均预测概率
for i in range(3):  # 遍历 3 个特征项
    if i == 2:  # 月份使用训练期实际类别
        feature_grid = np.unique(gam_features_train_array[:, i])  # 月份网格
    else:  # 连续变量限制在训练期1%至99%分位数
        grid_limits = np.quantile(gam_features_train_array[:, i], [0.01, 0.99])  # 稳健范围
        feature_grid = np.linspace(*grid_limits, 50)  # 连续网格
    average_probabilities = []  # 保存每个网格值的平均概率
    for grid_value in feature_grid:  # 逐个网格值计算完整预测器
        counterfactual_features = gam_features_train_array.copy()  # 保留其他变量联合分布
        counterfactual_features[:, i] = grid_value  # 只替换目标特征
        average_probabilities.append(linear_gam_model.predict_proba(counterfactual_features).mean())  # 平均概率
    partial_data.append((feature_grid, np.asarray(average_probabilities)))  # 保存PDP

基于上面计算好的偏依赖数据,图 7.17 分三个子图展示动量、波动率和月份对应的平均上涨概率:

# ── 绘制 pygam LogisticGAM 的偏依赖图(概率尺度)──
fig, axes = plt.subplots(1, 3, figsize=(18, 5))  # 1 行 3 列子图

# ── 将计算好的偏依赖数据绘制到子图 ──
for i, ax in enumerate(axes):  # 遍历 3 个子图坐标轴
    x_vals, pdep_prob = partial_data[i]  # 解包第 i 个特征的数据
    if i == 2:  # Month 是离散分类变量 → 用阶梯线 + 误差条展示
        ax.plot(x_vals, pdep_prob, 'b-', linewidth=3, drawstyle='steps-mid')  # 阶梯折线
        ax.scatter(x_vals, pdep_prob, color='b')  # 标出实际月份类别
    else:  # Momentum / Volatility 是连续变量 → 实线
        ax.plot(x_vals, pdep_prob, 'b-', linewidth=3)  # 主效应曲线

    ax.set_xlabel(labels[i], fontsize=12)       # 设置横轴标签
    ax.set_ylabel('P(上涨)', fontsize=12)       # 纵轴: 上涨概率
    ax.set_title(titles[i], fontsize=14)        # 子图标题
    ax.grid(True, alpha=0.3)                    # 添加半透明网格线
    ax.set_ylim([0, 1])                         # 概率范围 [0, 1]
    ax.axhline(0.5, color='gray', linestyle='--')  # 50% 参考线——涨跌分界

plt.tight_layout()  # 自动调整子图间距
plt.show()  # 渲染并显示图表
pyGAM完整线性预测器下,动量、波动率和月份对应的平均上涨概率曲线。
图 7.17: 逻辑回归GAM的偏依赖图 (涨跌概率)

三张图给出完整模型在训练期经验分布上的平均预测概率。曲线形状仅刻画模型所学到的预测关联;是否稳定、能否推广,仍由后段测试期指标决定,不能凭图形宣称显著性、日历异象或经济价值。

偏依赖图揭示了各特征的非线性效应模式。接下来,我们进一步通过方差分析(ANOVA / Deviance Analysis)来严格检验 Month 变量应当以何种形式进入模型——是作为因子项、线性项,还是非线性平滑项。这种嵌套模型逐步比较的方法是 GAM 模型选择中的标准流程。

列表 7.17: GAM模型的方差分析
from pygam import l  # l: 线性项(用于方差分析比较)

# ── 三个嵌套模型逐步比较(ANOVA / 偏差分析)──
# 模型 1:Month 作为因子项 f(2)——仅月份阶梯效应
gam_model_factor_month = LogisticGAM(s(0) + s(1) + f(2))  # 前两个特征用平滑项,月份为因子项
gam_model_factor_month.fit(gam_features_train_array, direction_train_array)  # 只拟合训练期

# 模型 2:Month 作为线性项 l(2)——假设月份有单调线性趋势
gam_model_linear_dist = LogisticGAM(s(0) + s(1) + l(2))  # 构建GAM模型:前两个特征用平滑项,第三个为线性项
gam_model_linear_dist.fit(gam_features_train_array, direction_train_array)  # 只拟合训练期

# 模型 3:Month 作为平滑项 s(2)——允许月份有非线性曲面效应
gam_model_smooth_dist = LogisticGAM(s(0) + s(1) + s(2))  # 构建GAM模型:三个特征均使用平滑项
gam_model_smooth_dist.fit(gam_features_train_array, direction_train_array)  # 只拟合训练期

# ── 提取各模型的 Deviance(偏差)和有效自由度 ──
# Deviance 越小 → 拟合越好;EDF 越大 → 模型越复杂
deviance_factor_month = float(np.ravel(gam_model_factor_month.statistics_['deviance'])[0])  # 提取标量
# 完成模型/Pipeline构建
deviance_linear_dist = float(np.ravel(gam_model_linear_dist.statistics_['deviance'])[0])
# 完成模型/Pipeline构建
deviance_smooth_dist = float(np.ravel(gam_model_smooth_dist.statistics_['deviance'])[0])

dof_factor_month = float(np.ravel(gam_model_factor_month.statistics_['edof'])[0])  # 有效自由度标量
dof_linear_dist = float(np.ravel(gam_model_linear_dist.statistics_['edof'])[0])  # 完成模型/Pipeline构建
dof_smooth_dist = float(np.ravel(gam_model_smooth_dist.statistics_['edof'])[0])  # 完成模型/Pipeline构建

test_brier_factor = brier_score_loss(direction_test_array, gam_model_factor_month.predict_proba(gam_features_test_array))
test_brier_linear = brier_score_loss(direction_test_array, gam_model_linear_dist.predict_proba(gam_features_test_array))
test_brier_smooth = brier_score_loss(direction_test_array, gam_model_smooth_dist.predict_proba(gam_features_test_array))

接下来只打印训练期偏差与有效自由度作为诊断。因子项、线性项和平滑项并不构成具有正自由度差的简单嵌套序列,不能把偏差差机械代入 F 分布。

# ── 打印模型比较汇总表 ──
print('训练期模型诊断:')  # 打印标题
print('-' * 80)  # 分隔线
print(f'{"模型":<40} {"训练偏差":<15} {"有效自由度":<15} {"测试Brier":<15}')
print('-' * 80)  # 分隔线
print(f'{"月份为因子项":<40} {deviance_factor_month:<15.2f} {dof_factor_month:<15.2f} {test_brier_factor:<15.4f}')
print(f'{"月份为线性项":<40} {deviance_linear_dist:<15.2f} {dof_linear_dist:<15.2f} {test_brier_linear:<15.4f}')
print(f'{"月份为平滑项":<40} {deviance_smooth_dist:<15.2f} {dof_smooth_dist:<15.2f} {test_brier_smooth:<15.4f}')
训练期模型诊断:
--------------------------------------------------------------------------------
模型                                       训练偏差            有效自由度           测试Brier        
--------------------------------------------------------------------------------
月份为因子项                                   4004.77         32.10           0.2643         
月份为线性项                                   4048.02         22.27           0.2639         
月份为平滑项                                   4005.19         31.66           0.2639         

训练期偏差只能用于检查拟合,不能在这里产生有效的 F 检验或样本外排序。月份编码应依据其类别/周期语义预先确定,候选平滑度则应在训练期内部用时间验证选择,再到未来测试期报告对数损失或 Brier 分数。

Tip: GAMs的实践建议

GAMs在实际应用中有以下最佳实践:

  1. 模型选择
    • 从简单的可加模型开始
    • 使用偏依赖图检查每个变量的效应
    • 考虑重要的交互作用
    • 使用AIC、BIC或交叉验证选择模型
  2. 平滑参数选择
    • 独立或近似稳定样本可考虑广义交叉验证(Generalized Cross-Validation, GCV);时间任务优先使用前向验证
    • 通常4-10个自由度足够
    • 对于定性变量,使用虚拟变量
  3. 诊断
    • 检查残差图
    • 检查偏依赖图的合理性
    • 比较GAM与线性模型
    • 进行ANOVA检验
  4. 解释性
    • 使用偏依赖图解释每个变量的效应
    • 注意可加性假设的限制
    • 考虑效应大小和置信区间
  5. 计算效率
    • 对于大数据集,考虑使用回归样条而非平滑样条
    • 使用backfitting算法提高效率
    • 考虑并行计算

GAMs在许多实际应用中是一个很好的折衷方案,特别是在需要平衡灵活性和可解释性的情况下。

7.9 本章小结 (Chapter Summary)

本章我们探讨了超越线性关系的方法,主要内容包括:

  1. 多项式回归:通过添加多项式项扩展线性模型,简单但可能导致边界问题
  2. 阶梯函数:将连续变量分段,简单易解释但可能损失信息
  3. 回归样条:使用分段多项式在节点处平滑连接,灵活且稳定
  4. 平滑样条:通过最小化带平滑性惩罚的目标函数获得样条
  5. 局部回归:在每个点附近拟合局部模型,适应性强
  6. 广义加性模型(GAMs):扩展到多个预测变量,保持可加性和可解释性

方法比较

方法 灵活性 可解释性 计算成本 适用场景
多项式回归 简单非线性关系
阶梯函数 分段常数关系
回归样条 一般非线性关系
平滑样条 中小型数据集
局部回归 探索性分析
GAMs 中-高 多变量非线性关系

实践建议

  • 从简单的多项式回归或样条开始
  • 使用交叉验证选择模型复杂性
  • 检查拟合曲线的合理性,特别是在边界处
  • 对于多个预测变量,优先考虑GAMs
  • 始终与线性基线模型比较

在中国房地产市场的案例中,我们展示了这些方法如何捕捉房价与房龄之间的复杂非线性关系。这些技术同样适用于其他中国金融和经济问题,如股票收益率建模、消费者行为分析、经济增长预测等。

目标—评价证据:O7.1 对应练习 1—5、9—10;O7.2 对应练习 6、8;O7.3—O7.4 对应练习 7 与 M07;O7.5 对应练习 6—8。达标要求概念/理论题至少 70% 分,所有时间断言通过,模型比较显示实际有限指标与支持域,且不预写胜者。

7.10 理论来源与前沿

从线性走向非线性并不是放弃理论,而是把更丰富的函数空间纳入建模。多项式回归、阶梯函数与样条方法都可以理解为‘选择一组基函数’并在参数空间中做线性估计;平滑样条通过在拟合误差与曲线粗糙度之间加入惩罚项,给出一个可控的偏差-方差权衡。GAM 则把多维非线性拆成一维可加平滑项,使得模型在保持一定可解释性的同时具有较强灵活度。

近年来的前沿发展集中在:

  1. 可解释的非线性约束:引入单调性、凸性等形状约束,使模型更符合经济学直觉与业务逻辑。
  2. 核方法与可扩展近似:用核技巧刻画非线性,并用随机特征、Nyström 等方法提升大规模可用性。
  3. 与因果推断结合:在非线性混杂控制、异质性效应估计等问题上,非参数方法与因果识别设计相互补强。

7.11 项目里程碑:M07 非线性模型卡与支持域图

本里程碑落实 小节 6.3,并复用 a-share-drawdown-20d-v1、M02 唯一 manifest/hash 与 M05 审计折。键、标签、label_date、末端 NA 与严格 purge 条件不得重定义。输出:训练事件率基线与 LogisticGAM 的验证 Brier/AUC、候选平滑参数、支持域图和 registry 记录;M07 不访问测试标签。

时点与资源\(\max(\text{train.label_date})<\min(\text{validation.prediction_date})\);偏效应只在训练支持域解释;课堂目标 60 分钟。失败条件:manifest/hash 不匹配、把 GCV/UBRE 称为时间 CV、访问测试键、标签跨界、缺训练事件率基线或固定胜者。量规(20 分):合同/时点 6 分,LogisticGAM 与基线 5 分,支持域 5 分,复现/解释 4 分。

本章闭环:M07 把 LogisticGAM 的平滑度选择限定在训练键内,验证键只产生一次 AUC/Brier 和支持域证据;随后按冻结配方在完整开发键上重拟,以与 M06 完全相同的 registry schema 交给 M09。

7.12 练习

7.12.1 概念题

  1. 多项式回归、阶梯函数与回归样条都能表达非线性。比较三者的连续性、局部性和外推行为,并说明“监管阈值生效后风险率跳变”与“温度—故障率连续弯曲”分别更适合先检验哪类表示;选择仍须由冻结验证证据确认。 [核心|难度:1|时间:8分钟|分值:5|项目:无]

  2. 解释“结点(knot)”在回归样条中的作用。结点数目越多一定越好吗?为什么? [核心|难度:1|时间:8分钟|分值:5|项目:无]

  3. 在平滑样条中,惩罚项 \(\int (f''(t))^2 dt\) 的直觉含义是什么?它如何影响曲线形状? [核心|难度:2|时间:8分钟|分值:5|项目:无]

  4. 广义可加模型(GAM)为什么既“灵活”又“可解释”?它的主要结构性假设是什么? [核心|难度:2|时间:8分钟|分值:5|项目:无]

  5. 局部回归(LOESS/LOWESS)中的带宽(span)扮演什么角色?带宽与偏差-方差之间如何权衡? [核心|难度:2|时间:8分钟|分值:5|项目:无]

7.12.2 应用题

  1. 选取一家长三角上市公司,构造“交易活跃度(如换手率/成交额)与未来波动率”的关系图,并分别用: [拓展|难度:2|时间:45分钟|分值:15|项目:无]

    • 线性回归
    • 三次多项式回归
    • 自然样条(或 B 样条)

进行拟合。比较三种拟合的稳定性与可解释性。

  1. a-share-drawdown-20d-v1 的 M05 训练—验证键上拟合 LogisticGAM,至少包含一个平滑项与一个线性项;报告训练事件率基线和模型的验证 Brier、候选平滑参数、manifest 哈希及严格 purge 断言,禁止访问测试键。 [核心|难度:3|时间:60分钟|分值:20|项目:M07]

    • 模型在冻结验证键上的 AUC 与 Brier score
    • 平滑项的形状(单调、U 型等)与经济含义
  2. 在同一数据上,比较“单一二维平滑(核/LOESS)”与“两个一维平滑相加(GAM)”在泛化误差与解释性上的差异。 [拓展|难度:3|时间:40分钟|分值:15|项目:无]

7.12.3 理论题

  1. 说明平滑样条的估计可以写成线性平滑器 \(\hat y = S_\lambda y\)。解释 \(\mathrm{tr}(S_\lambda)\) 为什么可被视为有效自由度,并描述其随 \(\lambda\) 变化的规律。 [补救|难度:2|时间:18分钟|分值:10|项目:无]

  2. 在局部线性回归中,简述为什么它相比局部常数回归(Nadaraya–Watson)在边界处具有更小的偏差。 [补救|难度:2|时间:15分钟|分值:10|项目:无]

展开完整参考解答与评分键

7.13 练习参考解答

统一评分与验收:每题分值见元数据;概念/理论题按“定义 40%—机制或推导 40%—支持域边界 20%”给分。代码题要求有限 MSE/AUC/Brier,概率在 \([0,1]\),日期断言通过,数值复算容差 \(10^{-8}\)。重新估计测试期样条基、把 GCV/UBRE 称时间 CV、标签窗口跨界或预写模型胜者均判失败。

7.13.1 概念题参考解答

  1. 连续性、局部性与外推:多项式是全局且连续光滑的基,一个局部点会影响整体,边界外推可能不稳定;阶梯函数在预注册切点处允许不连续跳变,适合制度阈值但切点内保持常数;回归样条是局部支撑的连续分段多项式,自然样条边界外线性。监管阈值先检验阶梯表示,连续温度关系先检验样条;最终仍比较同一冻结验证集。

  2. 结点作用:结点让函数在不同区间拥有不同的局部形状(分段多项式),同时通过连续性约束保证整体平滑。结点越多,训练误差通常更小,但方差更大、过拟合风险上升,需要用 CV/信息准则选择复杂度。

  3. 惩罚项直觉\(\int (f''(t))^2dt\) 惩罚“曲率/弯折”,鼓励函数更接近直线。\(\lambda\) 越大,越不允许快速弯折,曲线更光滑;\(\lambda\to 0\) 时更贴合数据。

  4. GAM 的结构假设

    • 灵活:每个变量可以用一维平滑函数 \(f_j(x_j)\) 表达非线性。
    • 可解释:每个 \(f_j\) 可单独画出来解释边际效应。
    • 关键假设:可加性(缺少高阶交互,或需显式加入交互项/张量积平滑)。
  5. 带宽与偏差-方差:带宽大(span 大)意味着用更多邻域点,方差低但偏差高;带宽小意味着更局部,偏差低但方差高。时间任务应使用前向验证;广义交叉验证(GCV)是线性平滑器在稳定、近似独立设定下的快捷准则,不是时间感知切分,AIC 也属于不同的信息准则。

7.13.2 应用题参考解答(流程与模板)

  1. 三种非线性拟合对比: 建议使用时间切分或滚动窗口。关注:边界是否发散、局部是否过度波动。
# ── 练习 6 参考代码:三种非线性方法拟合对比 ──
import pandas as pd  # 导入 pandas 用于数据处理
import numpy as np  # 数值计算库(分位数节点计算等)
import statsmodels.api as sm  # 导入 statsmodels 用于统计建模
from patsy import dmatrix, build_design_matrices           # 生成并复用训练期设计矩阵
from sklearn.linear_model import LinearRegression          # 普通线性回归
from sklearn.preprocessing import PolynomialFeatures       # 多项式特征变换

# ── 1. 读入本地股价数据 ──
import os  # 导入 os 用于跨平台路径拼接
import platform  # 导入 platform 用于严格按操作系统选择数据根目录
# 根据操作系统设置数据根目录路径
DATA_ROOT = os.path.abspath(os.path.expanduser(os.environ['BOOK_DATA_DIR']))  # 从必需环境变量解析数据根目录
assert os.path.isdir(DATA_ROOT), f'BOOK_DATA_DIR 不存在: {DATA_ROOT}'  # 缺少挂载时立即失败
path = os.path.join(DATA_ROOT, 'stock/stock_price_pre_adjusted.h5')  # 拼接本地数据文件路径
stock_analysis_data = pd.read_hdf(  # 只读取目标股票和所需字段
    path, key='data', where='order_book_id == "002415.XSHE"',  # HDF表内筛选
    columns=['close', 'total_turnover']  # 限定分析列
)  # 完成选择性读取

# 若 order_book_id 在索引中则重置为普通列
if 'order_book_id' in stock_analysis_data.index.names:  # 获取索引
    stock_analysis_data = stock_analysis_data.reset_index()  # 重置DataFrame索引
# 统一日期列名为 trade_date,兼容不同数据源的命名差异
# 获取列名列表
if 'date' in stock_analysis_data.columns and 'trade_date' not in stock_analysis_data.columns:
    stock_analysis_data = stock_analysis_data.rename(columns={'date': 'trade_date'})  # 重命名列或索引
stock_analysis_data['trade_date'] = pd.to_datetime(stock_analysis_data['trade_date'])  # 将交易日期统一转换为 datetime,便于后续时间特征提取
# 筛选海康威视 & 按日期排序
# 创建数据副本避免修改原始数据
stock_analysis_data = stock_analysis_data.sort_values('trade_date').copy()  # HDF读取阶段已筛选股票
# ── 2. 构造特征与目标变量 ──
# Activity: 对成交金额取对数——衡量市场交投活跃度
stock_analysis_data['Activity'] = np.log(stock_analysis_data['total_turnover'] + 1)  # 计算自然对数
# Future_Vol: 未来 20 个交易日收益率的标准差(向前看 20 天)
# 计算百分比变化(收益率)
stock_analysis_data['Future_Vol'] = stock_analysis_data['close'].pct_change().rolling(20).std().shift(-20)
stock_analysis_data['Target_End_Date'] = stock_analysis_data['trade_date'].shift(-20)  # 保存20日标签窗口终点
stock_analysis_data = stock_analysis_data.dropna()  # 删除因滚动窗口产生的缺失值

# ── 3. 时间切分训练 / 测试集(前 80% 训练 / 后 20% 测试)──
train_size = int(len(stock_analysis_data) * 0.8)  # 计算元素数量
exercise_purge_days = 20  # 未来20日波动率标签要求20个交易日隔离带
data_train = stock_analysis_data.iloc[:train_size - exercise_purge_days].copy()  # 排除标签窗口跨越测试起点的训练行
data_test = stock_analysis_data.iloc[train_size:].copy()  # 较晚样本作为测试集
assert data_train['Target_End_Date'].max() < data_test['trade_date'].min()  # 审计练习标签窗口不跨界

activity_features_train = data_train[['Activity']].values   # shape=(n_train, 1)
volatility_target_train = data_train['Future_Vol'].values    # shape=(n_train,)
activity_features_test = data_test[['Activity']].values      # shape=(n_test, 1)

接下来,分别建立三种非线性拟合模型并对比测试集表现。

# ── (1) 线性回归 ──
linear_model_obj = LinearRegression().fit(activity_features_train, volatility_target_train)  # 拟合线性回归模型
pred_linear = linear_model_obj.predict(activity_features_test)  # 预测测试集

# ── (2) 三次多项式回归 ──
poly = PolynomialFeatures(degree=3)                          # 生成 x, x², x³ 三列特征
poly_model_obj = LinearRegression().fit(poly.fit_transform(activity_features_train), volatility_target_train)  # 拟合多项式回归
pred_poly = poly_model_obj.predict(poly.transform(activity_features_test))  # 预测测试集

# ── (3) 自然样条 (df=4) ──
# patsy cr(): 自然三次样条基函数——边界外为线性,内部为三次
natural_spline_features_train = dmatrix('cr(x, df=4)', {'x': activity_features_train.flatten()}, return_type='dataframe')  # 训练集样条基
natural_spline_features_test = build_design_matrices(  # 复用训练期节点、边界和列顺序
    [natural_spline_features_train.design_info],  # 训练期设计信息
    {'x': activity_features_test.flatten()}  # 测试期自变量
)[0]  # 取出唯一设计矩阵
spline_model_obj = sm.OLS(volatility_target_train, natural_spline_features_train).fit()  # 拟合 OLS 样条回归
pred_spline = spline_model_obj.predict(natural_spline_features_test)  # 预测测试集

# ── 4. 对比三种方法的测试集 MSE ──
print(f'Linear MSE: {np.mean((data_test["Future_Vol"] - pred_linear)**2):.6f}')  # 线性回归 MSE
print(f'Poly(3) MSE: {np.mean((data_test["Future_Vol"] - pred_poly)**2):.6f}')  # 多项式回归 MSE
print(f'Spline(4) MSE: {np.mean((data_test["Future_Vol"] - pred_spline)**2):.6f}')  # 自然样条 MSE
Linear MSE: 0.000099
Poly(3) MSE: 0.000102
Spline(4) MSE: 0.000100

三种模型都在同一未来测试期计算有限的 MSE;自然样条测试矩阵复用训练期 design_info,避免测试期重新估计节点或发生列错配。哪一种模型更好应以实际打印的测试误差为准,不能预先宣称非线性显著或不显著。

  1. LogisticGAM 预测冻结的 20 日最大回撤事件:从 fresh kernel 读取运行器提供的开发特征/开发标签,复用 M02 manifest/hash 与 M05 持久化折,只比较训练事件率基线和验证指标,不读取行情 HDF、测试特征或测试标签。

fresh-kernel 执行顺序:依次运行 lst-m07-01-manifest-setuplst-m07-02-development-inputslst-m07-03-authoritative-joinlst-m07-04-fold-contractlst-m07-05-fold-membershiplst-m07-06-inner-selectionlst-m07-07-validation-consumerfig-m07-support-domainlst-m07-09-refit-recordlst-m07-10-registry-consumer。前六块只建立并核验输入/选择状态;第七块开始消费冻结选择,最后两块才写授权的候选制品与 registry。

列表 7.18: M07步骤01:导入依赖并核验唯一manifest
import hashlib  # 为 manifest、键、registry 和产物生成内容哈希
import json  # 读写公共 JSON Lines registry
import os  # 读取课程运行器显式提供的隔离入口
import pickle  # 序列化冻结 LogisticGAM 产物
from pathlib import Path  # 解析运行器输入路径
import matplotlib.pyplot as plt  # 绘制真实开发样本的支持域与偏效应
import numpy as np  # 执行有限平滑网格与数值守卫
import pandas as pd  # 读取权限分离的开发制品
from pygam import LogisticGAM, s, l  # 显式区分平滑项与线性项
from sklearn.metrics import brier_score_loss, roc_auc_score  # 记录概率与排序证据
m07_manifest_path = Path(os.environ['BOOK_PROJECT_MANIFEST']).expanduser().resolve()  # 锁定M02唯一manifest
m07_manifest_hash = hashlib.sha256(m07_manifest_path.read_bytes()).hexdigest()  # 计算实际字节哈希
assert m07_manifest_hash == os.environ['BOOK_PROJECT_MANIFEST_SHA256'].lower()  # 拒绝切分漂移
m07_manifest = pd.read_csv(m07_manifest_path)  # 读取冻结键与分段
m07_manifest[['prediction_date', 'label_date']] = m07_manifest[['prediction_date', 'label_date']].apply(pd.to_datetime)  # 统一原点与权威终点日类型
assert not m07_manifest.duplicated(['order_book_id', 'prediction_date']).any()  # 要求公司日键唯一
列表 7.19: M07步骤02:读取权限分离的开发特征与标签
m07_feature_path = Path(os.environ['BOOK_PROJECT_DEVELOPMENT_FEATURES']).resolve()  # 只解析运行器提供的开发特征
m07_label_path = Path(os.environ['BOOK_PROJECT_DEVELOPMENT_LABELS']).resolve()  # 只解析运行器提供的开发标签
m07_features = ['momentum_20', 'volatility_20']  # 冻结 LogisticGAM 特征 schema
m07_key_columns = ['order_book_id', 'prediction_date']  # 冻结跨制品公司日键
m07_feature_frame = pd.read_csv(m07_feature_path, parse_dates=['prediction_date'])  # 读取不含结果代理的开发特征
m07_label_frame = pd.read_csv(m07_label_path, parse_dates=['prediction_date'])  # 读取只含开发事件的标签源
m07_forbidden_features = {'y', 'event', 'label_date', 'close', 'future_min', 'future_min_close', 'split'}  # 禁止未来窗口与保留字段进入特征源
assert set(m07_feature_frame.columns).issuperset(m07_key_columns + m07_features)  # 特征源必须覆盖冻结模型列
assert not m07_forbidden_features.intersection(m07_feature_frame.columns)  # 审计前命名空间不得持有结果代理
assert set(m07_label_frame.columns) == set(m07_key_columns + ['y'])  # 标签源只能携带开发键和事件
assert not m07_feature_frame.duplicated(m07_key_columns).any()  # 开发特征键必须唯一
assert not m07_label_frame.duplicated(m07_key_columns).any()  # 开发标签键必须唯一
m07_feature_frame = m07_feature_frame[m07_key_columns + m07_features].copy()  # 只保留本候选事前登记的特征
assert set(m07_label_frame['y'].dropna().unique()) == {0, 1}  # 标签必须是完整二元开发结果
列表 7.20: M07步骤03:按权威开发键连接并执行purge审计
def m07_key_hash(key_frame):  # 定义与公共 registry 一致的键哈希
    canonical_keys = key_frame[['order_book_id', 'prediction_date']].sort_values(['order_book_id', 'prediction_date']).copy()  # 固定行列顺序
    canonical_keys['prediction_date'] = canonical_keys['prediction_date'].dt.strftime('%Y-%m-%d')  # 固定日期序列化
    return hashlib.sha256(canonical_keys.to_csv(index=False, lineterminator='\n').encode()).hexdigest()  # 返回 SHA-256
m07_development_manifest = m07_manifest[m07_manifest['split'].isin(['train', 'validation'])].copy()  # 仅领取开发键
assert {'order_book_id', 'prediction_date', 'label_date', 'split'}.issubset(m07_manifest.columns)  # manifest 提供唯一权威 label_date
assert not any(column.endswith(('_x', '_y')) for column in m07_manifest.columns)  # 禁止上游列冲突后缀
m07_development_key_set = set(map(tuple, m07_development_manifest[m07_key_columns].to_numpy()))  # 冻结权威开发键集合
assert m07_development_key_set == set(map(tuple, m07_feature_frame[m07_key_columns].to_numpy())) == set(map(tuple, m07_label_frame[m07_key_columns].to_numpy()))  # 拒绝开发制品缺键或增键
m07_frame = m07_development_manifest.merge(m07_feature_frame, on=m07_key_columns, how='left', validate='one_to_one', indicator='feature_merge')  # 权威左连接开发特征
assert m07_frame['feature_merge'].eq('both').all()  # 要求全部开发键有特征
m07_frame = m07_frame.drop(columns='feature_merge').merge(m07_label_frame, on=m07_key_columns, how='left', validate='one_to_one', indicator='label_merge')  # 权威左连接开发标签
assert m07_frame['label_merge'].eq('both').all()  # 要求全部开发键有标签
m07_frame = m07_frame.rename(columns={'y': 'event'}).drop(columns='label_merge')  # 连接后统一章内事件名
assert not any(column.endswith(('_x', '_y')) for column in m07_frame.columns)  # 明确拒绝静默重名
assert not m07_frame[m07_features + ['event', 'label_date']].isna().any().any()  # 禁止静默丢弃开发键
m07_frame = m07_frame.sort_values(['prediction_date', 'order_book_id']).reset_index(drop=True)  # 固定面板顺序
m07_train = m07_frame[m07_frame['split'].eq('train')].copy().reset_index(drop=True)  # 冻结训练键
m07_validation = m07_frame[m07_frame['split'].eq('validation')].copy().reset_index(drop=True)  # 冻结验证键
assert m07_train['label_date'].max() < m07_validation['prediction_date'].min()  # 执行严格 purge
assert m07_train['event'].nunique() == m07_validation['event'].nunique() == 2  # 禁止单类拟合或 AUC
m07_development_hashes = {'train': m07_key_hash(m07_train), 'validation': m07_key_hash(m07_validation), 'combined': m07_key_hash(m07_frame)}  # 登记开发键哈希
列表 7.21: M07步骤04:核验M05持久化折合同与文件哈希
m07_policy_path = Path(os.environ['BOOK_PROJECT_FOLD_POLICY']).resolve()  # 解析 M05 持久化折政策
m07_fold_path = Path(os.environ['BOOK_PROJECT_FOLD_ARTIFACT']).resolve()  # 解析 M05 持久化折成员
m07_policy_bytes = m07_policy_path.read_bytes()  # 保留政策原始字节
m07_fold_bytes = m07_fold_path.read_bytes()  # 保留折成员原始字节
m07_policy_hash = hashlib.sha256(m07_policy_bytes).hexdigest()  # 计算实际政策哈希
m07_fold_hash = hashlib.sha256(m07_fold_bytes).hexdigest()  # 计算实际折成员哈希
assert m07_policy_hash == os.environ['BOOK_PROJECT_FOLD_POLICY_SHA256'].lower()  # 拒绝政策漂移
assert m07_fold_hash == os.environ['BOOK_PROJECT_FOLD_SHA256'].lower()  # 拒绝折成员漂移
m07_policy = json.loads(m07_policy_bytes)  # 解析已核验政策
m07_fold_artifact = json.loads(m07_fold_bytes)  # 解析已核验折成员
assert m07_policy['schema_version'] == '1' and m07_policy['contract_id'] == 'a-share-drawdown-20d-v1'  # 核对政策身份
assert m07_policy['manifest_sha256'] == m07_manifest_hash and m07_policy['source_split'] == 'train'  # 绑定唯一 manifest 训练段
assert m07_policy['key_columns'] == m07_key_columns and m07_policy['construction'] == 'expanding-date-blocks'  # 禁止键或构造漂移
assert m07_policy['purge_rule'] == 'fit.label_date < min(score.prediction_date)'  # 固定 purge 规则
assert m07_fold_artifact['schema_version'] == '1' and m07_fold_artifact['contract_id'] == m07_policy['contract_id']  # 核对折制品版本与任务身份
assert m07_fold_artifact['policy_sha256'] == m07_policy_hash and m07_fold_artifact['manifest_sha256'] == m07_manifest_hash  # 关闭制品关联链
列表 7.22: M07步骤05:从持久化键恢复训练内前向折
m07_train_key_index = {(row.order_book_id, row.prediction_date.strftime('%Y-%m-%d')): row.Index for row in m07_train.itertuples()}  # 将持久化键映射为训练行索引
m07_forward_folds = []  # 只保存从 M05 制品还原的索引对
for m07_fold_record in m07_fold_artifact['folds']:  # 逐折核验持久化成员
    m07_fold_core = {key: value for key, value in m07_fold_record.items() if key != 'fold_sha256'}  # 分离单折内容与登记哈希
    m07_actual_fold_hash = hashlib.sha256((json.dumps(m07_fold_core, ensure_ascii=False, sort_keys=True, separators=(',', ':')) + '\n').encode()).hexdigest()  # 复算单折规范哈希
    assert m07_actual_fold_hash == m07_fold_record['fold_sha256']  # 任一折内容被改动即失败
    m07_fit_keys = [tuple(key) for key in m07_fold_record['fit_keys']]  # 恢复拟合键
    m07_score_keys = [tuple(key) for key in m07_fold_record['score_keys']]  # 恢复评分键
    assert len(m07_fit_keys) == len(set(m07_fit_keys)) and len(m07_score_keys) == len(set(m07_score_keys))  # 折内键唯一
    assert set(m07_fit_keys).isdisjoint(m07_score_keys)  # 拟合与评分成员互斥
    assert set(m07_fit_keys + m07_score_keys).issubset(m07_train_key_index)  # 折键只能来自 manifest train
    m07_fit_index = np.array([m07_train_key_index[key] for key in m07_fit_keys], dtype=int)  # 恢复拟合索引
    m07_score_index = np.array([m07_train_key_index[key] for key in m07_score_keys], dtype=int)  # 恢复评分索引
    assert m07_train.iloc[m07_fit_index]['label_date'].max() < m07_train.iloc[m07_score_index]['prediction_date'].min()  # 用权威终点日复核 purge
    assert m07_fold_record['fit_label_date_max'] == m07_train.iloc[m07_fit_index]['label_date'].max().strftime('%Y-%m-%d')  # 核对拟合边界
    assert m07_fold_record['score_prediction_date_min'] == m07_train.iloc[m07_score_index]['prediction_date'].min().strftime('%Y-%m-%d')  # 核对评分边界
    assert m07_train.iloc[m07_fit_index]['event'].nunique() == m07_train.iloc[m07_score_index]['event'].nunique() == 2  # 禁止单类折
    m07_forward_folds.append((m07_fit_index, m07_score_index))  # 登记已核验 M05 折
列表 7.23: M07步骤06:仅用训练内前向折选择平滑参数
m07_candidates = np.logspace(-2, 3, 8)  # 事前登记有限平滑参数
m07_inner_brier = []  # 保存每个参数的训练内层损失
for candidate_lambda in m07_candidates:  # 显式执行前向验证而不调用 GCV/UBRE
    fold_brier = []  # 收集当前平滑度的三折证据
    for fit_index, score_index in m07_forward_folds:  # 仅在训练键内选择
        fold_model = LogisticGAM(s(0, lam=candidate_lambda) + l(1, lam=candidate_lambda))  # 固定平滑动量与线性波动率结构
        fold_model.fit(m07_train.iloc[fit_index][m07_features].to_numpy(), m07_train.iloc[fit_index]['event'].astype(int).to_numpy())  # 拟合早期训练行
        fold_probability = fold_model.predict_proba(m07_train.iloc[score_index][m07_features].to_numpy())  # 预测更晚内层块
        fold_brier.append(brier_score_loss(m07_train.iloc[score_index]['event'].astype(int), fold_probability))  # 记录当折 Brier
    m07_inner_brier.append(float(np.mean(fold_brier)))  # 汇总当前平滑度的时间折损失
m07_selected_lambda = float(m07_candidates[int(np.argmin(m07_inner_brier))])  # 只依据训练内层冻结平滑度
列表 7.24: M07步骤07:消费冻结参数并对验证键评分一次
m07_model = LogisticGAM(s(0, lam=m07_selected_lambda) + l(1, lam=m07_selected_lambda))  # 恢复冻结 LogisticGAM 配方
m07_model.fit(m07_train[m07_features].to_numpy(), m07_train['event'].astype(int).to_numpy())  # 仅在训练键上拟合
m07_validation_probability = m07_model.predict_proba(m07_validation[m07_features].to_numpy())  # 验证键只预测一次
m07_validation_metrics = {'auc': float(roc_auc_score(m07_validation['event'].astype(int), m07_validation_probability)), 'brier': float(brier_score_loss(m07_validation['event'].astype(int), m07_validation_probability)), 'inner_forward_brier': float(min(m07_inner_brier))}  # 登记完整验证证据
m07_baseline_brier = brier_score_loss(m07_validation['event'].astype(int), np.repeat(m07_train['event'].mean(), len(m07_validation)))  # 建立训练事件率基线
assert np.isfinite(list(m07_validation_metrics.values()) + [m07_baseline_brier]).all()  # 禁止非有限结果
print('baseline Brier=', m07_baseline_brier, 'lambda=', m07_selected_lambda, 'GAM=', m07_validation_metrics, 'manifest_sha256=', m07_manifest_hash)  # 现场输出而不预写胜者

图 7.18 同时展示训练样本真正覆盖的二维支持域与动量平滑项的偏效应。红色验证点落在训练支持矩形之外时,应把相应预测解释为外推,而不是把曲线形状写成普遍规律。

m07_momentum_bounds = m07_train['momentum_20'].quantile([0.01, 0.99]).to_numpy()  # 用训练分位数定义稳健支持区间
m07_volatility_bounds = m07_train['volatility_20'].quantile([0.01, 0.99]).to_numpy()  # 用训练分位数定义波动率支持区间
m07_is_outside_support = ~m07_validation['momentum_20'].between(*m07_momentum_bounds) | ~m07_validation['volatility_20'].between(*m07_volatility_bounds)  # 按训练中央矩形识别验证域外点
m07_outside_support_count = int(m07_is_outside_support.sum())  # 计算index要求的域外计数
print({'validation_outside_support': m07_outside_support_count, 'validation_rows': len(m07_validation)})  # 输出可审计支持域证据
m07_partial_grid = m07_model.generate_X_grid(term=0, n=100)  # 只在已拟合 GAM 的基函数网格上评价偏效应
m07_partial_effect, m07_partial_interval = m07_model.partial_dependence(term=0, X=m07_partial_grid, width=0.95)  # 计算实际拟合偏效应与区间
m07_figure, m07_axes = plt.subplots(1, 2, figsize=(11, 4.5))  # 创建支持域与偏效应并列画布
m07_axes[0].scatter(m07_train['momentum_20'], m07_train['volatility_20'], s=8, alpha=0.18, color='#2C3E50', label='训练')  # 绘制真实训练支持点
m07_axes[0].scatter(m07_validation['momentum_20'], m07_validation['volatility_20'], s=10, alpha=0.28, color='#E3120B', label='验证')  # 叠加冻结验证点识别外推
m07_axes[0].axvspan(m07_momentum_bounds[0], m07_momentum_bounds[1], color='#008080', alpha=0.08)  # 标出训练动量中央支持区
m07_axes[0].axhspan(m07_volatility_bounds[0], m07_volatility_bounds[1], color='#F0A700', alpha=0.08)  # 标出训练波动率中央支持区
m07_axes[0].set(xlabel='二十日动量', ylabel='二十日波动率', title='训练支持域与验证位置')  # 提供业务可读坐标与标题
m07_axes[0].legend()  # 解释训练与验证样本颜色
m07_axes[1].plot(m07_partial_grid[:, 0], m07_partial_effect, color='#008080', linewidth=2)  # 绘制实际拟合的动量偏效应
m07_axes[1].fill_between(m07_partial_grid[:, 0], m07_partial_interval[:, 0], m07_partial_interval[:, 1], color='#8E9EAA', alpha=0.25)  # 绘制点态置信带
m07_axes[1].set_xlim(m07_momentum_bounds)  # 将解释限制在训练中央支持区
m07_axes[1].set(xlabel='二十日动量', ylabel='对数几率偏效应', title='冻结平滑项的支持域内形状')  # 明确曲线尺度与边界
m07_figure.tight_layout()  # 避免中文标签重叠
plt.show()  # 输出可评分支持域证据
图 7.18
列表 7.25: M07步骤09:按冻结配方重拟并构造候选记录
m07_development_model = LogisticGAM(s(0, lam=m07_selected_lambda) + l(1, lam=m07_selected_lambda))  # 按冻结配方新建最终开发模型
m07_development_model.fit(m07_frame[m07_features].to_numpy(), m07_frame['event'].astype(int).to_numpy())  # 在训练加验证键上重拟
m07_artifact_bytes = pickle.dumps(m07_development_model, protocol=5)  # 序列化 M09 将读取的真实候选
m07_registry_path = Path(os.environ['BOOK_CANDIDATE_REGISTRY']).expanduser().resolve()  # 锁定公共 registry
m07_registry_path.parent.mkdir(parents=True, exist_ok=True)  # 确保运行器授权目录存在
m07_artifact_directory = m07_registry_path.parent / 'candidate_artifacts'  # 统一产物根目录
m07_artifact_directory.mkdir(parents=True, exist_ok=True)  # 创建产物目录
m07_artifact_path = m07_artifact_directory / 'm07-logistic_gam.pkl'  # 固定候选文件名
m07_artifact_path.write_bytes(m07_artifact_bytes)  # 持久化冻结模型
m07_record = {'contract_id': 'a-share-drawdown-20d-v1', 'milestone_id': 'M07', 'candidate_id': 'logistic_gam', 'implementation': 'pygam.pygam.LogisticGAM', 'feature_schema': {'columns': m07_features, 'dtype': 'float64', 'as_of': 'company-day-t-close'}, 'manifest_sha256': m07_manifest_hash, 'fold_policy_sha256': m07_policy_hash, 'fold_sha256': m07_fold_hash, 'development_key_hashes': m07_development_hashes, 'preprocessing': ['pyGAM internal spline basis fitted on development keys'], 'frozen_hyperparameters': {'lambda': m07_selected_lambda, 'terms': ['s(momentum_20)', 'l(volatility_20)']}, 'validation_metrics': m07_validation_metrics, 'support_domain': {'definition': 'training 0.01-0.99 marginal quantile rectangle', 'validation_outside_count': m07_outside_support_count}, 'refit_recipe': {'fit_splits': ['train', 'validation'], 'target': 'event', 'probability_rule': 'predict_proba'}, 'artifact_sha256': hashlib.sha256(m07_artifact_bytes).hexdigest(), 'artifact_file': str(m07_artifact_path.relative_to(m07_registry_path.parent)), 'status': 'ready', 'failure_reason': None}  # 填满公共 schema、M05折证据与域外计数
列表 7.26: M07步骤10:校验并幂等写入公共registry
m07_required_fields = {'contract_id', 'milestone_id', 'candidate_id', 'implementation', 'feature_schema', 'manifest_sha256', 'fold_policy_sha256', 'fold_sha256', 'development_key_hashes', 'preprocessing', 'frozen_hyperparameters', 'validation_metrics', 'refit_recipe', 'artifact_sha256', 'status', 'failure_reason'}  # 复用含 M05 折证据的公共 schema
assert m07_required_fields.issubset(m07_record)  # 缺任一必需字段即失败
m07_existing_records = [json.loads(line) for line in m07_registry_path.read_text(encoding='utf-8').splitlines() if line.strip()] if m07_registry_path.exists() else []  # 读取已存在候选
m07_preserved_records = [record for record in m07_existing_records if not (record.get('milestone_id') == 'M07' and record.get('candidate_id') == 'logistic_gam')]  # 移除本候选旧版本
m07_complete_registry = sorted(m07_preserved_records + [m07_record], key=lambda record: (record['milestone_id'], record['candidate_id']))  # 固定 registry 行序
m07_registry_text = ''.join(json.dumps(record, ensure_ascii=False, sort_keys=True, separators=(',', ':')) + '\n' for record in m07_complete_registry)  # 生成可复算 JSONL
m07_registry_path.write_text(m07_registry_text, encoding='utf-8')  # 幂等写回完整 registry
print('registry_sha256=', hashlib.sha256(m07_registry_text.encode()).hexdigest(), 'candidate_id=logistic_gam')  # 输出 registry 审计证据

若当前 pyGAM 版本没有 l,必须登记失败而不能把因子项冒充线性项。LogisticGAM 只在训练键内选平滑度,冻结验证键只记一次指标;registry 写权限与测试集不可见性仍由外部运行器执行。

  1. 二维平滑 vs 可加平滑:二维平滑可捕捉交互,但更吃样本、解释更难;GAM 更稳健且解释清晰,但若真实关系强交互可能欠拟合。用外推测试误差与可视化(二维等高线 vs 两条一维曲线)对比即可。

7.13.3 理论题参考解答(要点)

  1. 线性平滑器与自由度:平滑样条解是对 \(y\) 的线性变换,存在矩阵 \(S_\lambda\) 使得 \(\hat y=S_\lambda y\)\(\mathrm{tr}(S_\lambda)\) 衡量“输出对输入的敏感度总量”,等价于线性模型中 hat matrix 的迹,因此可视为有效自由度。\(\lambda\uparrow\) 时更平滑,\(\mathrm{tr}(S_\lambda)\) 下降;\(\lambda\downarrow 0\) 时接近插值,\(\mathrm{tr}(S_\lambda)\) 上升。

  2. 边界偏差更小的直觉:局部常数回归在边界附近因为邻域不对称导致一阶项无法抵消,从而产生较大偏差;局部线性回归显式拟合截距与斜率,能在边界处用线性项补偿不对称,故边界偏差通常更小。

7.14 章末学习闭环

逐项自检:能否区分多项式、阶梯、样条、局部回归与 GAM;能否说明 GCV 不等于时间 CV;能否检查支持域;能否冻结 LogisticGAM record。禁止在稀疏支持区强外推,也禁止把偏依赖当因果。常见误区是预写最优阶数、用训练 RSS 推荐模型。

不看正文回答:阶梯函数与回归样条在节点处的连续性有何不同?平滑参数控制什么?为何同时报告支持密度?迁移任务:为季节性需求风险设计非线性候选。下一章复用共同样本与候选制品,比较树集成。

展开检索答案 阶梯函数可在切点跳变,回归三次样条以基/约束保证到二阶导连续;平滑参数权衡拟合与曲率;低密度区域由少量观测支持,外推风险更高。

出口决策:阶梯/回归样条/平滑样条的结构区别、时间验证、支持域三项为 must-pass;三项全过、其余目标至少一项有证据且检索至少两题正确才进入第 8 章。结构错回 小节 7.3小节 7.5,用“补贴门槛”比较跳变和 C2 连续;验证错把新标签窗画入前向折;支持域错以“设备温度—故障率”标出域外点。每项复测 2 分后重入 小节 7.14

Hastie, Trevor, 和 Robert Tibshirani. 1986年. 《Generalized Additive Models》. Statistical Science 1 (3): 297~318. https://doi.org/10.1214/ss/1177013604.