3  线性回归 (Linear Regression)

3.1 本章学习契约

先修：掌握第 2 章的条件关联、MSE 和训练/测试边界；会基础矩阵运算和财务比率构造。

学习目标与达成标准

  1. 从 RSS 推导简单 OLS 和正规方程；标准是写出一阶条件和维度。
  2. 正确解释多元系数、区间和 \(p\) 值；标准是包含“其他变量不变”与非因果边界。
  3. 解释 \(R^2\)、RSE 和样本外损失的不同角色；标准是不用单一指标选模。
  4. 识别异方差、非线性、异常点和共线性；标准是为残差证据匹配合适补救。
  5. 解释交互项的条件边际关联；标准是能在陌生财务变量上写出并解释模型。

3.1.1 入口检查与补修

入口题：若 \(X\) 为 \(n\times p\)，\(Y\) 为 \(n\times1\)，\((X^TX)^{-1}X^TY\) 的维度是什么？系数显著是否自动意味因果？

作答后展开答案、门槛与补修

答案与门槛：维度为 \(p\times1\)；显著性不自动提供因果识别。两点都对才通过；否则补修矩阵乘法维度和第 2 章“预测与推断”后重测。

3.1.2 低风险检索、渐隐链与迁移

检索：闭卷写出正规方程，并用一句话说明 VIF 衡量什么。

渐隐链：先在给定设计矩阵上计算 OLS；再只给残差图选择诊断；最后在陌生的门店销量场景中自行定义交互项、估计对象和结论边界。

3.1.3 目标—评价证据映射

目标 学习活动 评价证据
1 正规方程检索 练习 7—8
2—4 系数解释与诊断渐隐链 练习 1—5
5 陌生情境交互迁移 练习 6、小节 3.8

3.2 引言 (Introduction)

线性回归(linear regression)是统计学习和应用统计学中最为基础、也是最广泛使用的监督学习方法之一。尽管现代机器学习领域涌现了诸如深度学习、梯度提升树等复杂模型，但线性回归凭借其强大的可解释性、理论完备性以及计算上的高效性，依然在高水平学术研究和工程实践中占据核心地位。

3.2.1 线性回归在经济金融领域的典型应用

线性回归在经济金融领域有着广泛而深入的应用，以下列举几个最具代表性的前沿应用场景，以帮助读者理解这一基础工具的强大威力。

应用一：资产定价与因子模型。 CAPM 时间序列回归将单个资产的超额收益对市场超额收益回归，斜率 \(\beta\) 衡量市场因子暴露。Fama–French 三因子的经验实现通常以资产超额收益为因变量，以市场、SMB 和 HML 因子收益为自变量做时间序列回归；截面定价是另一层检验，不应与时间序列因子暴露混为一谈。

应用二：宏观经济预测与政策评估。 泰勒规则是以通胀偏离目标和产出缺口描述政策利率的反应函数，不是“GDP增长、通胀和货币供应量预测利率”的泛称。政策效果若要做因果解释，还需要额外的识别设计。

应用三：公司财务分析与信用评级。 Altman Z-score 原型是线性判别分数，用若干财务比率的线性组合做破产分类；它不是用 OLS 直接估计违约概率的线性回归模型。

应用四：房地产估值与特征价格模型。 Rosen (1974) 提出的特征价格模型(Hedonic Pricing Model)用属性（面积、楼龄、区位、学区、交通便利性等）的条件价格关联描述房屋成交价。若要量化地铁开通的因果效应，还需事件时间、对照组与可辩护的识别设计，不能由普通横截面 OLS 直接得到。

以上这些应用场景清楚地表明：线性回归不仅是一种统计技术，更是理解经济金融世界运行规律的基本思维工具。 掌握了线性回归，就掌握了进入量化分析和实证研究领域的钥匙。本章将从最基础的简单线性回归开始，逐步构建起完整的多元回归分析框架。

Note: 线性回归的理论基石

高斯-马尔可夫定理(Gauss-Markov Theorem)为线性回归提供了严谨的理论支持。该定理指出：在经典线性模型假设(Classical Linear Model Assumptions)下，OLS估计量是所有线性无偏估计量中方差最小的，即它是最佳线性无偏估计量(Best Linear Unbiased Estimator, BLUE)。

具体假设包括：

  1. 线性性：模型相对于参数是线性的。
  2. 零均值：误差项的期望值为零，\(E(\epsilon|X) = 0\)。
  3. 同方差：所有观测值的误差项具有相同的方差，\(\text{Var}(\epsilon|X) = \sigma^2\)。
  4. 无自相关：不同观测值的误差项互不相关。
  5. 无完全共线性：预测变量之间不存在精确的线性关系。

当这些假设满足时，OLS 在线性无偏估计量类中方差最小。这不是对所有可能估计量的无条件“最高效”声明；若存在异方差或自相关，经典方差公式与 BLUE 结论都需修改。

假设—结果对照：线性参数化和无完全共线性使 OLS 解唯一；\(E(\epsilon\mid X)=0\) 给出条件无偏性；再加同方差与无自相关才得到 Gauss–Markov 的 BLUE 结果。有限样本的精确 \(t/F\) 分布还需条件正态性；大样本渐近推断则依赖适当的矩、依赖与抽样条件，并应根据异方差或时序相关使用稳健标准误。

在本书第2章中，我们讨论了如何使用统计学习方法来估计函数\(f\)，该函数将输入变量\(X\)映射到输出变量\(Y\)。线性回归假设\(f\)的形式是线性的：

\[ Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \cdots + \beta_p X_p + \epsilon \tag{3.1}\]

其中：

  • \(Y\)是响应变量（定量变量），如公司的营业收入或股票收益率。
  • \(X_1, X_2, \ldots, X_p\)是预测变量，如研发投入、宏观经济指标。
  • \(\beta_0\) 是截距项，\(\beta_1, \ldots, \beta_p\)是回归系数；在模型设定成立时，它们描述保持其他变量不变的条件均值变化。只有额外因果识别假设成立时才可解释为干预效应。
  • \(\epsilon\) 是误差项，捕捉了所有未包含在模型中的随机因素。

案例背景：长三角金融科技公司

在本章中，我们继续使用第2章引入的金融科技公司案例。该公司拥有200家制造业上市公司的数据，包括它们的：

  • 研发费用（R&D）：用于技术创新和产品开发的费用（万元）
  • 销售费用（Sales）：用于市场推广和销售团队的费用（万元）
  • 管理费用（Admin）：用于行政管理和运营的费用（万元）

我们的目标是描述这些费用与公司营业收入（Revenue）的条件关联，并建立预测模型。任何资源配置建议都需要另行评估干预成本与因果效果。

3.3 简单线性回归 (Simple Linear Regression)

简单线性回归只涉及一个预测变量\(X\)和一个响应变量\(Y\)，假设它们之间的关系为：

\[ Y = \beta_0 + \beta_1 X + \epsilon \tag{3.2}\]

其中：

  • \(\beta_0\)是截距(intercept)：当\(X = 0\)时\(Y\)的平均值
  • \(\beta_1\)是斜率(slope)：\(X\)每增加1单位时\(Y\)的平均变化量
  • \(\epsilon\)是误差项，表示\(Y\)中无法被\(X\)解释的部分

3.3.1 估计系数 (Estimating the Coefficients)

假设我们有\(n\)个观测\((x_1, y_1), (x_2, y_2), \ldots, (x_n, y_n)\)。我们的目标是估计\(\beta_0\)和\(\beta_1\)，使得拟合的直线\(\hat{y} = \hat{\beta}_0 + \hat{\beta}_1 x\)尽可能接近实际观测点。

最小二乘法(Ordinary Least Squares, OLS)通过最小化残差平方和(Residual Sum of Squares, RSS)来寻找最优参数：

\[ \text{RSS}(\beta_0, \beta_1) = \sum_{i=1}^{n} (y_i - \hat{y}_i)^2 = \sum_{i=1}^{n} (y_i - \beta_0 - \beta_1 x_i)^2 \tag{3.3}\]

3.3.1.1 数学推导：一阶条件 (First Order Conditions)

为了找到使 RSS 最小的 \(\beta_0\) 和 \(\beta_1\)，我们分别对这两个参数求偏导，并令其等于 0：

  1. 对 \(\beta_0\) 求导： \[ \frac{\partial \text{RSS}}{\partial \beta_0} = -2 \sum_{i=1}^{n} (y_i - \beta_0 - \beta_1 x_i) = 0 \] 由此可得：\(\sum y_i - n\hat{\beta}_0 - \hat{\beta}_1 \sum x_i = 0\)，即： \[ \hat{\beta}_0 = \bar{y} - \hat{\beta}_1 \bar{x} \tag{3.4}\]

  2. 对 \(\beta_1\) 求导： \[ \frac{\partial \text{RSS}}{\partial \beta_1} = -2 \sum_{i=1}^{n} (y_i - \beta_0 - \beta_1 x_i)x_i = 0 \] 代入 \(\hat{\beta}_0\) 的表达式： \[ \sum_{i=1}^{n} (y_i - (\bar{y} - \hat{\beta}_1 \bar{x}) - \hat{\beta}_1 x_i)x_i = 0 \] \[ \sum (y_i - \bar{y})x_i - \hat{\beta}_1 \sum (x_i - \bar{x})x_i = 0 \] 最终得到斜率的估计公式： \[ \hat{\beta}_1 = \frac{\sum_{i=1}^{n}(x_i - \bar{x})(y_i - \bar{y})}{\sum_{i=1}^{n}(x_i - \bar{x})^2} \tag{3.5}\]

其中，\(\bar{x} = \frac{1}{n}\sum x_i\) 和 \(\bar{y} = \frac{1}{n}\sum y_i\) 分别是样本均值。

案例应用：销售费用与营业收入的关联

让我们首先研究销售费用与公司营业收入的条件关联。我们将使用本地存储的 A 股上市公司财务报表数据进行演示。为了确保口径一致，我们选择 2023 年度的年报数据；这项横截面分析不识别销售费用的因果效应。

为了将最小二乘法的理论与现实金融数据分析结合起来，我们在接下来的代码中，利用 Python 的 pandas 库从本地读取了 A 股制造业上市公司的真实财务测算数据。在这段代码里，我们专门提取了“销售费用”与“营业收入”这两个核心字段，并在清洗掉异常极小值之后，对它们进行了对数转换（Log-Transformation）。之所以在这里使用对数变换，是因为真实的财务绝对值往往呈现非常强烈的偏斜分布（少数巨无霸公司的数据会严重压扁图表），转化后不仅图形更美观，而且根据计量经济学常识，双对数线性回归的拟合斜率 (\(\hat{\beta}_1\)) 能够直接代表并衡量两个变量之间的弹性。代码的后半部分调用 statsmodels 库执行了核心的 OLS 拟合，打印出了精确的回归分析检验报告，并利用 matplotlib 将散点图和那条贯穿其中的耀眼的红色 OLS 拟合线直观地展示出来。

import numpy as np  # 导入numpy用于数组运算与多项式拟合
import pandas as pd  # 导入pandas用于读取和处理h5数据
import matplotlib.pyplot as plt  # 导入matplotlib用于科学绑图
import statsmodels.api as sm  # 导入statsmodels用于OLS回归拟合
from scipy import stats  # 导入scipy.stats用于统计分布函数（备用）
import os  # 导入os模块用于跨平台路径处理
from pathlib import Path  # 用路径对象解析数据根目录

# 设置出版用中文字体
# 配置全局绑图参数：指定中文字体族与负号渲染
plt.rcParams['font.family'] = ['Source Han Serif SC']
plt.rcParams['axes.unicode_minus'] = False  # 解决负号显示为方块的问题

接下来，我们从本地路径读取 A 股上市公司财务报表数据，筛选 2023 年年报中的营业收入与销售费用字段，清洗异常值后进行对数变换，为后续简单线性回归分析做好准备。

# 按项目的跨平台数据根目录约定定位财务报表文件
BOOK_DATA_DIR = Path(os.environ['BOOK_DATA_DIR']).expanduser().resolve()  # 从必需环境变量解析数据根目录
assert BOOK_DATA_DIR.is_dir(), f'BOOK_DATA_DIR 不存在: {BOOK_DATA_DIR}'  # 在读取前验证目录
local_data_path = BOOK_DATA_DIR / 'stock/financial_statement.h5'  # 拼接财务报表路径

df_2023_report = pd.read_hdf(  # 从本地财务报表h5文件中选择性读取2023年年报所需字段，避免全量载入全部季度和全部列
    local_data_path,  # 指定财务报表h5文件路径
    where="quarter='2023q4'",  # 仅保留2023年第四季度年报数据
    columns=['quarter', 'revenue', 'selling_expense']  # 只读取本节简单线性回归分析需要的季度、营收和销售费用字段
).copy()  # 创建独立副本以安全执行后续筛选和对数变换

target_variables = ['revenue', 'selling_expense']  # 目标变量：营业收入和销售费用
df_analysis = df_2023_report[target_variables].dropna()  # 剔除含缺失值的记录

# 排除营收低于1000万或销售费用低于100万的公司（保证对数变换有效）
df_analysis = df_analysis[(df_analysis['revenue'] > 1e7) & (df_analysis['selling_expense'] > 1e6)]  # 过滤极端小值

df_analysis['log_revenue'] = np.log10(df_analysis['revenue'])  # 对营业收入取常用对数
df_analysis['log_sell_exp'] = np.log10(df_analysis['selling_expense'])  # 对销售费用取常用对数

# 当样本量超过250时随机抽样，保持图表简洁清晰
if len(df_analysis) > 250:  # 判断清洗后样本量是否超过可视化阈值
    df_plot_sample = df_analysis.sample(250, random_state=42)  # 固定种子抽样250个观测
else:  # 样本量不足250的情形
    df_plot_sample = df_analysis  # 全量使用所有样本

x_input_data = df_plot_sample['log_sell_exp']  # 自变量：对数销售费用
y_output_data = df_plot_sample['log_revenue']  # 因变量：对数营业收入
x_axis_label = 'Log10 销售费用'  # X轴标签文字
y_axis_label = 'Log10 营业收入'  # Y轴标签文字

以上代码完成了数据准备和对数变换。接下来我们使用OLS拟合简单线性回归模型，并绘制散点图和拟合直线，直观展示销售费用与营业收入之间的弹性关系。

# 构建含截距的自变量矩阵用于OLS回归
marketing_features_matrix = sm.add_constant(x_input_data)  # 为X添加常数列（截距项）
ols_model_fit = sm.OLS(y_output_data, marketing_features_matrix).fit()  # 最小二乘法拟合回归模型

estimated_beta0 = ols_model_fit.params[0]  # 提取截距估计值β₀
estimated_beta1 = ols_model_fit.params[1]  # 提取斜率估计值β₁（即弹性系数）
r_squared_value = ols_model_fit.rsquared  # 提取模型拟合优度R²

print('--- 简单线性回归分析报告 ---')  # 输出分析标题
print(f'估计截距 (Beta0): {estimated_beta0:.4f}')  # 输出截距
print(f'估计斜率 (Beta1): {estimated_beta1:.4f}')  # 输出斜率
print(f'拟合优度 (R-squared): {r_squared_value:.4f}')  # 输出R²
print(f'斜率 P 值: {ols_model_fit.pvalues[1]:.4e}')  # 输出斜率的假设检验p值
print(f'\n经济解释：在双对数模型下，斜率即为弹性。')  # 经济学解读
print(f'意味着销售费用每增加 1%，预期营业收入将增加约 {estimated_beta1:.2f}%。')  # 弹性含义
fig_simple_reg, ax_simple_reg = plt.subplots(figsize=(11, 7))  # 创建画布
# 绘制原始散点（每个点代表一家上市公司）
ax_simple_reg.scatter(x_input_data, y_output_data, alpha=0.5, s=50,  # 在子图中绑制散点图
                     color='#34495e', edgecolors='white', label='样本上市公司')  # 绘制散点图，灰蓝色(#34495e)标注各上市公司观测值

x_line_range = np.linspace(x_input_data.min(), x_input_data.max(), 100)  # 生成拟合线的X坐标序列
y_line_predicted = ols_model_fit.predict(sm.add_constant(x_line_range))  # 计算拟合线上的预测值
ax_simple_reg.plot(x_line_range, y_line_predicted, color='#e74c3c',  # 在子图中绑制折线图
                  linewidth=3, label=f'OLS 拟合线 (R²={r_squared_value:.2f})')  # 绘制红色拟合直线

ax_simple_reg.set_xlabel(x_axis_label, fontsize=12)  # 设置X轴标签
ax_simple_reg.set_ylabel(y_axis_label, fontsize=12)  # 设置Y轴标签
ax_simple_reg.set_title('中国上市制造企业：销售费用与营业收入的关联',
                       fontsize=14, fontweight='bold', pad=20)  # 设置图标题
ax_simple_reg.legend(loc='upper left', frameon=True)  # 添加左上角图例
ax_simple_reg.grid(True, linestyle='--', alpha=0.6)  # 添加虚线网格

plt.tight_layout()  # 自动调整布局
plt.show()  # 显示图表
2023 年年报销售费用与营业收入的双对数散点，叠加红色 OLS 拟合直线。
图 3.1: 销售费用与营业收入的散点图及最小二乘拟合线（对数尺度）。数据显示了明显的正相关关系：销售投入的增加通常伴随着营业收入的增长，斜率代表了两者的弹性关系。数据来源：本地财务报表数据库。

从 图 3.1 的结果可以看出：

  1. 斜率估计：在双对数条件均值模型成立时，\(\hat{\beta}_1\) 是样本中的条件关联弹性。例如估计值为 0.9 时，销售费用高 1% 的公司，其拟合营收平均高约 0.9%；这不是“增加费用会使营收增加”的干预结论。

  2. 统计显著性：若按所用标准误计算的 \(p\) 值很小，只能拒绝该模型中的零斜率假设。遗漏的公司规模、行业与反向决定关系都可能解释关联，因此不能据此断言更多推广导致更高收入。

3.3.2 评估系数估计的准确性 (Assessing the Accuracy of the Coefficient Estimates)

回忆简单线性回归模型：

\[ Y = \beta_0 + \beta_1 X + \epsilon \]

我们在估计\(\beta_0\)和\(\beta_1\)时使用的是有限样本，因此估计值\(\hat{\beta}_0\)和\(\hat{\beta}_1\)会与真实值有所偏差。为了评估估计的准确性，我们需要了解估计量的抽样分布(sampling distribution)。

在经典假设下（误差项独立同分布，\(E(\epsilon) = 0\)，\(\text{Var}(\epsilon) = \sigma^2\)），OLS估计量具有以下性质：

  1. 无偏性 (Unbiasedness)：\(E(\hat{\beta}_0) = \beta_0\)，\(E(\hat{\beta}_1) = \beta_1\)

  2. 方差 (Variance)： \[ \text{Var}(\hat{\beta}_0) = \sigma^2 \left[\frac{1}{n} + \frac{\bar{x}^2}{\sum_{i=1}^{n}(x_i - \bar{x})^2}\right] \tag{3.6}\] \[ \text{Var}(\hat{\beta}_1) = \frac{\sigma^2}{\sum_{i=1}^{n}(x_i - \bar{x})^2} \tag{3.7}\]

  3. 标准误差 (Standard Error)： \[ \text{SE}(\hat{\beta}_0) = \sqrt{\frac{\sigma^2}{n} + \frac{\sigma^2 \bar{x}^2}{\sum_{i=1}^{n}(x_i - \bar{x})^2}} \tag{3.8}\] \[ \text{SE}(\hat{\beta}_1) = \sqrt{\frac{\sigma^2}{\sum_{i=1}^{n}(x_i - \bar{x})^2}} \tag{3.9}\]

其中，\(\sigma^2\)通常使用残差标准误(Residual Standard Error, RSE)来估计：

\[ \hat{\sigma}^2 = \text{RSE}^2 = \frac{1}{n-2}\sum_{i=1}^{n}(y_i - \hat{y}_i)^2 = \frac{\text{RSS}}{n-2} \tag{3.10}\]

Tip: 为什么除以\(n-2\)而不是\(n\)？

在估计\(\sigma^2\)时，我们使用\(n-2\)而不是\(n\)作为分母，这是为了得到\(\sigma^2\)的无偏估计。原因如下：

  1. 我们在估计\(\hat{\beta}_0\)和\(\hat{\beta}_1\)时，使用了数据的两个自由度(degrees of freedom)。
  2. 这使得残差\((y_i - \hat{y}_i)\)的方差略微小于真实误差\(\epsilon_i\)的方差。
  3. 除以\(n-2\)可以补偿这种低估，得到无偏估计。

更一般地，对于\(p\)个参数的模型，我们除以\(n-p-1\)。这是贝塞尔校正(Bessel’s correction)。

假设检验 (Hypothesis Testing)

最常见的假设检验是检验斜率是否显著不为零：

  • 零假设 \(H_0: \beta_1 = 0\)（\(X\)和\(Y\)之间没有线性关系）
  • 备择假设 \(H_a: \beta_1 \neq 0\)（\(X\)和\(Y\)之间存在线性关系）

t统计量 (t-statistic)：

\[ t = \frac{\hat{\beta}_1 - 0}{\text{SE}(\hat{\beta}_1)} \tag{3.11}\]

在误差项条件正态且\(H_0\)为真时，\(t\)精确服从自由度为\(n-2\)的\(t\)分布。非正态情形下，下述结论通常是大样本近似，或需稳健/自助推断。

\[ p\text{-value} = P(|T_{n-2}| \geq |t|) \tag{3.12}\]

如果p值小于显著性水平（通常为0.05），我们拒绝\(H_0\)，认为\(\beta_1\)显著不为零。

置信区间 (Confidence Interval)

\(\beta_1\)的\(95\%\)置信区间为：

\[ \hat{\beta}_1 \pm 2 \cdot \text{SE}(\hat{\beta}_1) \tag{3.13}\]

更一般地，置信水平为\(1-\alpha\)的置信区间为：

\[ \hat{\beta}_1 \pm t_{\alpha/2, n-2} \cdot \text{SE}(\hat{\beta}_1) \tag{3.14}\]

其中，\(t_{\alpha/2, n-2}\)是自由度为\(n-2\)的\(t\)分布的上\(\alpha/2\)分位数。

当我们完成了线性模型的拟合之后，评估各变量系数的准确性和显著性是计量分析的必经之路。下面的代码段简明地演示了这一流程。在这里我们构建了一个严密的模拟数据集，假设有 200 个分公司的电视广告预算和随之产生的销售收入。通过将这些特征矩阵喂给 statsmodels.api.OLS 模型进行拟合后，只需要调用其自带的 .summary() 方法，Python 便会为你输出一份极其详尽的、堪比任何专业统计学术软件（如 Stata 或 SPSS）的标准回归分析摘要报告。这份报告中包含了我们上文推导的所有关键统计量——不仅有每一项的估计系数 \(\hat{\beta}\)，还有对应的标准误 (Std. Err.)、用来进行假设检验的 \(t\) 统计量、决定是否拒绝零假设的关键 \(P\) 值，以及估计区间的 95% 置信区间边界。

表 3.1: 简单线性回归（销售量 ~ 电视广告）的完整统计摘要 
import statsmodels.api as sm  # 导入statsmodels用于OLS回归分析
import numpy as np  # 导入numpy用于数值计算
import pandas as pd  # 导入pandas（备用数据处理）

np.random.seed(42)  # 设置随机种子保证结果可复现
market_count = 200  # 模拟200个区域市场的广告投放数据
tv_ad_budget = np.random.uniform(0, 300, market_count)  # 均匀分布模拟电视广告预算（0~300万元）
# DGP: 销量 = 5 + 0.04*广告 + ε，其中ε~N(0,4)
sales_revenue = 5 + 0.04 * tv_ad_budget + np.random.normal(0, 2, market_count)  # 真实线性关系+高斯噪声

ad_features_matrix = sm.add_constant(tv_ad_budget)  # 为自变量矩阵添加截距列（常数1）
sales_linear_model = sm.OLS(sales_revenue, ad_features_matrix).fit()  # 最小二乘法拟合简单线性回归模型

print(sales_linear_model.summary())  # 输出包含系数、标准误、t值、p值、置信区间的完整回归报告
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.765
Model:                            OLS   Adj. R-squared:                  0.763
Method:                 Least Squares   F-statistic:                     643.4
Date:                Thu, 13 Aug 2026   Prob (F-statistic):           4.04e-64
Time:                        14:11:54   Log-Likelihood:                -415.57
No. Observations:                 200   AIC:                             835.1
Df Residuals:                     198   BIC:                             841.7
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          5.2104      0.264     19.702      0.000       4.689       5.732
x1             0.0395      0.002     25.365      0.000       0.036       0.043
==============================================================================
Omnibus:                        7.028   Durbin-Watson:                   2.131
Prob(Omnibus):                  0.030   Jarque-Bera (JB):                9.199
Skew:                           0.231   Prob(JB):                       0.0101
Kurtosis:                       3.943   Cond. No.                         327.
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

表 3.1 提供了丰富的统计信息：

  • Coef.：系数估计值
  • Std.Err.：标准误差
  • t：t统计量
  • P>|t|：p值（两个尾部的概率）
  • [0.025 0.975]：95%置信区间

3.3.3 评估模型的准确性 (Assessing the Accuracy of the Model)

在估计回归系数后，我们需要评估整体模型的拟合优度。两个最重要的指标是残差标准误(RSE)和\(R^2\)统计量。

3.3.3.1 残差标准误 (Residual Standard Error)

RSE是误差项标准差\(\sigma\)的估计：

\[ \text{RSE} = \sqrt{\frac{1}{n-2}\sum_{i=1}^{n}(y_i - \hat{y}_i)^2} \tag{3.15}\]

RSE 是对误差项条件标准差 \(\sigma\) 的估计，其单位与 \(Y\) 相同。它不是平均绝对偏差（MAE），也不是平均残差（含截距 OLS 的平均残差为零）。

在我们的案例中，如果 RSE = 2.0，可说拟合残差的典型标准差尺度约为2千件；不能把它精确翻译成“平均偏差”。

3.3.3.2 \(R^2\)统计量 (\(R^2\) Statistic)

\(R^2\)统计量衡量的是响应变量的变异中可以被预测变量解释的比例：

\[ R^2 = \frac{\text{TSS} - \text{RSS}}{\text{TSS}} = 1 - \frac{\text{RSS}}{\text{TSS}} \tag{3.16}\]

其中：

  • \(\text{TSS} = \sum_{i=1}^{n}(y_i - \bar{y})^2\)是总平方和(Total Sum of Squares)，衡量\(Y\)的总变异
  • \(\text{RSS} = \sum_{i=1}^{n}(y_i - \hat{y}_i)^2\)是残差平方和(Residual Sum of Squares)，衡量未被模型解释的变异

\(R^2\)的取值范围是\([0, 1]\)：

  • \(R^2 = 1\)：模型完美拟合数据
  • \(R^2 = 0\)：模型无法解释任何变异（等同于只用均值预测）

在简单线性回归中，\(R^2\)等于\(X\)和\(Y\)的相关系数的平方：

\[ R^2 = \text{Cor}(X, Y)^2 \tag{3.17}\]

Clarification on \(R^2\)的局限性

虽然\(R^2\)是广泛使用的模型拟合指标，但它有一些重要的局限性：

  1. \(R^2\)不能判断模型是否正确：高\(R^2\)并不意味着模型设定正确，可能是遗漏变量问题或伪相关。

  2. \(R^2\)会随着预测变量数量增加而增加：在多元回归中，添加任何预测变量（即使是无用的变量）都会使\(R^2\)增加或保持不变。这就是为什么我们通常使用调整\(R^2\)(Adjusted \(R^2\))。

  3. \(R^2\)在非线性模型中的解释不同：在非线性模型中，\(R^2\)不一定代表”解释变异的比例”。

  4. 高\(R^2\)不等于因果性：\(R^2\)只是衡量相关性，不能推出因果关系。

因此，\(R^2\)应该与其他指标（如残差分析、交叉验证性能）结合使用，而不是单独依赖它来评估模型质量。

3.4 多元线性回归 (Multiple Linear Regression)

当有多个预测变量时，我们使用多元线性回归：

\[ Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \cdots + \beta_p X_p + \epsilon \tag{3.18}\]

其中：

  • \(X_1, X_2, \ldots, X_p\)是\(p\)个不同的预测变量
  • \(\beta_1, \beta_2, \ldots, \beta_p\)是相应的回归系数，表示在其他变量不变的情况下，\(X_j\)每增加1单位时\(Y\)的平均变化量

案例应用：三种费用对营业收入的综合影响

在金融科技案例中，我们通常需要同时考虑多种费用的协同作用。这里我们分析：研发费用、销售费用及管理费用对营业收入的综合贡献。

下面是一个明确标注的受控模拟，用于演示多元回归几何。标准化使系数表示“\(X_j\) 增加一个样本标准差时 \(Y\) 的条件关联变化”，但系数绝对值仍会受共线性、测量误差和变量编码影响，不能自动解释为预测重要性或因果贡献。

import numpy as np  # 导入numpy用于数值运算
import pandas as pd  # 导入pandas用于数据框操作
import matplotlib.pyplot as plt  # 导入matplotlib用于绑图
from mpl_toolkits.mplot3d import Axes3D  # 导入3D绑图引擎
from sklearn.linear_model import LinearRegression  # 导入线性回归模型
from sklearn.preprocessing import StandardScaler  # 导入标准化工具

np.random.seed(42)  # 设置随机种子保证可复现
n_company_samples = 250  # 模拟250家企业样本

# 模拟三种费用数据（单位：万元）
rd_exp_input = np.random.uniform(100, 1000, n_company_samples)  # 研发费用范围100~1000万
sales_exp_input = np.random.uniform(80, 800, n_company_samples)  # 销售费用范围80~800万
admin_exp_input = np.random.uniform(50, 500, n_company_samples)  # 管理费用范围50~500万

# 构造真实数据生成过程(DGP): Revenue = 200 + 1.5*RD + 2.2*Sales + 0.8*Admin + ε
true_revenue_values = (200 + 1.5 * rd_exp_input +  # 研发费用贡献
                       2.2 * sales_exp_input +  # 销售费用贡献（弹性最大）
                       0.8 * admin_exp_input)  # 管理费用贡献（弹性最小）
observed_revenue = true_revenue_values + np.random.normal(0, 150, n_company_samples)  # 加入随机噪声

模拟数据就绪后，我们使用 scikit-learn 对三种费用进行标准化处理和多元线性回归拟合，最终将结果可视化为3D回归平面。

# 整合为DataFrame
df_multi_finance = pd.DataFrame({  # 构建包含所有变量的数据框
    'RD_Exp': rd_exp_input,  # 研发费用列
    'Sales_Exp': sales_exp_input,  # 销售费用列
    'Admin_Exp': admin_exp_input,  # 管理费用列
    'Revenue': observed_revenue  # 营业收入列（响应变量）
})  # 完成四列数据框构建，用于后续标准化与建模

feature_columns = ['RD_Exp', 'Sales_Exp', 'Admin_Exp']  # 定义自变量列名列表
X_features = df_multi_finance[feature_columns]  # 提取自变量矩阵
y_target_obs = df_multi_finance['Revenue']  # 提取因变量向量

scaler_tool = StandardScaler()  # 实例化Z-score标准化器
X_features_scaled = scaler_tool.fit_transform(X_features)  # 对特征矩阵做标准化（均值=0，标准差=1）

multi_ols_reg = LinearRegression()  # 实例化线性回归模型
multi_ols_reg.fit(X_features_scaled, y_target_obs)  # 在标准化特征上拟合OLS模型

print('--- 多元线性回归分析结果 ---')  # 输出报告标题
print(f'回归截距 (Beta0): {multi_ols_reg.intercept_:.4f}')  # 输出截距估计值
for name, coef in zip(feature_columns, multi_ols_reg.coef_):  # 遍历输出各特征的标准化系数
    print(f'{name} 的标准化回归系数: {coef:.4f}')  # 报告条件关联尺度，不把它当作因果重要性
print(f'模型 R-squared: {multi_ols_reg.score(X_features_scaled, y_target_obs):.4f}')  # 输出R²
--- 多元线性回归分析结果 ---
回归截距 (Beta0): 2212.0553
RD_Exp 的标准化回归系数: 399.3895
Sales_Exp 的标准化回归系数: 476.7641
Admin_Exp 的标准化回归系数: 98.4139
模型 R-squared: 0.9460

上述结果来自事先设定系数的受控模拟，因此只用于验证标准化多元回归能否回收数据生成过程。系数大小和 \(R^2\) 由模拟参数决定，不能解释为长三角真实企业中销售、研发或管理费用的贡献排序，更不能据此提出预算干预建议。

下面我们将拟合后的多元线性回归关系以3D图形的方式进行可视化展示。将管理费用固定在其标准化均值（即0）处，通过研发费用与销售费用这两个维度构建回归超平面，散点代表各企业的实际观测值，颜色深浅映射营收高低。

# ---- 3D可视化（取研发与销售两个维度） ----
fig_3d_reg = plt.figure(figsize=(14, 10))  # 创建大尺寸画布
ax_3d_reg = fig_3d_reg.add_subplot(111, projection='3d')  # 添加3D子图

# 绘制三维散点（颜色按营收深浅编码，viridis色阶映射收入高低）
scatter_plot_3d = ax_3d_reg.scatter(X_features_scaled[:, 0], X_features_scaled[:, 1], y_target_obs,
                                  c=y_target_obs, cmap='viridis', s=60, alpha=0.7)  # viridis渐变色映射营收

x_grid_range = np.linspace(X_features_scaled[:, 0].min(), X_features_scaled[:, 0].max(), 20)  # 研发维度网格
y_grid_range = np.linspace(X_features_scaled[:, 1].min(), X_features_scaled[:, 1].max(), 20)  # 销售维度网格
X_surf, Y_surf = np.meshgrid(x_grid_range, y_grid_range)  # 构造二维网格矩阵

# 计算回归超平面上的预测值（管理费用取标准化均值0）
Z_surf_pred = (multi_ols_reg.intercept_ +  # 截距项
              multi_ols_reg.coef_[0] * X_surf +  # 研发费用的边际贡献
              multi_ols_reg.coef_[1] * Y_surf +  # 销售费用的边际贡献
              multi_ols_reg.coef_[2] * 0)  # 管理费用固定为标准化均值水平

ax_3d_reg.plot_surface(X_surf, Y_surf, Z_surf_pred, alpha=0.2, color='red')  # 绘制半透明红色回归平面

ax_3d_reg.set_xlabel('研发费用 (标准化)', fontsize=11)  # X轴标签
ax_3d_reg.set_ylabel('销售费用 (标准化)', fontsize=11)  # Y轴标签
ax_3d_reg.set_zlabel('营业收入 (万元)', fontsize=11)  # Z轴标签
ax_3d_reg.set_title('多元回归：研发与销售规模对营收的联合贡献', fontsize=14, pad=20)  # 图标题

plt.show()  # 显示3D回归可视化
三维散点以研发费用和销售费用为平面轴、营业收入为高度，半透明回归平面表示联合线性拟合。
图 3.2: 多元线性回归的三维可视化：研发费用与销售费用对营业收入的联合预测超平面。散点颜色深浅反映营收水平高低（已标准化处理）。

图 3.2 将多元线性回归的拟合结果以三维视角直观呈现。图中每个散点代表一家长三角制造业上市公司，其位置由标准化后的研发费用（X轴）和销售费用（Y轴）决定，Z轴为实际营业收入。半透明的红色平面即为 OLS 估计得到的回归超平面。可以观察到：绝大多数散点围绕回归平面两侧分布，没有出现系统性的偏离模式，说明线性假设在此场景下是合理的。同时，沿销售费用维度的平面上升梯度略大于研发费用维度，与前面标准化系数的比较结论相呼应。颜色较深（营收较高）的企业倾向于聚集在两个费用标准化值都较高的区域，进一步验证了费用投入规模与营收之间的正向关联。

3.4.1 估计回归系数：矩阵视角 (Estimating via Matrix Algebra)

在多元回归中，模型可以简洁地表示为：

\[ \mathbf{Y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\epsilon} \tag{3.19}\]

为了理解这一数学模型，我们将各个组件拆解如下：

  • \(\mathbf{Y}\) 为 \(n \times 1\) 的向量，包含所有观测响应。
  • \(\mathbf{X}\) 为 \(n \times (p+1)\) 的设计矩阵(Design Matrix)。其第一列全是 1（用于截距项），随后的每列代表一个预测变量。
  • \(\boldsymbol{\beta}\) 为需估计的 \((p+1) \times 1\) 系数向量。
  • \(\boldsymbol{\epsilon}\) 为 \(n \times 1\) 的随机误差向量。

3.4.1.1 核心推导：正规方程 (Normal Equations)

我们寻找向量 \(\hat{\boldsymbol{\beta}}\)，使得残差平方和(RSS)最小：

\[ \text{RSS}(\boldsymbol{\beta}) = (\mathbf{Y} - \mathbf{X}\boldsymbol{\beta})^T (\mathbf{Y} - \mathbf{X}\boldsymbol{\beta}) \]

展开该标量形式： \[ \text{RSS}(\boldsymbol{\beta}) = \mathbf{Y}^T \mathbf{Y} - 2\boldsymbol{\beta}^T \mathbf{X}^T \mathbf{Y} + \boldsymbol{\beta}^T \mathbf{X}^T \mathbf{X} \boldsymbol{\beta} \]

利用矩阵微积分对 \(\boldsymbol{\beta}\) 求梯度： \[ \nabla_{\boldsymbol{\beta}} \text{RSS} = -2 \mathbf{X}^T \mathbf{Y} + 2 \mathbf{X}^T \mathbf{X} \boldsymbol{\beta} \]

令梯度为零，得到著名的正规方程： \[ \mathbf{X}^T \mathbf{X} \hat{\boldsymbol{\beta}} = \mathbf{X}^T \mathbf{Y} \tag{3.20}\]

由此解得 OLS 的矩阵解析解： \[ \hat{\boldsymbol{\beta}} = (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{X}^T \mathbf{Y} \tag{3.21}\]

Warning: 矩阵可逆性与多重共线性

注意，公式 式 3.21 的前提是 \(\mathbf{X}^T \mathbf{X}\) 必须可逆（即非奇异）。如果预测变量之间存在完全共线性(Perfect Collinearity)，即某一列是其他列的线性组合，则 \(\mathbf{X}^T \mathbf{X}\) 不存在逆矩阵，模型无法唯一求解。在金融数据中，过多的冗余财务指标常会导致“近似多重共线性”，造成系数估计不稳。

3.4.2 一些重要问题 (Some Important Questions)

在拟合多元线性回归模型后，我们通常关注以下问题：

3.4.2.1 1. 至少有一个预测变量与响应变量相关吗？ (Is There a Relationship Between the Response and Predictors?)

这个问题的答案可以通过F检验 (F-test) 来回答。

假设：

  • \(H_0: \beta_1 = \beta_2 = \cdots = \beta_p = 0\)（所有预测变量都与\(Y\)无关）
  • \(H_a:\) 至少有一个\(\beta_j \neq 0\)

F统计量：

\[ F = \frac{(\text{TSS} - \text{RSS}) / p}{\text{RSS} / (n - p - 1)} \tag{3.22}\]

在\(H_0\)为真时，\(F\)服从\(F_{p, n-p-1}\)分布。如果F统计量的值很大（对应的p值很小），我们拒绝\(H_0\)。

下面用受控模拟演示联合 \(F\) 检验。由于非零广告系数已写入数据生成过程，小 \(p\) 值只说明检验在该设定下识别到至少一个非零条件关联；它不是现实广告投放效果的证据。

表 3.2: 多元线性回归的ANOVA表和F检验 
import statsmodels.api as sm  # 导入statsmodels回归引擎
import numpy as np  # 导入numpy用于数值运算
import pandas as pd  # 导入pandas用于DataFrame操作
from scipy import stats  # 导入scipy的统计模块（备用）

np.random.seed(42)  # 设置可复现的随机种子
sample_size = 200  # 模拟200个分公司
tv_ad_budget = np.random.uniform(0, 300, sample_size)  # 电视广告预算（万元）
online_ad_budget = np.random.uniform(0, 100, sample_size)  # 网络广告预算（万元）
newspaper_ad_budget = np.random.uniform(0, 100, sample_size)  # 报纸广告预算（万元）
# DGP: 销量 = 5 + 0.04*电视 + 0.03*网络 + 0.01*报纸 + ε
# 根据线性DGP加高斯噪声生成模拟销售量
sales_revenue = 5 + 0.04*tv_ad_budget + 0.03*online_ad_budget + 0.01*newspaper_ad_budget + np.random.normal(0, 2, sample_size)

数据模拟完成后，我们将三种广告渠道作为特征矩阵输入OLS模型，提取F统计量进行整体联合显著性检验。

# 构造含截距的特征矩阵（三种广告投入）
media_features_matrix = sm.add_constant(pd.DataFrame({  # 合并三列并自动添加截距
    '电视': tv_ad_budget,  # 电视广告预算列
    '网络': online_ad_budget,  # 网络广告预算列
    '报纸': newspaper_ad_budget  # 报纸广告预算列
}))  # 完成含截距的特征矩阵构建，用于后续F检验
media_sales_model = sm.OLS(sales_revenue, media_features_matrix).fit()  # OLS拟合多元回归模型

f_stat = media_sales_model.fvalue  # 提取模型整体的F统计量
f_pvalue = media_sales_model.f_pvalue  # 提取F检验对应的联合p值

print('F检验结果：')  # 输出标题
print('=' * 50)  # 分隔线
print(f'F统计量 = {f_stat:.4f}')  # 输出F统计量数值
print(f'p值 = {f_pvalue:.4e}')  # 输出p值（科学计数法格式）
print(f'\n结论：p值 < 0.05，拒绝零假设。')  # 给出统计学结论
print('在受控模拟中，至少一个预设斜率被检验识别为非零。')  # 限定结论范围
print(f'\nR² = {media_sales_model.rsquared:.4f}')  # 输出总体拟合优度
print(f'调整R² = {media_sales_model.rsquared_adj:.4f}')  # 输出考虑自由度惩罚后的调整R²
F检验结果：
==================================================
F统计量 = 203.1853
p值 = 6.8172e-60

结论：p值 < 0.05，拒绝零假设。
在受控模拟中，至少一个预设斜率被检验识别为非零。

R² = 0.7567
调整R² = 0.7530

运行结果应以代码现场输出为准。若联合 \(p\) 值很小，可拒绝“受控模拟中的所有斜率同时为零”，但不能把这一机械检验解释为现实广告渠道的影响、重要性或因果贡献。\(R^2\) 同样只是该模拟样本的拟合摘要。

3.4.2.2 2. 所有预测变量都重要吗？哪些变量重要？ (Do All Predictors Help? Which Variables Matter?)

要判断单个预测变量是否重要，我们可以对每个系数进行t检验：

\[ H_0: \beta_j = 0 \quad \text{vs} \quad H_a: \beta_j \neq 0 \]

t统计量为：

\[ t = \frac{\hat{\beta}_j}{\text{SE}(\hat{\beta}_j)} \tag{3.23}\]

Clarification on 多重共线性

当预测变量之间高度相关时，会出现多重共线性(multicollinearity)问题。这会导致：

  1. 系数估计不稳定：小样本变化可能导致系数估计大幅变化。
  2. 标准误增大：系数的标准误会变大，导致t统计量变小，即使变量实际重要，也可能无法拒绝零假设。
  3. 解释困难：系数的符号可能与直觉相反，或难以解释。

检测多重共线性的方法：

  • 计算预测变量之间的相关系数矩阵
  • 计算方差膨胀因子 (VIF)：\(\text{VIF} = \frac{1}{1 - R_j^2}\)，其中\(R_j^2\)是用其他预测变量预测\(X_j\)的\(R^2\)
  • 5 或 10 有时被用作经验警戒线，但不是通用检验阈值；VIF 应结合研究目的、变量定义和估计不确定性解释

处理多重共线性的方法：

  • 若目标是解释，先明确要保持的估计对象，再考虑重参数化并报告扩大的不确定性
  • 若目标是预测，可在训练内部比较正则化、降维或精简特征的样本外表现
  • 不应只因越过某个经验阈值就机械删除有理论含义的变量

3.4.2.3 3. 模型拟合数据的效果如何？ (How Well Does the Model Fit the Data?)

在多元回归中，我们使用\(R^2\)和调整\(R^2\)来评估模型拟合：

\[ R^2 = 1 - \frac{\text{RSS}}{\text{TSS}} \tag{3.24}\] \[ \text{Adjusted } R^2 = 1 - \frac{\text{RSS} / (n - p - 1)}{\text{TSS} / (n - 1)} \tag{3.25}\]

调整\(R^2\)惩罚了添加无用变量的行为，因此在比较不同模型时更有用。

3.4.2.4 4. 给定一组预测变量的值，如何预测响应变量？ (Given a Set of Predictor Values, What Response Value Should We Predict?)

对于新的观测\(\mathbf{x}_0 = (1, x_{01}, x_{02}, \ldots, x_{0p})^T\)，预测值为：

\[ \hat{y}_0 = \mathbf{x}_0^T \hat{\boldsymbol{\beta}} \tag{3.26}\]

预测的置信区间和预测区间（prediction interval）为：

  • 95%置信区间（针对平均响应）： \[ \hat{y}_0 \pm 2 \cdot \text{SE}(\hat{y}_0) \]

  • 95%预测区间（针对单个观测）： \[ \hat{y}_0 \pm 2 \cdot \sqrt{\text{Var}(\epsilon) + \text{Var}(\hat{y}_0)} \tag{3.27}\]

3.5 回归模型的其他考虑 (Other Considerations in the Regression Model)

3.5.1 定性预测变量 (Qualitative Predictors)

在实际的商业与经济研究中，预测变量并不总是定量的。在中国资本市场的研究中，最典型的定性变量包括：

  • 企业性质(Ownership)：国企 vs 民企。
  • 地区分布(Region)：长三角、大湾区 vs 其他地区。
  • 行业类别(Industry)：科技、金融、制造等。

3.5.1.1 虚拟变量 (Dummy Variables)

对于只有两个水平的定性变量，我们引入一个虚拟变量(dummy variable) \(D\)：

\[ D_i = \begin{cases} 1 & \text{若第 } i \text{ 个公司是国有企业} \\ 0 & \text{若第 } i \text{ 个公司是非国有企业} \end{cases} \]

将其放入模型： \[ y_i = \beta_0 + \beta_1 D_i + \epsilon_i \]

此时系数的经济含义非常直观：

  • \(\beta_0\) 是非国有企业的平均响应值。
  • \(\beta_0 + \beta_1\) 是国有企业的平均响应值。
  • \(\beta_1\) 代表了两者之间的平均差异指标。

受控模拟：虚拟变量系数的解释

下面用带有明确数据生成过程的模拟样本说明虚拟变量系数；本节不估计真实国企与民企的盈利差异。

代码把产权标签编码为 0/1，并在模拟中控制对数资产规模。拟合系数说明如何读取“给定模型和其他变量时的组间条件均值差”，而不是提供关于现实产权效率的证据。

表 3.3: 受控模拟中的虚拟变量回归：产权标签与ROE 
import pandas as pd  # 导入pandas用于构建数据框
import numpy as np  # 导入numpy用于数值计算
import statsmodels.api as sm  # 导入statsmodels进行OLS回归

np.random.seed(42)  # 固定随机种子
n_firms = 300  # 模拟样本量为300家企业

# 模拟企业资产规模（对数值），范围约为20~26，对应现实中的中型上市公司
log_asset_size = np.random.uniform(20, 26, n_firms)  # 生成均匀分布随机样本

# 模拟产权性质：0代表民营企业、1代表国有企业，比例7:3
is_state_owned = np.random.choice([0, 1], size=n_firms, p=[0.7, 0.3])  # 随机抽样

# 真实DGP: ROE = 0.02 + 0.005*(Size-23) - 0.015*SOE + ε
# 负号意味着在该模拟设定下，国企的平均ROE低于民企
# 根据DGP线性结构加高斯噪声生成模拟ROE值
firm_roe_values = 0.02 + 0.005 * (log_asset_size - 23) - 0.015 * is_state_owned + np.random.normal(0, 0.02, n_firms)

模拟数据生成完毕后，我们利用 statsmodels 的OLS函数来拟合虚拟变量回归模型，并解读国有产权标签对ROE的系数含义。

# 构建数据框
df_firm_data = pd.DataFrame({  # 整合三列变量到同一个数据框
    'LogSize': log_asset_size,  # 对数资产规模
    'IsSOE': is_state_owned,  # 国有企业虚拟变量
    'ROE': firm_roe_values  # 净资产收益率
})  # 完成数据框构建，用于后续虚拟变量回归建模

# 构造含截距的自变量矩阵
X_dummy_matrix = sm.add_constant(df_firm_data[['LogSize', 'IsSOE']])  # 自动添加常数列作为截距项

# OLS拟合
ownership_ols_model = sm.OLS(df_firm_data['ROE'], X_dummy_matrix).fit()  # 最小二乘法估计回归参数

print(ownership_ols_model.summary().tables[1])  # 输出系数估计表（含标准误、t值、p值及置信区间）

print('\n--- 结果解读 ---')  # 输出结果到控制台
soe_coefficient = ownership_ols_model.params['IsSOE']  # 提取国企虚拟变量的回归系数
print(f'国有企业(IsSOE)的系数为: {soe_coefficient:.4f}')  # 输出系数值
# 解读：控制资产规模后，国企ROE平均低于民企的幅度
# 格式化输出国企相对民企的ROE差距（百分点）
print('模拟样本中，在给定规模后，产权标签对应的条件均值差为 {:.2f} 个百分点。'.format(soe_coefficient*100))
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const         -0.0796      0.015     -5.178      0.000      -0.110      -0.049
LogSize        0.0043      0.001      6.498      0.000       0.003       0.006
IsSOE         -0.0179      0.003     -7.126      0.000      -0.023      -0.013
==============================================================================

--- 结果解读 ---
国有企业(IsSOE)的系数为: -0.0179
模拟样本中，在给定规模后，产权标签对应的条件均值差为 -1.79 个百分点。

若拟合结果接近代码设定的系数，只能说明估计程序在这一受控模拟中工作正常。显著性、符号和数值都由数据生成过程预设，不能被称为实证发现，也不能支持规模经济、产权效率或治理机制的现实解释。真实研究还需要可追溯数据、明确的抽样框架与因果识别设计。

3.5.2 线性模型的扩展 (Extensions of the Linear Model)

标准的线性模型基于两个极强的假设：可加性(Additivity)和线性性(Linearity)。

3.5.2.1 交互效应与协同作用 (Interaction Effects)

可加性假设认为一个预测变量对响应的影响与各其他变量的水平无关。但在商业决策中，这种假设往往过于简单。

例如，一个公司的研发投入 (\(X_{RD}\)) 和行业集中度 (\(X_{HHI}\)) 可能会产生相互影响：低竞争行业的研发投入产出比可能远高于高竞争行业。这种协同作用可以通过交互项(Interaction Term)来建模：

\[ Y = \beta_0 + \beta_1 X_{RD} + \beta_2 X_{HHI} + \beta_3 (X_{RD} \times X_{HHI}) + \epsilon \tag{3.28}\]

3.5.2.2 非线性关系：多项式回归 (Non-Linearity)

线性假设认为 \(X\) 每变化一单位，\(Y\) 的变化是恒定的。然而，许多金融关系呈现 U 型、倒 U 型或随 \(X\) 改变的边际斜率。

案例：公司规模与创新效率

著名的“熊彼特假设”认为大企业更具创新优势，但当企业规模过大时，组织官僚化可能导致创新效率下降。这可以用二次项来捕捉：

\[ \text{Innovation} = \beta_0 + \beta_1 \text{Size} + \beta_2 \text{Size}^2 + \epsilon \tag{3.29}\]

若 \(\beta_2<0\)，拟合二次函数是凹的；其候选转折点为 \(x^*=-\beta_1/(2\beta_2)\)。只有当 \(x^*\) 位于事前定义的可行域与观测支持域内，且转折点及两侧斜率的不确定性允许这种解释时，才可把样本曲线称为支持域内的倒 U 型；否则只能报告凹性，不能宣称存在可行“最优规模”。观察回归也不把这个转折点识别为因果最优政策。

下面用受控模拟检验多项式管道是否能回收一个已知的凹二次数据生成过程。真实条件均值设为 \(5+0.09x-0.0002x^2\)，所以理论转折点为 225 万元；固定种子实现的广告预算支持域约为 \([1.657,296.066]\) 万元，转折点两侧分别有 149 与 51 个观测。因而本例可以讨论支持域内先增后减的形状，但它只验证方法，不是现实广告因果证据。

import numpy as np  # 导入numpy用于数组运算和随机数生成
import matplotlib.pyplot as plt  # 导入matplotlib用于绑图
import pandas as pd  # 整理训练与验证证据表
from sklearn.preprocessing import PolynomialFeatures  # 导入多项式特征生成器
from sklearn.preprocessing import StandardScaler  # 在训练期缩放原始预算以改善幂基数值条件
from sklearn.linear_model import LinearRegression  # 导入线性回归模型
from sklearn.pipeline import Pipeline  # 导入管道工具用于串联预处理与建模
from sklearn.metrics import mean_squared_error  # 分别计算训练与验证均方误差
from sklearn.model_selection import train_test_split  # 用固定切分构造未参与拟合的验证证据

np.random.seed(42)  # 设置随机种子确保可复现
observation_count = 200  # 设定模拟样本量为200个观测
tv_ad_budget = np.random.uniform(0, 300, observation_count)  # 模拟电视广告预算，范围0~300万
turning_point_budget = 225.0  # 由已知二次DGP的一阶条件得到支持域内转折点
true_sales_revenue = 5 + 0.09 * tv_ad_budget - 0.0002 * tv_ad_budget**2  # 构造先增后减的凹二次条件均值
sales_revenue_obs = true_sales_revenue + np.random.normal(0, 2, observation_count)  # 添加高斯噪声模拟观测误差
budget_train, budget_validation, sales_train, sales_validation = train_test_split(tv_ad_budget.reshape(-1, 1), sales_revenue_obs, test_size=0.3, random_state=42)  # 固定留出30%验证样本
assert tv_ad_budget.min() < turning_point_budget < tv_ad_budget.max()  # 拒绝把支持域外顶点解释为样本内转折
assert np.sum(tv_ad_budget < turning_point_budget) > 0 and np.sum(tv_ad_budget > turning_point_budget) > 0  # 要求转折点两侧都有证据

数据准备完成后，表 3.4 在同一固定训练—验证切分上比较三种复杂度。缩放器只在训练样本中拟合；十次模型是否过拟合由验证损失而不是曲线是否难看来判断。

degrees = [1, 2, 10]  # 分别尝试线性、二次和十次多项式
colors = ['steelblue', 'darkgreen', 'crimson']  # 各阶次对应的曲线颜色
titles = ['线性模型 (Linear)', '二次多项式 (Quadratic)', '10次多项式 (Degree 10)']  # 子图标题
polynomial_models = {}  # 保存只在训练样本拟合的候选模型
polynomial_metric_records = []  # 收集训练损失、验证损失与条件数

for degree in degrees:  # 在同一切分上逐一拟合预先声明的阶数
    polynomial_model = Pipeline([  # 创建sklearn Pipeline串联预处理与建模
        ('scale', StandardScaler()),  # 只用训练预算估计中心与尺度
        ('poly', PolynomialFeatures(degree=degree, include_bias=False)),  # 将标准化预算扩展为幂基
        ('linear', LinearRegression())  # 对扩展后的特征矩阵做OLS拟合
    ])  # 完成Pipeline管道组装
    polynomial_model.fit(budget_train, sales_train)  # 禁止验证样本参与缩放与拟合
    polynomial_models[degree] = polynomial_model  # 保存冻结候选供同一支持域绘图
    transformed_train = polynomial_model[:-1].transform(budget_train)  # 取得训练设计矩阵用于条件数诊断
    design_condition_number = np.linalg.cond(np.column_stack([np.ones(len(transformed_train)), transformed_train]))  # 量化幂基的数值敏感性
    polynomial_metric_records.append({'degree': degree, 'train_mse': mean_squared_error(sales_train, polynomial_model.predict(budget_train)), 'validation_mse': mean_squared_error(sales_validation, polynomial_model.predict(budget_validation)), 'condition_number': design_condition_number})  # 登记可复算证据

polynomial_validation_evidence = pd.DataFrame(polynomial_metric_records).round(4)  # 形成实际输出表而不预写胜者
polynomial_validation_evidence  # 显示固定种子证据
表 3.4: 固定种子多项式回归的训练、验证与数值条件证据
degree train_mse validation_mse condition_number
0 1 4.9719 5.4079 1.0000
1 2 3.5801 4.0633 2.8806
2 10 3.3807 4.5165 4984.3221

图 3.3 把拟合限制在实现样本的支持域，并标出已知转折点；虚线是真实模拟条件均值，不是模型可访问的训练标签。

x_grid_range = np.linspace(tv_ad_budget.min(), tv_ad_budget.max(), 300).reshape(-1, 1)  # 只在实现样本支持域内绘制
true_sales_grid = 5 + 0.09 * x_grid_range[:, 0] - 0.0002 * x_grid_range[:, 0]**2  # 计算已知DGP供视觉校验
fig, axes = plt.subplots(1, 3, figsize=(18, 5))  # 创建1行3列的画布

for axis, degree, color, title in zip(axes, degrees, colors, titles):  # 为每个冻结候选绘制独立子图
    y_predicted_values = polynomial_models[degree].predict(x_grid_range)  # 在共同支持域网格预测
    axis.scatter(budget_train[:, 0], sales_train, alpha=0.45, s=45, color='gray', label='训练样本')  # 展示用于拟合的观测
    axis.scatter(budget_validation[:, 0], sales_validation, alpha=0.65, s=45, color='#F0A700', label='验证样本')  # 展示未参与拟合的观测
    axis.plot(x_grid_range[:, 0], true_sales_grid, color='#2C3E50', linestyle='--', linewidth=1.8, label='真实条件均值')  # 绘制已知DGP
    axis.plot(x_grid_range[:, 0], y_predicted_values, color=color, linewidth=2.5, label=title)  # 绘制实际拟合曲线
    axis.axvline(turning_point_budget, color='#008080', linestyle=':', linewidth=2, label='DGP转折点')  # 标记支持域内顶点
    axis.set_xlabel('电视广告预算（万元）', fontsize=11)  # 设置预算轴标签
    axis.set_ylabel('销售量（千件）', fontsize=11)  # 设置响应轴标签
    axis.set_title(title, fontsize=12, fontweight='bold')  # 标明候选复杂度
    axis.legend(fontsize=9)  # 解释样本、真实曲线与拟合曲线
    axis.grid(True, alpha=0.3)  # 添加半透明网格

plt.tight_layout()  # 自动调整子图间距
plt.show()  # 显示三图对比
三幅并排图分别展示训练点、验证点、真实凹二次曲线和候选拟合；竖虚线标出支持域内225万元转折点。
图 3.3: 固定种子凹二次数据生成过程上的线性、二次与十次多项式拟合

当前固定种子下，表 3.4 给出的二次模型训练/验证 MSE 约为 3.5801/4.0633，十次模型约为 3.3807/4.5165：十次模型降低训练损失却提高验证损失，因此这里的“过拟合”有样本外证据。十次幂基的条件数也远高于二次模型，提示数值敏感性；但数值条件差与验证过拟合是两种诊断，前者不能单独证明后者。图 3.3 中二次拟合在 225 万元附近转为下降，与 DGP 和两侧样本覆盖一致；支持域之外不作形状外推。

3.5.3 潜在问题与诊断 (Potential Problems and Diagnostics)

在将线性回归应用于金融场景（如多因子选股或财务预测）时，经常会遇到破坏 OLS 假设的情形：

  1. 非线性关系 (Non-linearity)：响应变量与预测变量的真实关系并非线性。
    • 诊断：绘制残差图（Residual Plot）。若残差呈现明显模式（如抛物线），则需考虑多项式项。
  2. 误差项的相关性 (Correlation of Errors)：观测值之间不独立（常见于时间序列数据）。
    • 诊断：计算 Durbin-Watson 统计量。相关性会导致标准误被低估，t 统计量虚高。
  3. 异方差性 (Heteroscedasticity)：模型误差的方差不是常数。
    • 诊断：漏斗形的残差图。
    • 解决：使用加权最小二乘（WLS）或稳健标准误（Robust Standard Errors）。
  4. 异常值与高杠杆点 (Outliers and High Leverage Points)：
    • 诊断：计算学生化残差（Studentized Residuals）和 Cook 距离。
  5. 多重共线性 (Collinearity)：预测变量之间高度相关。

    • 诊断：计算方差膨胀因子(VIF)。 \[ \text{VIF}(\hat{\beta}_j) = \frac{1}{1 - R_j^2} \] 其中 \(R_j^2\) 是将 \(X_j\) 对其他所有预测变量进行回归得到的拟合优度。

    • 判定：较大的 VIF 表明该系数在当前设计矩阵中难以精确分离，但 5 或 10 只是经验警戒线。是否重参数化、正则化或保留变量，应由估计对象与预测目标决定，而非机械删列。

3.6 本章小结 (Chapter Summary)

本章详细介绍了线性回归，这是统计学习中最基础也最重要的方法之一：

  1. 简单线性回归：使用一个预测变量来预测响应变量，通过最小二乘法估计截距和斜率。

  2. 多元线性回归：使用多个预测变量，可以分离每个变量的独立效应。

  3. 模型评估：

    • 使用RSE衡量预测精度
    • 使用\(R^2\)衡量拟合优度
    • 使用F检验判断整体显著性
    • 使用t检验判断单个变量的显著性

  4. 模型扩展：

    • 虚拟变量处理定性预测变量
    • 交互项建模变量间的交互效应
    • 多项式项和变换建模非线性关系

  5. 潜在问题：非线性、多重共线性、异方差、异常值等，需要诊断和适当处理。

线性回归是许多更复杂方法的基础，掌握它对于理解后续章节中的高级方法至关重要。关于经典线性模型、统计学习视角及其与更灵活方法的联系，可参见 James 等 (2023年) 与 Hastie 等 (2009年)。

3.7 理论来源与前沿

线性回归的理论基础来自最小二乘与正态线性模型：在误差独立同分布且同方差时，OLS 具有最佳线性无偏性（Gauss-Markov 定理）；在正态误差假设下，OLS 与极大似然估计一致，并能导出精确的 t/F 检验与区间估计。现代统计学习强调在线性回归中系统处理三类工程问题：诊断（残差与杠杆点）、稳健性（厚尾/异方差）与可解释性（系数含义与变量选择）。

研究前沿方面，线性回归并未‘过时’，而是以新的形式持续活跃：

  1. 高维与稀疏建模：Lasso/Elastic Net 使 \(p\gg n\) 仍可估计，并把变量选择纳入优化。
  2. 稳健回归与分位数回归：在异常值与异方差普遍存在的金融数据中更可靠。
  3. 因果推断中的回归：在差分、匹配、工具变量等识别设计下，回归承担‘控制混杂’与‘估计处理效应’的角色。

3.8 贯穿项目里程碑 M03

本里程碑与 小节 6.3 的 M03 一致，预计 60 分钟，输入是已验收的 M02 a-share-drawdown-20d-v1 契约与同一逻辑快照。线性回归在这里作为二元标签的线性概率基线，不把系数解释为因果效应。

契约项 可审计要求
输入 M02 的同一键、20 日标签、label_date、特征可得日和训练/验证成员；不得创建新标签或访问最终测试
输出 训练事件率/多数类基线、线性概率模型的验证 Brier、系数表及条件关联解释；概率超出 \([0,1]\) 的数量须报告
时点与冻结 特征必须在 prediction_date 收盘时可得；训练拟合、验证比较，最终测试仍封存
失败条件 改变任务 ID、键、窗口或阈值；用目标构成期未来价格作特征；全样本拟合预处理；把系数称为因果；用验证/测试标签定义基线
量规 M 方法与验证 40%，C 结论边界 35%，E 表达与图表 25%

3.9 习题

3.9.1 概念题

  1. [核心｜难度:1｜时间:10分钟｜分值:10｜项目:无] BLUE 与 OLS：为什么说 OLS 在高斯-马尔可夫假设下是 BLUE？如果误差项不符合正态分布，OLS 是否还是无偏的？

  2. [核心｜难度:2｜时间:15分钟｜分值:15｜项目:无] 系数解释：在多元回归 \(Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \epsilon\) 中，如何解释 \(\beta_1\)？如果 \(X_1\) 和 \(X_2\) 高度正相关，\(\beta_1\) 的估计值可能会出现什么现象？

  3. [核心｜难度:1｜时间:10分钟｜分值:10｜项目:无] 模型评估：解释 \(R^2\) 与调整 \(R^2\) 的区别。为什么在解释性研究中，我们不能仅仅依靠 \(R^2\) 来选择模型？

  4. [核心｜难度:2｜时间:15分钟｜分值:15｜项目:无] 共线性诊断：什么是方差膨胀因子（VIF）？如果你在分析公司财务数据时发现“资产规模”和“营业成本”的 VIF 极高，你会如何处理？

3.9.2 应用题

  1. [核心｜难度:3｜时间:60分钟｜分值:40｜项目:M03] 同任务线性概率基线：复用 M02 的 a-share-drawdown-20d-v1 键、标签和训练/验证成员，使用当日及以前可得特征拟合线性概率模型。提交训练事件率与训练多数类基线、验证 Brier、超出 \([0,1]\) 的预测数、系数表和明确的非因果解释；不得访问最终测试。

  2. [拓展｜难度:3｜时间:35分钟｜分值:25｜项目:无] 交互项分析：研究“公司规模”与“研发投入”对“营收增长”是否存在交互效应。即：大企业的研发投入是否比较小企业更有利于营收增长？

3.9.3 理论题

  1. [拓展｜难度:2｜时间:15分钟｜分值:15｜项目:无] 参数推导：证明在简单线性回归中，拟合线一定经过样本中心点 \((\bar{x}, \bar{y})\)。

  2. [拓展｜难度:3｜时间:25分钟｜分值:20｜项目:无] 方差推导：给定 \(\hat{\beta} = (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{X}^T \mathbf{Y}\)，假设 \(\text{Var}(\boldsymbol{\epsilon} | \mathbf{X}) = \sigma^2 \mathbf{I}\)，证明 \(\text{Var}(\hat{\boldsymbol{\beta}} | \mathbf{X}) = \sigma^2 (\mathbf{X}^T \mathbf{X})^{-1}\)。

3.10 练习参考解答

展开第 3 章完整解答与评分键

统一评分与验收说明：按题面元数据分值计分；推导题中一阶条件、矩阵维度、交叉项消去和最终式均为独立评分点，数据题须交付系数/区间、VIF 与诊断。浮点结果相对容差 \(10^{-6}\)；系数表允许因软件版本导致最后一位显示差异，但样本、变量口径与稳健标准误必须一致。常见失败包括机械按 VIF 删列、把条件关联当因果、分母为零/无穷值未清理，以及交互项缺主效应。

3.10.1 概念题解答

  1. BLUE 指最佳线性无偏估计。无偏性只需 \(E[\epsilon|X]=0\) 即可满足，与正态性无关。正态分布假设主要用于小样本下的 \(t\) 检验和区间估计。

  2. \(\beta_1\) 表示在模型设定成立时，保持 \(X_2\) 不变，\(X_1\) 每增加 1 单位对应的 \(Y\) 条件均值变化。它不是自动成立的因果效应。若高度共线性，\(\beta_1\) 的标准误可能增大，且估计值会对样本扰动更敏感。

  3. 调整 \(R^2\) 对自由度进行了惩罚。\(R^2\) 总是随变量增加而不减，可能奖励无效复杂度。解释性研究还必须检查研究设计、模型假设、效应量及其不确定性；不能只按显著性或 \(R^2\) 选模。

  4. VIF 衡量在其他条件相同时，设计矩阵共线性对应的系数方差膨胀。先核对“资产规模”和“营业成本”的定义与研究对象；若做解释，应保留必要控制并报告不确定性，或用有经济含义的重参数化；若做预测，可在训练内部比较岭回归、降维或精简特征。不得仅因跨过经验阈值就机械删变量。

3.10.2 应用题参考解答（代码模板）

  1. 同任务线性概率基线：

答案可在 fresh kernel 顺序运行，不依赖正文状态。必须显式提供 BOOK_PROJECT_ARTIFACT_DIR、BOOK_PROJECT_DEVELOPMENT_FEATURES 与 BOOK_PROJECT_DEVELOPMENT_FEATURES_SHA256；project-development-v1.csv 的模式是键、label_date、split、y，特征文件的模式只能是相同键加数值特征。最终 test 键和标签均不载入。

import hashlib  # 核验全部上游与外部特征字节
import json  # 读取M02机器合同
import os  # 读取fresh-kernel所需显式入口
from pathlib import Path  # 使用跨平台路径
import numpy as np  # 导入数值与有限性工具
import pandas as pd  # 读取冻结制品
from sklearn.linear_model import LinearRegression  # 导入线性概率模型
from sklearn.metrics import brier_score_loss  # 导入概率平方损失

artifact_dir = Path(os.environ['BOOK_PROJECT_ARTIFACT_DIR']).resolve()  # 定位M02唯一制品桶
manifest_path = artifact_dir / 'project-manifest-v1.csv'  # 定位权威全成员表
development_path = artifact_dir / 'project-development-v1.csv'  # 定位M02开发标签
contract_path = artifact_dir / 'project-contract-v1.json'  # 定位总体与哈希合同
feature_path = Path(os.environ['BOOK_PROJECT_DEVELOPMENT_FEATURES']).resolve()  # 定位仅开发特征
assert all(path.is_file() for path in [manifest_path, development_path, contract_path, feature_path])  # 输入缺失时关闭流程
contract = json.loads(contract_path.read_text(encoding='utf-8'))  # 读取M02登记哈希
expected_hashes = contract['artifact_hashes']  # 取得权威制品哈希字段
assert hashlib.sha256(manifest_path.read_bytes()).hexdigest() == expected_hashes['manifest_sha256']  # 拒绝manifest变化
assert hashlib.sha256(development_path.read_bytes()).hexdigest() == expected_hashes['development_sha256']  # 拒绝开发标签变化
assert hashlib.sha256(feature_path.read_bytes()).hexdigest() == os.environ['BOOK_PROJECT_DEVELOPMENT_FEATURES_SHA256'].lower()  # 绑定真实特征文件
manifest = pd.read_csv(manifest_path, parse_dates=['prediction_date', 'label_date'])  # 读取权威全成员表
project_labels = pd.read_csv(development_path, parse_dates=['prediction_date', 'label_date'])  # 读取M02开发标签
project_features = pd.read_csv(feature_path, parse_dates=['prediction_date'])  # 读取预测时点可得特征
key_columns = ['order_book_id', 'prediction_date']  # 冻结同一项目键
reserved_columns = {'label_date', 'y', 'split', 'event', 'future_min_close', 'target', 'outcome'}  # 禁止特征源携带结果或成员字段
assert reserved_columns.isdisjoint(set(project_features.columns))  # 字段冲突在合并前失败
assert not project_features.duplicated(key_columns).any()  # 特征键必须唯一
development_manifest = manifest[manifest['split'].isin(['train', 'validation'])].copy()  # 以manifest为权威开发左表
label_columns = key_columns + ['label_date', 'split', 'y']  # 冻结开发标签模式
assert set(project_labels.columns) == set(label_columns)  # 拒绝标签制品模式漂移
authoritative_panel = development_manifest.merge(project_labels[key_columns + ['y']], on=key_columns, how='left', validate='one_to_one', indicator='label_merge')  # 审计标签精确覆盖
assert authoritative_panel['label_merge'].eq('both').all() and len(authoritative_panel) == len(development_manifest)  # 禁止标签缺键缩样
assert set(map(tuple, project_labels[key_columns].to_numpy())) == set(map(tuple, development_manifest[key_columns].to_numpy()))  # 禁止额外标签键
authoritative_panel = authoritative_panel.drop(columns='label_merge')  # 移除标签审计字段
project_panel = authoritative_panel.merge(project_features, on=key_columns, how='left', validate='one_to_one', indicator='feature_merge')  # 权威左连接特征
assert project_panel['feature_merge'].eq('both').all() and len(project_panel) == len(development_manifest)  # 禁止特征缺键缩样
assert set(map(tuple, project_features[key_columns].to_numpy())) == set(map(tuple, development_manifest[key_columns].to_numpy()))  # 禁止额外特征键
project_panel = project_panel.drop(columns='feature_merge')  # 移除特征审计字段
feature_columns = [column for column in project_features.columns if column not in key_columns]  # 冻结无保留冲突的特征列
assert feature_columns and np.isfinite(project_panel[feature_columns].to_numpy(dtype=float)).all()  # 拒绝空特征与非有限输入
key_bytes = project_panel[key_columns].to_csv(index=False, date_format='%Y-%m-%d', lineterminator='\n').encode('utf-8')  # 规范化实际开发键
assert hashlib.sha256(key_bytes).hexdigest() == hashlib.sha256(development_manifest[key_columns].to_csv(index=False, date_format='%Y-%m-%d', lineterminator='\n').encode('utf-8')).hexdigest()  # 精确核对排序、行与键哈希
print({'development_rows': len(project_panel), 'development_key_sha256': hashlib.sha256(key_bytes).hexdigest(), 'features_sha256': hashlib.sha256(feature_path.read_bytes()).hexdigest()})  # 输出血缘证据
train_panel = project_panel[project_panel['split'] == 'train'].copy()  # 只取得训练成员
validation_panel = project_panel[project_panel['split'] == 'validation'].copy()  # 只取得验证成员
assert train_panel['label_date'].max() < validation_panel['prediction_date'].min()  # 核验20日purge边界
assert train_panel['y'].nunique() == 2 and validation_panel['y'].nunique() == 2  # 概率评价要求双类别
training_event_rate = train_panel['y'].mean()  # 只用训练标签计算事件率
training_majority = int(training_event_rate >= 0.5)  # 只用训练标签锁定多数类
linear_probability_model = LinearRegression().fit(train_panel[feature_columns], train_panel['y'])  # 拟合LPM
validation_probability_raw = linear_probability_model.predict(validation_panel[feature_columns])  # 生成验证预测
validation_probability = np.clip(validation_probability_raw, 0, 1)  # 为Brier报告截断概率并保留越界计数
linear_brier = brier_score_loss(validation_panel['y'], validation_probability)  # 计算验证Brier
baseline_probability = np.repeat(training_event_rate, len(validation_panel))  # 构造训练事件率概率基线
baseline_brier = brier_score_loss(validation_panel['y'], baseline_probability)  # 计算基线Brier
coefficient_table = pd.DataFrame({'term': ['intercept'] + feature_columns, 'coefficient': [linear_probability_model.intercept_] + linear_probability_model.coef_.tolist()})  # 生成具名系数表
assert coefficient_table['term'].is_unique and np.isfinite(coefficient_table['coefficient']).all()  # 拒绝匿名、重复或非有限系数
print({'training_event_rate': training_event_rate, 'training_majority': training_majority})  # 输出训练基线
print({'lpm_brier': linear_brier, 'baseline_brier': baseline_brier, 'outside_unit_interval': int(((validation_probability_raw < 0) | (validation_probability_raw > 1)).sum())})  # 输出验证证据
print(coefficient_table)  # 交付题面要求的具名有限系数表

系数只描述冻结特征与未来 20 日回撤事件概率的条件关联，不是交易、因果或市场机制效应。若 LPM 验证 Brier 未低于训练事件率基线，应报告“无增量概率预测证据”，不能改用测试挑选结论。

  1. 交互项分析：

本题要求检验”公司规模”与”研发投入”对”营收增长”是否存在交互效应。我们构建如下模型：

\[ \text{RevenueGrowth} = \beta_0 + \beta_1 \text{Size} + \beta_2 \text{RD} + \beta_3 (\text{Size} \times \text{RD}) + \epsilon \]

若 \(\beta_3\) 显著不为零，则说明研发费用率与下一期营收增长的条件关联随公司规模而变。下面使用本地真实年报面板，以 \(t\) 年已公布的研发费用率和规模预测 \(t+1\) 年营收增长，并使用 HC3 异方差稳健标准误。

import os  # 导入环境变量接口以显式解析数据入口
from pathlib import Path  # 导入跨平台路径对象
import numpy as np  # 本题独立完成有限值与对数变换
import pandas as pd  # 本题独立读取财务面板
import statsmodels.api as sm  # 本题独立拟合HC3稳健OLS

BOOK_DATA_DIR = Path(os.environ['BOOK_DATA_DIR']).expanduser().resolve()  # 要求调用者显式提供数据根目录
local_data_path = BOOK_DATA_DIR / 'stock/financial_statement.h5'  # 明确声明本题使用的财务面板文件
interaction_source = pd.read_hdf(  # 选择性读取年报面板所需字段
    local_data_path, where="quarter>='2018q4' & quarter<='2023q4'",  # 限定可比年报时段
    columns=['order_book_id', 'quarter', 'revenue', 'r_n_d', 'total_assets']  # 仅读取信息时点可得变量
).copy()  # 创建面板变量副本
interaction_source = interaction_source[interaction_source['quarter'].str.endswith('q4')].copy()  # 只保留年报
interaction_source['fiscal_year'] = interaction_source['quarter'].str[:4].astype(int)  # 提取排序用会计年度
interaction_source = interaction_source.sort_values(['order_book_id', 'fiscal_year'])  # 按公司与年度排序
interaction_source['future_revenue'] = interaction_source.groupby('order_book_id')['revenue'].shift(-1)  # 下一年营收作为未来目标
interaction_source['future_year'] = interaction_source.groupby('order_book_id')['fiscal_year'].shift(-1)  # 记录目标年度以核验连续性
interaction_source = interaction_source[interaction_source['future_year'] == interaction_source['fiscal_year'] + 1].copy()  # 拒绝跨年缺口错配
interaction_source['revenue_growth'] = interaction_source['future_revenue'] / interaction_source['revenue'] - 1  # 构造下一年营收增长
interaction_source['log_assets'] = np.log(interaction_source['total_assets'])  # 用对数总资产衡量规模
interaction_source['rd_ratio'] = interaction_source['r_n_d'] / interaction_source['revenue']  # 构造当期研发费用率
interaction_columns = ['revenue_growth', 'log_assets', 'rd_ratio']  # 定义完整案例变量契约
interaction_panel = interaction_source[interaction_columns].replace([np.inf, -np.inf], np.nan).dropna().copy()  # 清理不可计算观测
interaction_panel = interaction_panel[interaction_panel['revenue_growth'].between(-1, 3)]  # 剔除明显口径异常增长
interaction_panel = interaction_panel[interaction_panel['rd_ratio'].between(0, 0.5)]  # 限定具有经济含义的研发费用率
interaction_panel['centered_log_assets'] = interaction_panel['log_assets'] - interaction_panel['log_assets'].mean()  # 中心化便于解释主效应
interaction_panel['size_by_rd'] = interaction_panel['centered_log_assets'] * interaction_panel['rd_ratio']  # 构造规模与研发的交互项
interaction_features = ['centered_log_assets', 'rd_ratio', 'size_by_rd']  # 定义主效应和交互项
interaction_design_matrix = sm.add_constant(interaction_panel[interaction_features])  # 添加截距
interaction_ols_model = sm.OLS(interaction_panel['revenue_growth'], interaction_design_matrix).fit(cov_type='HC3')  # 拟合异方差稳健 OLS
interaction_coef = interaction_ols_model.params['size_by_rd']  # 提取交互项条件关联估计
interaction_pvalue = interaction_ols_model.pvalues['size_by_rd']  # 提取 HC3 稳健 p 值
print(interaction_ols_model.summary().tables[1])  # 输出可复算系数、稳健标准误与区间
print(f'规模×研发系数: {interaction_coef:.4f}; HC3 p值: {interaction_pvalue:.4g}')  # 报告不预设显著性的结果
=======================================================================================
                          coef    std err          z      P>|z|      [0.025      0.975]
---------------------------------------------------------------------------------------
const                   0.0970      0.003     29.417      0.000       0.091       0.103
centered_log_assets    -0.0073      0.002     -3.789      0.000      -0.011      -0.004
rd_ratio                0.3816      0.055      6.920      0.000       0.274       0.490
size_by_rd             -0.3579      0.038     -9.432      0.000      -0.432      -0.284
=======================================================================================
规模×研发系数: -0.3579; HC3 p值: 4.034e-21

交互项的边际含义是 \(\partial E(\text{RevenueGrowth}\mid X)/\partial\text{RD}=\hat\beta_2+\hat\beta_3\text{CenteredSize}\)。无论 \(p\) 值如何，这里都只是观察性条件关联；研发预算可能与未观测的管理质量、行业机会同时变动，因此不作因果解释。

3.10.3 理论题解答

  1. 证明中心点： 由一阶条件 \(\frac{\partial RSS}{\partial \beta_0} = -2\sum(y_i - \beta_0 - \beta_1 x_i) = 0\)。 展开得：\(\sum y_i - n\beta_0 - \beta_1 \sum x_i = 0\)。 两边除以 \(n\)：\(\bar{y} - \beta_0 - \beta_1 \bar{x} = 0 \Rightarrow \bar{y} = \beta_0 + \beta_1 \bar{x}\)。 证毕。

  2. 方差证明： \(\hat{\beta} = (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{X}^T (\mathbf{X}\boldsymbol{\beta} + \boldsymbol{\epsilon}) = \boldsymbol{\beta} + (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{X}^T \boldsymbol{\epsilon}\)。 \(\text{Var}(\hat{\beta} | \mathbf{X}) = \text{Var}((\mathbf{X}^T \mathbf{X})^{-1} \mathbf{X}^T \boldsymbol{\epsilon} | \mathbf{X})\)。 \(= [(\mathbf{X}^T \mathbf{X})^{-1} \mathbf{X}^T] \text{Var}(\boldsymbol{\epsilon} | \mathbf{X}) [(\mathbf{X}^T \mathbf{X})^{-1} \mathbf{X}^T]^T\)。 \(= (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{X}^T (\sigma^2 \mathbf{I}) \mathbf{X} (\mathbf{X}^T \mathbf{X})^{-1}\)。 \(= \sigma^2 (\mathbf{X}^T \mathbf{X})^{-1} (\mathbf{X}^T \mathbf{X}) (\mathbf{X}^T \mathbf{X})^{-1} = \sigma^2 (\mathbf{X}^T \mathbf{X})^{-1}\)。 证毕。

3.11 章末学习闭环

逐项自检：能否解释多元系数的条件含义；能否检查残差、杠杆与 VIF；能否区分 OLS 推断条件与因果识别；能否核验 M02 hash 并报告 LPM 基线。禁止从观测回归自动推出因果，也禁止全样本预处理。常见误区是机械按 VIF 删列、把高 \(R^2\) 当有效外推。

不看正文回答：OLS 的条件均值假设是什么？异方差首先影响什么？LPM 概率越界如何报告？迁移任务：为营销响应构造只作概率基线的 LPM。下一章复用线性预测子与概率评价，转入分类模型。

展开检索答案 核心条件是给定设计下误差条件均值为零；异方差使经典标准误失效但不自动造成系数偏误；保留原始越界计数并明确 Brier 的裁剪规则。

出口决策：MP3-1、MP3-2、MP3-3 均为不可补偿的必过项；任一项未通过，即使总分达标也不得进入第 4 章。三项均通过、其余目标至少一项有证据且闭卷检索至少两题正确，才完成 小节 3.11。

  • MP3-1 条件系数：未能写出“保持其他变量不变”的条件含义者，回到 式 3.28，并以练习 6 的交互边际式作为锚点；补交一份逐项系数解释和条件边际关联证据，再完成一张变量名、量纲和交互基准均不同的新题。新题 2 分，条件对象 1 分、交互边际式 1 分；满 2 分才重新进入本出口判断。
  • MP3-2 训练期预处理与开发切分：在全样本拟合缩放、分箱或多项式基，或让验证/测试参与选择者，回到 表 3.4，并以“训练拟合—验证比较—测试封存”的代码顺序作为锚点；补交带成员与时点断言的训练/验证证据，再完成一份使用不同切分日期和不同预处理器的新代码题。新题 2 分，训练期拟合 1 分、验证只评分且测试封存 1 分；满 2 分才重新进入本出口判断。
  • MP3-3 非因果边界：把观察性系数、转折点或交互项称为干预效应者，回到 式 3.29，并以练习 5 的系数边界说明作为锚点；补交“估计对象—可能混杂—允许结论”三栏证据，再完成一份陌生行业观察回归的新解释题。新题 2 分，明确条件关联 1 分、拒绝未经识别的因果/政策结论 1 分；满 2 分才重新进入本出口判断。
Hastie, Trevor, Robert Tibshirani, 和 Jerome Friedman. 2009年. The Elements of Statistical Learning: Data Mining, Inference, and Prediction. 2 版. Springer. https://doi.org/10.1007/978-0-387-84858-7.
James, Gareth, Daniela Witten, Trevor Hastie, Robert Tibshirani, 和 Jonathan Taylor. 2023年. An Introduction to Statistical Learning: With Applications in Python. Springer. https://doi.org/10.1007/978-3-031-38747-0.