← 返回模块
3.2.3.4beta 可读 · 未来付费内容校验中内容版本 2026-05-27

面向 quant 的 scipy.optimize 与 scipy.linalg

3.2.3 · SciPy 与统计工具 · 编程

周一上午十点，浦东一家中型私募的研究台。PM 把 3.2.2 L5 那张已经稳定跑通的 tear-sheet 推过来，篮子是 {600519.SH, 000001.SZ, 600036.SH}、252 个交易日的 NumPy 收益矩阵 returns，形状 (252,3)(252, 3)。「我现在不要单只票的 alpha 也不要 Sharpe——给我四个数：第一，这只 3 票篮子的​​最小方差​​长仓权重；第二，顶端主成分占多少方差，看篮子风险是不是被一个共因子吃掉；第三，按 510300.SH 月度平值 call 的隐含波动率（implied volatility, IV）算一张 Black-Scholes 模型理论价；第四，同一张 call 用数值积分跑一遍，跟闭式价对账。」四个问题，四种调用：scipy.optimize.minimize、scipy.linalg.eigh、scipy.integrate.quad、外加一段闭式 BS。本课把这三个 SciPy 子包压成一节，并在收尾把四件事拧成一函数。

一、scipy.optimize.minimize + SLSQP：约束最小方差组合

问题陈述。给定 (N,N)(N, N) 协方差矩阵 Σ\Sigma（从 3.2.2 L5 的 df.cov() 直接喂进来，年化后 Sigma = np.cov(returns, rowvar=False) * 252），在长仓全额投资的两个约束下最小化组合方差：

min⁡w∈RN  w⊤Σws.t.1⊤w=1,w≥0.\min_{w \in \mathbb{R}^N} \; w^\top \Sigma w \quad \text{s.t.} \quad \mathbf{1}^\top w = 1, \quad w \ge 0.

这是一道​​二次规划​​（quadratic program）——目标光滑、等式约束加不等式约束都齐。SciPy 的接口是 scipy.optimize.minimize，求解器选 ​SLSQP​​（Sequential Least Squares Programming，序列二次规划）。一行规则：SLSQP handles smooth objective with equality + inequality constraints （最小方差长仓正中靶心）；L-BFGS-B 只吃 bounds 的拟牛顿（只有上下界时用它）；Nelder-Mead 无导数（仅在目标非光滑且维数小才考虑）；trust-constr 处理一般非线性约束（约束本身非线性时用，比如波动率预算依赖 ww 的二次项）。

代码反射。等式约束写成 {'type': 'eq', 'fun': ...} 字典；非负约束直接用 bounds；初值用等权 1/N：

import numpy as np
from scipy.optimize import minimize
N = Sigma.shape[0]
w0 = np.full(N, 1.0 / N)
objective = lambda w: w @ Sigma @ w
constraints = [{'type': 'eq', 'fun': lambda w: w.sum() - 1.0}]
bounds = [(0.0, 1.0)] * N
result = minimize(objective, x0=w0, method='SLSQP', bounds=bounds, constraints=constraints, options={'ftol': 1e-10, 'maxiter': 200})
assert result.success

result.x 是最优权重向量，result.fun 是达到的组合方差，result.message 在收敛失败时给一条人类可读的提示。三个可行性自查在课后练习里强制：assert result.success 必须为真、np.isclose(result.x.sum(), 1.0) 必须成立、(result.x >= -1e-9).all() 允许一个量级 10−910^{-9} 的数值松弛。ftol=1e-10 是为了让 SLSQP 在二次目标上收得足够紧；maxiter=200 在三票篮子上完全够用，仅在病态协方差或高维情形下需要放大。在 A 股 3 票篮子（茅台、平安、招行）上跑出来的示意权重大约 [0.35, 0.30, 0.35]（​​仅作教学示意，非投资建议​​）；之所以三只票权重接近均匀，是因为协方差结构本身决定了最优解会沿主对角分散——一旦把篮子换成 10 只同业银行股，最小方差会大幅集中到方差最低的两三只票上，这是协方差矩阵主导优化的最直观信号。

最小方差只是 Markowitz ​均值方差优化​​框架最简单的一个落点：把目标改成 μ⊤w−λ2w⊤Σw\mu^\top w - \tfrac{\lambda}{2} w^\top \Sigma w、加入预期收益估计 μ\mu、风险厌恶 λ\lambda、行业中性、单票上限、交易成本、跟踪误差等约束之后，整套方法论留给 4.4.x。值得一句话点名的兄弟接口：scipy.optimize.differential_evolution / basinhopping / dual_annealing 是全局优化器；scipy.optimize.root / brentq / newton 是根求解；scipy.optimize.least_squares 是带 bounds 与解析 Jacobian 的非线性最小二乘——本课只点名、不动手，详见 SciPy 官方用户指南 optimize 章。

二、scipy.linalg.eigh：对称特征分解与 PCA 顶端方差占比

协方差矩阵是实对称且​​半正定​​——三个性质合在一起意味着所有​​特征值​​实数非负，​​特征向量​​可取正交单位向量。SciPy 的对应函数是 scipy.linalg.eigh，专门吃对称 / Hermitian 矩阵。规则一句：For symmetric / Hermitian matrices use scipy.linalg.eigh, which exploits symmetry, guarantees real eigenvalues, sorts ascending, and is ~2× faster and more numerically stable than np.linalg.eig。np.linalg.eig 不利用对称性、不排序、还可能返回复数虚部数值噪声；在协方差矩阵上选 eigh。

from scipy.linalg import eigh
eigvals, eigvecs = eigh(Sigma)
total_var = eigvals.sum()
pc_shares = eigvals / total_var
top_pc_share = pc_shares[-1]

签名契约：eigvals 升序排列的实数特征值数组，eigvecs 是 (N,N)(N, N) 矩阵且​​列向量​​为对应特征向量。重建一次以确认：

assert np.allclose(Sigma, eigvecs @ np.diag(eigvals) @ eigvecs.T)

PCA 读法。组合总方差等于 Σ\Sigma 的迹（trace），也等于特征值之和：eigvals.sum() == np.trace(Sigma) 这条恒等式可以用 np.allclose 验证。第 kk 个​​主成分分析​​（principal component analysis, PCA）方向的方差占比即 eigvals[k] / eigvals.sum()；由于 eigh 升序排列，顶端 PC 在 eigvals[-1]。在三票 A 股篮子（茅台、平安、招行——后两者同属申万一级行业「银行」）上典型值在 0.76 左右，反映篮子风险被一个准市场因子主导，仅作示意。一条经验阈值：顶端 PC 占比超过 70% 通常意味着篮子风险高度集中在一个共同因子（A 股宽基里多半是市场因子）；低于 50% 才能算分散。生产研究里用顶端 PC 占比做日度监控，是发现「篮子悄悄变得越来越像一个 beta」的最便宜信号。

边界条件。​​生产级 PCA 的标准做法不是 eigh(cov)，而是对中心化数据矩阵直接做 SVD（奇异值分解）​​：np.linalg.svd 或 scipy.linalg.svd 跳过显式构造协方差，对近奇异矩阵在数值上更稳——本节用 eigh(Sigma) 是因为它一行讲得清，方法论与因子模型构造（决定保留多少 PC、用 varimax / promax 旋转到可解释因子、与 Fama-French / Carhart 风格因子模型对接）留给 2.4.x 与 4.3.x。scipy.linalg 里另外两个 quant 反射常用的接口是 cho_solve（基于 Cholesky 分解的对称正定线性系统求解）与 solve_triangular（三角线性系统）——一并点名。

三、scipy.integrate.quad：Black-Scholes 期望积分与闭式对账

第三件事，按 510300.SH 平值近月 call 算价。在风险中性测度下，标的服从对数正态分布（lognormal），欧式 call 期望支付折现即为：

C(S0,K,r,σ,T)=e−rT∫0∞max⁡(ST−K,0)⋅fST(ST) dST,ST∼Lognormal ⁣(ln⁡S0+(r−12σ2)T,  σ2T).C(S_0, K, r, \sigma, T) = e^{-rT} \int_0^{\infty} \max(S_T - K, 0) \cdot f_{S_T}(S_T) \, dS_T, \quad S_T \sim \mathrm{Lognormal}\!\left(\ln S_0 + \left(r - \tfrac{1}{2}\sigma^2\right)T, \; \sigma^2 T\right).

数值积分一行交付——scipy.integrate.quad(f, a, b) 是自适应高斯-勒让德求积，对一维光滑被积函数稳。示意参数取 510300.SH ETF 期权场景：S_0 = 5.0、K = 5.0（at-the-money）、r = 0.018（CCDC 短端国债收益率 2024 年水平）、sigma = 0.20（沪深300 ETF 期权常见 IV）、T = 30/252（约一个月）。注意 A 股个股期权尚未普遍上市，工作样例落在 510300.SH 这类已上市的 ETF 期权上。

from scipy.integrate import quad
from scipy.stats import lognorm
def integrand(S_T, K): return np.maximum(S_T - K, 0.0) * lognorm.pdf(S_T, s=sigma*np.sqrt(T), scale=S_0*np.exp((r - 0.5*sigma**2)*T))
value, abserr = quad(integrand, 0, np.inf, args=(K,), limit=200)
call_price = np.exp(-r * T) * value

quad 返回 (value, abserr)：积分值与绝对误差估计——abserr 是 QUADPACK 内部对积分值绝对误差的估计上界，对一维光滑被积函数而言通常远低于我们后面给 np.isclose 的 rtol=1e-4，所以拿来当数值断言的护栏足够稳。args=(K,) 把行权价以位置参数形式转给被积函数；limit=200 给一个不至于在默认上界提前退出的子区间预算；np.inf 直接当作积分上界传入，QUADPACK 内部会做无穷区间到有限区间的变量替换，不需要手工截断。需要一次对一个数组的积分时换 scipy.integrate.quad_vec；金融里通常无需——闭式公式或解析展开多半已经存在。scipy.integrate 的其它子接口 dblquad、nquad（高维数值积分）与 solve_ivp（常微分方程数值解）属于一并点名、不动手的层级。

闭式价做横向对账。Black-Scholes 公式把这个一维积分压成两次正态 CDF 调用：

C(S0,K,r,σ,T)=S0 Φ(d1)−Ke−rT Φ(d2),d1,2=ln⁡(S0/K)+(r±12σ2)TσT.C(S_0, K, r, \sigma, T) = S_0 \, \Phi(d_1) - K e^{-rT} \, \Phi(d_2), \quad d_{1,2} = \frac{\ln(S_0 / K) + (r \pm \tfrac{1}{2}\sigma^2)T}{\sigma \sqrt{T}}.

def bs_call(S_0, K, r, sigma, T):
    from scipy.stats import norm
    d1 = (np.log(S_0 / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)
    return S_0 * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2)
closed_form = bs_call(S_0, K, r, sigma, T)
assert np.isclose(call_price, closed_form, rtol=1e-4)

公式推导（delta hedging、风险中性论证、PDE 解法）见 1.4.3 与 2.7.x；本课只把 quad 当作一行配方，并以闭式价做数值断言。隐含波动率反解（scipy.optimize.brentq 求 σ\sigma 使 BS(sigma) - market_price == 0）、波动率曲面构造、局部 / 随机波动率模型留给 1.4.x。

四、Capstone：四件套压成一函数

把三段串起来。给定 252 天 (N,3)(N, 3) 的收益矩阵与一组期权参数，一次性返回 min_var_weights、min_var_variance、top_pc_variance_share、call_price 四个键。课后练习要求你交付这一个函数，并在 pytest 里加形状与范围自查（权重和为 1、非负；顶端 PC 占比落在 [0,1][0, 1]；call 价为正、上界 S0S_0 来自无套利）。

Exercise

实现函数 build_basket_analytics(returns: np.ndarray, S_0: float, K: float, r: float, sigma_implied: float, T: float) -> dict，返回 4 键字典：min_var_weights（SLSQP 解 w⊤Σww^\top \Sigma w 最小化，约束 sum(w) == 1 且 0 <= w <= 1，其中 Sigma = np.cov(returns, rowvar=False) * 252）、min_var_variance（result.fun）、top_pc_variance_share（scipy.linalg.eigh 的最大特征值除以特征值之和）、call_price（用 scipy.integrate.quad 对风险中性对数正态被积函数积分，参数 s = sigma_implied * np.sqrt(T)、scale = S_0 * np.exp((r - 0.5 * sigma_implied**2) * T)，再乘以折现因子 np.exp(-r * T)）。注意 sigma_implied 是期权的隐含波动率，与 Sigma 蕴含的历史回报波动率是两个不同的量。

提示
提示一：SLSQP 需要 constraints=[{"type": "eq", "fun": lambda w: w.sum() - 1.0}] 与 bounds=[(0.0, 1.0)] * N；初值 w0 = np.full(N, 1/N)。scipy.optimize.minimize 的 options 建议带 ftol=1e-10。
提示
提示二：eigh 返回升序 eigvals；顶端 PC 占比是 eigvals[-1] / eigvals.sum()。对 quad，积分 max(S_T - K, 0) * lognorm.pdf(...) 从 0 到 np.inf，再用 exp(-r * T) 折现。

五、收尾：工具与方法论的边界

到这里你能在 252 天 A 股篮子上跑出最小方差权重、读出顶端主成分占比、用 quad 给 510300.SH 平值 call 数值定价并与 Black-Scholes 闭式价在 10−410^{-4} 相对容差内对账。两条边界请先记住。​​到 Track 4 的边界​​：最小方差是组合构造里最简单的一处落点；均值方差优化、风险平价（risk parity）、Black-Litterman、稳健优化（robust optimization）、带交易成本的优化、因子组合等方法论住在 4.4.x，PCA 因子模型构造住在 4.3.x——本课只搭工具。​​到性能的边界​​：SLSQP 对 NN 几百量级毫无压力；上千维或更复杂的凸约束就该上 cvxpy、clarabel、osqp 这类现代凸 / 锥规划求解器，性能调优是 3.3.1 的事。最后一条全局优化（differential_evolution、basinhopping、dual_annealing）、根求解（brentq、newton）、least_squares 与高维积分（dblquad、nquad）、ODE 求解（solve_ivp）这些 SciPy 兄弟接口在此点名、不动手。

3.2.3 这一模块到此收束：你已经把 scipy.stats 的分布对象、假设检验与 bootstrap、linregress + curve_fit 的回归与曲线拟合、以及今天这节的 optimize + linalg + integrate 三件套全部压进反射。下一模块 3.2.4 Synthetic Data & APIs 解决一个上游问题——前面所有工具吃进去的 returns 矩阵从哪儿来？当真实行情数据要授权、要清洗、要拼接时，如何用 SciPy 的分布对象与可控随机种子构造一只统计性质可控的合成篮子，并把它包装成可在线上 / 线下复用的数据接口。