用 SciPy 做回归与曲线拟合

周三午后，浦东陆家嘴一家中型私募的研究台上，PM 把一张 252 天的样本期跑出来推过来：「600519.SH 对沪深300 ETF（510300.SH）的 beta 我刚才用 np.linalg.lstsq 解出来是 0.91——但 0.91 离 1 到底有多远？是抽样噪音里飘出来的一格，还是这只票就比沪深300 系统性低 beta？」3.2.1 L3 教过你三种从设计矩阵 $X$ 解出 $(\alpha, \beta)$ 的等价方法：np.linalg.lstsq、正规方程 np.linalg.solve(X.T @ X, X.T @ y)、以及 Cholesky 分解。三种都只返回​​点估计​​就停了。研究台、tear-sheet、面试题里的下一句话永远是「斜率的​​标准误​​是多少？β 是否显著偏离 1？」——这一步要的是​​普通最小二乘​​（ordinary least squares, OLS）的标准误与 t 统计量，np.linalg.lstsq 不出。scipy.stats.linregress 出。

一、scipy.stats.linregress 作为一行 OLS 配方

把 252 天的市场代理收益与个股收益分别叫 r_mkt（沪深300 ETF 510300.SH）与 r_stock（贵州茅台 600519.SH），都是日简单收益（daily simple return）。SciPy 的反射动作只有三行：

import scipy.stats as st
result = st.linregress(r_mkt, r_stock)
alpha, beta, r, p, beta_se, alpha_se = result.intercept, result.slope, result.rvalue, result.pvalue, result.stderr, result.intercept_stderr

返回值是 LinregressResult 命名元组（named tuple），六个字段一字排开：slope 是回归斜率，即 $\beta$，写作 $y = \beta x + \alpha$；intercept 是截距 $\alpha$；rvalue 是 Pearson 相关系数，其平方即决定系数 R²；pvalue 是对零假设 $H_0: \beta = 0$ 的双侧 Wald 检验 p 值；stderr 是斜率的标准误；intercept_stderr 是截距的标准误（SciPy 1.6 起加入，现版默认可用）。承重等价关系是 t 统计量的定义 $t = \text{slope} / \text{stderr}$；pvalue 即 t 分布（自由度 $N - 2$）落在 $\pm|t|$ 之外的尾部面积。一个示意数：在 2024 年某段 252 个交易日上跑出 beta ≈ 0.91、beta_se ≈ 0.06，rvalue ≈ 0.78——rvalue 的平方告诉你沪深300 解释了样本里 60% 出头的 600519.SH 日收益方差。

二、检验 H₀: β = 1.0 与 95% 置信区间

回到 PM 的问题。研究台真正想问的零假设不是 $H_0: \beta = 0$（那对一只蓝筹股没有意义），而是 $H_0: \beta = 1.0$——即「这只票是否与沪深300 一比一同步移动」。检验的 t 统计量改用 1.0 作为基准：

t_one = (beta - 1.0) / beta_se
p_one = 2 * (1 - st.t.cdf(abs(t_one), df=252-2))
ci_low, ci_high = beta - 1.96 * beta_se, beta + 1.96 * beta_se

自由度 df = 252 - 2 = 250：252 个观测减去两个被估参数（截距与斜率）。p_one 把 t 分布在 $\pm|t_{\text{one}}|$ 之外的累积面积乘 2，即双侧 p 值；若 p_one < 0.05，在 5% 水平拒绝「β 等于 1」。95% 置信区间 [beta - 1.96 * beta_se, beta + 1.96 * beta_se] 走的是高斯近似（Gaussian approximation），大样本下渐近有效；小样本要把 1.96 换成 st.t.ppf(0.975, df=N-2)。把示意数代进去：t_one ≈ (0.91 - 1.0) / 0.06 ≈ -1.5，CI 大约是 [0.79, 1.03]——区间盖住了 1.0，单凭这一条样本你写不下「显著低于市场 beta」。这才是 PM 真正要带去风控会议的句子。

三、linregress 与 np.linalg.lstsq 对照，以及 statsmodels 毕业路径

承接 3.2.1 L3 的对照规则：对一个特征加一个截距的 OLS，st.linregress(x, y) 与 np.linalg.lstsq(np.column_stack([np.ones_like(x), x]), y) 解出的 (intercept, slope) 点估计是逐位相同的；只有 linregress 同时返回 rvalue、pvalue、stderr、intercept_stderr。np.linalg.lstsq 的优势在于直接吃形如 $(N, p)$ 的多特征设计矩阵——一旦特征超过一个，linregress 不可用。三句决策规则：单特征加截距、需要标准误用 linregress；多特征、只要系数用 np.linalg.lstsq；多特征还要标准误、t 统计量、F 检验、稳健协方差，毕业到 statsmodels。

OLS 的标准误公式承重一个假设：残差 i.i.d.、零均值、方差恒定（同方差）。日收益数据上两个假设都被打破：波动有时间聚集（异方差），rolling 残差有自相关。研究台的标准修法是异方差稳健（heteroskedasticity-consistent，Eicker-Huber-White）标准误，或对时间序列的 Newey-West / HAC 自相关稳健标准误。毕业规则一句：当你需要异方差稳健（HC0/HC1/HC2/HC3）、自相关稳健（Newey-West / HAC）标准误、多特征回归加 F 检验、GLM、或面板回归时，毕业到 statsmodels.api.OLS(y, X).fit().get_robustcov_results(cov_type="HC1").summary()。statsmodels 是另一个独立的库，本课只点名 API、不动手；何时切过去由残差诊断决定。

四、scipy.optimize.curve_fit 与 Nelson-Siegel 利率曲线

研究台第二常用的回归不是线性 OLS，而是把一条​​收益率曲线​​（yield curve）拟合到一个参数化函数上。SciPy 的工具是 scipy.optimize.curve_fit，签名 popt, pcov = curve_fit(f, xdata, ydata, p0=initial_guess, maxfev=10000)，其中 f 必须写成 f(x, *params) -> y。

工作例子。5 个期限点的国债收益率曲线，tau = [1, 3, 5, 7, 10]（年），从中央国债登记结算有限责任公司（CCDC）每日发布的中债到期收益率曲线上采样——2024 年年中的近似水平 observed_yields = [0.018, 0.020, 0.022, 0.024, 0.025]（即 1.8% / 2.0% / 2.2% / 2.4% / 2.5%）。拟合一段 ​Nelson-Siegel 模型​​单段形式：

三个参数 $(b_0, b_1, \lambda)$ 分别承担长期水平、短端斜率与衰减时间尺度。代码：

import numpy as np
from scipy.optimize import curve_fit
def ns(tau, b0, b1, lam): return b0 + b1 * (1 - np.exp(-tau / lam)) / (tau / lam)
popt, pcov = curve_fit(ns, tau, observed_yields, p0=[0.03, -0.01, 2.0], maxfev=10000)
sigma_params = np.sqrt(np.diag(pcov))

popt 是拟合后的参数向量（NumPy 数组）；pcov 是参数的协方差矩阵——对角线即各参数的方差，开根号即标准误。初值 p0=[0.03, -0.01, 2.0] 不能省：curve_fit 默认 method='lm' 走 Levenberg-Marquardt 这一类局部优化器，初值离谱直接发散；maxfev=10000 给一个不至于在默认 100 次评估上提前退出的迭代预算。若传入 bounds，方法自动落到 trust-region reflective（method='trf'）。拟合完做一次自查：

curve_residuals = observed_yields - ns(tau, *popt)
rss = np.sum(curve_residuals**2)

rss（residual sum of squares）越小越好；五个观测点拟合三个参数的情况下，你期望 rss 落在 $10^{-6}$ 数量级，因为这条曲线本身就贴近 Nelson-Siegel 的形状。生产级利率曲线工程（中债综合财富指数曲线、Nelson-Siegel-Svensson 扩展、平滑样条无套利曲线）见 1.2.x 固定收益方向；本课只把 Nelson-Siegel 当作 curve_fit 的教学样例，不当作建模方法论。

五、残差诊断与 OLS 的承重假设

任何回归在引用「标准误」之前都要做两件诊断。其一，残差对拟合值的散点图：肉眼看是否有「漏斗形」扩散，提示异方差；本课只描述这一步，不动 matplotlib（可视化不是承重）。其二，残差的正态性检验：调用 L2 教过的 Shapiro-Wilk。

residuals = r_stock - (alpha + beta * r_mkt)
shapiro_stat, shapiro_p = st.shapiro(residuals)

零假设是「残差服从​​正态分布​​（Gaussian distribution，又称高斯分布）」。若 shapiro_p < 0.05，拒绝正态——此时 linregress 报的 stderr 不再严格有效，正确做法是用 L2 的 scipy.stats.bootstrap 套一个 lambda x, y: st.linregress(x, y).slope 跑出 95% 的非参数置信区间。OLS 标准误的承重等价关系：「残差 i.i.d.、同方差、近似正态」⇨「linregress 报的 stderr 严格有效」；任一项被打破，要么换稳健标准误（毕业到 statsmodels），要么换非参数 CI（bootstrap）。OLS 估计为何无偏、标准误公式为何长那样、Gauss-Markov 定理等推导见 2.2.2 linear-regression-ols-and-geometry；本课把 linregress 与 curve_fit 当成两条可以一行交付的配方。

下一课接力

到这里你能在两条日收益序列上交付一组带标准误的 OLS 估计、检验 β 是否显著偏离 1.0、给出 95% CI；再用 curve_fit 把一条 5 点 CCDC 利率曲线拟合到 3 参数 Nelson-Siegel 上、并报出参数标准误。残差正态性诊断与 statsmodels 毕业 API 也都在台面上。L4 把视线收回到组合层：给定一只 $(N, N)$ 的协方差矩阵（3.2.2 L5 的 df.cov() 直接喂进来），用 scipy.optimize.minimize 的 SLSQP 解出最小方差组合的权重；同一只矩阵丢进 scipy.linalg.eigh 拿到对称特征分解，是 2.4 主成分分析的入口；再用 scipy.integrate.quad 把 Black-Scholes 期望积分跑一遍，与闭式价做一次数值对账——三件套压成一节课。

练习