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

7.1 导读

经济金融关系很少在整个支持域内严格线性,但直接提高模型复杂度也会增加方差、外推风险与解释难度。本章依次介绍多项式、阶梯函数、回归样条、平滑样条、局部回归和广义加性模型,并以真实 A 股行情为例比较它们与线性基线的样本外表现。方法选择取决于数据支持、验证误差、外推需求与解释目标,而不是固定的复杂度排序。

7.2 学习目标

完成本章后,读者应能够:

  1. 区分多项式、阶梯函数、回归样条、平滑样条和局部回归的基函数与支持域。
  2. 在同一前向时间切分上比较线性与非线性模型,并复用训练期的样条设计信息。
  3. 区分节点数、平滑参数与有效自由度,以及 GCV/UBRE 与显式时间验证。
  4. 解释自然样条的边界约束和局部回归的维度限制。
  5. 把 GAM 的偏依赖限定为模型预测形状,并结合支持密度、业务损失和样本外证据作判断。
  6. 区分固定基设计、数据驱动选择与惩罚平滑的推断条件,并为时间数据选择匹配依赖结构的不确定性方法。

7.2.1 一周核心路线与进阶证据

一周核心路线分成两条证据链。连续响应路线先用多项式和阶梯函数理解“把输入变成基函数”,再掌握自然样条的节点、自由度与边界线性约束,并由练习 5 在共同日期前向折上比较线性、三次多项式与自然样条的价格水平 MSE。二元响应路线先按 列表 7.18列表 7.19列表 7.20列表 7.21 构造共享对象,再按 列表 7.22列表 7.23表 7.2表 7.3表 7.4图 7.18 运行唯一一套 SplineTransformer + LogisticRegression 逻辑样条 GAM,在固定训练—验证边界上与训练事件率和线性 Logit 比较。两条路线的响应、损失和切分用途不同,绝不组成同一“四模型”排名。标为 [拓展] 或“阅读示例”的平滑样条、LOWESS、连续响应 Ridge GAM 与第二套 pygam 实现,第一次学习均可跳过。

进入第 8 章前的最低证据是:能解释同一输入为何可对应不同基函数设计;能说明自然样条的边界约束为何不等于安全外推;能在连续价格响应的共同前向折报告线性、三次多项式和自然样条的逐折 MSE 与时间支持域;再在二元回撤响应的固定训练—验证键上报告事件率、线性 Logit 与锁定逻辑样条 GAM 的 Brier、AUROC、四特征支持域及模型内预测形状。连续路线失败时返回 列表 7.24表 7.5,二元路线失败时返回 列表 7.22 和其后核心标签;不要求先完成任何拓展实现,也不得跨两种响应比较分数。

到目前为止,本书主要关注线性模型。线性模型易于描述和实现,在解释和推断方面具有优势,但线性条件均值可能无法表示数据支持域内的弯曲或分段关系。在 章节 6 中,岭回归、Lasso 和主成分回归通过约束线性模型降低估计方差;本章则改变基函数,以容纳可验证的非线性形状。

在本章中,我们将放宽线性假设,同时尽可能保持可解释性。多项式回归和阶梯函数直接扩展线性模型的基函数设计;样条、局部回归和广义加性模型则提供更细的局部形状控制。

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

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

回归样条(Regression Splines)通过在不同区间使用低次多项式扩展线性模型和阶梯函数。各段在区域边界或节点(knots)处满足连续性约束;增加节点可以提高局部灵活性,也会增加估计方差,因此节点与自由度需要用训练期验证选择。

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

局部回归(Local Regression)与样条类似,但邻域可以重叠,权重通常随目标点距离连续变化。

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

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

7.3 多项式回归 (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\),然后可用线性最小二乘的计算工具估计系数。但“计算形式是线性的”不等于经典推断自动有效:只有次数和基函数预先固定、估计未受惩罚,并且误差结构与所用协方差估计相匹配时,才能采用相应的标准误、区间或检验;具体边界见 小节 7.5.1

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

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

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

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

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

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

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

列表 7.1: 导入必要的库与数据
# 准备多项式拟合、时间切分、区间估计和曲线比较所需接口
import numpy as np              # 构造多项式设计、日期网格与有限指标
import pandas as pd             # 对齐行情日期、特征、响应与折叠结果
import matplotlib.pyplot as plt # 比较拟合形状、验证误差和支持域
import statsmodels.api as sm    # 估计 OLS/GLM 并计算模型内区间
from sklearn.preprocessing import PolynomialFeatures  # 在每个拟合流程内生成预定次数的幂基
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  # 拼接文件路径并为后续独立代码块登记统一数据根
BOOK_DATA_DIR = os.path.abspath('/home/ubuntu/r2_data_mount/data')  # 明文定义在线教材的BOOK_DATA_DIR绝对路径
DATA_ROOT = BOOK_DATA_DIR  # 保留本章后续代码使用的数据根名称
os.environ['BOOK_DATA_DIR'] = 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')  # 定位海康威视多项式案例的前复权行情
assert os.path.isfile(path), f'缺少前复权行情文件: {path};请检查 BOOK_DATA_DIR'  # 在HDF读取前报告确切缺失输入
haikang_order_book_id = '002415.XSHE'  # 锁定与修订前完全相同的海康威视统计对象
haikang_data = pd.read_hdf(  # 利用HDF表的data_columns在存储层筛选公司并限制载入字段
    path,  # 读取已核验存在的全市场前复权行情表
    where=f'order_book_id == {haikang_order_book_id!r}',  # 只载入海康威视全部可用交易日,保持原样本支持域
    columns=['close']  # 只读取响应构造所需收盘价,日期和公司键由表索引保留
).reset_index()  # 将HDF多重索引恢复为公司与日期列

# 统一日期列名为 trade_date,兼容不同数据源的命名差异
# 获取列名列表
if 'date' in haikang_data.columns and 'trade_date' not in haikang_data.columns:  # 仅在本地表使用date列名时执行统一命名
    haikang_data = haikang_data.rename(columns={'date': 'trade_date'})  # 重命名日期列以保持后续估计接口不变
haikang_data['trade_date'] = pd.to_datetime(haikang_data['trade_date'])  # 将交易日期统一转换为datetime以确保日期轴正确显示

# 存储层筛选已完成;这里只按交易日升序排列
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: 海康威视股价数据摘要
# 汇总时间支持与价格尺度,为后续基函数图提供读数基准
# 包含均值、标准差、四分位数等,帮助了解价格分布
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

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

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

下面用同一海康威视价格序列比较 1、3、5、10、15 和 20 次多项式的样本内拟合。次数升高会降低或保持训练误差,却可能在支持域边缘产生较大振荡。这一图形用于展示高次全局多项式的数值与边界风险;它不能仅凭曲线形状判定真实价格过程,也不能替代前向验证。

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

# 只在观测日期支持域内设置五百个评价点以比较边界形状
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))  # 用同一画布尺度并列比较六种全局次数

# 在共同样本上逐一拟合预定次数并保存样本内形状诊断
for i, degree in enumerate(polynomial_degrees, 1):  # 在共同样本上逐一估计预定次数
    # Pipeline 将 "升维 → 回归" 两步封装:
    #   PolynomialFeatures(d) 把 X 变成 [1, X, X², ..., X^d]
    #   LinearRegression() 对升维后的矩阵做最小二乘
    current_model = Pipeline([  # 将当前次数升维与 OLS 绑定,避免步骤口径漂移
        ('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)  # 计算当前模型在训练集上的均方误差

    # 把当前次数的样本内形状放入对应面板,供后文检查边界振荡
    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  # 以扩张窗口估计各多项式次数的未来块误差

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

# 对每个候选多项式次数执行交叉验证
for degree in cv_polynomial_degrees:  # 在同一前向折上逐一评价预先声明的次数
    # 构建多项式回归 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

# 将共同扩张折的误差按次数对齐,读取验证最低点
plt.figure(figsize=(10, 6))  # 展示多项式次数与前向验证 MSE 的对应关系
plt.plot(cv_polynomial_degrees, cross_validation_scores, 'bo-', linewidth=2)  # 蓝色圆点折线:各次数的 CV MSE
plt.xlabel('多项式次数', fontsize=14)  # 明确横轴是候选模型复杂度
plt.ylabel('CV MSE (TS-Split)', fontsize=14)  # 明确纵轴是前向折预测损失
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.3.2 方差分析

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

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

# ── 逐次拟合多项式(1 次 → 5 次),记录残差平方和 RSS 与自由度 ──
max_polynomial_degree = 5  # 限定逐级 F 诊断比较到五次项
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('多项式次数的方差分析:')  # 标明下表只比较嵌套模型的训练期 RSS
print('-' * 80)  # 分隔嵌套模型诊断的标题与列名
print(f'{"模型":<20} {"残差自由度":<15} {"RSS (x1e6)":<15} {"F统计量":<15} {"p值":<15}')  # 明示嵌套比较所需的自由度、拟合改善与尾部概率
print('-' * 80)  # 分隔列名与逐级 F 诊断结果

# 逐级比较相邻两个模型的拟合改善程度
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 / 当前残差自由度)
    # 将相邻模型的 RSS 改善按新增自由度与高阶模型残差方差标准化
    f_statistic = max(0.0, (rss_difference / dof_difference) / (residual_sum_squares[i] / degrees_of_freedom_values[i]))
    # p 值:在 F 分布下,观测到的 F 统计量以上的尾部面积
    # 在对应分子、分母自由度的 F 分布下计算右尾概率
    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 计算经典参照分布下的尾部面积。每一行只有在次数预先固定、误差独立同方差且相应分布条件成立时,才是新增一阶多项式项的有效有限样本 F 检验。这里的同一价格序列既有时序依赖,五个次数又被连续查看,因此打印出的 \(p\) 值只是演示经典计算的未校正诊断量,不具备本案例中的名义显著性解释,也不是未来预测证据。若研究目标是预先指定系数的关联推断,可在条件合适时使用 HAC 协方差;若目标是选中次数后的整条曲线不确定性,则需按 小节 7.5.1 所述重做选择并采用块重抽样或有效的选择后方法。复杂度预测仍由后文前向验证选择。

# 用五次多项式说明模型内均值区间,不把区间当样本外误差
degree = 5  # 选择的多项式次数

# 构建 sklearn Pipeline:先生成多项式特征,再做线性回归
poly_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%点态均值区间')  # 标明该区间未校正时间依赖与复杂度选择
plt.xlabel('时间', fontsize=14)  # 明确横轴对应交易日期
plt.ylabel('股价 (元)', fontsize=14)  # 保留响应变量的人民币单位
plt.title('海康威视股价的多项式回归拟合', fontsize=16)  # 标明条件均值模型与案例对象
plt.legend(fontsize=12)  # 区分观测价格、条件均值与模型内置信带
plt.grid(True, alpha=0.3)  # 辅助读取拟合曲线与点态区间的日期位置
plt.show()  # 输出多项式拟合及其模型内不确定性带
海康威视价格散点、五次多项式拟合曲线及其置信带随交易日期变化。
图 7.3: 五次多项式拟合与经典点态均值区间(时序推断条件未满足)

图 7.3 的蓝色曲线展示五次多项式拟合,浅蓝色阴影是把次数视为固定并假设独立同方差误差时得到的 95% 经典点态均值区间。它只演示软件在该工作模型下计算什么,不代表当前股价序列具有 95% 覆盖率:时序相关、异方差、查看多个次数后的选择以及同时观察整条曲线都会破坏这一解释。边界处区间变宽只能说明该工作模型的杠杆与估计不确定性增加;该带既不是未来价格预测区间,也不是同时置信带。若要对当前时间序列作形状推断,应采用与选择过程及依赖结构匹配的 HAC 或块重抽样方案,并明确支持域。

7.3.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] 之间

# 获取固定设计与独立观测工作假设下的经典点态区间
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='预测概率')  # 把二项GLM条件概率映射回观测日期支持域
# 填充区域
plt.fill_between(date_grid, predicted_confidence_lower, predicted_confidence_upper,
                 alpha=0.2, color='blue', label='经典95%点态区间')  # 标明区间未校正时间依赖与样本阈值选择

plt.xlabel('交易日期', fontsize=14)  # 明确概率形状对应的样本日期
plt.ylabel(f'P(股价 ≥ {high_price_threshold:.2f}元)', fontsize=14)  # 明确纵轴事件定义与价格阈值
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.5.1 处理选择与时间依赖。该图只描述样本内条件模式,不支持均值回复、趋势机制或因果效应。

7.4 阶梯函数 (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.4.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  # 固定等频时间箱的哑变量列顺序
# 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)  # 明确横轴为交易日期
plt.ylabel('股价 (元)', fontsize=14)  # 保留分段响应的人民币单位
plt.title('海康威视股价的阶梯函数拟合 (分段常数)', fontsize=16)  # 标明面板展示分段常数条件均值
plt.legend(fontsize=12)  # 区分观测路径与分箱后的条件均值
plt.grid(True, alpha=0.3)  # 辅助比较分段均值与箱边界位置
plt.show()  # 输出阶梯函数的分段常数拟合形状
海康威视价格散点叠加按时间箱分段的阶梯状预测线。
图 7.5: 阶梯函数回归拟合

箱数决定了阶梯函数的灵活度,应在预先给定的候选集合中用训练期交叉验证选择。下面比较 5、10、15、20 和 30 个箱子,并用 sklearn.preprocessing.KBinsDiscretizer 保证验证数据沿用训练折学到的分箱边界。图中最低的验证 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  # 特征:时间索引

    # 对每个交叉验证折进行训练和验证
    # 在每个共同前向折内重新学习分箱边界,避免验证信息进入阈值
    for train_indices, validation_indices in kfold_iterator.split(time_index_features):
        # ── 分割训练集和验证集 ──
        # 保留二维形状以满足折内离散器和线性模型接口
        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 随分箱数量的变化;最低点只是当前候选集内的运行时选择。

# 将共同扩张折的误差按分箱数对齐,避免用训练拟合选阈值
plt.figure(figsize=(10, 6))  # 展示阶梯函数箱数与前向验证 MSE 的对应关系
plt.plot(candidate_bin_counts, step_cv_mse_scores, 'bo-', linewidth=2, markersize=8)  # 蓝色圆点折线
plt.xlabel('分箱数量 (Bins)', fontsize=14)  # 明确横轴是阶梯函数复杂度
plt.ylabel('CV MSE', fontsize=14)  # 明确纵轴是前向折预测损失
plt.title('阶梯函数:箱子数量选择', fontsize=16)  # 标明分箱数由时间验证比较
plt.grid(True, alpha=0.3)  # 辅助读取候选箱数之间的验证误差差距
plt.show()  # 输出候选分箱数对应的前向验证误差
折线以箱数为横轴、前向验证MSE为纵轴,展示五个候选分箱复杂度的误差。
图 7.6: 阶梯函数候选箱数的前向验证均方误差

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

7.5 基函数 (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),\ldots,b_K(x_i)\) 为预测变量的线性模型,并用最小二乘估计系数。至于标准误、置信区间和整体 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.1 基函数模型的推断边界

基函数模型常同时承担“预测未来 \(Y\)”与“描述条件均值形状 \(E(Y\mid X=x)\)”两项任务。前者主要依靠严格样本外损失评价;后者若要报告标准误、区间或显著性,还必须先登记推断对象、基函数选择过程和误差依赖结构。以下四种情形不能混用:

  1. 固定、未惩罚的基设计。 若次数、节点和基函数在查看本次响应之前已指定,则条件于设计矩阵可以使用 OLS 估计。经典同方差标准误以及精确 \(t/F\) 分布还要求独立同方差误差和相应正态条件;异方差下应采用适当的稳健协方差,并把结果理解为满足相应大样本条件时的近似推断。
  2. 同一数据驱动的复杂度选择。 若先用响应选择次数、节点、自由度或候选变换,再把选中设计当作预先固定,朴素标准误、置信区间和 F 检验会遗漏选择不确定性。交叉验证可以支持预测复杂度选择,却不会自动修复选择后推断。若形状推断是主要目标,可预先指定设计,或使用样本分割、有效的选择后推断程序与外部复现样本,并清楚说明条件于何种选择事件。
  3. 惩罚样条与惩罚式 GAM。 惩罚会收缩系数,复杂度由平滑参数和有效自由度共同决定;普通未惩罚 OLS 的精确 F 分布不能机械套用。软件给出的平滑项检验或区间通常是依赖惩罚选择与模型假设的近似量,应采用与惩罚估计匹配的协方差或后验近似,必要时用完整重选平滑参数的重抽样评估不确定性。
  4. 金融时间序列。 异方差、自相关和结构变化使独立误差区间失效。在预先固定的低维参数目标下,可在条件适合时使用 HAC 协方差;对整条拟合曲线或数据驱动平滑过程,可使用保留局部依赖的块重抽样,并在每次重抽样内重做全部选择。块长度、平稳性近似和支持域必须说明。无论采用哪种区间,时间趋势关联都不自动识别因果效应,也不应外推到观测支持域之外。

因此,本章的前向验证回答候选函数能否改善未来预测;它不把选中的形状转换为未经校正的显著性结论。后续图中的经典区间若用于教学展示,将明确标注其条件和在当前价格序列中的局限。

7.6 回归样条 (Regression Splines)

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

7.6.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.6.2 约束与样条

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

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

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

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

7.6.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 \] 设观测支持区间为 \([a,b]\),区间内部有 \(K\) 个节点。自然样条要求在 \(a\) 左侧和 \(b\) 右侧保持线性,等价于施加 \(f''(a)=f''(b)=0\)。从普通三次样条的 \(K+4\) 维空间出发,这两个独立线性约束把维数减少为 \(K+2\)

在截断幂表示中,把边界外二阶导数设为零后,可将约束写成对三次项系数的两个独立线性关系。不要把左右两侧展开时出现的多个代数等式重复计数:三次样条已经满足函数及前两阶导数在内部节点连续,新增的自然边界条件只有两个。因此,若 \(K\) 明确表示内部节点数,自然三次样条的函数空间维数为 \(K+2\)(包括截距与线性趋势);软件中的 df 还可能因是否显式包含截距而采用不同列数约定,必须查看设计矩阵而不能仅凭节点数猜测。

7.6.4 回归样条的Python实现

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

列表 7.7: 导入样条相关库
from patsy import dmatrix                    # 生成训练期节点固定的三次样条基
from sklearn.linear_model import LinearRegression  # 在给定样条基上估计线性系数
import statsmodels.api as sm                  # 对给定样条设计计算 OLS 推断量

回归三次样条是在节点之间使用分段三次多项式,但各段并非独立估计。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)',
                   # 用时间索引数值构造预设节点的 B 样条基
                   {'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()  # 在预设节点的 B 样条基上估计训练拟合

模型已经拟合完成,接下来在 图 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)  # 明确横轴为交易日期
plt.ylabel('股价 (元)', fontsize=14)  # 保留样条响应的人民币单位
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() 在两端施加"线性"约束,防止边界发散
# 用时间索引数值构造带线性边界约束的自然样条基
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)  # 明确横轴为交易日期
plt.ylabel('股价 (元)', fontsize=14)  # 保留自然样条响应的人民币单位
plt.title('自然样条回归拟合', fontsize=16)  # 标明曲线施加自然边界约束
plt.legend(fontsize=12)  # 区分观测路径与自然边界约束下的拟合
plt.grid(True, alpha=0.3)  # 辅助读取自然边界约束下的拟合形状
plt.show()  # 输出自然样条在当前支持域内的拟合形状
海康威视价格散点叠加自然样条曲线,边界附近曲线趋于线性。
图 7.8: 自然样条回归拟合

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

7.6.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):  # 在共同前向折内估计当前自由度的未来误差
        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},  # 只用当前训练折时间支持域拟合自然样条基
                         return_type='dataframe')  # 指定返回DataFrame格式
        features_val = build_design_matrices([features_train.design_info], {'Time': data_validation['Time'].values})[0]  # 复用训练期节点、边界与列顺序

        # 在训练集上拟合 OLS
        # 沿用训练折 design_info 变换后续验证期,禁止重估节点
        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 随自由度的变化,并从实际输出中选择当前候选集内的自由度。

# 将共同扩张折的误差按自由度对齐,读取训练期选择证据
plt.figure(figsize=(10, 6))  # 展示自然样条自由度与前向验证 MSE 的对应关系
plt.plot(candidate_degrees_of_freedom, spline_cv_mse_scores, 'bo-', linewidth=2, markersize=8)  # 蓝色圆点折线
plt.xlabel('自由度', fontsize=14)  # 明确横轴为自然样条候选复杂度
plt.ylabel('交叉验证均方误差', fontsize=14)  # 明确纵轴为前向折预测损失
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}')  # 报告当前候选与切分下的选择
print(f'对应的最小交叉验证MSE: {minimum_cv_mse:.4f}')  # 报告选择所依据的前向验证损失
折线以自然样条自由度为横轴、前向验证MSE为纵轴,并以竖虚线标出当次最低误差对应的自由度。
图 7.9: 自然样条候选自由度的前向验证均方误差
最优自由度: 4
对应的最小交叉验证MSE: 129.3479

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

7.6.6 与多项式回归的比较

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

下面在同一数据与复杂度设置下比较 15 次全局多项式和 15 自由度自然样条。这个设定主要用于观察边界行为;最终选择必须依据同一验证切分下的损失,而不是曲线看起来是否平滑。 下图只展示当前训练样本中的两种拟合形状。若 15 次全局多项式在边界振荡而自然样条更平稳,这与自然样条的边界约束一致;但图形本身不证明它捕捉了真实信号,也不证明样条普遍胜出。模型选择必须读取同一前向折误差、离散度、支持域和基线;若样条未改善验证损失,应明确报告“无增量证据”。

列表 7.11: 比较自然样条与多项式回归
# ── 对手 A:拟合 15 次多项式 ──
polynomial_model_15deg = 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},  # 在同一时间支持域构造自然样条对照
                       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()  # 汇总各折负MSE作为统一比较准则

# 对样条需手动循环(statsmodels 不兼容 sklearn cross_val_score)
fold_mse_natural_15dof = []  # 保存自然样条在各前向折的MSE以评估波动
# 手动循环 5 折交叉验证,计算自然样条的 CV MSE
for train_indices, validation_indices in kfold_iterator.split(haikang_data):  # 在同一前向折比较两类基函数的未来误差
    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},  # 仅用当前训练折估计样条设计信息
                     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}')  # 报告全局高次基的前向误差
print(f'15自由度自然样条的交叉验证MSE: {cv_mse_natural_15dof:.4f}')  # 报告自然边界基的前向误差
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)  # 明确横轴为共同交易日索引
axes[0].set_ylabel('股价 (元)', fontsize=14)  # 保留响应的人民币单位
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)  # 保持右图使用相同交易日索引
axes[1].set_ylabel('股价 (元)', fontsize=14)  # 保持两图响应尺度一致
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.7 平滑样条 (Smoothing Splines)

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

7.7.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.7.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 最小的候选值。平滑样条的线性平滑矩阵允许用一次全样本拟合的残差与杠杆值计算 LOOCV,而不必显式重拟合 \(n\) 次:

\[ \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(x_i)\) 与平滑矩阵对角元,就能计算每个留一残差。章节 5式 5.2 对最小二乘线性回归给出同类恒等式;对回归样条或其他线性基函数模型,它避免了逐个观测重新拟合的计算。

7.7.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  # 以单一平滑项实现惩罚样条机制比较

# ── 准备特征矩阵和目标向量 ──
time_features = haikang_data[['Time']].values  # shape=(n,1)
price_target = haikang_data['Price'].values     # 取每个交易日价格作为长度为 n 的一维连续响应向量

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

plt.figure(figsize=(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)  # 明确横轴为交易日期
plt.ylabel('股价 (元)', fontsize=14)  # 保留平滑响应的人民币单位
plt.title('不同平滑参数的平滑样条拟合', fontsize=16)  # 标明面板比较惩罚强度与曲率
plt.legend(fontsize=11)  # 对齐各曲线的lambda与现场有效自由度
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])  # 读取 GCV 在预设网格内选中的惩罚强度
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))  # 为GCV选中曲线与观测路径保留共同坐标轴
plt.plot(haikang_data['trade_date'], price_target, 'k-', alpha=0.3, linewidth=1, label='观测股价')  # 提供读取GCV拟合形状的价格参照

optimal_predicted_price = smoothing_gam.predict(time_grid)  # 最优模型网格预测
plt.plot(date_grid, optimal_predicted_price, 'b-', linewidth=3,  # 把GCV选中条件均值映射回观测日期
         label=f'GCV选中拟合 (λ={optimal_lambda:.2f})')  # 图例标注准则选中的 λ 值

plt.xlabel('时间', fontsize=14)  # 明确横轴为交易日期
plt.ylabel('股价', fontsize=14)  # 明确纵轴为价格水平
plt.title('通过GCV准则选择的平滑样条', fontsize=16)  # 标明曲线只反映 GCV 候选选择
plt.legend(fontsize=12)  # 报告GCV在预设网格内选中的lambda
plt.grid(True, alpha=0.3)  # 辅助读取 GCV 所选曲线与观测路径的偏离
plt.show()  # 输出GCV准则选中平滑器的训练拟合形状
  0% (0 of 30) |                         | Elapsed Time: 0:00:00 ETA:  --:--:--

 10% (3 of 30) |##                       | Elapsed Time: 0:00:00 ETA:   0:00:00

 23% (7 of 30) |#####                    | Elapsed Time: 0:00:00 ETA:   0:00:00

 36% (11 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

 76% (23 of 30) |##################      | Elapsed Time: 0:00:00 ETA:   0:00:00

 90% (27 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. 计算复杂度
    • 回归样条:若设计矩阵有 \(K\) 列,稠密最小二乘形成交叉乘积约需 \(O(nK^2)\)、直接求解约需 \(O(K^3)\);具体成本还取决于 QR、迭代或稀疏求解器
    • 平滑样条:以每个训练点为节点不等于计算必然为 \(O(n^2)\)\(O(n^3)\);朴素稠密矩阵会达到二次存储和三次求解成本,利用带状结构或专门算法则可显著降低,某些实现接近关于 \(n\) 的线性成本
  4. 灵活性与规模
    • 回归样条:较小的 \(K\) 可控制存储与估计成本,但大样本适用性仍取决于基维数、矩阵条件数和求解器
    • 平滑样条:复杂度由节点表示、惩罚矩阵结构和求解器共同决定,不能仅凭方法名称判定大样本可行性

方法选择条件

  • 回归样条用较少节点控制计算量;是否适合较大样本还取决于基函数数目与实现。
  • 平滑样条把复杂度选择转化为 \(\lambda\) 的选择,但仍须指定候选准则并核对数据依赖结构。
  • 线性平滑器的 LOOCV 可由杠杆值恒等式计算;时间预测仍应使用保持顺序的验证。
  • 两者是否产生相近结果取决于节点、惩罚、支持域与评价样本,不能预先假定。

本段拓展实现到此结束;一周核心路线可直接前往 小节 7.12,从 列表 7.18 开始沿正文依赖链运行。

7.8 局部回归 (Local Regression)

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

7.8.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.8.2 [拓展] LOWESS 历史平滑阅读示例

本阅读示例使用 statsmodelslowess(Locally Weighted Scatterplot Smoothing)展示邻域宽度如何改变历史平滑形状;它不属于一周核心路线,也不提供未来预测证据。

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

from statsmodels.nonparametric.smoothers_lowess import lowess  # 在各目标时点按邻域权重估计局部直线

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

plt.figure(figsize=(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)  # 明确横轴为历史交易日期
plt.ylabel('股价 (元)', fontsize=14)  # 保留局部平滑响应的人民币单位
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.8.2.1 [拓展] 局部窗口的训练 RSS 调参演示

下面三块代码仅展示为何训练 RSS 会机械偏向较小窗口,以及为什么预测窗口必须改由训练期前向验证选择;它不属于一周核心路线。

列表 7.12: 局部回归不同 frac 参数的残差平方和比较
from statsmodels.nonparametric.smoothers_lowess import lowess  # 对每个候选窗口复算训练样本局部拟合
import numpy as np  # 汇总候选窗口的训练残差平方和

# 事前固定七个窗口,仅演示训练RSS为何不能承担预测选参
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:  # 逐一计算预定窗口的样本内 RSS 诊断
    # 执行 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))  # 展示训练 RSS 随窗口缩小的机械变化
 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))  # 展示 LOESS 窗口比例与训练 RSS 的诊断关系
plt.plot(rss_comparison_df['frac'], rss_comparison_df['RSS'], 'o-', linewidth=2, markersize=8)  # 比较窗口比例与样本内残差的机械关系
plt.xlabel('frac (窗口比例)', fontsize=14)  # 明确横轴为局部邻域占比
plt.ylabel('残差平方和 (RSS)', fontsize=14)  # 明确纵轴仅为样本内拟合损失
plt.title('局部回归 frac 参数选择', fontsize=16)  # 图标题
plt.grid(True, alpha=0.3)  # 辅助读取各窗口比例对应的训练RSS差距

# 标注最优 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),  # 避免文字遮挡相邻窗口的训练 RSS
             arrowprops=dict(arrowstyle='->', color='red'), fontsize=12, color='red')  # 指向样本内最低点而不暗示预测最优
plt.show()  # 输出训练RSS诊断,不把最低点解释为预测最优窗口
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\)

不适用场景

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

本节 LOWESS 阅读示例到此结束;二元核心路线可直接前往 小节 7.12,从 列表 7.18 开始沿正文依赖链运行。

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

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

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

7.9.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.3小节 7.8 中,我们讨论了许多拟合单变量函数的方法。GAM 把这些单变量方法作为可加预测器的构建块,使每一项的函数形式可以分别指定与检查。

以自然样条为例,考虑使用自然样条对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.6 所讨论,自然样条可以使用适当选择的基函数来构造。因此,整个模型只是对样条基变量和虚拟变量的大型回归,所有这些都打包在一个大型回归矩阵中。

7.9.2 GAMs的优缺点

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

优点

  • GAM 允许为每个 \(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.9.3 [拓展阅读] 连续响应的 Ridge GAM 实现

scikit-learn 可把 SplineTransformer 生成的单变量样条基与线性末端组合成可加预测器。下面的连续收益示例用于阅读 Ridge 末端与偏依赖,不属于一周核心路线;二元核心路线从 列表 7.18 开始,不依赖本拓展对象。

我们将构建如下金融因子模型: \[ \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)
# 准备本拓展连续响应示例的基函数、类别编码和正则线性末端
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))  # 并列承载三个连续特征的边际预测形状

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

# 关注的三个特征名称列表
selected_gam_features = ['Momentum', 'Volatility', 'Month']  # 保持偏依赖面板与拟合特征顺序一致

# 一次性生成 3 张偏依赖图(average = 平均 ICE 曲线)
partial_dependence_display = PartialDependenceDisplay.from_estimator(  # 对训练期联合经验分布平均其他特征的预测
    ridge_gam_pipeline,      # 已拟合的 GAM 管道
    gam_features_train,      # 只用训练期经验分布确定网格并做边际平均
    selected_gam_features,   # 要展示的特征
    kind='average',          # 只画均值线,不画每条 ICE 曲线
    ax=ax,  # 把三个特征的偏效应绑定到同一出版画布
    **partial_dependence_params,  # 展开参数字典传入函数
)  # 生成基于训练期经验分布平均的三个偏效应面板

# 标明面板显示模型预测形状,并为三个标题保留可读间距
# 在偏效应面板上方标明预测量而非因果效应
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 把各特征的非线性基与线性末端组合成可加预测器,便于把变换、折内预处理和评价放进同一管道。可加结构提高了单变量形状的可读性,但不自动包含交互,也不保证优于线性基线。本段连续响应阅读示例到此结束;二元核心路线可直接前往 小节 7.12,并从 列表 7.18 开始。

7.9.4 [拓展阅读] 同期行情样本的逻辑 GAM

GAM 也可以用于定性响应。下面沿用本节行情样本预测未来五日收益是否为正,用来观察把 Ridge 末端替换为逻辑链接后的概率输出;它不是一周核心实现。

\[ \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))  # 并列承载逻辑GAM的三个概率偏效应

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)  # 标出概率尺度的决策参考而非真值边界

plt.show()  # 输出逻辑GAM在概率尺度上的边际预测形状
三个面板显示动量、波动率和月份与模型预测上涨概率的平均偏依赖及0.5参考线。
图 7.16: 逻辑回归GAM的偏依赖图 (涨跌概率)

逻辑回归 GAM 的偏依赖图展示各特征与模型预测上涨概率之间的边际预测关联,灰色虚线是 0.5 决策参考。倒 U 形、单调抑制或月份差异都只能在图中确实出现时描述,并须结合测试期 Brier、基线和样本密度判断稳定性。本段阅读示例到此结束;二元回撤核心路线可直接前往 小节 7.12,并从 列表 7.18 开始。

7.9.5 [拓展] 第二套 pygam 实现与完整预测器偏依赖

本节用于比较另一套平滑项语法与完整预测器偏依赖,不属于一周核心路线;核心进阶不要求安装或运行 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 要求 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 = ['动量对上涨概率的影响',      # 标明第一项是模型中的动量预测形状
          '波动率对上涨概率的影响',    # 标明第二项是模型中的波动率预测形状
          '季节性对上涨概率的影响']    # 标明第三项是模型中的月份预测形状

partial_data = []  # 存储每个特征的网格与平均预测概率
for i in range(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 分三个子图展示动量、波动率和月份对应的平均上涨概率:

# 在概率尺度汇总第二套GAM的单项预测形状
fig, axes = plt.subplots(1, 3, figsize=(18, 5))  # 1 行 3 列子图

# 将每个特征的现场概率区间放到独立面板以检查支持差异
for i, ax in enumerate(axes):  # 把三个可加项放入共同概率尺度的面板
    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])                         # 保持所有面板使用合法概率尺度
    ax.axhline(0.5, color='gray', linestyle='--')  # 50% 参考线——涨跌分界

plt.tight_layout()  # 防止三个概率偏效应面板的标签互相遮挡
plt.show()  # 输出第二套GAM实现的完整预测器边际概率形状
pyGAM完整线性预测器下,动量、波动率和月份对应的平均上涨概率曲线。
图 7.17: 逻辑回归GAM的偏依赖图 (涨跌概率)

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

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

列表 7.17: GAM模型的方差分析
from pygam import 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])  # 汇总月份因子编码的训练偏差
deviance_linear_dist = float(np.ravel(gam_model_linear_dist.statistics_['deviance'])[0])  # 汇总月份线性编码的训练偏差
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])  # 量化月份线性编码的训练期复杂度
dof_smooth_dist = float(np.ravel(gam_model_smooth_dist.statistics_['edof'])[0])  # 量化月份平滑编码的训练期复杂度

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 分数。本节第二套实现到此结束;二元核心路线返回 小节 7.12,并从 列表 7.18 开始。

GAM 的使用条件

以下检查项用于限定 GAM 的解释与选择范围:

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

当可加结构近似合理、训练支持充分且前向验证优于基线时,GAM 可以作为灵活性与可解释性之间的候选方案。

7.10 本章小结 (Chapter Summary)

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

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

方法比较

方法 灵活性 可解释性 典型成本驱动 适用场景
多项式回归 样本量、次数与矩阵条件数 简单非线性关系
阶梯函数 样本量与分箱数 分段常数关系
回归样条 样本量、基维数 \(K\) 与求解器 一般非线性关系
平滑样条 节点表示、惩罚矩阵结构与求解器 需要自动平滑选择的关系
局部回归 查询点数、邻域大小与维数 探索性分析
GAMs 总基维数、响应族与迭代求解器 多变量非线性关系

实践建议

  • 从简单的多项式回归或样条开始
  • 使用交叉验证选择模型复杂性
  • 检查拟合曲线的合理性,特别是在边界处
  • 多个预测变量的非线性主要可由可加项表示时,把 GAM 纳入候选
  • 始终与线性基线模型比较

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

一周核心进阶证据由练习 1—5 与 小节 7.12 共同提供:练习 5 在连续价格响应上报告线性、三次多项式与自然样条的共同前向折 MSE 和时间支持域;小节 7.12 则在二元回撤响应上,从 列表 7.22 开始依次报告事件率、线性 Logit 与唯一逻辑样条 GAM 的固定验证 Brier、AUROC、四特征支持域和模型内预测形状。两条路线不得跨 estimand 排名。练习 6—9 以及标为 [拓展] 或“阅读示例”的实现可后读,不作为进入第 8 章的前提;任何答案都不得预写非线性方法的胜者。

7.11 理论来源与前沿

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

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

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

7.12 综合案例:逻辑样条能否改善未来回撤概率

本章的独立行情案例保留连续收益、波动率和多种非线性方法,用于完整学习基函数与 GAM。贯穿项目不改变估计对象,并把第 3 章的线性概率基线升级为合法概率模型。为保证从干净内核运行,下面在本章内定义并调用与 小节 2.7 相同的 BOOK_DATA_DIR 构造器;真实峰谷最大回撤标签、日期边界和四项返回对象均保持一致,不读取外部辅助脚本。

先定义前向最大回撤函数。对路径 \(P_0,\ldots,P_{20}\),先计算运行峰值 \(H_j=\max_{0\le k\le j}P_k\),再取 \(\min_j(P_j/H_j-1)\)。因此价格路径 \(100\to120\to105\) 的最大回撤是 \(105/120-1=-12.5\%\),即使终点仍高于预测日价格,也必须记为事件。

列表 7.18: 二元核心路线:峰谷最大回撤标签函数与夹具
展开共享回撤标签函数
import os  # 读取共享行情样本的数据根配置
from pathlib import Path  # 核验后复权行情文件并构造可移植路径
import numpy as np  # 计算滚动特征、峰谷回撤与有限概率指标
import pandas as pd  # 按公司—预测日对齐特征、标签和时间分段

def calculate_forward_max_drawdown(price_values, horizon_days=20):  # 计算含预测日在内的前向峰谷最大回撤
    price_array = np.asarray(price_values, dtype=float)  # 转为连续数值数组以构造滑动路径
    window_size = horizon_days + 1  # 路径包含t日与其后horizon_days个交易日
    drawdowns = np.full(price_array.size, np.nan, dtype=float)  # 末端窗口不完整时保持未知
    if price_array.size < window_size:  # 保护短序列而不虚构标签
        return drawdowns  # 返回全缺失结果表示没有完整路径
    price_windows = np.lib.stride_tricks.sliding_window_view(price_array, window_size)  # 构造每个预测日的完整价格路径
    running_peaks = np.maximum.accumulate(price_windows, axis=1)  # 逐路径记录每一时点之前的最高价
    drawdowns[:price_windows.shape[0]] = np.min(price_windows / running_peaks - 1, axis=1)  # 取峰值后最深谷值跌幅
    return drawdowns  # 返回与原价格序列等长的前向最大回撤

m02_fixture_drawdown = calculate_forward_max_drawdown([100, 120, 105], horizon_days=2)[0]  # 用先涨后跌路径区分峰谷回撤与起点收益
assert np.isclose(m02_fixture_drawdown, -0.125)  # 核对峰值120到谷值105的跌幅为百分之十二点五
assert m02_fixture_drawdown <= -0.10  # 核对该路径必须被判为百分之十回撤事件

下面把数据读取、预测日特征和峰谷标签封装为一个准备函数。固定公司名单是教学研究边界,不是事后按收益挑出的“优胜者”;return_1d 等特征均以 \(t\) 日结束,只有标签计算向未来展开。

列表 7.19: 二元核心路线:共享行情特征与标签准备函数
展开共享行情与标签准备函数
def prepare_m02_price_panel(book_data_dir):  # 从统一数据根构造未切分的公司—预测日样本
    data_root = Path(book_data_dir).expanduser().resolve()  # 同时接受环境变量字符串或Path对象
    assert data_root.is_dir(), f'BOOK_DATA_DIR 不存在: {data_root}'  # 在读取前验证数据挂载
    price_path = data_root / 'stock/stock_price_post_adjusted.h5'  # 指向后复权日行情
    assert price_path.is_file(), f'缺少后复权行情文件: {price_path};请检查 BOOK_DATA_DIR'  # 在read_hdf前给出明确缺失文件
    company_ids = ('600104.XSHG', '002415.XSHE', '600276.XSHG', '002230.XSHE')  # 固定四家长三角教学公司
    price_frames = [pd.read_hdf(price_path, where=f'order_book_id={company_id!r}', columns=['close', 'volume']).reset_index() for company_id in company_ids]  # 在存储层选择所需公司和字段
    price_panel = pd.concat(price_frames, ignore_index=True)  # 合并为公司—交易日面板
    price_panel['prediction_date'] = pd.to_datetime(price_panel['date'])  # 把交易日定义为收盘后预测原点
    price_panel = price_panel.loc[price_panel['prediction_date'].between('2012-01-01', '2024-12-31')].copy()  # 固定研究观察期
    price_panel = price_panel.sort_values(['order_book_id', 'prediction_date']).reset_index(drop=True)  # 保证窗口只在公司内按时间推进
    grouped_prices = price_panel.groupby('order_book_id', sort=False)  # 建立公司内窗口边界
    price_panel['return_1d'] = grouped_prices['close'].transform(lambda prices: prices.pct_change(fill_method=None))  # 计算截至预测日的一日收益
    price_panel['momentum_5d'] = grouped_prices['close'].transform(lambda prices: prices.pct_change(5, fill_method=None))  # 计算截至预测日的五日动量
    price_panel['volatility_20d'] = grouped_prices['return_1d'].transform(lambda returns: returns.rolling(20, min_periods=20).std())  # 计算过去二十日波动率
    volume_mean_20d = grouped_prices['volume'].transform(lambda volume: volume.rolling(20, min_periods=20).mean())  # 估计预测日前成交量常态
    price_panel['volume_ratio_20d'] = price_panel['volume'] / volume_mean_20d  # 表示当日成交量相对过去均值的活跃度
    feature_names = ['return_1d', 'momentum_5d', 'volatility_20d', 'volume_ratio_20d']  # 固定跨章特征名称与顺序
    price_panel['max_drawdown_20d'] = grouped_prices['close'].transform(lambda prices: calculate_forward_max_drawdown(prices, horizon_days=20))  # 计算t至t+20真实峰谷最大回撤
    price_panel['label_end_date'] = grouped_prices['prediction_date'].shift(-20)  # 记录未来第二十个公司交易日
    price_panel['drawdown_event_20d'] = price_panel['max_drawdown_20d'].le(-0.10).where(price_panel['max_drawdown_20d'].notna()).astype('Int64')  # 完整路径才生成二元事件
    price_panel = price_panel.replace([np.inf, -np.inf], np.nan).dropna(subset=feature_names + ['max_drawdown_20d', 'drawdown_event_20d', 'label_end_date']).copy()  # 排除无效特征和未完成路径
    return price_panel, feature_names  # 把未切分样本交给公共构造函数

公共函数负责固定切分、边界断言和四项跨章交付。任何章节在干净内核中都可以用 build_m02_shared_sample(BOOK_DATA_DIR) 重建同一对象,不需要依赖前一章的内存状态。

列表 7.20: 二元核心路线:共享时间切分与四项返回对象
展开共享样本切分函数
def build_m02_shared_sample(book_data_dir):  # 返回第2至第9章唯一的共享分析对象
    price_panel, feature_names = prepare_m02_price_panel(book_data_dir)  # 从统一数据根重建真实峰谷标签
    validation_start = pd.Timestamp('2020-01-01')  # 固定验证段起点
    test_start = pd.Timestamp('2022-01-01')  # 固定封存测试段起点
    study_end = pd.Timestamp('2024-12-31')  # 固定研究终点
    train_mask = (price_panel['prediction_date'] < validation_start) & (price_panel['label_end_date'] < validation_start)  # 清除伸入验证期的训练标签
    validation_mask = price_panel['prediction_date'].between(validation_start, test_start, inclusive='left') & (price_panel['label_end_date'] < test_start)  # 清除伸入测试期的验证标签
    test_mask = price_panel['prediction_date'].between(test_start, study_end, inclusive='both')  # 标记只供第九章最终评价的日期
    price_panel['split'] = np.select([train_mask, validation_mask, test_mask], ['train', 'validation', 'test'], default='unused')  # 赋予互斥时间段
    shared_columns = ['order_book_id', 'prediction_date', 'label_end_date', *feature_names, 'max_drawdown_20d', 'drawdown_event_20d', 'split']  # 复用第2章权威字段顺序
    shared_sample = price_panel.loc[price_panel['split'] != 'unused', shared_columns].copy()  # 排除边界隔离行并形成共享对象
    assert list(shared_sample.columns) == shared_columns  # 阻止同名对象发生字段或顺序漂移
    assert not shared_sample.duplicated(['order_book_id', 'prediction_date']).any()  # 验证公司—预测日键唯一
    assert np.isfinite(shared_sample['max_drawdown_20d']).all() and shared_sample['max_drawdown_20d'].le(0).all()  # 验证连续目标有限且非正
    assert shared_sample['drawdown_event_20d'].astype(bool).eq(shared_sample['max_drawdown_20d'].le(-0.10)).all()  # 验证二元标签由连续目标唯一派生
    assert set(shared_sample['drawdown_event_20d'].unique()) == {0, 1}  # 验证二元标签包含两类
    assert shared_sample.loc[shared_sample['split'] == 'train', 'label_end_date'].max() < validation_start  # 验证训练标签不触及验证期
    assert shared_sample.loc[shared_sample['split'] == 'validation', 'label_end_date'].max() < test_start  # 验证验证标签不触及测试期
    training_event_rate = float(shared_sample.loc[shared_sample['split'] == 'train', 'drawdown_event_20d'].mean())  # 只用训练标签估计概率基线
    test_keys = shared_sample.loc[shared_sample['split'] == 'test', ['order_book_id', 'prediction_date']].copy()  # 固定第九章最终比较键
    return shared_sample, feature_names, training_event_rate, test_keys  # 返回约定的四项跨章交付
列表 7.21: 二元核心路线:从BOOK_DATA_DIR构造共享样本
book_data_dir_value = os.environ.get('BOOK_DATA_DIR')  # 安全读取统一数据根配置
assert book_data_dir_value, '请先设置 BOOK_DATA_DIR,使其指向包含 stock/ 子目录的数据根'  # 缺失时说明修复方法
BOOK_DATA_DIR = Path(book_data_dir_value).expanduser().resolve()  # 解析统一数据根
m02_shared_sample, m02_features, m02_training_event_rate, m02_test_keys = build_m02_shared_sample(BOOK_DATA_DIR)  # 在当前干净内核构造共享对象

本章后续只读取训练段与验证段;测试键虽已固定,但测试标签和成绩继续封存。干净内核按 列表 7.18列表 7.19列表 7.20列表 7.21列表 7.22 的顺序开始,再运行后续核心标签;若出现未定义对象,应回到这条精确依赖链,不应转去运行任何 pygam 拓展块。

比较对象是线性逻辑回归与逻辑样条 GAM。后者对每个连续特征分别生成三次 B 样条基,再在线性预测器中相加,所以能够表达单变量弯曲,但不会自动包含特征交互。节点数只按验证 Brier 在预先声明的有限集合中选择;测试段保持封存。

列表 7.22: 逻辑样条核心路线的固定训练—验证准备
from sklearn.pipeline import make_pipeline  # 把核心案例的折内变换与概率模型绑定
from sklearn.preprocessing import SplineTransformer, StandardScaler  # 在训练段生成样条基并估计尺度
from sklearn.linear_model import LogisticRegression  # 估计线性基线与逻辑样条的事件概率
from sklearn.metrics import brier_score_loss, roc_auc_score  # 比较概率误差与排序能力
import matplotlib.pyplot as plt  # 绘制锁定GAM的四特征支持域内条件切片
m07_expected_columns = ['order_book_id', 'prediction_date', 'label_end_date', *m02_features, 'max_drawdown_20d', 'drawdown_event_20d', 'split']  # 锁定第2章连续与二元目标的有序模式
assert list(m02_shared_sample.columns) == m07_expected_columns  # 阻止字段、目标或顺序悄然变化
m07_train_sample = m02_shared_sample.loc[m02_shared_sample['split'] == 'train'].copy()  # 取得固定训练键
m07_validation_sample = m02_shared_sample.loc[m02_shared_sample['split'] == 'validation'].copy()  # 取得固定验证键
m07_train_features = m07_train_sample[m02_features]  # 按第2章顺序提取训练特征
m07_validation_features = m07_validation_sample[m02_features]  # 按相同顺序提取验证特征
m07_train_target = m07_train_sample['drawdown_event_20d'].astype(int)  # 使用同一二元回撤标签
m07_validation_target = m07_validation_sample['drawdown_event_20d'].astype(int)  # 保持验证估计对象不变
m07_development_keys = pd.concat([m07_train_sample, m07_validation_sample])[['order_book_id', 'prediction_date']]  # 汇总本章允许访问的开发键
assert m07_development_keys.merge(m02_test_keys, on=['order_book_id', 'prediction_date']).empty  # 证明核心路线没有打开测试键
m07_linear_logit = make_pipeline(StandardScaler(), LogisticRegression(C=1.0, max_iter=2000, random_state=42))  # 建立可校准的线性逻辑基线
m07_linear_logit.fit(m07_train_features, m07_train_target)  # 只用训练段估计预处理与系数
m07_linear_probability = m07_linear_logit.predict_proba(m07_validation_features)[:, 1]  # 输出验证期合法概率

下面只在训练—验证开发边界内选择节点数。每个候选管道都会在训练段重新拟合样条基,因此验证日期不会参与节点位置或标准化估计。

列表 7.23: 一周核心:逻辑样条 GAM 的训练期拟合与验证选择
m07_candidate_knots = [3, 5, 7]  # 预先声明有限的样条复杂度集合
m07_candidate_models = []  # 保存每个已拟合候选供锁定选择
m07_candidate_records = []  # 保存验证证据而不预写胜者
for knot_count in m07_candidate_knots:  # 逐一比较预先声明的节点数
    candidate_model = make_pipeline(SplineTransformer(n_knots=knot_count, degree=3, include_bias=False), StandardScaler(), LogisticRegression(C=1.0, max_iter=2000, random_state=42))  # 构造可加逻辑样条
    candidate_model.fit(m07_train_features, m07_train_target)  # 只在训练段学习基函数和系数
    candidate_probability = candidate_model.predict_proba(m07_validation_features)[:, 1]  # 生成同一验证键概率
    candidate_brier = brier_score_loss(m07_validation_target, candidate_probability)  # 用概率误差选择复杂度
    m07_candidate_models.append(candidate_model)  # 保存候选以避免选择后重复窥视验证集
    m07_candidate_records.append({'n_knots': knot_count, 'validation_brier': candidate_brier, 'validation_auc': roc_auc_score(m07_validation_target, candidate_probability)})  # 记录概率与排序证据
m07_candidate_table = pd.DataFrame(m07_candidate_records)  # 形成可读验证表
m07_selected_index = int(m07_candidate_table['validation_brier'].argmin())  # 按预先指定的Brier准则锁定候选
m07_spline_gam = m07_candidate_models[m07_selected_index]  # 保存后续可重拟合的锁定逻辑样条
m07_spline_probability = m07_spline_gam.predict_proba(m07_validation_features)[:, 1]  # 读取锁定候选的验证概率

动态结果按四步阅读:先看 表 7.2 的节点数与验证 Brier,确认选择确实遵循预设准则;再由 表 7.3 把锁定样条与线性逻辑、训练事件率放在同一验证键比较;随后核对 表 7.4 的越界比例;最后读取 图 7.18 的支持域内模型预测形状。Brier 较低而 AUROC 无变化,可能只是校准改善;反之则可能只是排序改善。任何一次验证胜出都不是未来稳定性保证。

m07_candidate_table.assign(selected=lambda frame: frame.index == m07_selected_index)  # 同表标出按验证Brier锁定的节点候选
表 7.2: 逻辑样条节点候选的验证期选择证据
n_knots validation_brier validation_auc selected
0 3 0.204605 0.694188 False
1 5 0.204503 0.693817 True
2 7 0.205515 0.692703 False
m07_event_rate_probability = np.repeat(m02_training_event_rate, len(m07_validation_target))  # 构造同一训练事件率概率基线
m07_model_table = pd.DataFrame({'model': ['training_event_rate', 'linear_logit', 'selected_spline_gam'], 'validation_brier': [brier_score_loss(m07_validation_target, m07_event_rate_probability), brier_score_loss(m07_validation_target, m07_linear_probability), brier_score_loss(m07_validation_target, m07_spline_probability)], 'validation_auc': [0.5, roc_auc_score(m07_validation_target, m07_linear_probability), roc_auc_score(m07_validation_target, m07_spline_probability)]})  # 在相同键上比较三种概率预测
print('锁定节点数:', m07_candidate_knots[m07_selected_index])  # 明示被锁定的复杂度
m07_model_table  # 以可渲染表格呈现基线差、排序和概率误差
锁定节点数: 5
表 7.3: 第7章线性逻辑与逻辑样条GAM的验证期动态比较
model validation_brier validation_auc
0 training_event_rate 0.229024 0.500000
1 linear_logit 0.207203 0.691998
2 selected_spline_gam 0.204503 0.693817

下面逐项核对四个特征的训练范围,以及固定验证期落在该范围之外的比例。越界比例不是模型成绩,却直接标记样条外推风险;即使比例为零,稀疏的联合支持仍可能限制解释。

m07_support_records = []  # 累积每个共享特征的边际支持证据
for feature_name in m02_features:  # 对四项固定输入执行相同范围核对
    feature_minimum = float(m07_train_features[feature_name].min())  # 只用训练段确定支持下界
    feature_maximum = float(m07_train_features[feature_name].max())  # 只用训练段确定支持上界
    outside_support = ~m07_validation_features[feature_name].between(feature_minimum, feature_maximum)  # 标记验证值是否需要边际外推
    m07_support_records.append({'feature': feature_name, 'train_min': feature_minimum, 'train_max': feature_maximum, 'validation_outside_rate': float(outside_support.mean())})  # 保存范围与越界率
m07_support_table = pd.DataFrame(m07_support_records)  # 形成四特征支持域审计表
m07_support_table  # 渲染训练范围和验证越界比例
表 7.4: 逻辑样条四特征的训练支持域与验证越界比例
feature train_min train_max validation_outside_rate
0 return_1d -0.100177 0.100249 0.0
1 momentum_5d -0.316233 0.433268 0.0
2 volatility_20d 0.005162 0.072038 0.0
3 volume_ratio_20d 0.000000 20.000000 0.0

图 7.18 在每个特征的训练期 5%—95% 分位区间内逐一改变该特征,其他特征固定为训练中位数。它展示的是锁定模型的一维条件切片,不是边际因果效应,也不揭示未纳入的交互。

m07_reference_values = m07_train_features.median().to_dict()  # 以训练中位数定义共同条件切片基准
m07_shape_figure, m07_shape_axes = plt.subplots(2, 2, figsize=(12, 8))  # 为四项共享特征建立同尺度诊断面板
for shape_axis, feature_name in zip(m07_shape_axes.ravel(), m02_features):  # 逐项绘制锁定模型的支持域内条件形状
    shape_lower, shape_upper = m07_train_features[feature_name].quantile([0.05, 0.95])  # 避开训练支持极端稀疏尾部
    shape_grid = np.linspace(shape_lower, shape_upper, 100)  # 在训练期中央支持域生成评价网格
    shape_frame = pd.DataFrame([m07_reference_values] * len(shape_grid))  # 其他特征保持训练中位数
    shape_frame[feature_name] = shape_grid  # 只改变当前特征以形成条件切片
    shape_probability = m07_spline_gam.predict_proba(shape_frame[m02_features])[:, 1]  # 用锁定模型生成事件概率
    shape_axis.plot(shape_grid, shape_probability)  # 展示支持域内预测形状而非因果效应
    shape_axis.set(xlabel=feature_name, ylabel='预测回撤概率', title=f'{feature_name} 条件切片')  # 标明特征和概率尺度
m07_shape_figure.tight_layout()  # 防止四面板标题与坐标标签遮挡
plt.show()  # 输出锁定GAM的四项支持域内形状证据
四个面板分别在训练期百分之五至百分之九十五分位范围内改变单一特征,纵轴为锁定逻辑样条模型的回撤事件预测概率。
图 7.18: 锁定逻辑样条GAM在四项训练支持域内的预测形状

只有当锁定样条在验证 Brier 上改善线性逻辑和事件率基线、AUROC 不出现实质退化、验证越界比例可接受且支持域内形状可解释时,才有理由把它交给第 9 章做一次最终测试。条件切片仍是模型预测形状,不是某个行情变量对回撤的因果效应;第 7 章不得访问测试标签来补救不理想的验证结果。

7.13 练习

7.13.1 概念题

  1. 区分多项式次数、回归样条节点数、平滑样条参数 \(\lambda\) 与有效自由度。 [核心|难度:1|分值:5|任务:独立]

  2. 为什么自然三次样条在边界外为线性?若有 \(K\) 个内部节点,其函数空间维数是多少? [核心|难度:2|分值:8|任务:独立]

  3. 比较阶梯函数、局部回归与 GAM 的适用场景和主要局限。 [核心|难度:2|分值:8|任务:独立]

  4. 解释为什么 GCV 最优不等于时间前向验证最优,并说明用同一数据选自由度后,普通 OLS 的 F 检验或点态区间为何不能直接作名义推断。 [核心|难度:2|分值:6|任务:独立]

7.13.2 应用题

  1. 使用海康威视真实行情,在同一前向切分上比较线性、三次多项式和自然样条,并报告逐折 MSE 与支持域。 [核心|难度:3|分值:20|任务:独立]

  2. 为“杠杆率与下一期现金流风险”设计 GAM,说明每个平滑项、交互需求、信息时点和解释边界。 [拓展|难度:3|分值:15|任务:综合案例]

  3. 比较二维平滑 \(f(x_1,x_2)\) 与可加平滑 \(f_1(x_1)+f_2(x_2)\) 的泛化误差和解释性。 [拓展|难度:3|分值:15|任务:结对]

7.13.3 理论题

  1. 推导普通三次样条与自然三次样条的维数。 [拓展|难度:3|分值:15|任务:推导]

  2. 说明平滑样条目标函数在 \(\lambda\to0\)\(\lambda\to\infty\) 时的极限行为。 [拓展|难度:3|分值:15|任务:推导]

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

7.14 练习参考解答

7.14.1 概念题参考解答

  1. 多项式次数决定全局幂基最高阶;回归样条节点数决定局部拼接位置;\(\lambda\) 是平滑样条粗糙度惩罚强度;有效自由度是平滑矩阵对响应的实际适应程度,通常随 \(\lambda\) 增大而下降。它们相关但不能互换。

  2. 自然约束为支持区间两端的二阶导数为零,因而边界外只保留截距与线性趋势。普通三次样条若有 \(K\) 个内部节点,维数为 \(K+4\);两个独立自然边界约束使维数降为 \(K+2\)。软件 df 是否含截距要以设计矩阵列数为准。

  3. 阶梯函数适合业务本身按区间决策且阈值可解释的场景,但边界不连续;局部回归适合低维、支持密集的平滑关系,但计算成本高且受维度灾难影响;GAM 以可加结构兼顾非线性与解释性,却会漏掉未显式加入的交互。三者都需同一验证设计和支持域检查。

  4. GCV 近似留一风险,通常建立在观测可交换和所用平滑矩阵的条件上;时间预测中未来观测不能参与训练,序列相关和制度漂移也会改变误差结构。因此 GCV 可用于训练期候选搜索,却不能替代按日期推进的前向验证。若同一响应又用于选择自由度,选中设计的普通 OLS F 检验或点态区间会遗漏选择不确定性;惩罚平滑还改变有效自由度,不能机械套用未惩罚 OLS 分布。固定低维参数可在条件适合时用 HAC 协方差;整条曲线或重选平滑度的过程宜用保留依赖且每次重做选择的块重抽样,或采用有效的选择后推断。

7.14.2 应用题参考解答

  1. 完整量规(20 分)

    • 数据(4 分):读取本地海康威视后复权收盘价,把未来日期上的价格水平作为连续验证响应;只在各训练折拟合标准化和基函数。
    • 切分(4 分):使用至少三个前向窗口,三个候选共享完全相同的训练与验证日期。
    • 模型(4 分):线性为基线;三次多项式固定次数;自然样条事前固定为 6 个自由度;多项式与样条设计均在折内拟合。
    • 输出(4 分):至少报告 foldmodelmse、训练日期支持域和验证越界比例,并汇总模型的 MSE 均值与离散度。
    • 结论(4 分):若样条改善不稳定或误差区间重叠,结论为“无稳定增量”;即使改善,也只适用于该数据窗口,不说明价格路径存在结构性规律。

    可运行核心答案。 这里的 estimand 是“只根据时间趋势拟合后来日期的后复权价格水平”,损失是价格平方误差。它不预测第 2—9 章的回撤事件,因此其 MSE 不能与 表 7.3 的 Brier 或 AUROC 排名。

列表 7.24: 练习5:后复权价格水平与共同前向折数据
import os  # 从学生显式配置读取练习所需数据根
from pathlib import Path  # 在读取前核验后复权行情文件
import numpy as np  # 构造交易时点索引并计算逐折平方误差
import pandas as pd  # 整理海康威视日期、价格与证据表
from patsy import dmatrix, build_design_matrices  # 在各训练折拟合并复用自然样条设计
from sklearn.linear_model import LinearRegression  # 估计共同折上的连续价格条件均值
from sklearn.model_selection import TimeSeriesSplit  # 生成保持日期顺序的扩张窗口折
from sklearn.pipeline import Pipeline  # 把折内基变换与线性末端绑定
from sklearn.preprocessing import PolynomialFeatures, StandardScaler  # 生成固定三次幂基并只用训练折估计尺度
m07_ex5_root_text = os.environ.get('BOOK_DATA_DIR')  # 区分环境变量缺失与文件缺失
assert m07_ex5_root_text, '请先设置 BOOK_DATA_DIR,使其指向包含 stock/ 子目录的数据根'  # 给出可执行的数据配置提示
m07_ex5_root = Path(m07_ex5_root_text).expanduser().resolve()  # 规范化跨平台数据根
assert m07_ex5_root.is_dir(), f'BOOK_DATA_DIR 不存在或不是目录: {m07_ex5_root}'  # 在拼接文件前验证目录
m07_ex5_path = m07_ex5_root / 'stock' / 'stock_price_post_adjusted.h5'  # 定位题面要求的后复权行情
assert m07_ex5_path.is_file(), f'缺少后复权行情文件: {m07_ex5_path};请检查 BOOK_DATA_DIR'  # 在HDF读取前报告确切缺口
m07_ex5_prices = pd.read_hdf(m07_ex5_path, where="order_book_id='002415.XSHE'", columns=['close']).reset_index()  # 在存储层只读取海康威视收盘价
assert {'date', 'close'}.issubset(m07_ex5_prices.columns), '后复权行情缺少 date 或 close 字段'  # 核对日期与连续响应身份
m07_ex5_prices['date'] = pd.to_datetime(m07_ex5_prices['date'], errors='coerce')  # 把交易日统一为可排序日期
m07_ex5_prices['close'] = pd.to_numeric(m07_ex5_prices['close'], errors='coerce')  # 把后复权收盘价统一为数值响应
m07_ex5_prices = m07_ex5_prices.loc[m07_ex5_prices['date'].notna() & np.isfinite(m07_ex5_prices['close']) & m07_ex5_prices['close'].gt(0)].sort_values('date').reset_index(drop=True)  # 形成有效日期—价格序列
m07_ex5_prices['time_index'] = np.arange(len(m07_ex5_prices), dtype=float)  # 用有序交易时点作为三种趋势模型的共同输入
assert len(m07_ex5_prices) >= 200, '海康威视有效价格少于200条,无法形成三个前向折'  # 保护前向评价的最低样本支持

表 7.5 中三个候选共享完全相同的外层日期折。Pipeline.fit() 只在当前训练折学习多项式设计与尺度;自然样条则把训练折的 design_info 原样用于后来验证日期。

m07_ex5_splitter = TimeSeriesSplit(n_splits=3, gap=1)  # 预先固定三个扩张窗口和一交易时段间隔
m07_ex5_records = []  # 累积每折每模型的损失与支持证据
for fold_number, (fit_indices, score_indices) in enumerate(m07_ex5_splitter.split(m07_ex5_prices), start=1):  # 按日期顺序推进共同外层折
    fit_frame = m07_ex5_prices.iloc[fit_indices].copy()  # 取得当前较早训练日期
    score_frame = m07_ex5_prices.iloc[score_indices].copy()  # 取得当前后来验证日期
    fit_x = fit_frame[['time_index']]  # 构造当前折单一时间输入
    score_x = score_frame[['time_index']]  # 保持验证列名与顺序一致
    fit_y = fit_frame['close'].to_numpy()  # 取得当前折训练价格响应
    score_y = score_frame['close'].to_numpy()  # 取得当前折验证价格响应
    candidate_models = {'linear': Pipeline([('scale', StandardScaler()), ('model', LinearRegression())]), 'cubic_polynomial': Pipeline([('basis', PolynomialFeatures(degree=3, include_bias=False)), ('scale', StandardScaler()), ('model', LinearRegression())])}  # 声明折内线性与三次幂基管道
    support_text = f"{fit_frame['date'].min().date()}--{fit_frame['date'].max().date()}"  # 保存当前训练日期支持区间
    outside_rate = float((~score_frame['time_index'].between(fit_frame['time_index'].min(), fit_frame['time_index'].max())).mean())  # 量化后来日期位于训练时间支持外的比例
    for model_name, candidate_model in candidate_models.items():  # 在完全相同的日期折比较两种参数化
        candidate_model.fit(fit_x, fit_y)  # 只用当前训练折拟合尺度、基函数和系数
        candidate_prediction = candidate_model.predict(score_x)  # 对当前后来日期生成价格预测
        m07_ex5_records.append({'fold': fold_number, 'model': model_name, 'mse': float(np.mean((score_y - candidate_prediction) ** 2)), 'support': support_text, 'validation_outside_rate': outside_rate})  # 保存逐折MSE与支持证据
    spline_fit = dmatrix('cr(time_index, df=6)', fit_frame, return_type='dataframe')  # 只用当前训练折确定自然样条边界与基
    spline_score = build_design_matrices([spline_fit.design_info], score_frame)[0]  # 复用训练折设计转换后来验证日期
    spline_model = LinearRegression(fit_intercept=False).fit(spline_fit.to_numpy(), fit_y)  # 以无名称矩阵在训练折自然样条基上估计系数
    spline_prediction = spline_model.predict(np.asarray(spline_score))  # 用同类无名称矩阵对验证折生成自然边界预测
    m07_ex5_records.append({'fold': fold_number, 'model': 'natural_spline_df6', 'mse': float(np.mean((score_y - spline_prediction) ** 2)), 'support': support_text, 'validation_outside_rate': outside_rate})  # 保存自然样条同折证据
m07_ex5_fold_table = pd.DataFrame(m07_ex5_records).sort_values(['fold', 'model'])  # 对齐三折三模型的审计顺序
m07_ex5_fold_table  # 渲染至少含fold、model、mse和support的核心结果
表 7.5: 练习5:线性、三次多项式与自然样条的共同前向折MSE和支持域
fold model mse support validation_outside_rate
1 1 cubic_polynomial 1.969245e+05 2010-05-28--2014-04-28 1.0
0 1 linear 4.036792e+04 2010-05-28--2014-04-28 1.0
2 1 natural_spline_df6 1.261129e+05 2010-05-28--2014-04-28 1.0
4 2 cubic_polynomial 1.636529e+06 2010-05-28--2018-03-14 1.0
3 2 linear 7.497180e+04 2010-05-28--2018-03-14 1.0
5 2 natural_spline_df6 8.595579e+05 2010-05-28--2018-03-14 1.0
7 3 cubic_polynomial 1.023564e+06 2010-05-28--2022-02-08 1.0
6 3 linear 2.149537e+05 2010-05-28--2022-02-08 1.0
8 3 natural_spline_df6 1.576648e+06 2010-05-28--2022-02-08 1.0
m07_ex5_summary = m07_ex5_fold_table.groupby('model')['mse'].agg(['mean', 'std', 'count']).reset_index()  # 汇总三折误差水平与离散度
m07_ex5_summary  # 渲染条件化选择所需的均值、标准差与折数
表 7.6: 练习5:三种连续响应模型的前向MSE汇总
model mean std count
0 cubic_polynomial 952339.026680 722440.333780 3
1 linear 110097.808747 92441.481375 3
2 natural_spline_df6 854106.429320 725283.131163 3

先逐行读取 表 7.5 的日期支持与误差,再用 表 7.6 比较三折均值和离散度。所有验证日期都晚于相应训练日期,因此时间输入的 validation_outside_rate 通常为 1;这正是趋势外推任务的风险证据,而不是代码失败。只有当某个非线性候选在三个共同折上均稳定降低 MSE,且边界预测没有数值发散,才可说它在这一连续价格任务上提供增量;否则选择线性基线。无论结果如何,都不能据此声称二元回撤概率的逻辑样条胜出或失败。

  1. 可设 \[\operatorname{logit}P(Y_{i,t+1}=1) =\beta_0+f_1(\text{leverage}_{it})+f_2(\text{size}_{it}) +f_3(\text{cashflow}_{it}),\] 其中 \(Y_{i,t+1}\) 表示下一报告期经营现金流为负。财务特征按实际披露日对齐,只使用预测时点已知信息。若理论上认为杠杆与规模存在非可加交互,应把二维平滑列为单独候选,而不是从一维偏依赖推断交互。报告线性 Logit 基线、前向 Brier/AUROC、每个平滑项的支持密度和置信带;曲线表示模型内条件预测形状,不是提高杠杆的因果效应。

  2. 示范答案与评分量规(15 分)

    • 在相同前向折上拟合两个模型,并用相同损失比较(4 分)。
    • 二维平滑能表示交互,但基函数数量增长更快,需要更密集的联合支持;可加模型方差较低,每个分量更易解释,却无法表示交互(4 分)。
    • 报告逐折误差差、二维支持热图和一维分量图;若联合支持稀疏,二维曲面的局部形状不可解释(4 分)。
    • 结论按验证证据条件化:只有误差改善稳定且发生在有数据支持的区域,才保留二维模型;否则选择可加模型(3 分)。

7.14.3 理论题参考解答

  1. 无约束时,\(K+1\) 个区间各有四个三次多项式系数,共 \(4(K+1)\) 个。每个内部节点要求函数值、一阶导数和二阶导数连续,提供 \(3K\) 个独立约束,因此普通三次样条维数为 \[4(K+1)-3K=K+4.\] 自然样条再施加左右边界二阶导数为零的两个独立约束,所以维数为 \[K+4-2=K+2.\]

  2. 平滑样条最小化 \[\sum_i[y_i-g(x_i)]^2+\lambda\int[g''(t)]^2dt.\]\(\lambda\to0\) 时,粗糙度几乎不受罚,拟合趋向穿过训练点的自然三次插值样条,有效自由度增大、方差可能很高;当 \(\lambda\to\infty\) 时,有限目标要求 \(g''(t)\to0\),所以 \(g\) 趋向最小二乘直线。有限样本中的最佳 \(\lambda\) 取决于验证损失,不由这两个极限直接决定。

7.15 章末回顾

本章的方法没有无条件赢家。低次多项式适合简单全局形状,自然样条适合受控的局部弯曲,局部回归依赖低维密集支持,GAM 依赖可加结构。节点、平滑参数和自由度必须区分;任何经验范围都只是候选网格起点,最终选择应服从时间验证、外推边界、解释需求和业务损失。预测验证不自动提供形状推断:固定设计、选择后设计、惩罚估计与时间依赖各自需要匹配的不确定性方法。

无提示检索

  1. 自然样条在边界外施加了什么约束,它为什么仍不消除外推风险?
  2. GCV 与前向验证回答的问题有何不同,金融预测为何以后者作为部署证据?
  3. GAM 的单变量偏效应图为什么不能识别交互,更不能作因果解释?
展开检索反馈与定向返回

若第 1 题不确定,返回“回归样条与自然样条”;若第 2 题不确定,返回“平滑样条”并对照本章时间验证;若第 3 题不确定,返回“广义加性模型”的可加性与识别边界。

下一章将让树模型通过递归分裂学习非线性与交互,同时保持相同的前向验证边界。

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