元信息：Ruey S. Tsay《Analysis of Financial Time Series》（3rd ed., Wiley）精读笔记，负责范围 PDF 第 457–707 页（第 8 章 8.6.3 节中段起，至第 12 章末、全书索引及尾部空白页）。

第 8 章多元时间序列分析及其应用（Multivariate Time Series Analysis and Its Applications）（接上一块）

8.6.3 协整检验：Johansen 检验（PDF p.457）（接上一块）

上一块已介绍误差修正模型（ECM）：$\Delta x_t=\mu_t+\Pi x_{t-1}+\sum_{i=1}^{p-1}\Phi_i^*\Delta x_{t-i}+a_t$，其中 $\Pi=\alpha\beta'$，$\mathrm{Rank}(\Pi)=m$ 即协整向量个数。本页接续检验部分。

偏调整（partial adjustment）：先把 $\Delta x_t$ 与 $x_{t-1}$ 分别对确定项 $d_t$ 及滞后差分 $\Delta x_{t-i}$（$i=1,\dots,p-1$）做多元线性回归，得到残差 $\hat u_t$（调整后的 $\Delta x_t$）和 $\hat v_t$（调整后的 $x_{t-1}$）。检验方程化为 $\hat u_t=\Pi\hat v_t+a_t .$
在正态假设下，$\Pi$ 秩的似然比检验可通过 $\hat u_t$ 与 $\hat v_t$ 的典型相关分析（canonical correlation analysis）完成；这里的典型相关是 $\Delta x_t$ 与 $x_{t-1}$ 在剔除 $d_t$ 与 $\Delta x_{t-i}$ 影响后的偏典型相关。$\{\hat\lambda_i\}$（降序排列）为 $\hat u_t,\hat v_t$ 之间的平方典型相关系数。
迹检验（trace test）：$H_0:\mathrm{Rank}(\Pi)=m$ vs $H_a:\mathrm{Rank}(\Pi)>m$，Johansen (1988) 的 LR 统计量 $LR_{tr}(m)=-(T-p)\sum_{i=m+1}^{k}\ln(1-\hat\lambda_i). \tag{8.42}$ 若秩为 $m$，则 $i>m$ 的 $\hat\lambda_i$ 应很小，统计量也小。由于存在单位根，渐近分布不是卡方，而是标准布朗运动的泛函，临界值须模拟得到。
最大特征值检验（maximum eigenvalue test）：序贯检验 $H_0:\mathrm{Rank}=m$ vs $H_a:\mathrm{Rank}=m+1$， $LR_{\max}(m)=-(T-p)\ln(1-\hat\lambda_{m+1}),$ 同样为非标准分布，临界值靠模拟。
使用方式：从 $m=0$ 开始依次检验，第一个不被拒绝的 $m$ 即协整秩估计。

8.6.4 协整 VAR 模型的预测（PDF p.457）

拟合好的 ECM 在给定估计参数的条件下先预测差分序列 $\Delta x_t$，再累加得到 $x_t$ 的预测。与普通 VAR 预测的区别：ECM 预测强制施加了协整关系（长期均衡约束），因而长期预测中各分量保持协整联系。

8.6.5 实例：美国 3 个月与 6 个月国库券利率（PDF p.458–462）

数据：1958-12-12 至 2004-08-06 的周度二级市场 3 个月（tb3m）与 6 个月（tb6m）国库券利率，2383 个观测，来自圣路易斯联储。图 8.12 显示两者高度同步。软件为 S-Plus（VAR、coint、VECM）。
单变量 ADF 检验（AR(3)）：统计量 −2.34 与 −2.33，p 值约 0.16，不能拒绝单位根。对 $x_t=(tb3m_t,tb6m_t)'$ 用 BIC 选出 VAR(3)。
确定项选受限常数（restricted constant, trend='rc'），因为没有理由认为利率有漂移；lags = p-1 = 2。
检验结果：特征值 0.0322、0.0023。迹统计量 $H(0)$：83.27（5% 临界值 19.96，1% 为 24.60），$H(1)$：5.49（临界 9.24/12.97）；最大特征值统计量 $H(0)$：77.78（15.67/20.20），$H(1)$：5.49。结论：恰有一个协整向量。
VECM 极大似然估计：协整向量标准化为 $(1,-1.0124)$（标准误 0.0086），截距 0.2254（标准误 0.0545）。调整系数（loading）$\alpha=(-0.0949,-0.0211)'$，t 值 −4.76、−1.18。拟合模型近似为 $\Delta x_t=\begin{bmatrix}-0.09\\-0.02\end{bmatrix}(w_{t-1}+0.23)+\begin{bmatrix}0.05&0.27\\-0.04&0.32\end{bmatrix}\Delta x_{t-1}+\begin{bmatrix}-0.21&0.25\\-0.03&0.10\end{bmatrix}\Delta x_{t-2}+a_t,$ 其中平稳序列 $w_t\approx tb3m_t-tb6m_t$，均值约 −0.225；残差标准差 0.20、0.18；$R^2$ 分别 0.108、0.091。
图 8.13 协整残差：1980 年代初高利率高波动时期出现大残差。
预测（原点 2004-08-06）：差分序列 1 步预测 (−0.0378, −0.0642)，标准误 (0.2009, 0.1807)；10 步预测 (−0.2276, −0.1314)，标准误 (0.846, 0.816)。水平预测（levels=T）：1 步 (1.4501, 1.7057)，10 步 (1.4722, 1.7078)。图 8.14、8.15 展示，95% 逐点区间因单位根非平稳而很宽、信息量有限。
备注：R 包 urca 的 ca.jo 可做 Johansen 检验（见配对交易一节）。

8.7 门限协整与套利（Threshold Cointegration and Arbitrage）（PDF p.462–466）

目标：用多元时间序列方法识别指数期现套利机会，并说明第 4 章的单变量非线性模型可结合协整推广到多元。

持有成本模型（cost-of-carry model）：设 $f_{t,\ell}$ 为到期日 $\ell$ 的 S&P 500 指数期货对数价格，$s_t$ 为现货成分股对数价格， $f_{t,\ell}-s_t=(r_{t,\ell}-q_{t,\ell})(\ell-t)+z_t^*, \tag{8.43}$ $r_{t,\ell}$ 为无风险利率，$q_{t,\ell}$ 为股息率，$\ell-t$ 为剩余期限。$z_t^*$ 必须单位根平稳，否则存在持续套利机会。套利交易：当对数价差偏离持有成本过大时，同时买（卖空）现货、卖（买）期货。只有当 $|z_t^*|$ 超过由交易成本及其他经济、风险因素决定的某一水平时套利才有利可图。
调整利率与股息后 $f$ 与 $s$ 协整，协整向量 $(1,-1)$，协整序列为 $z_t^*$。因此应对收益 $r_t=(\Delta f_t,\Delta s_t)'$ 用误差修正形式建模。

8.7.1 多元门限模型（Multivariate Threshold Model）（PDF p.464）

套利交易本身改变市场动态，故模型随是否存在套利而切换：

r_t=\begin{cases}c_1+\sum_{i=1}^p\Phi_i^{(1)}r_{t-i}+\beta_1 z_{t-1}+a_t^{(1)}, & z_{t-1}\le\gamma_1,\\ c_2+\sum_{i=1}^p\Phi_i^{(2)}r_{t-i}+\beta_2 z_{t-1}+a_t^{(2)}, & \gamma_1<z_{t-1}\le\gamma_2,\\ c_3+\sum_{i=1}^p\Phi_i^{(3)}r_{t-i}+\beta_3 z_{t-1}+a_t^{(3)}, & \gamma_2<z_{t-1},\end{cases}\tag{8.44}

其中 $z_t=100z_t^*$（因 $z_t^*$ 数值很小），$\gamma_1<0<\gamma_2$ 为门限，$z_{t-1}$ 为门限变量（一般可在 $d\in\{1,\dots,d_0\}$ 中选 $z_{t-d}$），各区制的 $a_t^{(i)}$ 为相互独立的二维白噪声。

这是第 4 章 TAR 模型与式 (8.33) ECM 的共同推广，称三区制多元门限模型。
直观：只有 $|z_t|$ 大时套利才盈利，故套利只发生在区制 1 和 3；中间区制由正常市场力量主导，两价格近似随机游走、不受协整约束，计量上 $\beta_2$ 应不显著。这一现象称门限协整（threshold cointegration）（Balke & Fomby, 1997）。

8.7.2 数据（PDF p.465）

1993 年 5 月 S&P 500 指数及其在 CME 交易的 6 月期货合约的日内成交数据（Forbes, Kalb & Kofman, 1999），构造分钟级双变量价格序列，共 7060 个观测。为避免异常值影响，将 10 个极端值（两侧各 5 个）替换为相邻两点均值；不考虑条件异方差。图 8.16 为期货、现货 1 分钟对数收益及 $z_t$ 的时序图。

8.7.3 估计（PDF p.465–467）

完整设定包括选门限变量、区制数及各区制阶数 $p$（参见 Tsay 1998）。门限可用 AIC 或残差平方和等信息准则估计。
设定 $p=8$，$d\in\{1,2,3,4\}$，$\gamma_1\in[-0.15,-0.02]$，$\gamma_2\in[0.025,0.145]$，各区间 300 个格点做网格搜索；AIC 选 $z_{t-1}$ 为门限变量，$\hat\gamma_1=-0.0226$，$\hat\gamma_2=0.0377$。三区制样本量分别 2234、2410、2408。
表 8.8 主要结论：
1. 中间区制 $z_{t-1}$ 系数（$\Delta f_t$ 方程 −0.00010，t=−0.30；$\Delta s_t$ 方程 0.00012，t=0.86）在 5% 水平不显著，证实无套利机会时两者无协整；区制 1、3 中 $\Delta s_t$ 方程的 $z_{t-1}$ 系数高度显著（t=10.47、9.75）。
2. $\Delta f_t$ 在三个区制都负依赖于 $\Delta f_{t-1}$（−0.085、−0.039、−0.041），与第 5 章的买卖价反弹（bid–ask bounce）一致。
3. 期货滞后收益比现货滞后收益更有信息量（$\Delta s_t$ 方程中 $\Delta f_{t-1},\dots,\Delta f_{t-7}$ 多数 t 值很大），因期货流动性更好——即期货领先现货（价格发现）。

8.8 配对交易（Pairs Trading）（PDF p.466–476）

配对交易是一种市场中性（market-neutral）策略。本节聚焦统计套利配对交易，利用协整与 ECM 思想（延伸阅读：Vidyamurthy 2004；Pole 2007）。核心思想是相对定价：依据套利定价理论（APT），风险特征相近的两只股票价格应大致相同；若出现差距，可能一只高估一只低估，卖高买低，等待错误定价修正。真实价格并不重要，重要的是两观测价格（适当缩放后）应一致；两者（缩放后）的差距称价差（spread），价差越大，错误定价幅度和潜在利润越大。

8.8.1 理论框架（PDF p.466–468）

令 $p_{it}=\ln P_{it}$，假设其为随机游走 $p_{it}=p_{i,t-1}+r_{it}$。若两股票风险因子相近，由 APT 收益相近，价格受共同成分驱动而协整：存在 $w_t=p_{1t}-\gamma p_{2t}$ 单位根平稳、均值回复。两价格满足误差修正形式 $\begin{bmatrix}p_{1t}-p_{1,t-1}\\p_{2t}-p_{2,t-1}\end{bmatrix}=\begin{bmatrix}\alpha_1\\\alpha_2\end{bmatrix}(w_{t-1}-\mu_w)+\begin{bmatrix}\epsilon_{1t}\\\epsilon_{2t}\end{bmatrix},\tag{8.45}$ $\mu_w=E(w_t)$。参数 $\gamma,\mu_w,\alpha_1,\alpha_2$ 可用 ML 或 LS 估计；$w_t$ 称两对数价格的价差。
含义：收益依赖于上期对长期均衡的偏离 $w_{t-1}-\mu_w$；实践中 $\alpha_1,\alpha_2$ 应异号，表示向均衡回复。
组合收益：做多 1 股股票 1、做空 $\gamma$ 股股票 2，第 $i$ 期收益 $r_{p,t+i}=(p_{1,t+i}-p_{1t})-\gamma(p_{2,t+i}-p_{2t})=w_{t+i}-w_t,$ 即组合收益等于价差的增量，与 $\mu_w$ 无关。（注：严格说这是对数价格意义下的近似"收益"。）

8.8.2 交易策略（PDF p.468–469）

思路：利用价差围绕均衡值 $\mu_w$ 的振荡（均值回复）交易：偏离大时建仓，回归时平仓。
设 $\eta$ 为一次配对交易的成本（交易费、保证金利率、两只股票的买卖价差等），$\Delta$ 为目标偏离，在 $2\Delta>\eta$ 条件下：
- 当 $w_t=\mu_w-\Delta$ 时，买 1 股股票 1、卖空 $\gamma$ 股股票 2；
- 当 $w_{t+i}=\mu_w+\Delta$ 时平仓。组合收益 $w_{t+i}-w_t=2\Delta$，净利润 $2\Delta-\eta>0$。只要 $\Delta$ 相对 $w_t$ 标准差不太大，入场点就能出现；均值回复保证出场点会出现。
讨论：若 $\Delta>\eta$，也可在 $w_{t+i}=\mu_w$ 时平仓，净利 $\Delta-\eta$，交易更频繁、成本更高，但持仓期更短。反方向（$w_t=\mu_w+\Delta$）则卖空股票 1、买入 $\gamma$ 股股票 2。$\eta$ 是交易门槛。

8.8.3 简单示例：BHP 与 VALE（PDF p.469–476）

数据：纽交所交易的 BHP Billiton（澳大利亚，自然资源）与 Vale（巴西，金属采矿），同属自然资源行业、风险因子相似。Yahoo Finance 日度复权收盘价，2002-07-01 至 2006-03-31（946 个观测）。图 8.17 显示两者共同运动。$p_{1t}$=BHP，$p_{2t}$=VALE。
最小二乘法（Engle–Granger 两步思路）：回归 $p_{1t}=1.823+0.717p_{2t}+\hat w_t$，$\sigma_w=0.044$，$R^2=0.9899$。残差图 8.18 显示围绕 0 在固定范围内波动，ACF 指数衰减。对 $\hat w_t$ 拟合 AR(2)：$(1-0.805B-0.122B^2)\hat w_t=a_t$，$\sigma_a=0.018$，可分解为 $(1-0.935B)(1-0.130B)$，故平稳；ADF 检验统计量 −6.04，p 值 0.01，拒绝单位根。
- 注意（教材可补）：对估计残差做 ADF 时应使用 Engle–Granger 临界值而非普通 DF 临界值，原书未展开。
极大似然法（Johansen）：信息准则选 VAR(1)（R 中 ar 选 2 阶，对应 K=2），受限常数下检验：特征值 0.0415、0.0082；迹统计量 $H(0)$=47.74（临界 19.96/24.60），$H(1)$=7.77（9.24/12.97）；最大特征值统计量 $H(0)$=39.97（15.67/20.20），$H(1)$=7.77。确认协整秩为 1。
VECM 估计：协整向量 $(1,-0.7177)$（标准误 0.0112），截距 −1.8144；loading $\alpha=(-0.0671,0.0263)'$，t 值 −4.65、1.57；拟合模型 $\Delta x_t=\begin{bmatrix}-0.067\\0.026\end{bmatrix}(w_{t-1}-1.81)+\begin{bmatrix}-0.11&0.07\\0.07&0.04\end{bmatrix}\Delta x_{t-1}+a_t,$ 残差标准差 0.019、0.022。价差 $w_t=p_{1t}-0.718p_{2t}$，均值 1.81，与 LS 结果非常接近；$\hat\gamma=0.718$；如预期 $\alpha_1<0$、$\alpha_2>0$。
交易策略：$w_t$ 标准差 0.044，取 $\Delta=0.045$（略大于 1 个标准差），正态假设下偏离至少 $\Delta$ 的概率约 30%。图 8.19 画出 $\mu_w$、$\mu_w\pm0.045$ 三条线，价差多次在上下边界间穿越，交易机会多；每次配对交易对数收益 $2\Delta=0.09$。更真实的做法应在样本外实施。
选对问题：识别协整股票对是关键，应选风险因子相近的股票，用金融理论指导筛选。
R 演示：library(urca)；lm(bhp~vale) 得截距 1.822648、斜率 0.716664（残差标准误 0.04421，944 自由度）；arima(wt,order=c(2,0,0),include.mean=F) 得 ar1=0.8051、ar2=0.1219；polyroot 求特征根，1/Mod(x) 得 0.935、0.130；ca.jo(xt, ecdet="const", type='trace', K=2, spec='transitory') 得迹统计量 r=0：47.77，r≤1：7.78（10%/5%/1% 临界值 7.52/9.24/12.97）；type='eigen' 得 r=0：40.00。第一协整向量 (1, −0.7177, −1.8285)，loading (−0.0673, 0.0255)。

附录 A 向量与矩阵复习（Review of Vectors and Matrices）（PDF p.476–483）

不给证明，参见 Graybill (1969)。

基本定义：$m\times n$ 实矩阵 $A=[a_{ij}]$；行维数、列维数；对角元 $a_{ii}$；$m\times1$ 为列向量，$1\times n$ 为行向量（"向量"默认指列向量）；方阵、对角阵、单位阵 $I_m$；转置 $A'$，$a'_{ij}=a_{ji}$，$(A')'=A$；$A'=A$ 为对称阵。
基本运算：加减（同维）、数乘、乘法 $AC=[\sum_{v=1}^n a_{iv}c_{vj}]_{m\times q}$（需 $n=p$，称可相乘 conformable）。例：$\begin{bmatrix}2&1\\1&1\end{bmatrix}\begin{bmatrix}1&2&3\\-1&2&-4\end{bmatrix}=\begin{bmatrix}1&6&2\\0&4&-1\end{bmatrix}$。规则：$(AC)'=C'A'$；一般 $AC\ne CA$。
逆、迹、特征值：$A$ 非奇异（可逆）指存在唯一 $C$ 使 $AC=CA=I$，记 $A^{-1}$。迹 $\mathrm{tr}(A)=\sum a_{ii}$，满足 $\mathrm{tr}(A+C)=\mathrm{tr}A+\mathrm{tr}C$，$\mathrm{tr}(A)=\mathrm{tr}(A')$，$\mathrm{tr}(AC)=\mathrm{tr}(CA)$。右特征值/特征向量对：$Ab=\lambda b$；$m$ 个特征值，实矩阵的复特征值成共轭对；$A$ 非奇异当且仅当特征值全非零；$\mathrm{tr}(A)=\sum\lambda_i$，$|A|=\prod\lambda_i$。秩：$AA'$ 非零特征值个数。非奇异时 $(A^{-1})'=(A')^{-1}$。

（附录 A 续，PDF p.478–480）

正定矩阵（positive-definite matrix）：$A$ 对称且全部特征值为正；等价地，对任意非零 $b$ 有 $b'Ab>0$。性质：特征值实且正；谱分解（spectral decomposition） $A=P\Lambda P'$，$\Lambda$ 为特征值对角阵（常记 $\lambda_1\ge\cdots\ge\lambda_m$），$P$ 由单位特征向量 $e_i$（$Ae_i=\lambda_ie_i$，$e_i'e_i=1$）组成，特征值互异时 $e_i'e_j=0$，$P$ 为正交矩阵。
- 例：$\Sigma=\begin{bmatrix}2&1\\1&2\end{bmatrix}$，$\Sigma(1,1)'=3(1,1)'$，$\Sigma(1,-1)'=(1,-1)'$，特征值 3 和 1，单位特征向量 $(1/\sqrt2,1/\sqrt2)'$、$(1/\sqrt2,-1/\sqrt2)'$，且 $P'\Sigma P=\mathrm{diag}(3,1)$。
LDL 与 Cholesky 分解：对称阵 $A$ 存在单位下三角 $L$ 与对角阵 $G$ 使 $A=LGL'$（Strang 1980 第 1 章）；若 $A$ 正定则 $G$ 对角元为正，此时 $A=(L\sqrt G)(L\sqrt G)'$，称 Cholesky 分解（$L\sqrt G$ 仍为下三角，开方逐元素进行）。由此可对角化：$L^{-1}A(L^{-1})'=G$；$L^{-1}$ 仍为单位下三角。
- 例：上面的 $\Sigma$ 有 $L=\begin{bmatrix}1&0\\0.5&1\end{bmatrix}$，$G=\mathrm{diag}(2,1.5)$，$L^{-1}=\begin{bmatrix}1&0\\-0.5&1\end{bmatrix}$，$L^{-1}\Sigma(L^{-1})'=G$。（这正是第 8 章 VAR 结构形式/正交化冲击响应所用的分解。）
向量化与 Kronecker 积：$A=[a_1,\dots,a_n]$，$\mathrm{vec}(A)=(a_1',\dots,a_n')'$ 为 $mn\times1$ 向量（原文下标误写为 $a_m$）。Kronecker 积 $A\otimes C=[a_{ij}C]$，维数 $mp\times nq$。
- 例：$A=\begin{bmatrix}2&1\\-1&3\end{bmatrix}$，$C=\begin{bmatrix}4&-1&3\\-2&5&2\end{bmatrix}$，$\mathrm{vec}(A)=(2,-1,1,3)'$，$\mathrm{vec}(C)=(4,-2,-1,5,3,2)'$， $A\otimes C=\begin{bmatrix}8&-2&6&4&-1&3\\-4&10&4&-2&5&2\\-4&1&-3&12&-3&9\\2&-5&-2&-6&15&6\end{bmatrix}$。
- 性质：(1) 一般 $A\otimes C\ne C\otimes A$；(2) $(A\otimes C)'=A'\otimes C'$；(3) $A\otimes(C+D)=A\otimes C+A\otimes D$；(4) $(A\otimes C)(F\otimes G)=(AF)\otimes(CG)$；(5) $(A\otimes C)^{-1}=A^{-1}\otimes C^{-1}$；(6) 方阵 $\mathrm{tr}(A\otimes C)=\mathrm{tr}(A)\mathrm{tr}(C)$；(7) $\mathrm{vec}(A+C)=\mathrm{vec}A+\mathrm{vec}C$；(8) $\mathrm{vec}(ABC)=(C'\otimes A)\mathrm{vec}(B)$；(9) $\mathrm{tr}(AC)=\mathrm{vec}(C')'\mathrm{vec}(A)=\mathrm{vec}(A')'\mathrm{vec}(C)$；(10) $\mathrm{tr}(ABC)=\mathrm{vec}(A')'(C'\otimes I)\mathrm{vec}(B)=\mathrm{vec}(A')'(I\otimes B)\mathrm{vec}(C)=\cdots$（共六种循环写法）。
半向量化 vech：对 $k\times k$ 对称阵只堆叠主对角线及以下元素，$\mathrm{vech}(A)$ 维数 $k(k+1)/2$；$k=3$ 时 $\mathrm{vech}(A)=(a_{11},a_{21},a_{31},a_{22},a_{32},a_{33})'$。在多元波动率模型（第 10 章 VEC 模型）中大量使用。

附录 B 多元正态分布（Multivariate Normal Distributions）（PDF p.480–481）

密度：$x\sim N_k(\mu,\Sigma)$， $f(x|\mu,\Sigma)=\frac{1}{(2\pi)^{k/2}|\Sigma|^{1/2}}\exp\Big[-\tfrac12(x-\mu)'\Sigma^{-1}(x-\mu)\Big].\tag{8.47}$
二元情形：$\Sigma^{-1}=\frac{1}{\sigma_{11}\sigma_{22}-\sigma_{12}^2}\begin{bmatrix}\sigma_{22}&-\sigma_{12}\\-\sigma_{12}&\sigma_{11}\end{bmatrix}$，$\rho=\sigma_{12}/(\sigma_1\sigma_2)$，$|\Sigma|=\sigma_{11}\sigma_{22}(1-\rho^2)$， $f(x_1,x_2)=\frac{1}{2\pi\sigma_1\sigma_2\sqrt{1-\rho^2}}\exp\Big[-\frac{Q}{2(1-\rho^2)}\Big],\quad Q=\Big(\tfrac{x_1-\mu_1}{\sigma_1}\Big)^2+\Big(\tfrac{x_2-\mu_2}{\sigma_2}\Big)^2-2\rho\Big(\tfrac{x_1-\mu_1}{\sigma_1}\Big)\Big(\tfrac{x_2-\mu_2}{\sigma_2}\Big).$
分块 $x=(x_1',x_2')'$，$x_1$ 为 $p$ 维。性质：
1. $c'x\sim N(c'\mu,c'\Sigma c)$；反之，若任意非零 $c$ 使 $c'x$ 为一元正态，则 $x$ 多元正态。
2. 边际分布 $x_i\sim N_{k_i}(\mu_i,\Sigma_{ii})$。
3. $\Sigma_{12}=0$ 当且仅当 $x_1,x_2$ 独立（正态特有）。
4. $(x-\mu)'\Sigma^{-1}(x-\mu)\sim\chi^2_k$（原文写作 $m$ 自由度，应为 $k$）。
5. 条件分布：$(x_1|x_2=b)\sim N_p\big[\mu_1+\Sigma_{12}\Sigma_{22}^{-1}(b-\mu_2),\ \Sigma_{11}-\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}\big]$。这是正态假设下时间序列预测与递推最小二乘（以及 Kalman 滤波）的基础。

附录 C 部分 SCA 命令（PDF p.481–482）

列出例 8.6（月度 1 年与 3 年期国债利率对数序列）所用 SCA 命令：读数据、取对数、miden 用 AR 拟合 1–8 阶定阶、mtsm 设定 VARMA(2,1) 模型 (i-p1*b-p2*b**2)series = c+(i-t1*b)noise、mestim 初步估计，再用 p1(2,1)=0、cp1(2,1)=1 等把不显著参数固定为 0，method exact 精确似然重估并保存残差，最后对残差做 miden 检验。

第 8 章习题（PDF p.482–484）

8.1：Merck、J&J、GE、GM、Ford 与价值加权指数 1960–2008 月度对数收益（m-mrk2vw.txt）：样本均值/协方差/相关阵；多元混成检验 $H_0:\rho_1=\cdots=\rho_6=0$；领先–滞后关系。
8.2：1953-04 至 2009-10 月度 1 年与 10 年期国债恒定期限利率（679 个观测）变化序列：建二元 AR 并转成结构形式；建二元 MA 并比较。
8.3：对数序列建 VARMA。
8.4：以利差 $s_t=r_{10,t}-r_{1,t}$ 为门限变量检验门限协整，并建多元门限模型。
8.5：季节模型 $x_t-\Phi_4x_{t-4}=\phi_0+a_t$：均值与协方差、弱平稳充要条件、证明 $\Gamma_\ell=\Phi_4\Gamma_{\ell-4}$。
8.6：$x_t=a_t-\Theta_4a_{t-4}$ 求 $\Gamma_0,\dots,\Gamma_5$。
8.7：1953-04 至 2004-03 的 1 年与 3 年期利率水平：识别 VAR、6 期冲击响应、1–12 步预测；受限常数下 5% 水平协整检验；建 ECM 并预测；比较 VAR 与 ECM 预测。

参考文献（PDF p.484–485）：Balke & Fomby (1997)、Engle & Granger (1987)、Johansen (1988, 1995)、Reinsel & Ahn (1992)、Tsay (1998)、Vidyamurthy (2004)、Pole (2007)、Zivot & Wang (2003) 等。PDF p.486 为空白页。

第 8 章（后半）本章要点

Johansen 检验基于 $\Delta x_t$ 与 $x_{t-1}$ 的偏典型相关，迹统计量与最大特征值统计量渐近分布非标准，需查模拟临界值；确定项设定（无常数、受限常数、非受限常数等）影响临界值。
ECM 预测强制协整约束；水平预测区间随步长扩大（单位根），差分预测较快收敛到均值。
门限协整：协整（误差修正）效应只在偏离超过交易成本的区制中起作用，中间区制近似随机游走；多元门限模型用网格搜索 + AIC 选门限。
配对交易：协整对数价格的价差 $w_t=p_{1t}-\gamma p_{2t}$ 均值回复；做多 1 股/做空 $\gamma$ 股的组合收益即价差增量；以 $\mu_w\pm\Delta$ 为开平仓边界，净利 $2\Delta-\eta$。
附录给出后续章节常用的矩阵工具：谱分解、Cholesky、vec/vech、Kronecker 积，以及多元正态的条件分布公式。

与量化交易的关联

统计套利/配对交易：8.8 节就是配对交易的标准模板——Engle–Granger 两步法（OLS 求对冲比 $\gamma$ + 残差 ADF）或 Johansen 检验 + VECM 求协整向量与 loading；用价差标准差设定开平仓阈值。实盘需补充：样本外检验、滚动重估 $\gamma$、协整关系断裂的止损、交易成本 $\eta$ 的实测、卖空约束。loading $\alpha$ 的大小反映均值回复速度（半衰期 $\approx\ln 0.5/\ln(1+\alpha_1-\gamma\alpha_2)$ 一类的推算可由教材补充）。
期现套利与基差交易：8.7 节的持有成本模型与门限协整直接对应股指期货期现套利、基差回归策略；门限 $\gamma_1,\gamma_2$ 即"无套利区间"，可用于确定开仓阈值；期货领先现货的发现可用于日内价格发现与执行择时。
利率曲线交易：3 月/6 月国库券利率协整示例对应期限利差（曲线）均值回复交易。
风险建模与组合优化：Cholesky 分解用于协方差矩阵模拟（蒙特卡洛 VaR）、正交化冲击；多元正态条件分布用于条件预测、情景分析与 Kalman 滤波。

第 9 章主成分分析与因子模型（Principal Component Analysis and Factor Models）（PDF p.487 起）

本章引言（PDF p.487–488）

多资产组合收益同时、动态地依赖许多经济金融变量，第 8 章的高维模型复杂难用，本章讨论降维方法。**主成分分析（PCA）是最常用的降维方法；观测收益常有相似特征，说明受共同因子（common factors）**驱动。三类因子模型（Connor 1995；Campbell, Lo & MacKinlay 1997）：

宏观经济因子模型（macroeconomic factor models）：用 GDP 增长率、利率、通胀、失业率等可观测宏观变量作因子，线性回归估计。
基本面因子模型（fundamental factor models）：用公司规模、账面/市场价值、行业分类等公司属性构造因子。
统计因子模型（statistical factor models）：因子不可观测（潜变量），从收益序列中估计。

章节结构：9.1 一般因子模型；9.2 宏观因子模型；9.3 基本面因子模型；9.4 PCA；9.5 正交因子模型（含因子旋转和估计）；9.6 渐近主成分分析。参考 Alexander (2001)、Zivot & Wang (2003)。

9.1 因子模型（A Factor Model）（PDF p.488–489）

$k$ 个资产、$T$ 期，一般形式 $r_{it}=\alpha_i+\beta_{i1}f_{1t}+\cdots+\beta_{im}f_{mt}+\epsilon_{it},\quad t=1,\dots,T;\ i=1,\dots,k,\tag{9.1}$ $\alpha_i$ 截距，$f_{jt}$ 为 $m$ 个共同因子，$\beta_{ij}$ 为资产 $i$ 对第 $j$ 个因子的因子载荷（factor loading），$\epsilon_{it}$ 为特质因子（specific factor）。
假设：$f_t$ 为 $m$ 维平稳过程，$E(f_t)=\mu_f$，$\mathrm{Cov}(f_t)=\Sigma_f$；$E(\epsilon_{it})=0$；$\mathrm{Cov}(f_{jt},\epsilon_{is})=0$；$\mathrm{Cov}(\epsilon_{it},\epsilon_{js})=\sigma_i^2$（$i=j,t=s$），否则为 0。即共同因子与特质因子不相关，特质因子彼此不相关；共同因子之间不必不相关。
资产数 $k$ 可能大于期数 $T$（见 9.6 节）。因子分析通常假设因子（从而 $r_t$）无序列相关；若收益有序列相关，先用第 8 章模型去除。
矩阵形式：$r_{it}=\alpha_i+\beta_if_t+\epsilon_{it}$（$\beta_i$ 为行向量）；截面形式 $r_t=\alpha+\beta f_t+\epsilon_t,\quad \mathrm{Cov}(\epsilon_t)=D=\mathrm{diag}\{\sigma_1^2,\dots,\sigma_k^2\},\tag{9.2}$ $\mathrm{Cov}(r_t)=\beta\Sigma_f\beta'+D .$ 因子可观测时 (9.2) 是截面回归形式。
时间序列形式（第 $i$ 个资产）：$R_i=\alpha_i1_T+F\beta_i'+E_i$（9.3），$F$ 为 $T\times m$，第 $t$ 行为 $f_t'$，$\mathrm{Cov}(E_i)=\sigma_i^2I_T$。
合并形式：$r_t=\xi g_t+\epsilon_t$，$g_t=(1,f_t')'$，$\xi=[\alpha,\beta]$（$k\times(m+1)$），堆叠为 $R=G\xi'+E,\tag{9.4}$ $R$ 为 $T\times k$，$G$ 为 $T\times(m+1)$。因子可观测时为**多元线性回归（MLR）**的特例（一般 MLR 的误差协方差不必对角）。

9.2 宏观经济因子模型（Macroeconometric Factor Models）（PDF p.490–496）

因子可观测，对 (9.4) 用最小二乘： $\hat\xi'=\begin{bmatrix}\hat\alpha'\\\hat\beta'\end{bmatrix}=(G'G)^{-1}(G'R),\qquad \hat E=R-G\hat\xi',$ $\hat D=\mathrm{diag}\{\hat\sigma_1^2,\dots,\hat\sigma_k^2\}$，$\hat\sigma_i^2$ 为 $\hat E'\hat E/(T-m-1)$ 的第 $(i,i)$ 元；第 $i$ 个资产 $R_i^2=1-[\hat E'\hat E]_{ii}/[R'R]_{ii}$（$R$ 已去均值时）。
注意：LS 估计未施加特质因子互不相关的约束，一般非有效；施加正交约束计算繁琐，常被忽略。可检查 $\hat E'\hat E/(T-m-1)$ 的非对角元是否接近 0 来检验模型充分性。

9.2.1 单因子模型：市场模型（PDF p.490–494）

市场模型（market model）（Sharpe 1970）： $r_{it}=\alpha_i+\beta_ir_{mt}+\epsilon_{it},\tag{9.5}$ $r_{it}$、$r_{mt}$ 为超额收益，$\beta_i$ 即股票 β。
数据：13 只股票（AA、AGE、CAT、F、FDX、GM、HPQ、KMB、MEL、NYT、PG、TRB、TXN）月度超额收益（百分比），S&P 500 为市场，3 个月国库券为无风险利率，1990-01 至 2003-12，$k=13$，$T=168$。表 9.1 给出样本均值（标准差），如 TXN 2.19(13.8)，SP5 0.42(4.33)。
S-Plus/R 实现：xmtx=cbind(1, 市场)，xit.hat=solve(xmtx,rtn)（最小二乘），E.hat、D.hat=diag(crossprod(E.hat)/(168-2))、r.square。
结果（$\hat\beta$, $\hat\sigma_i$, $R^2$）：AA 1.292/7.694/0.347；AGE 1.514/7.808/0.415；CAT 0.941/7.725/0.219；F 1.219/8.241/0.292；FDX 0.805/8.854/0.135；GM 1.046/8.130/0.238；HPQ 1.628/9.469/0.358；KMB 0.550/6.070/0.134；MEL 1.123/6.120/0.388；NYT 0.771/6.590/0.205；PG 0.469/6.459/0.090；TRB 0.718/7.215/0.157；TXN 1.796/11.474/0.316。金融股（AGE、MEL）与高科技股（HPQ、TXN）β 与 $R^2$ 较高；KMB、PG 较低。$R^2$ 介于 0.09–0.41，市场解释不到一半的个股波动。
模型隐含协方差 $\hat\Sigma_r=\hat\sigma_m^2\hat\beta\hat\beta'+\hat D$，隐含相关阵元素多在 0.1–0.4 且结构单调；样本相关阵则有 Cor(AA,CAT)=0.6、Cor(F,GM)=0.6、Cor(HPQ,TXN)=0.6 等行业内高相关，单因子模型捕捉不到。
**全局最小方差组合（GMVP）**用于比较模型协方差与样本协方差： $\min_\omega\ \omega'\Sigma\omega\quad\text{s.t.}\ \omega'1=1\ \Rightarrow\ \omega=\frac{\Sigma^{-1}1}{1'\Sigma^{-1}1}.$ 模型 GMVP 权重：AA 0.0117, AGE −0.0306, CAT 0.0792, F 0.0225, FDX 0.0802, GM 0.0533, HPQ −0.0354, KMB 0.2503, MEL 0.0703, NYT 0.1539, PG 0.2434, TRB 0.1400, TXN −0.0388。样本 GMVP：AA −0.0073, AGE −0.0085, CAT 0.0866, F −0.0232, FDX 0.0943, GM 0.0916, HPQ 0.0345, KMB 0.2296, MEL 0.0495, NYT 0.1790, PG 0.2651, TRB 0.0168, TXN −0.0080。两者对 TRB 权重差异明显（0.14 vs 0.017），但都重仓 KMB、NYT、PG（低 β 低波动股）。
残差相关阵检验特质因子不相关假设：存在较大值，如 Cor(CAT,AA)=0.45、Cor(GM,F)=0.48，说明单因子不足，遗漏了行业因子。

9.2.2 多因子模型（PDF p.494–496）

Chen, Roll & Ross (1986)：用宏观变量的**意外变化（surprises）**作因子，即去除动态依赖后的残差，可用 VAR 模型残差获得。
示例：CPI（城市消费者全项目，1982–84=100）与 16 岁以上平民就业人数（CE16，千人），均季调，1975-01 至 2003-12（用更长样本求意外序列），取对数差分得增长率（%）。BIC 选 VAR(3)（BIC 值：ar(1) 36.99、ar(2) 38.09、ar(3) 28.23、ar(4) 46.24…），取 1990–2003 的 VAR(3) 残差作两个因子，对 13 只股票超额收益回归（自由度 168−3）。
结果（图 9.2）：所有股票对 CPI 意外增长的 β 均为负（约 −14 到 −2），合理；但 $R^2$ 全都很低（< 0.06），两个宏观变量解释力很弱。模型隐含相关阵接近单位阵，拟合差；残差相关阵与原始收益相关阵几乎相同。
教学要点：宏观因子是可观测的、经济含义明确，但对个股月度收益的解释力往往很低。

9.3 基本面因子模型（Fundamental Factor Models）（PDF p.496 起）

用可观测的资产特定基本面（行业分类、市值、账面价值、成长/价值风格）构造解释超额收益的共同因子。两种方法：

BARRA 方法（Bar Rosenberg 创立，Grinold & Kahn 2000）：把观测到的基本面当作因子 β（载荷），在每个时点用截面回归估计因子实现值 $f_t$；β 不随时间变（本书设定），$f_t$ 随时间变化。
Fama–French 方法（Fama & French 1992）：基于基本面构造对冲组合（hedge portfolio），其收益作为因子实现值 $f_{jt}$。

9.3.1 BARRA 因子模型（PDF p.497–）

假设超额收益已去均值，每个时点 (9.2) 化为 $\tilde r_t=\beta f_t+\epsilon_t,\tag{9.6}$ $\beta$ 已知，是 $k$ 个观测、$m$ 个未知数的多元线性回归（$m<k$ 可估）。因 $\mathrm{Cov}(\epsilon_t)=D$ 非齐方差，应用加权最小二乘（WLS）： $\hat f_t=(\beta'D^{-1}\beta)^{-1}(\beta'D^{-1}\tilde r_t).\tag{9.7}$
$D$ 未知，两步法：
1. 每个 $t$ 做 OLS：$\hat f_{t,o}=(\beta'\beta)^{-1}(\beta'\tilde r_t)$（一致但非有效），残差 $\hat\epsilon_{t,o}=\tilde r_t-\beta\hat f_{t,o}$；因残差协方差不随时间变，合并全部时点：$\hat D_o=\mathrm{diag}\{\frac{1}{T-1}\sum_t\hat\epsilon_{t,o}\hat\epsilon_{t,o}'\}$。
2. 代入得 **GLS（可行 WLS）**估计 $\hat f_{t,g}=(\beta'\hat D_o^{-1}\beta)^{-1}(\beta'\hat D_o^{-1}\tilde r_t),\tag{9.8}$ 残差 $\hat\epsilon_{t,g}=\tilde r_t-\beta\hat f_{t,g}$，$\hat D_g=\mathrm{diag}\{\frac1{T-1}\sum_t\hat\epsilon_{t,g}\hat\epsilon_{t,g}'\}$。
因子协方差 $\hat\Sigma_f=\frac1{T-1}\sum_t(\hat f_{t,g}-\bar f_g)(\hat f_{t,g}-\bar f_g)'$，收益协方差 $\widehat{\mathrm{Cov}}(r_t)=\beta\hat\Sigma_f\beta'+\hat D_g$。

行业因子模型（Industry Factor Model）示例（PDF p.498–）：

10 只股票月度超额收益，1990-01 至 2003-12：金融（AGE 1.36(10.2)、C 2.08(9.60)、MWD 1.87(11.2)、MER 2.08(10.4)）、高科技（DELL 4.82(16.4)、HPQ 1.37(11.8)、IBM 1.06(9.47)）、其他（AA 1.09(9.49)、CAT 1.23(8.71)、PG 1.08(6.75)）。
三个因子分别代表三个行业，β 为行业指示变量： $\tilde r_{it}=\beta_{i1}f_{1t}+\beta_{i2}f_{2t}+\beta_{i3}f_{3t}+\epsilon_{it},\quad \beta_{ij}=\begin{cases}1&\text{资产 }i\text{ 属于行业 }j\\0&\text{否则}\end{cases}\tag{9.9, 9.10}$ 如 IBM 的 $\beta_i=(0,1,0)'$，Alcoa 为 $(0,0,1)'$。
因 β 是哑变量，OLS 估计就是各时点行业内平均超额收益：$\hat f_{t,o}=\big(\frac{AGE_t+C_t+MWD_t+MER_t}{4},\frac{DELL_t+HPQ_t+IBM_t}{3},\frac{AA_t+CAT_t+PG_t}{3}\big)'$；特质因子是个股超额收益与行业均值之差。再据此估计 $D$ 做 GLS。
代码：读 m-barra-9003.txt，去均值，构造 fin/tech/oth 哑变量矩阵 ind.dum，计算样本相关阵。（续）
样本相关阵（PDF p.499–500）：金融股之间 0.6–0.8（MWD–MER 0.8），HPQ–IBM 0.5，AA–CAT 0.6，PG 与其他股相关都很低（0.0–0.3）。
计算：F.hat.o = solve(crossprod(ind.dum)) %*% t(ind.dum) %*% rtn.rm；E.hat.o；diagD.hat.o = rowVars(E.hat.o)；GLS 权矩阵 Hmtx = (β'D^{-1}β)^{-1}β'D^{-1}，F.hat.g = Hmtx %*% rtn.rm。
GLS 权重（$H'$，PDF p.501）：金融因子对 AGE、C、MWD、MER 权重 0.187、0.255、0.259、0.300；科技因子对 DELL、HPQ、IBM 为 0.227、0.402、0.371；其他因子对 AA、CAT、PG 为 0.332、0.432、0.236。即 GLS 对特质波动大的股票（如 DELL）赋较小权重，不再是简单等权。
模型隐含相关阵：$\beta\hat\Sigma_f\beta'+\hat D_g$ 标准化后，行业内相关高于样本值（金融 0.7–0.8；HPQ–IBM 0.7；CAT–PG 样本仅 0.1，模型为 0.6），说明行业哑变量模型把同行业股票强行视为高度相关，可能高估行业内相关。图 9.3 为三个行业因子实现值时序。

因子模拟组合（Factor Mimicking Portfolio）（PDF p.501–502）：

单因子 BARRA 模型下，WLS 估计有组合解释：求解 $\min_\omega\ \tfrac12\omega'D\omega\quad\text{s.t.}\ \omega'\beta=1,$ 解为 $\omega'=(\beta'D^{-1}\beta)^{-1}\beta'D^{-1}$，于是 $\hat f_t=\omega'r_t$ 就是该组合的收益——在对因子暴露为 1 的组合中特质风险最小的那个。若再归一化 $\sum\omega_i=1$，称因子模拟组合；多因子时可逐因子套用。
备注：实践中超额收益样本均值常与 0 无显著差异，拟合 BARRA 模型前可不必去均值。

9.3.2 Fama–French 方法（PDF p.502–503）

对给定基本面（如账面市值比 B/M），两步确定因子实现值：先按该基本面对资产排序，再构造对冲组合——做多排序最高的一组、做空最低的一组（原文写"top quintile (1/3)"，实际 Fama–French 用的是三分位分组），该对冲组合在 $t$ 的收益即该基本面因子实现值。对每个基本面重复；得到 $\{f_t\}$ 后，再对每个资产做时间序列回归估计 β。
三因子：(a) 市场超额收益；(b) SMB（small minus big，小盘相对大盘）；(c) HML（high minus low，价值相对成长）。按市值排序定义规模，按账面权益/市场权益比定义价值（高 B/M）与成长股。
备注：不同因子模型里"因子"概念不同。FF 三个因子是三个金融基本面，若把它们线性组合成一个新属性也可称单因子模型（因为是线性模型），所以谈论因子个数要谨慎；统计因子模型中的因子个数则较明确。
对比（教材可强调）：BARRA 是"已知 β，截面回归求 $f_t$"；FF 是"构造组合得到 $f_t$，时间序列回归求 β"；宏观因子模型是"已知 $f_t$，时间序列回归求 β"。

9.4 主成分分析（Principal Component Analysis）（PDF p.503–508）

协方差（相关）结构是多元时间序列研究重点，在组合选择中尤为重要。PCA 用少数线性组合解释 $k$ 维向量 $r$ 的协方差结构 $\Sigma_r$，可用来研究 $k$ 个资产收益变动的主要来源。

9.4.1 PCA 理论（PDF p.503–505）

可用协方差阵 $\Sigma_r$ 或相关阵 $\rho_r$；相关阵是标准化向量 $r^*=S^{-1}r$ 的协方差阵（$S$ 为标准差对角阵），理论上以协方差阵讨论。
$y_i=w_i'r=\sum_jw_{ij}r_j$；若 $r$ 为简单收益，$y_i$ 即权重 $w_{ij}$ 的组合收益。规范化 $w_i'w_i=1$。 $\mathrm{Var}(y_i)=w_i'\Sigma_rw_i,\qquad \mathrm{Cov}(y_i,y_j)=w_i'\Sigma_rw_j.\tag{9.11, 9.12}$
定义：第 1 主成分 $y_1=w_1'r$ 在 $w_1'w_1=1$ 下最大化 $\mathrm{Var}(y_1)$；第 $i$ 主成分在 $w_i'w_i=1$ 且 $\mathrm{Cov}(y_i,y_j)=0$（$j<i$）下最大化方差。
Result 9.1：设 $\Sigma_r$ 的特征对 $(\lambda_i,e_i)$，$\lambda_1\ge\cdots\ge\lambda_k\ge0$，则第 $i$ 主成分 $y_i=e_i'r$，$\mathrm{Var}(y_i)=\lambda_i$，$\mathrm{Cov}(y_i,y_j)=0$（$i\ne j$）。特征值相等时特征向量不唯一。且
\[\sum_{i=1}^k\mathrm{Var}(r_i)=\mathrm{tr}(\Sigma_r)=\sum_i\lambda_i=\sum_i\mathrm{Var}(y_i).\tag{9.13}\]
第 $i$ 主成分解释总方差比例 $\lambda_i/\sum_j\lambda_j$，前 $i$ 个累计比例 $\sum_{j\le i}\lambda_j/\sum_j\lambda_j$；实践中取小的 $i$ 使累计比例足够大。用相关阵时 $\mathrm{tr}(\rho_r)=k$，比例为 $\lambda_i/k$。
- 推导思路（教材可补）：拉格朗日乘子法，$\max w'\Sigma w-\lambda(w'w-1)$ 得 $\Sigma w=\lambda w$，方差即 $\lambda$，故取最大特征值；正交约束逐次给出后续特征向量。
副产品：零特征值意味着分量间存在精确线性关系。若 $\lambda_k=0$，则 $y_k=\sum_je_{kj}r_j$ 为常数，只有 $k-1$ 个独立随机量，可降维。

9.4.2 经验 PCA（Empirical PCA）（PDF p.505–508）

弱平稳下用样本矩估计： $\hat\Sigma_r=\frac1{T-1}\sum_{t=1}^T(r_t-\bar r)(r_t-\bar r)',\quad \bar r=\frac1T\sum r_t,\tag{9.14}$ $\hat\rho_r=\hat S^{-1}\hat\Sigma_r\hat S^{-1},\tag{9.15}$ 再做对称阵特征分解。R/S-Plus：princomp；FinMetrics：mfactor。
例 9.1：IBM、HPQ、INTC、JPM、BAC 月度对数收益（含股息，%），1990-01 至 2008-12，228 个观测。样本均值 (0.70, 0.99, 1.20, 0.82, 0.41)'。样本协方差阵与相关阵（下三角；方差依次为 74.64、112.22、146.50、106.04、91.83）：
\[\hat\Sigma_r=\begin{bmatrix}74.64\\42.28&112.22\\48.03&70.45&146.50\\30.10&42.42&44.59&106.04\\21.07&26.30&29.24&67.45&91.83\end{bmatrix},\quad \hat\rho_r=\begin{bmatrix}1\\0.46&1\\0.46&0.55&1\\0.34&0.39&0.36&1\\0.25&0.26&0.25&0.68&1\end{bmatrix}.\]
表 9.3：
- 协方差阵 PCA 特征值 284.17、112.93、57.43、46.81、29.87；比例 0.535、0.213、0.108、0.088、0.056；累计 0.535、0.748、0.856、0.944、1.000。第一特征向量 (0.330, 0.483, 0.581, 0.448, 0.347)'，第二 (0.139, 0.279, 0.478, −0.550, −0.610)'。
- 相关阵 PCA 特征值 2.607、1.072、0.569、0.451、0.301；比例 0.522、0.214、0.114、0.090、0.060；累计 0.522、0.736、0.850、0.940、1.000。
- $\hat\lambda_1=2.608$，$e_1=(0.428,0.460,0.451,0.479,0.416)'$——近似等权组合，代表市场成分；$\hat\lambda_2=1.072$，$e_2=(0.341,0.356,0.385,-0.469,-0.623)'$——科技股与金融股之差，代表行业成分。两者共解释约 74% 的总变异。
碎石图（scree plot）：$\hat\lambda_i$ 对 $i$ 的图，找"肘部"（之后特征值较小且大小相近）定主成分个数；图 9.5 中例 9.1 与例 9.3 都取 2 个。除非 $\lambda_j=0$（$j>i$），取前 $i$ 个主成分只是近似。
备注：princomp 输出的"standard deviation"是特征值的平方根。命令：princomp(rtn)、summary、$loadings、screeplot、princomp(rtn, cor=T)。
常见误区：协方差阵 PCA 受高波动资产主导，量纲不同时应使用相关阵；主成分的符号任意（特征向量乘 −1 仍有效）。

9.5 统计因子分析（Statistical Factor Analysis）（PDF p.509–518）

动机："维数灾难（curse of dimensionality）"——参数随阶数与维数急剧增长；多元数据常有共同结构。统计因子分析从数据中识别少数可解释协方差/相关阵大部分变异的因子。
传统因子分析假设无序列相关：周频及更高频金融数据常违背，月度等低频收益基本合理；违背时先用参数模型去除线性动态依赖，再对残差做因子分析。
模型：$r_t$ 弱平稳，均值 $\mu$，协方差 $\Sigma_r$， $r_t-\mu=\beta f_t+\epsilon_t,\tag{9.16}$ $\beta=[\beta_{ij}]_{k\times m}$ 为载荷矩阵，$m<k$。关键特征：因子与载荷都不可观测，因此 (9.16) 不是多元线性回归。
正交因子模型（orthogonal factor model）假设：(1) $E(f_t)=0$，$\mathrm{Cov}(f_t)=I_m$；(2) $E(\epsilon_t)=0$，$\mathrm{Cov}(\epsilon_t)=D$ 对角；(3) $f_t$ 与 $\epsilon_t$ 独立。于是 $\Sigma_r=\beta\beta'+D,\tag{9.17}\qquad \mathrm{Cov}(r_t,f_t)=\beta.\tag{9.18}$ 即 $\mathrm{Var}(r_{it})=\sum_j\beta_{ij}^2+\sigma_i^2$，$\mathrm{Cov}(r_{it},r_{jt})=\sum_l\beta_{il}\beta_{jl}$，$\mathrm{Cov}(r_{it},f_{jt})=\beta_{ij}$。$c_i^2=\sum_j\beta_{ij}^2$ 称共同度（communality），$\sigma_i^2$ 称唯一性/特殊方差（uniqueness / specific variance），$\mathrm{Var}(r_{it})=c_i^2+\sigma_i^2$。
并非所有协方差阵都有正交因子表示；且表示不唯一：对任意 $m\times m$ 正交阵 $P$，$\beta^*=\beta P$，$f_t^*=P'f_t$ 仍满足所有假设（$\mathrm{Cov}(f_t^*)=P'P=I$）。不唯一性既是缺点（载荷含义任意）也是优点（可旋转以获得好解释）；$f_t^*=P'f_t$ 是 $m$ 维空间中的旋转。

9.5.1 估计（PDF p.510–512）

两种方法：

主成分法（Principal Component Method）：不需正态假设、不需预设因子数，协方差阵与相关阵都适用，但结果是近似。设样本协方差阵特征对 $(\hat\lambda_i,\hat e_i)$，取 $m$ 个因子， $\hat\beta=\big[\sqrt{\hat\lambda_1}\hat e_1\,|\,\sqrt{\hat\lambda_2}\hat e_2\,|\cdots|\,\sqrt{\hat\lambda_m}\hat e_m\big].\tag{9.19}$ 特殊方差 $\hat\sigma_i^2=\hat\sigma_{ii,r}-\sum_{j=1}^m\hat\beta_{ij}^2$，共同度 $\hat c_i^2=\sum_j\hat\beta_{ij}^2$。近似误差矩阵 $\hat\Sigma_r-(\hat\beta\hat\beta'+\hat D)$ 的元素平方和 $\le\hat\lambda_{m+1}^2+\cdots+\hat\lambda_k^2$，即误差受被忽略特征值平方和控制。增加 $m$ 时已有载荷不变。
极大似然法（Maximum-Likelihood Method）：若 $f_t,\epsilon_t$ 联合正态，$r_t\sim N(\mu,\beta\beta'+D)$，在唯一性约束 $\beta'D^{-1}\beta=\Delta$（对角阵）下极大化似然，$\mu$ 用样本均值（详见 Johnson & Wichern 2007）。须事先给定 $m$。用修正 LR 检验 $m$ 因子模型的充分性： $LR(m)=-\Big[T-1-\tfrac16(2k+5)-\tfrac23m\Big]\Big[\ln|\hat\Sigma_r|-\ln|\hat\beta\hat\beta'+\hat D|\Big],\tag{9.20}$ 原假设下渐近 $\chi^2$，自由度 $\tfrac12[(k-m)^2-k-m]$。

9.5.2 因子旋转（Factor Rotation）（PDF p.512）

正交变换下 $\beta\beta'+D=\beta^*\beta^{*\prime}+D$，共同度和特殊方差不变，故可寻找 $P$ 使因子易解释；旋转有无穷多种。
Varimax 准则（Kaiser 1958）：令 $\tilde\beta^*_{ij}=\beta^*_{ij}/c_i$（按共同度平方根缩放），选正交阵 $P$ 最大化 $V=\frac1k\sum_{j=1}^m\Big[\sum_{i=1}^k(\tilde\beta^*_{ij})^4-\frac1k\Big(\sum_{i=1}^k\tilde\beta^{*2}_{ij}\Big)^2\Big].$ 含义：每个因子列上载荷平方的"方差"尽量大，使每列载荷分化为大值组与可忽略组，便于解释。旋转是辅助解释工具，有时有用有时无信息；还有很多其他准则（如 quartimax）。

9.5.3 应用（PDF p.512–518）

先用多元混成统计量检验无序列相关假设；若有，建 VARMA 去动态依赖后对残差做因子分析。很多收益序列的线性模型残差相关阵与原数据很接近，此时动态依赖影响可忽略。
例 9.2（五只股票，同例 9.1）：$Q_5(1)=39.99$、$Q_5(5)=160.60$、$Q_5(10)=293.04$，对应 $\chi^2_{25},\chi^2_{125},\chi^2_{250}$ 的 p 值 0.029、0.017、0.032——有轻微序列相关但 1% 水平不显著，忽略。基于相关阵 ML 估计、2 因子（表 9.4）：
- 未旋转载荷 $f_1$：IBM 0.327、HPQ 0.348、INTC 0.337、JPM 0.734、BAC 0.960；$f_2$：0.530、0.669、0.647、0.186、−0.111。
- Varimax 旋转后 $f_1^*$：0.593、0.733、0.709、0.358、0.124；$f_2^*$：0.189、0.177、0.171、0.667、0.958。
- 共同度 $1-\sigma_i^2$：0.387、0.568、0.531、0.573、0.934。方差贡献：未旋转 1.801、1.193（比例 0.360、0.239）；旋转后 1.535、1.459（0.307、0.292）；合计 2.994（0.599）。
- 结论：两因子解释约 60% 变异；旋转后科技股（IBM、HPQ、INTC）重载 $f_1^*$，金融股（JPM、BAC）重载 $f_2^*$，区分了行业；varimax 改变了两因子的次序；IBM 特殊方差较大（共同度仅 0.387），有自身特征值得进一步研究。
例 9.3（美国债券指数月度对数收益，期限 30、20、10、5、1 年，1942-01 至 1999-12，696 个观测；数据见例 8.2）：数据有序列相关，但拟合 VARMA(2,1) 前后同期相关阵几乎不变（如原始 30Y–1Y 相关 0.63，残差 0.66），故直接分析。表 9.5（2 因子，相关阵）：
- 主成分法：$f_1$ 载荷 0.952、0.954、0.956、0.955、0.800；$f_2$ 载荷 0.253、0.240、0.140、−0.142、−0.585；旋转后 $f_1^*$ 0.927、0.922、0.866、0.704、0.325，$f_2^*$ 0.333、0.345、0.429、0.660、0.936；共同度 0.970、0.968、0.934、0.931、0.982；方差比例 0.856、0.101（合计 0.957）。
- ML 法：$f_1$ 0.849、0.857、0.896、1.000、0.813；$f_2$ −0.513、−0.486、−0.303、0.000、0.123；共同度 0.985、0.970、0.895、1.000、0.675；比例 0.784、0.121（合计 0.905）。
- 解释（主成分法）：$f_1$ 五个序列载荷近似相等，代表整体债券收益（水平）；$f_2$ 载荷随期限单调递增且和约为 0，代表期限效应，即长债（≥10 年）与短债之差（斜率）。旋转后 $f_1^*$ 载荷与期限成正比，$f_2^*$ 与期限成反比（长端因子/短端因子）。两种方法前两因子均解释 90% 以上，特殊方差很小。
例 9.4（表 9.2 的 10 只股票月度超额收益，1990–2003，R/S-Plus factanal）：
- 2 因子 ML 被 LR 检验拒绝：$LR(2)=72.96$，$\chi^2_{26}$，p≈0。
- 3 因子：$\chi^2=26.48$，自由度 18，p=0.0892，5% 水平可接受。SS loadings 2.635、1.825、1.326；比例 0.264、0.183、0.133；累计 0.579。Uniqueness：AGE 0.479、C 0.341、MWD 0.201、MER 0.216、DELL 0.690、HPQ 0.346、IBM 0.638、AA 0.417、CAT 0.000、PG 0.885。
- 载荷：Factor1 金融股高（AGE 0.678、C 0.739、MWD 0.817、MER 0.819）；Factor2 科技股和 AA（DELL 0.547、HPQ 0.771、IBM 0.515、AA 0.546）；Factor3 主要是 CAT（0.970）与 AA（0.497），代表其余工业部门。图 9.6 为载荷图。
- rotate(stat.fac, rotation='quartimax') 做 quartimax 旋转，结果类似（MWD 0.844、HPQ 0.753、CAT 0.931 等）；predict(stat.fac, type='weighted.ls') 得因子实现值。
- 模型拟合相关阵（fitted）比 9.3.1 的行业因子模型更接近样本相关阵（如 CAT–PG 0.1，与样本一致）；也可用 GMVP 比较。
- 注意：CAT 的 uniqueness 为 0 属于 Heywood 情形（ML 估计中特殊方差估到边界），教材可提醒。

9.6 渐近主成分分析（Asymptotic Principal Component Analysis）（PDF p.518–）

前述 PCA 假设 $k<T$。针对 $T$ 小、$k$ 大，Connor & Korajczyk (1986, 1988) 提出 APCA：依赖 $k\to\infty$ 的渐近结果，对 $T\times T$ 矩阵 $\hat\Omega_T=\frac1k(R-1_T\bar r')(R-1_T\bar r')'$ 做特征分解（$R$ 为 $T\times k$ 收益矩阵，$\bar r_i=1_T'R_i/T$）。$k\to\infty$ 时这等价于传统统计因子分析：因子 $f_t$ 的 APCA 估计就是 $\hat\Omega_T$ 前 $m$ 个特征向量，令 $F_t$ 为由它们组成的 $m\times T$ 矩阵，$f_t$ 是其第 $t$ 列。
仿 BARRA 的精化步骤（Connor & Korajczyk 1988）：
1. 用 $\hat\Omega_T$ 得 $f_t$ 初始估计；
2. 对每个资产做时间序列 OLS $r_{it}=\alpha_i+\beta_if_t+\epsilon_{it}$，得残差方差 $\hat\sigma_i^2$；
3. $\hat D=\mathrm{diag}\{\hat\sigma_1^2,\dots,\hat\sigma_k^2\}$，缩放收益 $R^*=R\hat D^{-1/2}$；
4. 计算 $\hat\Omega^*=\frac1k(R^*-1_T\bar r_*')(R^*-1_T\bar r_*')'$ 并特征分解，得精化的 $f_t$。
直观：$T\times T$ 矩阵与 $k\times k$ 矩阵 $R'R$ 有相同非零特征值（SVD 对偶），当 $k\gg T$ 时分解 $T\times T$ 矩阵计算量小得多。

9.6.1 因子个数的选择（PDF p.519）

Connor & Korajczyk (1993)：若 $m$ 是正确因子数，从 $m$ 增到 $m+1$ 时特质误差的截面方差不应显著下降。
Bai & Ng (2002) 信息准则：$\hat\Omega_T$ 的特征分解求解最小二乘问题 $\min_{\alpha,\beta,f_t}\frac1{kT}\sum_{i=1}^k\sum_{t=1}^T(r_{it}-\alpha_i-\beta_if_t)^2 .$ 用 APCA 得到的 $m$ 维 $f_t$ 对每个资产回归得残差方差 $\hat\sigma_i^2(m)$，截面平均 $\hat\sigma^2(m)=\frac1k\sum_i\hat\sigma_i^2(m)$， $C_{p1}(m)=\hat\sigma^2(m)+m\hat\sigma^2(M)\Big(\frac{k+T}{kT}\Big)\ln\Big(\frac{kT}{k+T}\Big),$ $C_{p2}(m)=\hat\sigma^2(m)+m\hat\sigma^2(M)\Big(\frac{k+T}{kT}\Big)\ln(P_{kT}^2),\quad P_{kT}=\min(\sqrt k,\sqrt T),$ $M$ 为预设最大因子数，在 $0\le m\le M$ 中取最小者。两准则可能选出不同个数。（续）

9.6.2 实例（PDF p.520–523）

40 只股票月度简单收益，2001-01 至 2003-12，$k=40$，$T=36$（$k>T$）。股票为 2004 年 9 月某日 NASDAQ（INTC、MSFT、SUNW、CSCO、AMAT、ORCL、SIRI、COCO、CORV、SUPG、YHOO、JDSU、QCOM、CIEN、DELL、ERTS、EBAY、ADCT、AAPL、JNPR）与 NYSE（LU、PFE、NT、BAC、BSX、GE、TXN、XOM、FRX、Q、F、TWX、C、MOT、JPM、TYC、HPQ、NOK、WMT、AMD）成交最活跃者（表 9.6）。S-Plus 命令 mfactor。
选因子数：Connor–Korajczyk 法（k='ck', max.k=10, sig=0.05）选 $m=1$，此时载荷中位数 0.629，回归 $R^2$ 中位数 0.487（0.090–0.831）；Bai–Ng 法（k='bn'）选 $m=6$（$C_{p1}$、$C_{p2}$ 结果不同，取较小者）。
用 $m=6$ 做 APCA：F.1 载荷中位数 0.561（0.048–2.222），其余因子载荷正负都有；回归 $R^2$ 中位数 0.695，均值 0.651（0.219–0.999）。碎石图（图 9.7）累计方差比例：0.446、0.648、0.779、0.829、0.865、0.894、0.909，6 个因子解释约 89.4%。图 9.8 为 6 个因子收益时序。命令 screeplot.mfactor、fplot(factors(apca))。
教学提示：两种选因子数方法结论差异很大（1 vs 6），$T=36$ 时尤其不稳定，实务中需结合经济解释与样本外表现。

第 9 章习题（PDF p.521–523）

9.1：13 只股票（AA、AXP、CAT、DE、F、FDX、HPQ、IBM、JNJ、KMB、MMM、PG、WFC）与 S&P 500 的 1990–2008 月度超额收益，做市场模型，求 $\beta_i,\sigma_i^2,R^2$。
9.2：Merck、J&J、GE、GM、Ford、价值加权指数 1960–2008：协方差阵/相关阵 PCA；统计因子分析，确定因子数，用主成分法和 ML 法估计载荷。
9.3：10 只股票（制药 ABT、LLY、MRK、PFE；汽车 F、GM；石油 BP、CVX、RD、XOM）与 SP5 的 1990–2003 月度超额收益：单因子市场模型、β 和 $R^2$ 图，用 GMVP 比较模型与样本协方差。
9.4：同数据做 BARRA 行业因子模型，画三因子实现值并评价。
9.5：PCA 与碎石图，确定并解释共同因子个数。
9.6：统计因子分析，5% 水平下因子个数、载荷图与解释。
9.7：用联邦基金利率与工业生产指数（1954-07 至 2003-12，VAR 求意外序列）做宏观因子模型。

参考文献（PDF p.523–524）：Alexander (2001)、Bai & Ng (2002)、Campbell, Lo & MacKinlay (1997)、Chen, Roll & Ross (1986)、Connor (1995)、Connor & Korajczyk (1986, 1988, 1993)、Fama & French (1992)、Grinold & Kahn (2000)、Johnson & Wichern (2007)、Kaiser (1958)、Sharpe (1970)、Zivot & Wang (2003)。

第 9 章本章要点

一般因子模型 $r_t=\alpha+\beta f_t+\epsilon_t$，协方差结构 $\Sigma_r=\beta\Sigma_f\beta'+D$，把 $k(k+1)/2$ 个参数压缩为 $km+m(m+1)/2+k$ 个。
三类因子模型的"已知/未知"不同：宏观（$f_t$ 已知，时序回归求 β）、BARRA（β 已知，截面 WLS/GLS 求 $f_t$）、Fama–French（由排序对冲组合得 $f_t$，再时序回归求 β）、统计因子（两者都未知，PCA/ML 估计）。
单因子市场模型 $R^2$ 仅 0.09–0.41，残差存在行业相关；宏观意外因子解释力很弱；行业因子模型可能高估行业内相关；3 因子统计模型最贴近样本相关阵。
GMVP 权重 $\Sigma^{-1}1/(1'\Sigma^{-1}1)$ 是比较协方差估计的实用工具；单因子 BARRA 的 WLS 因子估计等价于"单位暴露、最小特质风险"的因子模拟组合。
PCA：主成分为协方差（相关）阵特征向量，方差为特征值；用累计方差比例和碎石图定个数；股票第一主成分≈市场，债券前两主成分≈水平与斜率。
正交因子模型的旋转不确定性；主成分法与 ML 法估计；LR 检验因子数；varimax 旋转。
$k>T$ 时用 APCA（分解 $T\times T$ 矩阵）；因子数用 Connor–Korajczyk 检验或 Bai–Ng 信息准则。

与量化交易的关联

风险模型（最直接）：BARRA 方法即商业多因子风险模型（Barra USE/CNE、Axioma）的原型——行业哑变量 + 风格暴露作为 β，每期截面 WLS 回归得到因子收益，再估 $\Sigma_f$ 与特质方差 $D$，组合风险 $=w'(\beta\Sigma_f\beta'+D)w$。A 股多因子风险模型（行业 + 市值、动量、波动等风格）完全沿用此框架；实务中 WLS 权重常取市值平方根而非 $1/\sigma_i^2$。
因子研究：Fama–French 排序对冲组合法是构造因子收益（SMB、HML 及各类异象因子）的标准做法；截面回归法得到的因子收益与因子模拟组合是"纯因子收益"的计算方式，可用于因子绩效归因与 IC 之外的因子检验。
组合优化：因子模型协方差比样本协方差更稳定，尤其 $k$ 接近或大于 $T$ 时样本阵奇异/病态；GMVP 比较展示了协方差估计对最优权重的敏感性。
统计套利与残差收益：PCA/APCA 提取的统计因子可用于构建市场中性组合，残差（特质收益）做均值回复交易（Avellaneda–Lee 式 PCA 统计套利）；Bai–Ng 准则用于确定剔除的因子个数。
固定收益：债券收益的水平/斜率主成分即利率曲线 PCA，用于久期/曲线风险对冲（DV01 按主成分分解）。
宏观因子：意外序列（VAR 残差）作为因子的方法可用于宏观择时、资产配置中的宏观敏感度分析，但需注意解释力通常很低。

第 10 章多元波动率模型及其应用（Multivariate Volatility Models and Their Applications）（PDF p.525 起）

本章引言（PDF p.525–526）

将第 3 章单变量波动率模型推广到多元，多元波动率指多资产收益的条件协方差矩阵。用途：组合选择、资产配置、多资产头寸的 VaR。
记 $r_t=\mu_t+a_t$，$\mu_t=E(r_t|F_{t-1})$，$a_t$ 为新息。均值方程通常用带外生变量的 VARMA： $\mu_t=\Upsilon x_t+\sum_{i=1}^p\Phi_ir_{t-i}-\sum_{i=1}^q\Theta_ia_{t-i},\tag{10.1}$ $x_t$ 为 $m$ 维外生变量（$x_{1t}=1$），$\Upsilon$ 为 $k\times m$。
$\Sigma_t=\mathrm{Cov}(a_t|F_{t-1})$ 为 $k\times k$ 正定阵，多元波动率建模即刻画 $\Sigma_t$ 的时间演化。维数灾难：$\Sigma_t$ 有 $k(k+1)/2$ 个量（5 维即 15 个）。本章介绍相对简单、实用的模型，特别是允许时变相关系数的模型（可估计时变市场 β）。
结构：10.1 指数加权估计（基准）；10.2 多元 GARCH 推广；10.3 两种 $\Sigma_t$ 重参数化（Cholesky 尤其有用）；10.4 二元 GARCH；10.5 高维；10.6 降维；10.7 应用；10.8 多元 Student-t 分布。

10.1 指数加权估计（Exponentially Weighted Estimate）（PDF p.526–530）

等权估计 $\hat\Sigma=\frac1{t-1}\sum_{j=1}^{t-1}a_ja_j'$（$a_j$ 均值为零）。为体现时变且强调近期，用指数平滑： $\hat\Sigma_t=\frac{1-\lambda}{1-\lambda^{t-1}}\sum_{j=1}^{t-1}\lambda^{j-1}a_{t-j}a_{t-j}',\quad 0<\lambda<1,\tag{10.2}$ 权重和为 1。$t$ 足够大（$\lambda^{t-1}\approx0$）时递推为 $\hat\Sigma_t=(1-\lambda)a_{t-1}a_{t-1}'+\lambda\hat\Sigma_{t-1},$ 称 EWMA（指数加权移动平均）协方差估计。
给定 $\lambda$ 与初值 $\hat\Sigma_1$ 可递推。若 $a_t=r_t-\mu_t(\Theta)$ 服从 $N(0,\Sigma_t)$，可用 ML 联合估计 $\lambda$ 与 $\Theta$： $\ln L(\Theta,\lambda)\propto-\frac12\sum_{t=1}^T\ln|\Sigma_t|-\frac12\sum_{t=1}^T(r_t-\mu_t)'\Sigma_t^{-1}(r_t-\mu_t)$ （原文漏写 ln），用 $\hat\Sigma_t$ 代替 $\Sigma_t$ 递推计算。
例 10.1：恒生指数与日经 225 日对数收益（%），2006-01-04 至 2008-12-30，仅取两市同时开市日，713 个观测（Yahoo Finance）。图 10.1 显示金融危机影响。单变量 GARCH(1,1)：
\[r_{1t}=0.109+a_{1t},\quad \sigma_{1t}^2=0.038+0.143a_{1,t-1}^2+0.855\sigma_{1,t-1}^2,\tag{10.3}\]

\[r_{2t}=0.003+a_{2t},\quad \sigma_{2t}^2=0.044+0.127a_{2,t-1}^2+0.861\sigma_{2,t-1}^2.\tag{10.4}\]
除日经均值常数外均 5% 显著；标准化残差及其平方的 Ljung–Box 无不足；两波动率方程接近 IGARCH(1,1)（次贷危机推高波动）；图 10.2 显示 2008 年波动显著升高。
- 二元 EWMA（FinMetrics mgarch(formula.mean=~arma(0,0), formula.var=~ewma1)）：C(1)=0.0824（t=2.67），C(2)=−0.0068（不显著），ALPHA=0.0695（标准误 0.0049），故 $\hat\lambda=1-0.0695\approx0.9305$，处于实践常见范围（RiskMetrics 日频取 0.94）。图 10.3：EWMA 波动率比单变量 GARCH 更平滑，但形态相似。

10.2 若干多元 GARCH 模型（Some Multivariate GARCH Models）（PDF p.530–536）

综述见 Bauwens, Laurent & Rombouts (2004)。

10.2.1 对角向量化（DVEC）模型（PDF p.530–533）

Bollerslev, Engle & Wooldridge (1988) 推广 EWMA： $\Sigma_t=A_0+\sum_{i=1}^mA_i\odot(a_{t-i}a_{t-i}')+\sum_{j=1}^sB_j\odot\Sigma_{t-j},\tag{10.5}$ $A_i,B_j$ 为对称阵，$\odot$ 为 Hadamard 积（逐元素乘）。称 DVEC($m,s$)。二元 DVEC(1,1)（只写下三角）： $\sigma_{11,t}=A_{11,0}+A_{11,1}a_{1,t-1}^2+B_{11,1}\sigma_{11,t-1},$ $\sigma_{21,t}=A_{21,0}+A_{21,1}a_{1,t-1}a_{2,t-1}+B_{21,1}\sigma_{21,t-1},$ $\sigma_{22,t}=A_{22,0}+A_{22,1}a_{2,t-1}^2+B_{22,1}\sigma_{22,t-1}.$ 每个元素只依赖自身滞后值和对应乘积项，即逐元素 GARCH(1,1)。优点简单；缺点：不保证 $\Sigma_t$ 正定，且不允许波动率之间的动态依赖（溢出）。EWMA 是其特例（$A_i=(1-\lambda)J$，$B=\lambda J$，$A_0=0$）。
例 10.2：辉瑞（PFE）与默克（MRK）月度简单收益（含股息），1965-01 至 2008-12，528 个观测。$Q(10)$=10.48(0.40) 与 11.42(0.33)，无显著序列相关，均值方程仅含常数。mgarch(rtn~1, ~dvec(1,1)) 估计：C(1)=0.0135、C(2)=0.0131；除 A(1,1)（p=0.056）外均 5% 显著。拟合模型： $\sigma_{11,t}=0.00075+0.071a_{1,t-1}^2+0.786\sigma_{11,t-1},$ $\sigma_{21,t}=0.00008+0.025a_{1,t-1}a_{2,t-1}+0.950\sigma_{21,t-1},$ $\sigma_{22,t}=0.00008+0.041a_{2,t-1}^2+0.945\sigma_{22,t-1}.$ 个体检验：PFE 标准化残差 $Q(12)=9.53(0.66)$，平方残差 $Q(12)=22.08(0.037)$（原文正文写 12.35(0.42)，与输出不符——输出中 12.349 是 MRK 标准化残差的统计量）；MRK 平方残差 6.44(0.89)。更有信息的检验是对二元标准化残差及其平方用多元 Q 统计量（Li 2004）。输出对象含 sigma.t（波动率）、R.t（相关）。图 10.5：时变相关在 0.37–0.83 之间。

10.2.2 BEKK 模型（PDF p.533–536）

Engle & Kroner (1995) 的 Baba–Engle–Kraft–Kroner 模型保证正定： $\Sigma_t=AA'+\sum_{i=1}^mA_i(a_{t-i}a_{t-i}')A_i'+\sum_{j=1}^sB_j\Sigma_{t-j}B_j',\tag{10.6}$ $A$ 下三角，$A_i,B_j$ 为 $k\times k$。只要 $AA'$ 正定，$\Sigma_t$ 几乎必然正定；允许波动率间动态依赖。缺点：(1) $A_i,B_j$ 参数对滞后波动/冲击无直接解释；(2) 参数个数 $k^2(m+s)+k(k+1)/2$ 随 $m,s$ 迅速增加；经验上许多估计不显著。
例 10.3：同 PFE/MRK 数据，mgarch(rtn~1, ~bekk(1,1))。估计：C=(0.0133, 0.0127)；$A$：A(1,1)=0.025、A(2,1)=0.013、A(2,2)=3×10⁻⁶（不显著）；ARCH 矩阵 $\begin{bmatrix}0.213&0.063\\0.100&0.182\end{bmatrix}$（非对角两元 p=0.168、0.405 不显著）；GARCH 矩阵 $\begin{bmatrix}0.909&-0.008\\-0.059&0.982\end{bmatrix}$（非对角不显著）。拟合方程 $\Sigma_t=\begin{bmatrix}0.025&0\\0.013&3\times10^{-6}\end{bmatrix}\begin{bmatrix}0.025&0.013\\0&3\times10^{-6}\end{bmatrix}+\begin{bmatrix}0.213&0.063\\0.100&0.182\end{bmatrix}a_{t-1}a_{t-1}'\begin{bmatrix}0.213&0.100\\0.063&0.182\end{bmatrix}+\begin{bmatrix}0.901&-0.008\\-0.059&0.982\end{bmatrix}\Sigma_{t-1}\begin{bmatrix}0.901&-0.059\\-0.008&0.982\end{bmatrix}.$ Ljung–Box：PFE 标准化残差 9.47(0.66)，平方 21.55(0.043)；MRK 11.59(0.48)，平方 9.19(0.69)，个体检验未显示不足。图 10.6：与 DVEC 相比，BEKK 的时变相关波动更大。BEKK 常含不显著参数，且需做矩阵乘法才能解读。

10.3 重参数化（Reparameterization）（PDF p.536–）

利用 $\Sigma_t$ 的对称性进行重参数化，介绍两种。

10.3.1 使用相关系数（Use of Correlations）（PDF p.536–537）

$\Sigma_t\equiv[\sigma_{ij,t}]=D_t\rho_tD_t,\tag{10.7}$

$\rho_t$ 为条件相关阵，$D_t=\mathrm{diag}\{\sqrt{\sigma_{11,t}},\dots,\sqrt{\sigma_{kk,t}}\}$。$\Sigma_t$ 的演化由条件方差 $\sigma_{ii,t}$ 和 $\rho_{ij,t}$（$j<i$）决定。定义 $k(k+1)/2$ 维向量

\Xi_t=(\sigma_{11,t},\dots,\sigma_{kk,t},\varrho_t')',\tag{10.8}

$\varrho_t=(\rho_{21,t},\dots,\rho_{k1,t}|\rho_{32,t},\dots,\rho_{k2,t}|\cdots|\rho_{k,k-1,t})'$（按列堆叠下三角）。$k=2$：$\Xi_t=(\sigma_{11,t},\sigma_{22,t},\rho_{21,t})'$（10.9）；$k=3$：$\Xi_t=(\sigma_{11,t},\sigma_{22,t},\sigma_{33,t},\rho_{21,t},\rho_{31,t},\rho_{32,t})'$（10.10）。

二元正态条件密度 $f(a_{1t},a_{2t}|\Xi_t)=\frac{1}{2\pi\sqrt{\sigma_{11,t}\sigma_{22,t}(1-\rho_{21,t}^2)}}\exp\big[-\frac{Q}{2(1-\rho_{21,t}^2)}\big]$，$Q=\frac{a_{1t}^2}{\sigma_{11,t}}+\frac{a_{2t}^2}{\sigma_{22,t}}-\frac{2\rho_{21,t}a_{1t}a_{2t}}{\sqrt{\sigma_{11,t}\sigma_{22,t}}}$；对数密度 $\ell=-\frac12\Big\{\ln[\sigma_{11,t}\sigma_{22,t}(1-\rho_{21,t}^2)]+\frac{1}{1-\rho_{21,t}^2}\Big(\frac{a_{1t}^2}{\sigma_{11,t}}+\frac{a_{2t}^2}{\sigma_{22,t}}-\frac{2\rho_{21,t}a_{1t}a_{2t}}{\sqrt{\sigma_{11,t}\sigma_{22,t}}}\Big)\Big\}.\tag{10.11}$
优点：直接建模协方差和相关；缺点：$k\ge3$ 时似然复杂；需约束极大化保证正定，$k$ 大时约束复杂。

10.3.2 Cholesky 分解（PDF p.537–）

优点：估计时无需正定约束（Pourahmadi 1999）；是正交变换，似然极其简单。
$\Sigma_t$ 正定，存在单位下三角 $L_t$ 与正对角阵 $G_t$： $\Sigma_t=L_tG_tL_t'.\tag{10.12}$
二元情形：$L_t=\begin{bmatrix}1&0\\q_{21,t}&1\end{bmatrix}$，$G_t=\mathrm{diag}(g_{11,t},g_{22,t})$，展开得 $\sigma_{11,t}=g_{11,t},\ \sigma_{21,t}=q_{21,t}g_{11,t},\ \sigma_{22,t}=g_{22,t}+q_{21,t}^2g_{11,t},\tag{10.13}$ $g_{11,t}=\sigma_{11,t},\ q_{21,t}=\frac{\sigma_{21,t}}{\sigma_{11,t}},\ g_{22,t}=\sigma_{22,t}-\frac{\sigma_{21,t}^2}{\sigma_{11,t}}.\tag{10.14}$ 对照回归 $a_{2t}=\beta a_{1t}+b_{2t}$（10.15）：$\beta=\sigma_{21,t}/\sigma_{11,t}$，$\mathrm{Var}(b_{2t})=\sigma_{22,t}-\sigma_{21,t}^2/\sigma_{11,t}$，$b_{2t}\perp a_{1t}$。故 $g_{11,t}$ 是 $a_{1t}$ 的方差，$g_{22,t}$ 是回归残差方差，$q_{21,t}$ 是回归系数 β。Cholesky 分解等价于正交变换 $b_{1t}=a_{1t}$，$b_{2t}=a_{2t}-q_{21,t}a_{1t}$，$\mathrm{Cov}(b_t)$ 对角。
三维情形：$\Sigma_t$ 元素与 $(g,q)$ 关系 $\sigma_{11}=g_{11}$，$\sigma_{21}=q_{21}g_{11}$，$\sigma_{22}=q_{21}^2g_{11}+g_{22}$，$\sigma_{31}=q_{31}g_{11}$，$\sigma_{32}=q_{31}q_{21}g_{11}+q_{32}g_{22}$，$\sigma_{33}=q_{31}^2g_{11}+q_{32}^2g_{22}+g_{33}$；反解 $q_{31}=\sigma_{31}/\sigma_{11}$，$q_{32}=\frac{1}{g_{22}}\big(\sigma_{32}-\frac{\sigma_{31}\sigma_{21}}{\sigma_{11}}\big)$，$g_{33}=\sigma_{33}-q_{31}^2g_{11}-q_{32}^2g_{22}$（均带下标 $t$）。它们正是逐次正交化回归 $b_{1t}=a_{1t}$，$a_{2t}=\beta_{21}b_{1t}+b_{2t}$，$a_{3t}=\beta_{31}b_{1t}+\beta_{32}b_{2t}+b_{3t}$ 的系数与残差方差：$q_{ij,t}=\beta_{ij}$，$g_{ii,t}=\mathrm{Var}(b_{it})$，$b_{it}\perp b_{jt}$。
一般情形：$b_{1t}=a_{1t}$，对 $1<i\le k$ 递归回归 $a_{it}=q_{i1,t}b_{1t}+\cdots+q_{i(i-1),t}b_{(i-1)t}+b_{it},\tag{10.16}$ $b_t=L_t^{-1}a_t,\quad a_t=L_tb_t,\tag{10.17}$ $\mathrm{Cov}(b_t)=L_t^{-1}\Sigma_t(L_t^{-1})'=G_t$。参数向量 $\Xi_t=(g_{11,t},\dots,g_{kk,t},q_{21,t},q_{31,t},q_{32,t},\dots,q_{k1,t},\dots,q_{k(k-1),t})',\tag{10.18}$ 仍为 $k(k+1)/2$ 维。
似然简化：$|L_t|=1$， $|\Sigma_t|=|G_t|=\prod_{i=1}^kg_{ii,t};\tag{10.19}$ 若 $a_t|F_{t-1}\sim N(0,\Sigma_t)$，则 $b_t\sim N(0,G_t)$， $\ell(a_t,\Sigma_t)=\ell(b_t,\Sigma_t)=-\frac12\sum_{i=1}^k\Big[\ln(g_{ii,t})+\frac{b_{it}^2}{g_{ii,t}}\Big].\tag{10.20}$
优点：(1) 只要 $g_{ii,t}>0$ 就正定，建模 $\ln g_{ii,t}$ 即可自动满足；(2) 参数有清晰的回归解释（正交化冲击的系数与残差方差）；(3) 相关系数 $\rho_{21,t}=q_{21,t}\sqrt{\sigma_{11,t}}/\sqrt{\sigma_{22,t}}$，即使 $q_{21,t}=c$ 为常数，只要方差比时变，相关仍时变——这是与相关系数参数化的主要区别。
由 (10.16) 与正交性，$\sigma_{ii,t}=\mathrm{Var}(a_{it}|F_{t-1})=\sum_{v=1}^iq_{iv,t}^2g_{vv,t}$（$q_{ii,t}=1$）。（续）协方差 $\sigma_{ij,t}=\sum_{v=1}^jq_{iv,t}q_{jv,t}g_{vv,t}$（$j<i$），$q_{vv,t}=1$。这就是 Cholesky 参数化下的 $\Sigma_t$。

10.4 二元收益的 GARCH 模型（GARCH Models for Bivariate Returns）（PDF p.541–551）

以 GARCH 为例（其他单变量波动率模型可同法推广）。多元 GARCH 用"确定性方程（exact equations）"（不含随机冲击）描述 $k(k+1)/2$ 维 $\Xi_t$ 的演化；即使 $k=2$（$\Xi_t$ 三维）也可能复杂，故常加限制。

10.4.1 常相关模型（Constant-Correlation Models）（PDF p.541–545）

Bollerslev (1990)：设 $\rho_{21,t}=\rho_{21}$ 不变（$|\rho_{21}|<1$），只需对 $\Xi_t^*=(\sigma_{11,t},\sigma_{22,t})'$ 建两条方程。GARCH(1,1)： $\Xi_t^*=\alpha_0+\alpha_1a_{t-1}^2+\beta_1\Xi_{t-1}^*,\tag{10.21}$ $\begin{bmatrix}\sigma_{11,t}\\\sigma_{22,t}\end{bmatrix}=\begin{bmatrix}\alpha_{10}\\\alpha_{20}\end{bmatrix}+\begin{bmatrix}\alpha_{11}&\alpha_{12}\\\alpha_{21}&\alpha_{22}\end{bmatrix}\begin{bmatrix}a_{1,t-1}^2\\a_{2,t-1}^2\end{bmatrix}+\begin{bmatrix}\beta_{11}&\beta_{12}\\\beta_{21}&\beta_{22}\end{bmatrix}\begin{bmatrix}\sigma_{11,t-1}\\\sigma_{22,t-1}\end{bmatrix},\tag{10.22}$ $a_{t-1}^2=(a_{1,t-1}^2,a_{2,t-1}^2)'$，$\alpha_{i0}>0$，$\alpha_1,\beta_1$ 非负。
令 $\eta_t=a_t^2-\Xi_t^*$，得 $a_t^2=\alpha_0+(\alpha_1+\beta_1)a_{t-1}^2+\eta_t-\beta_1\eta_{t-1}$，即 $a_t^2$ 的二元 ARMA(1,1)（单变量 GARCH(1,1) 的直接推广）。性质：
1. 若 $\alpha_1+\beta_1$ 的特征值都在 (0,1)，则 $a_t^2$ 弱平稳，无条件协方差正定；无条件方差 $(\sigma_1^2,\sigma_2^2)'=(I-\alpha_1-\beta_1)^{-1}\alpha_0$（原文写 $\phi_0$），协方差 $\rho_{21}\sigma_1\sigma_2$。
2. $\alpha_{12}=\beta_{12}=0$ 时 $a_{1t}$ 的波动率不依赖 $a_{2t}$ 的过去波动；$\alpha_{21}=\beta_{21}=0$ 同理。
3. $\alpha_1,\beta_1$ 都对角时退化为两个单变量 GARCH(1,1)，波动率无动态关联。
4. 预测：$\Xi_h^*(1)=\alpha_0+\alpha_1a_h^2+\beta_1\Xi_h^*$，$\Xi_h^*(\ell)=\alpha_0+(\alpha_1+\beta_1)\Xi_h^*(\ell-1)$（$\ell>1$）；协方差预测 $\hat\rho_{21}[\sigma_{11,h}(\ell)\sigma_{22,h}(\ell)]^{0.5}$。
例 10.4（恒生/日经日收益，同例 10.1）：均值 $r_{1t}=0.101+a_{1t}$，$r_{2t}=0.002+a_{2t}$（标准误 0.050、0.048）；波动率（对角） $\sigma_{11,t}=0.079+0.145a_{1,t-1}^2+0.833\sigma_{11,t-1},\quad \sigma_{22,t}=0.054+0.105a_{2,t-1}^2+0.875\sigma_{22,t-1},\tag{10.23}$ 常相关 0.668。标准化残差 $\tilde a_{it}=a_{it}/\sqrt{\sigma_{ii,t}}$ 的多元 Ljung–Box：$Q_2(4)=17.29(0.37)$，$Q_2(12)=48.21(0.46)$（自由度 16、48），拟合合理。称二元对角常相关模型：两市场波动率无动态关联但同期相关；实际中可能存在波动溢出（volatility spillover），此类模型未必合适。S-Plus：mgarch(rtn~1, ~ccc(1,1))。图 10.7。
例 10.5（IBM 与 S&P 500 月度对数收益 %，1926-01 至 1999-12）：常相关 GARCH(1,1)。均值 $r_{1t}=1.351+0.072r_{1,t-1}+0.055r_{1,t-2}-0.119r_{2,t-2}+a_{1t}$，$r_{2t}=0.703+a_{2t}$。波动率 $\begin{bmatrix}\sigma_{11,t}\\\sigma_{22,t}\end{bmatrix}=\begin{bmatrix}2.98\\2.09\end{bmatrix}+\begin{bmatrix}0.079&0\\0.042&0.045\end{bmatrix}\begin{bmatrix}a_{1,t-1}^2\\a_{2,t-1}^2\end{bmatrix}+\begin{bmatrix}0.873&-0.031\\-0.066&0.913\end{bmatrix}\begin{bmatrix}\sigma_{11,t-1}\\\sigma_{22,t-1}\end{bmatrix},\tag{10.24}$ 常相关 0.614（标准误 0.020）。$Q_2(4)=16.77(0.21)$、$Q_2(8)=32.40(0.30)$（自由度因均值方程含 3 个滞后项调整为 13、29），平方残差 $Q_2^*(4)=18.00(0.16)$、$Q_2^*(8)=39.09(0.10)$，5% 水平无问题。该模型显示两收益波动率之间存在反馈关系。

10.4.2 时变相关模型（Time-Varying Correlation Models）（PDF p.545–551）

常相关的缺陷：实际相关会变。IBM 与 S&P 500 用 120 个月滚动窗口的样本相关（图 10.8）随时间变化且近年下降（IBM 市值排名变化）。Tse (2000) 提出检验常相关的 LM 统计量。
方法一：直接建模相关系数。因相关为正，用 logistic 变换： $\rho_{21,t}=\frac{\exp(q_t)}{1+\exp(q_t)},\quad q_t=\varpi_0+\varpi_1\rho_{21,t-1}+\varpi_2\frac{a_{1,t-1}a_{2,t-1}}{\sqrt{\sigma_{11,t-1}\sigma_{22,t-1}}},\tag{10.25}$ 称相关系数的 GARCH(1,1) 模型；$\varpi_1=\varpi_2=0$ 时退化为常相关。相关为负可加负号；符号未知时用 Fisher 变换 $q_t=\ln\frac{1+\rho_{21,t}}{1-\rho_{21,t}}$，$\rho_{21,t}=\frac{\exp(q_t)-1}{\exp(q_t)+1}$，对 $q_t$ 建 GARCH 型方程。
例 10.5 续（方法一）：联合估计得均值 $r_{1t}=1.318+0.076r_{1,t-1}-0.068r_{2,t-2}+a_{1t}$，$r_{2t}=0.673+a_{2t}$；
\[\begin{bmatrix}\sigma_{11,t}\\\sigma_{22,t}\end{bmatrix}=\begin{bmatrix}2.80\\1.71\end{bmatrix}+\begin{bmatrix}0.084&0\\0.037&0.054\end{bmatrix}a_{t-1}^2+\begin{bmatrix}0.864&-0.020\\-0.058&0.914\end{bmatrix}\Xi_{t-1}^*,\tag{10.26}\]

\[q_t=-2.024+3.983\rho_{t-1}+0.088\frac{a_{1,t-1}a_{2,t-1}}{\sqrt{\sigma_{11,t-1}\sigma_{22,t-1}}}\tag{10.27}\]
（标准误 0.050、0.090、0.019，高度显著）。$Q_2(4)=20.57(0.11)$、$Q_2(8)=36.08(0.21)$；$Q_2^*(4)=16.69(0.27)$、$Q_2^*(8)=36.71(0.19)$。比较：均值与波动方程接近常相关模型；拟合相关（图 10.9）随时间波动、近年变小；平均 0.612 ≈ 常相关 0.614；用 $t=4$ 到 888 的观测，最大对数似然常相关 −3691.21，时变相关 −3679.64，有显著改进。
- 1 步预测（原点 $h=888$）：常相关模型 $a_{1,888}=3.075$，$a_{2,888}=4.931$，$\sigma_{11,888}=77.91$，$\sigma_{22,888}=21.19$，$\hat\Sigma_{888}(1)=\begin{bmatrix}71.09&21.83\\21.83&17.79\end{bmatrix}$；时变相关模型 $a_{1}=3.287$，$a_2=4.950$，$\sigma_{11}=83.35$，$\sigma_{22}=28.56$，$\rho_{888}=0.546$，$\hat\Sigma_{888}(1)=\begin{bmatrix}75.15&23.48\\23.48&24.70\end{bmatrix}$，相关预测 0.545。
方法二：Cholesky 分解。二元时 $\Xi_t=(g_{11,t},g_{22,t},q_{21,t})'$，一个简单 GARCH(1,1) 型模型： $g_{11,t}=\alpha_{10}+\alpha_{11}b_{1,t-1}^2+\beta_{11}g_{11,t-1},$ $q_{21,t}=\gamma_0+\gamma_1q_{21,t-1}+\gamma_2a_{2,t-1},\tag{10.28}$ $g_{22,t}=\alpha_{20}+\alpha_{21}b_{1,t-1}^2+\alpha_{22}b_{2,t-1}^2+\beta_{21}g_{11,t-1}+\beta_{22}g_{22,t-1},$ $b_{1t}=a_{1t}$，$b_{2t}=a_{2t}-q_{21,t}a_{1t}$。$b_{1t}$ 为单变量 GARCH(1,1)，$b_{2t}$ 为二元 GARCH(1,1)，$q_{21,t}$ 自相关并以 $a_{2,t-1}$ 为额外解释变量；似然用 (10.20)，$k=2$。
例 10.5 续（方法二）：均值 $r_{1t}=1.364+0.075r_{1,t-1}-0.058r_{2,t-2}+a_{1t}$，$r_{2t}=0.643+a_{2t}$；
\[g_{11,t}=3.714+0.113b_{1,t-1}^2+0.804g_{11,t-1},\]

\[q_{21,t}=0.0029+0.9915q_{21,t-1}-0.0041a_{2,t-1},\tag{10.29}\]

\[g_{22,t}=1.023+0.021b_{1,t-1}^2+0.052b_{2,t-1}^2-0.040g_{11,t-1}+0.937g_{22,t-1},\]
全部 1% 显著。时变相关
\[\rho_t=\frac{\sigma_{21,t}}{\sqrt{\sigma_{11,t}\sigma_{22,t}}}=\frac{q_{21,t}\sqrt{g_{11,t}}}{\sqrt{g_{22,t}+q_{21,t}^2g_{11,t}}}.\tag{10.30}\]
$Q_2(4)=19.77(0.14)$、$Q_2(8)=34.22(0.27)$；$Q_2^*(4)=15.34(0.36)$、$Q_2^*(8)=31.87(0.37)$，模型充分；相关有很强的动态依赖（0.9915）。图 10.10：相关比图 10.9 更平滑，确认下降趋势且近年更小。两种时变相关模型最大似然都约 −3672，拟合相近。
- Cholesky 法优点：无需正定约束（再对 $g_{ii,t}$ 取对数则整个模型无约束）；似然简单；$q_{ij,t},g_{ii,t}$ 可解释。缺点：结果依赖 $a_t$ 分量的排序（$a_{1t}$ 不变换），而理论上排序不应影响波动率，推断稍复杂。
- 1 步预测 $\hat\Sigma_{888}(1)=\begin{bmatrix}73.45&7.34\\7.34&17.87\end{bmatrix}$，相关 0.203，远小于前两模型；方差预测相近。教学提示：不同多元波动模型对相关的预测可差异很大，直接影响组合风险与对冲比。

10.4.3 动态相关模型（Dynamic Correlation Models, DCC）（PDF p.551–557）

基于 (10.7) $\Sigma_t=D_t\rho_tD_t$，对 $\rho_t$ 建简约模型，统称动态条件相关（DCC）模型。

Tse & Tsui (2002)： $\rho_t=(1-\theta_1-\theta_2)\bar\rho+\theta_1\rho_{t-1}+\theta_2\psi_{t-1},$ $\theta_1,\theta_2$ 为标量，$\bar\rho$ 为单位对角正定阵，$\psi_{t-1}$ 为用 $t-m,\dots,t-1$ 期冲击（标准化）计算的样本相关阵（$m$ 预设）。通常 $0\le\theta_i<1$、$\theta_1+\theta_2<1$ 保证 $\rho_t$ 正定。$\bar\rho$ 可取收益样本相关阵，则相关方程只有两个参数；$\bar\rho$ 和 $m$ 的选择需仔细研究。
Engle (2002)： $\rho_t=J_tQ_tJ_t,\quad J_t=\mathrm{diag}\{q_{11,t}^{-1/2},\dots,q_{kk,t}^{-1/2}\},$ $Q_t=(1-\theta_1-\theta_2)\bar Q+\theta_1\epsilon_{t-1}\epsilon_{t-1}'+\theta_2Q_{t-1},$ $\epsilon_{it}=a_{it}/\sqrt{\sigma_{ii,t}}$ 为标准化新息，$\bar Q$ 为 $\epsilon_t$ 的无条件协方差阵，$\theta_1,\theta_2\ge0$，$0<\theta_1+\theta_2<1$。$J_t$ 是归一化矩阵，保证结果为相关阵。实务上分两步估计：先逐个单变量 GARCH，再估相关方程（教材可补充）。
共同缺点：$\theta_1,\theta_2$ 为标量，所有相关对动态相同，$k$ 大时难以自圆其说。
Tsay (2006) 推广：(1) 标准化新息服从多元 Student-t（式 10.42）；(2) 边际波动率含杠杆效应： $D_t^2=\Lambda_0+\Lambda_1D_{t-1}^2+\Lambda_2A_{t-1}^2+\Lambda_3L_{t-1}^2,\tag{10.31}$ $\Lambda_i=\mathrm{diag}\{\ell_{1i},\dots,\ell_{ki}\}$，$A_{t-1}=\mathrm{diag}\{a_{1,t-1},\dots,a_{k,t-1}\}$，$L_{t-1}$ 对角元 $L_{i,t-1}=a_{i,t-1}$（若 $a_{i,t-1}<0$）否则 0。约束 $0\le\sum_{j=1}^3\ell_{ij}<1$，$\ell_{i0}>0$，$\ell_{ji}\ge0$ 保证波动率存在；$\Lambda_3=0$ 无杠杆。相关方程 $\rho_t=(1-\theta_1-\theta_2)\hat\rho+\theta_1\psi_{t-1}+\theta_2\rho_{t-1},\tag{10.32}$ $\hat\rho$ 为收益样本相关阵，$\theta_i\ge0$，$\theta_1+\theta_2<1$。
例 10.6：1999-01 至 2004-12 美元兑欧元、日元日汇率（圣路易斯联储正午即期）及 IBM、Dell 日股价（CRSP），简单收益（%），剔除任一市场休市日，1496 个观测，$r_t$=(USEU, JPUS, IBM, DELL)。表 10.1：均值 0.0091、−0.0059、0.0066、0.0028（≈0）；标准差 0.6469、0.6626、5.4280、10.1954（原表数值如此，股票明显更大；疑为方差或单位问题，按原文记录）；偏度 0.034、−0.167、−0.053、−0.038；超额峰度 2.71、2.03、6.22、3.31（均厚尾）；Q(12)=12.5、6.4、24.1、24.1。多元 Ljung–Box $Q(3)=59.12$（p=0.13），$Q(5)=106.44$（p=0.03），股票有轻微序列相关，简单起见均值方程用样本均值。
- 表 10.2（$L_{\max}$ 为最大似然值，$v$ 为多元 t 自由度）：(a) 全模型 $L_{\max}=-9175.80$：四个资产的 $(\Lambda_0,\Lambda_1,\Lambda_2)$ 分别为 USEU (0.0041, 0.9701, 0.0214)、JPUS (0.0088, 0.9515, 0.0281)、IBM (0.0071, 0.9636, 0.0326)、DELL (0.0150, 0.9531, 0.0461)；$(v,\theta_1,\theta_2)=(7.8729, 0.9808, 0.0137)$。注意原表把 0.98 记为 $\theta_1$，但按 (10.32) 的写法 $\theta_1$ 乘的是 $\psi_{t-1}$；从数值看 0.98 应是滞后相关阵 $\rho_{t-1}$ 的系数（高持续性），0.014 是局部样本相关 $\psi_{t-1}$ 的系数，编写时宜统一记号；(b) 限制 $\Lambda_1=\lambda I$，$L_{\max}=-9176.62$；(c) 最终限制模型（$\Lambda_0=(\lambda_1,\lambda_1,\lambda_3,\lambda_4)$，$\Lambda_1=\lambda I$，$\Lambda_2=(b_1,b_1,b_2,b_2)$），$L_{\max}=-9177.44$：$\lambda_1=0.0067$、$\lambda_3=0.0061$（不显著）、$\lambda_4=0.0148$（不显著），$\lambda=0.9603$，$b_1=0.0248$，$b_2=0.0372$，$v=7.918$，$\theta=(0.9809,0.0137)$；(d) 含杠杆效应 $L_{\max}=-9169.04$，$v=8.45$。
- 模型检验：对 $\hat\epsilon_t=\hat\Sigma_t^{-1/2}a_t$（对称平方根）做多元 Ljung–Box：(a) $Q(10)=167.79(0.32)$，平方 110.19(1.00)；(b) 168.59(0.31)、109.93(1.00)；(c) 168.50(0.31)、111.75(1.00)；都充分。
- 结论：(1) LR 检验不能拒绝最终限制模型，仅 9 个参数刻画四维收益的时变相关与波动；(2) 两只股票 $\Lambda_0$ 常数不显著，GARCH 参数和 $0.0372+0.9603=0.9975\approx1$，呈 IGARCH；汇率有非零常数和高持续性；(3) 与 69 日（约一个季度）滚动窗口估计比较（图 10.12、10.13）：形态相似，但滚动估计对大冲击反应更慢，模型波动率升降更快；模型相关更平滑。
- (d) 杠杆效应仅对股票显著且呈 IGARCH 形式：$\Lambda_3=\mathrm{diag}\{0,0,1-0.96-0.0241,1-0.96-0.0286\}=\mathrm{diag}\{0,0,0.0159,0.0114\}$。量虽小但显著：(b) 与 (d) 比较 LR=15.16，$\chi^2_2$ 下 p=0.0005。

10.5 高维波动率模型（Higher Dimensional Volatility Models）（PDF p.557–）

利用 Cholesky 分解的序贯性：$a_{it}\to b_{it}$ 的正交变换只涉及 $b_{jt}$（$j<i$），10.4 的模型是嵌套的（$g_{ii,t}$ 只依赖 $j<i$ 的量）。序贯建模步骤：
1. 选一个最关心的市场指数或股票，建单变量波动率模型；
2. 加入第二个序列，对其冲击做正交变换，建二元模型（以第 1 步估计为初值）；
3. 加入第三个序列，同样正交变换建三维模型（以二元估计为初值）；
4. 继续直到包含所有序列。每步做模型检验。经验表明可大幅简化高维建模并显著减少计算时间。
例 10.7：S&P 500 指数、Cisco、Intel 日对数收益（%），1991-01-02 至 1999-12-31，2275 个观测，排序 $r_t=(SP5,CSCO,INTC)'$。样本均值 (0.066, 0.257, 0.156)，标准差 (0.875, 2.853, 2.464)，相关 0.52、0.50、0.47。$Q_3(1)=26.20$、$Q_3(4)=79.73$、$Q_3(8)=123.68$（自由度 9、36、72，p≈0），有序列相关。表 10.3 交叉相关阵：(a) 指数不依赖 Cisco、Intel 过去收益；(b) Cisco 有自相关且依赖指数滞后 2、5 期；(c) Intel 依赖指数滞后 1、5 期——大盘股收益受市场过去表现影响，反之不然（与第 8 章 IBM 结果一致）。
- 第 1 步（S&P 500）： $r_{1t}=0.078+0.042r_{1,t-1}-0.062r_{1,t-3}-0.048r_{1,t-4}-0.052r_{1,t-5}+a_{1t},\quad \sigma_{11,t}=0.013+0.092a_{1,t-1}^2+0.894\sigma_{11,t-1}.\tag{10.33}$ $Q(10)=7.38(0.69)$，平方 3.14(0.98)。
- 第 2 步（加 Cisco）：均值 $r_{1t}=0.065-0.046r_{1,t-3}+a_{1t}$，$r_{2t}=0.325+0.195r_{1,t-2}-0.091r_{2,t-2}+a_{2t}$（10.34）； $g_{11,t}=0.006+0.051b_{1,t-1}^2+0.943g_{11,t-1},\ q_{21,t}=0.331+0.790q_{21,t-1}-0.041a_{2,t-1},\ g_{22,t}=0.177+0.082b_{2,t-1}^2+0.890g_{22,t-1}.\tag{10.35}$ 二元 Ljung–Box 无问题；$r_{1t}$ 的边际模型与单变量模型差别很小。
- 第 3 步（加 Intel）：均值 $r_{1t}=0.065-0.043r_{1,t-3}+a_{1t}$，$r_{2t}=0.326+0.201r_{1,t-2}-0.089r_{2,t-1}+a_{2t}$，$r_{3t}=0.192-0.264r_{1,t-1}+0.059r_{3,t-1}+a_{3t}$（10.36）； $\begin{aligned}g_{11,t}&=0.006+0.050b_{1,t-1}^2+0.943g_{11,t-1},\\ q_{21,t}&=0.277+0.824q_{21,t-1}-0.035a_{2,t-1},\\ g_{22,t}&=0.178+0.082b_{2,t-1}^2+0.889g_{22,t-1},\\ q_{31,t}&=0.039+0.973q_{31,t-1}+0.010a_{3,t-1},\\ q_{32,t}&=0.006+0.981q_{32,t-1}+0.004a_{2,t-1},\\ g_{33,t}&=1.188+0.053b_{3,t-1}^2+0.687g_{33,t-1}-0.019g_{22,t-1},\end{aligned}\tag{10.37}$ $b_{3t}=a_{3t}-q_{31,t}b_{1t}-q_{32,t}b_{2t}$；标准误见表 10.4（如 $g_{33,t}$ 常数 0.407）；除 $q_{32,t}$ 常数外均 5% 显著。标准化残差 $Q_3(4)=34.48(0.31)$、$Q_3(8)=60.42(0.70)$（自由度调整为 31、67）；平方 $Q_3^*(4)=28.71(0.58)$、$Q_3^*(8)=52.00(0.91)$，充分。
- 特征：(1) 本质上是时变相关 GARCH(1,1)；(2) 指数波动率不依赖 Cisco、Intel 的过去波动；(3) 经 Cholesky 逆变换，Cisco、Intel 的波动率依赖市场过去波动；(4) $q_{ij,t}$ 持续性高（AR(1) 系数大）。图 10.15：指数波动远小于个股，近年指数波动上升而 Cisco 并非如此。图 10.16：收益波动大时相关系数上升（续）（PDF p.561–563 续）：与国际股市指数的实证研究一致——金融危机期间市场间相关上升（波动–相关正向关系）。
模型 (10.37) 由两组方程组成：条件方差（$g_{ii,t}$）与相关量（$q_{ij,t}$，$i>j$）。此数据集中 AR(1) 型方程可能足够。令 $v_t=(\ln g_{11,t},\ln g_{22,t},\ln g_{33,t})'$，$q_t=(q_{21,t},q_{31,t},q_{32,t})'$，可用确定性滞后 1 模型 $v_t=c_1+\beta_1v_{t-1}$，$q_t=c_2+\beta_2q_{t-1}$；若加入噪声项 $v_t=c_1+\beta_1v_{t-1}+e_{1t}$，$q_t=c_2+\beta_2q_{t-1}+e_{2t}$，就得到简单的多元随机波动率（SV）模型。Chib, Nardari & Shephard (1999) 用 MCMC 研究高维 SV（允许时变相关但较受限）；另见 Harvey, Ruiz & Shephard (1994)；第 12 章讨论 MCMC。

10.6 因子–波动率模型（Factor–Volatility Models）（PDF p.563–566）

用因子模型简化多元波动率动态。"共同因子"可由实质性知识或经验方法确定；简单办法是对 $a_t=r_t-\mu_t$ 做 PCA。三步：(1) 选解释 $a_t$ 大部分变异的前几个主成分；(2) 对所选主成分建波动率模型；(3) 把各 $a_{it}$ 的波动率与主成分波动率联系起来。目标：降维同时保持多元波动率的准确近似。
例 10.8（IBM 与 S&P 500 月度对数收益，同例 10.5）：用例 8.4 的二元 AR(3) 得新息 $a_t$，协方差 PCA 特征值 63.373、13.489，第一主成分解释 82.2% 的广义方差，$x_t=0.797a_{1t}+0.604a_{2t}$。因序列相关弱，直接对 $r_t$ 做 PCA：特征值 63.625、13.513，第一主成分解释约 82.5%，$x_t=0.796r_{1t}+0.605r_{2t}$——均值方程对 PCA 影响可忽略，取后者为共同因子（图 10.17(a)）。
- 因子的单变量高斯 GARCH： $x_t=1.317+0.096x_{t-1}+a_t,\quad \sigma_t^2=3.834+0.110a_{t-1}^2+0.825\sigma_{t-1}^2,\tag{10.38}$ 全部 1% 显著，Ljung–Box 无问题（图 10.17(b) 为拟合波动率）。
- 以 $\sigma_t^2$ 为共同波动率因子：均值 $r_{1t}=1.140+0.079r_{1,t-1}+0.067r_{1,t-2}-0.122r_{2,t-2}+a_{1t}$，$r_{2t}=0.537+a_{2t}$； $\begin{bmatrix}\sigma_{11,t}\\\sigma_{22,t}\end{bmatrix}=\begin{bmatrix}19.08\\-5.62\end{bmatrix}+\begin{bmatrix}0.098&0\\0&0\end{bmatrix}\begin{bmatrix}a_{1,t-1}^2\\a_{2,t-1}^2\end{bmatrix}+\begin{bmatrix}0.333\\0.596\end{bmatrix}\sigma_t^2,\tag{10.39}$ $\rho_t=\frac{\exp(q_t)}{1+\exp(q_t)},\quad q_t=-2.098+4.120\rho_{t-1}+0.078\frac{a_{1,t-1}a_{2,t-1}}{\sqrt{\sigma_{11,t-1}\sigma_{22,t-1}}}.\tag{10.40}$ （标准误：常数 3.70、2.36；ARCH 0.044；因子载荷 0.076、0.050；相关方程 0.025、0.038、0.015。）
- 检验：$Q_2(4)=15.37(0.29)$、$Q_2(8)=34.24(0.23)$ 无序列相关；但平方 $Q_2^*(4)=20.25(0.09)$、$Q_2^*(8)=61.95(0.0004)$——高阶滞后的条件异方差未处理好，因单一因子只解释约 82.5% 的广义方差。
- 与时变相关模型 (10.26)–(10.27) 比较：相关方程基本相同；因子模型波动方程参数更少；是合理近似。
备注：此处用两步估计（先建因子波动模型，再把估计波动率当已知），简单但未必有效；若共同因子已知（如 $x_t=0.796r_{1t}+0.605r_{2t}$，原文此处误写为 0.769），可对 (10.38)–(10.40) 联合估计，更有效。

10.7 应用：多资产 VaR（Application）（PDF p.566–568）

投资者各持有 Cisco 与 Intel 多头 100 万美元，用 1991-01-02 至 1999-12-31 日对数收益（%）建模，以数据末端（$t=2275$）的 1 步预测和 5% 临界值计算日 VaR。由第 7 章，组合 VaR $\mathrm{VaR}=\sqrt{\mathrm{VaR}_1^2+\mathrm{VaR}_2^2+2\rho\,\mathrm{VaR}_1\mathrm{VaR}_2}.$ 收益为百分比，分位数除以 100；$r_{1t}$=Cisco，$r_{2t}$=Intel；所有估计 5% 显著、模型均通过 Ljung–Box 检验。
单变量模型 + 样本相关： $r_{1t}=0.380+0.034r_{1,t-1}-0.061r_{1,t-2}-0.055r_{1,t-3}+a_{1t}$，$\sigma_{1t}^2=0.599+0.117a_{1,t-1}^2+0.814\sigma_{1,t-1}^2$；$r_{2t}=0.187+a_{2t}$，$\sigma_{2t}^2=0.310+0.032a_{2,t-1}^2+0.918\sigma_{2,t-1}^2$；样本相关 0.473。预测 $\hat r_1=0.626$，$\hat\sigma_1^2=4.152$，$\hat r_2=0.187$，$\hat\sigma_2^2=6.087$。5% 分位数 $q_1=0.626-1.65\sqrt{4.152}=-2.736$，$q_2=0.187-1.65\sqrt{6.087}=-3.884$；$\mathrm{VaR}_1=\$27{,}360$，$\mathrm{VaR}_2=\$38{,}840$，总 VaR=$57,117。
常相关二元 GARCH(1,1)（对角）：$r_{1t}=0.385+0.038r_{1,t-1}-0.060r_{1,t-2}-0.047r_{1,t-3}+a_{1t}$，$r_{2t}=0.222+a_{2t}$，$\sigma_{11,t}=0.624+0.110a_{1,t-1}^2+0.816\sigma_{11,t-1}$，$\sigma_{22,t}=0.664+0.038a_{2,t-1}^2+0.853\sigma_{22,t-1}$，$\hat\rho=0.475$。预测 $\hat r_1=0.373$，$\hat\sigma_1^2=4.287$，$\hat r_2=0.222$，$\hat\sigma_2^2=5.706$；$\mathrm{VaR}_1=\$30{,}432$，$\mathrm{VaR}_2=\$37{,}195$，总 VaR=$58,180。
时变相关（Cholesky）：$r_{1t}=0.355+0.039r_{1,t-1}-0.057r_{1,t-2}-0.038r_{1,t-3}+a_{1t}$，$r_{2t}=0.206+a_{2t}$；$g_{11,t}=0.420+0.091b_{1,t-1}^2+0.858g_{11,t-1}$，$q_{21,t}=0.123+0.689q_{21,t-1}-0.014a_{2,t-1}$，$g_{22,t}=0.080+0.013b_{2,t-1}^2+0.971g_{22,t-1}$。预测 $\hat r_1=0.352$，$\hat r_2=0.206$，$\hat g_{11}=4.252$，$\hat q_{21}=0.421$，$\hat g_{22}=5.594$，故 $\hat\sigma_1^2=4.252$，$\hat\sigma_{21}=1.791$，$\hat\sigma_2^2=6.348$，$\hat\rho=0.345$；$\mathrm{VaR}_1=\$30{,}504$，$\mathrm{VaR}_2=\$39{,}512$，总 VaR=$57,648。
结论：三种方法结果相近，单变量最低、常相关最高，差距约 1100 美元；时变相关模型居中。（注意：此处 VaR 公式在均值非零时只是近似，教材可说明。）

10.8 多元 t 分布（Multivariate t Distribution）（PDF p.568–569）

多元高斯新息可能无法刻画收益峰度，用多元 Student-t。$k$ 维 $x$，自由度 $v$，$\mu=0$，$\Sigma=I$： $f(x|v)=\frac{\Gamma[(v+k)/2]}{(\pi v)^{k/2}\Gamma(v/2)}(1+v^{-1}x'x)^{-(v+k)/2}\tag{10.41}$ （Mardia, Kent & Bibby 1979, p.57）。各分量方差 $v/(v-2)$，定义标准化 $\epsilon_t=\sqrt{(v-2)/v}\,x$： $f(\epsilon_t|v)=\frac{\Gamma[(v+k)/2]}{[\pi(v-2)]^{k/2}\Gamma(v/2)}[1+(v-2)^{-1}\epsilon_t'\epsilon_t]^{-(v+k)/2}.\tag{10.42}$
波动率建模中 $a_t=\Sigma_t^{1/2}\epsilon_t$： $f(a_t|v,\Sigma_t)=\frac{\Gamma[(v+k)/2]}{[\pi(v-2)]^{k/2}\Gamma(v/2)|\Sigma_t|^{1/2}}[1+(v-2)^{-1}a_t'\Sigma_t^{-1}a_t]^{-(v+k)/2}.$ 用 Cholesky 分解（$a_t=L_tb_t$）： $f(b_t|v,L_t,G_t)=\frac{\Gamma[(v+k)/2]}{[\pi(v-2)]^{k/2}\Gamma(v/2)\prod_jg_{jj,t}^{1/2}}\Big[1+(v-2)^{-1}\sum_{j=1}^k\frac{b_{jt}^2}{g_{jj,t}}\Big]^{-(v+k)/2},$ 无需矩阵求逆，条件似然易算。注意：多元 t 下 $b_{jt}$ 不相关但不独立，似然不能拆成单变量乘积（与正态不同）。

附录：估计说明（Appendix: Some Remarks on Estimation）（PDF p.569–574）

多元 ARMA 用 SCA；多元波动率用 S-Plus FinMetrics、RATS 或 Matlab。给出 RATS 程序：
1. 例 10.5 对角常相关 AR(2)–GARCH(1,1)：nonlin 声明参数，frml 定义 $a_{1t},a_{2t}$、gvar1/gvar2 方差递推、二元正态对数似然 gln（含 $\ln(1-\rho^2)$ 项），maximize(method=bhhh, recursive, iterations=150)，再算标准化残差及其平方的 cor(qstats)，打印最后几期冲击和方差用于预测。数据 m-ibmspln.txt，888 个观测。
2. 例 10.5 时变相关模型：增加 rh1 = q0 + q1*rho(t-1) + q2*a1t*a2t/sqrt(h1*h2)，rh = exp(rh1)/(1+exp(rh1))，方差方程含交叉项 f1*h2(t-1)、f11*h1(t-1)、d11*a1t(t-1)**2。
3. 例 10.5 Cholesky 时变相关：qt = t0 + t1*q(t-1) + t2*a2t(t-1)，bt = a2t - qt*a1t，似然 $-\frac12[\ln h_1+\ln h_2+a_{1t}^2/h_1+b_t^2/h_2]$；$r_{2t}$ 方差 fv2 = v2 + qt^2*v1，相关 rhohat = qt*sqrt(v1/fv2)。
4. 例 10.7 三维 Cholesky 模型（d-cscointc.txt，2275 个观测，序贯估计得到初值）：b1t = a3t - q2t*a1t - q3t*bt，fv3 = v3 + q2t^2 v1 + q3t^2 v2，相关 rho21 = q1t*sqrt(v1/fv2)、rho31 = q2t*sqrt(v1/fv3)、rho32 = (q2t*q1t*v1 + q3t*v2)/sqrt(fv2*fv3)。

第 10 章习题（PDF p.574–575）

10.1：IBM、HPQ、S&P 指数 1962–2008 月度对数收益，EWMA 多元波动率，估计 $\lambda$ 并作图。
10.2：IBM 与 HPQ 拟合 DVEC(1,1)，检验并画波动率与时变相关。
10.3：S&P 与 HPQ 建 BEKK 模型。
10.4：三序列常相关模型并检验。
10.5：GE 与 S&P 500（1926–2008）常相关 GARCH，检验并给出 2008-12 原点的协方差 1 步预测。
10.6：GE、IBM、S&P 三维 DCC 模型（$\rho$ 取样本相关阵）。
10.7：GE 与 S&P（1926–1999）logistic 时变相关 GARCH 与 1 步预测。
10.8：同数据用 Cholesky 时变相关 GARCH，与 10.7 比较。
10.9：三维 Cholesky 时变波动率模型，$t=888$ 处 1 步预测。
10.10：持有 Dell 50 万与 Cisco 100 万美元多头，1990-02-20 至 1999-12-31 日数据，用 10.7 节三种方法计算 5% 日 VaR 并比较。

参考文献（PDF p.575）：Bauwens, Laurent & Rombouts (2004)、Bollerslev (1990)、Bollerslev, Engle & Wooldridge (1988)、Engle (2002)、Engle & Kroner (1995)、Pourahmadi (1999)、Tsay (2006)、Tse (2000)、Tse & Tsui (2002) 等。PDF p.576 为空白页。

第 10 章本章要点

多元波动率即条件协方差阵 $\Sigma_t$，难点是维数灾难与正定约束。
EWMA：$\hat\Sigma_t=(1-\lambda)a_{t-1}a_{t-1}'+\lambda\hat\Sigma_{t-1}$，一个参数，可作基准。
DVEC：逐元素 GARCH，简单但不保证正定、无溢出；BEKK：保证正定、允许溢出，但参数多且难解释。
两种重参数化：相关系数（$D_t\rho_tD_t$）直接但需约束；Cholesky（$L_tG_tL_t'$）等价于逐次回归正交化，无约束、似然可分解，但依赖变量排序。
常相关（CCC）模型、logistic/Fisher 变换的时变相关模型、Cholesky 时变相关模型、DCC（Tse–Tsui、Engle、Tsay 2006 含 t 分布与杠杆）。
高维时利用 Cholesky 嵌套性序贯建模；也可用 PCA 因子–波动率模型降维。
实证规律：波动大时相关上升；个股波动受市场过去波动影响，反之不然；股票日波动近 IGARCH，且有杠杆效应。
多元 VaR：$\sqrt{\mathrm{VaR}_1^2+\mathrm{VaR}_2^2+2\rho\mathrm{VaR}_1\mathrm{VaR}_2}$，相关估计方法会影响结果。

与量化交易的关联

风险管理：多元 GARCH/DCC 用于组合 VaR/ES 的条件协方差预测；EWMA（RiskMetrics，$\lambda\approx0.94$ 日频）是业界最常用的协方差估计；危机期间相关上升意味着分散化失效，压力测试应使用条件相关而非长期均值。
组合优化与动态资产配置：均值–方差/风险平价/最小方差组合需要 $\Sigma_t$ 预测，DCC 是标准选择；高维时用因子–波动率模型（因子 GARCH + 特质方差）或 Cholesky 序贯法控制参数数。
对冲与时变 β：Cholesky 中的 $q_{21,t}=\sigma_{21,t}/\sigma_{11,t}$ 正是时变回归系数（如个股对市场的时变 β、期货最优对冲比），可直接用于动态对冲与 β 中性调整。
配对/统计套利：时变相关与条件协方差可用于动态调整价差标准差和开仓阈值。
工程实现：Cholesky 参数化与多元 t 的无矩阵求逆似然便于实现与数值稳定；参数排序依赖性需注意（通常把市场指数放在第一位）。

第 11 章状态空间模型与卡尔曼滤波（State-Space Models and Kalman Filter）（PDF p.577 起）

本章引言（PDF p.577）

状态空间模型为时间序列分析提供灵活框架，尤其便于简化极大似然估计和处理缺失值。本章讨论状态空间模型与 ARIMA 的关系、Kalman 滤波算法、各种平滑方法及应用（已实现波动率、时变系数市场模型、公司季度 EPS）。参考书：Durbin & Koopman (2001)、Kim & Nelson (1999，经济应用与区制转换)、Anderson & Moore (1979，工程与最优控制)、Chan (2002)、Shumway & Stoffer (2000)、Hamilton (1994)、Harvey (1993)、West & Harrison (1997，贝叶斯预测)、Kitagawa & Gersch (1996，平滑先验)。11.4 节推导记号繁重，初读可跳过。

11.1 局部趋势模型（Local Trend Model）（PDF p.578–）

模型： $y_t=\mu_t+e_t,\quad e_t\sim N(0,\sigma_e^2),\tag{11.1}$ $\mu_{t+1}=\mu_t+\eta_t,\quad \eta_t\sim N(0,\sigma_\eta^2),\tag{11.2}$ $\{e_t\},\{\eta_t\}$ 为独立高斯白噪声，初值 $\mu_1$ 给定或服从已知分布且与噪声独立。$\mu_t$ 为随机游走，称趋势（trend），不可直接观测；$y_t$ 是带观测噪声的观测值；$y_t$ 的动态依赖完全来自 $\mu_t$。
应用：已实现波动率——$\mu_t$ 为真实对数波动率（随机游走演化），$y_t$ 为对数已实现波动率（由高频数据构造，受市场微观结构噪声影响），$\sigma_e$ 衡量微观结构噪声大小。
术语：这是线性高斯状态空间模型的特例。$\mu_t$ 为状态（state）；(11.1) 为观测方程（observation equation），$e_t$ 为测量误差；(11.2) 为状态方程/状态转移方程（state transition equation），$\eta_t$ 为状态新息。Durbin & Koopman 称局部水平模型（local-level model），是 Harvey (1993) 结构时间序列模型的简单情形。

与 ARIMA 的关系（PDF p.578–579）：

$\sigma_e=0$ 时 $y_t=\mu_t$ 为 ARIMA(0,1,0)。$\sigma_e>0$ 时 $y_t$ 为 ARIMA(0,1,1)：$(1-B)y_t=(1-\theta B)a_t$（11.3）。
推导：$\mu_{t+1}=\frac{1}{1-B}\eta_t$，故 $y_t=\frac{1}{1-B}\eta_{t-1}+e_t$，乘 $(1-B)$ 得 $w_t=(1-B)y_t=\eta_{t-1}+e_t-e_{t-1}$。$w_t$ 高斯，$\mathrm{Var}(w_t)=2\sigma_e^2+\sigma_\eta^2$，$\mathrm{Cov}(w_t,w_{t-1})=-\sigma_e^2$，更高阶为 0，故为 MA(1)。匹配矩： $(1+\theta^2)\sigma_a^2=2\sigma_e^2+\sigma_\eta^2,\tag{11.4}$ $\theta\sigma_a^2=\sigma_e^2.\tag{11.5}$ 两式相除得 $\theta$ 的二次方程，取 $|\theta|<1$ 的根，再求 $\sigma_a^2$。这就是第 2 章的简单指数平滑模型。
反过来，$\theta>0$ 的 ARIMA(0,1,1) 可解出 $\sigma_e^2,\sigma_\eta^2$ 得局部趋势模型；$\theta<0$ 时仍可写成无观测误差的状态空间形式。ARIMA 可以多种方式转为状态空间模型。仅凭数据，选 ARIMA 还是状态空间形式并不关键，取决于分析目的、实质问题与经验。
例 11.1：Alcoa 股票日内已实现波动率，2003-01-02 至 2004-05-07，340 个观测；日已实现波动率为日内 10 分钟对数收益（%）平方和（不含隔夜收益和第一个 10 分钟收益），分析其对数（数据来自 NYSE TAQ，图 11.1）。ARIMA 拟合 $(1-B)y_t=(1-0.858B)a_t,\quad\hat\sigma_a=0.5184,\tag{11.6}$ $\hat\theta$ 标准误 0.028；残差 $Q(12)=12.4$（p=0.33），平方残差 $Q(12)=8.2$（p=0.77），无 ARCH 效应。因 $\hat\theta>0$ 可转为局部趋势模型：MLE $\hat\sigma_\eta=0.0735$，$\hat\sigma_e=0.4803$——测量误差方差远大于状态新息，证实日内高频收益受测量误差影响。由 (11.6) 与 (11.4)(11.5) 换算得 $\sigma_e=0.480$，$\sigma_\eta=0.0736$，与 MLE 接近。

11.1.1 统计推断（Statistical Inference）（PDF p.581–582）

设 $F_t=\{y_1,\dots,y_t\}$，模型参数已知。三类推断：
- 滤波（filtering）：给定 $F_t$ 恢复 $\mu_t$（去除测量误差）；
- 预测（prediction）：给定 $F_t$ 预测 $\mu_{t+h}$ 或 $y_{t+h}$（$h>0$）；
- 平滑（smoothing）：给定 $F_T$（$T>t$）估计 $\mu_t$。类比读手写便条：滤波是根据已读内容辨认当前字，预测是猜下一个字，平滑是读完全文后再辨认某个字。
记号：$\mu_{t|j}=E(\mu_t|F_j)$，$\Sigma_{t|j}=\mathrm{Var}(\mu_t|F_j)$，$y_{t|j}=E(y_t|F_j)$；1 步预测误差 $v_t=y_t-y_{t|t-1}$，方差 $V_t=\mathrm{Var}(v_t|F_{t-1})=\mathrm{Var}(v_t)$（与 $F_{t-1}$ 独立）。由 (11.1)，$y_{t|t-1}=\mu_{t|t-1}$， $v_t=y_t-\mu_{t|t-1},\tag{11.7}$ $V_t=\Sigma_{t|t-1}+\sigma_e^2.\tag{11.8}$ 且 $E(v_t)=0$，$\mathrm{Cov}(v_t,y_j)=0$（$j<t$）：预测误差（新息）与过去观测不相关（正态下独立）。$F_t=\{F_{t-1},y_t\}=\{F_{t-1},v_t\}$。
定理 11.1（多元正态条件分布）：$x,y,z$ 联合正态，$\Sigma_{ww}$ 非奇异，且 $\Sigma_{yz}=0$，则
1. $E(x|y)=\mu_x+\Sigma_{xy}\Sigma_{yy}^{-1}(y-\mu_y)$；
2. $\mathrm{Var}(x|y)=\Sigma_{xx}-\Sigma_{xy}\Sigma_{yy}^{-1}\Sigma_{yx}$；
3. $E(x|y,z)=E(x|y)+\Sigma_{xz}\Sigma_{zz}^{-1}(z-\mu_z)$；
4. $\mathrm{Var}(x|y,z)=\mathrm{Var}(x|y)-\Sigma_{xz}\Sigma_{zz}^{-1}\Sigma_{zx}$。（第 3、4 条说明新增一个与已有信息不相关的变量时，条件均值和方差可"增量更新"——Kalman 滤波的核心。）

11.1.2 Kalman 滤波（Kalman Filter）（PDF p.582–）

目标：已知 $\mu_t|F_{t-1}$ 的分布与新数据 $y_t$，递推得到 $\mu_t|F_t$。因 $F_t=\{F_{t-1},v_t\}$，只需 $(\mu_t,v_t)'|F_{t-1}$ 的联合分布。
关键协方差： $\mathrm{Cov}(\mu_t,v_t|F_{t-1})=E[\mu_t(\mu_t+e_t-\mu_{t|t-1})|F_{t-1}]=E[(\mu_t-\mu_{t|t-1})^2|F_{t-1}]=\Sigma_{t|t-1},\tag{11.9}$ （用到 $E[\mu_{t|t-1}(\mu_t-\mu_{t|t-1})|F_{t-1}]=0$）。于是 $\begin{bmatrix}\mu_t\\v_t\end{bmatrix}\Big|F_{t-1}\sim N\left(\begin{bmatrix}\mu_{t|t-1}\\0\end{bmatrix},\begin{bmatrix}\Sigma_{t|t-1}&\Sigma_{t|t-1}\\\Sigma_{t|t-1}&V_t\end{bmatrix}\right).$
由定理 11.1，更新步： $\mu_{t|t}=\mu_{t|t-1}+\frac{\Sigma_{t|t-1}}{V_t}v_t=\mu_{t|t-1}+K_tv_t,\tag{11.10}$ $\Sigma_{t|t}=\Sigma_{t|t-1}-\frac{\Sigma_{t|t-1}^2}{V_t}=\Sigma_{t|t-1}(1-K_t),\tag{11.11}$ $K_t=\Sigma_{t|t-1}/V_t$ 称 Kalman 增益（Kalman gain），即 $\mu_t$ 对 $v_t$ 的回归系数，决定新冲击 $v_t$ 对状态估计的贡献。
预测步（由 11.2）： $\mu_{t+1|t}=\mu_{t|t},\tag{11.12}$ $\Sigma_{t+1|t}=\Sigma_{t|t}+\sigma_\eta^2.\tag{11.13}$ 观测到 $y_{t+1}$ 后重复，即 Kalman (1960) 算法。（续）
局部趋势模型的 Kalman 滤波汇总（初值 $\mu_1\sim N(\mu_{1|0},\Sigma_{1|0})$）：
\[\begin{aligned}v_t&=y_t-\mu_{t|t-1},\quad V_t=\Sigma_{t|t-1}+\sigma_e^2,\quad K_t=\Sigma_{t|t-1}/V_t,\\ \mu_{t+1|t}&=\mu_{t|t-1}+K_tv_t,\quad \Sigma_{t+1|t}=\Sigma_{t|t-1}(1-K_t)+\sigma_\eta^2,\quad t=1,\dots,T.\end{aligned}\tag{11.14}\]
Kalman 滤波推导方法很多，此处用定理 11.1 最简。初值选择见 11.1.6；$\sigma_e,\sigma_\eta$ 用 ML 估计，似然也靠 Kalman 滤波计算（11.1.7）。
- 稳态直观：$K_t$ 收敛到常数 $K$ 时，$\mu_{t+1|t}=(1-K)\mu_{t|t-1}+Ky_t$ 就是平滑常数为 $K$ 的指数平滑，与 ARIMA(0,1,1) 等价（$K=1-\theta$）。
例 11.1 续：取 $\Sigma_{1|0}=\infty$（扩散初始化）、$\mu_{1|0}=0$，图 11.2(a) 滤波状态 $\mu_{t|t}$ 比原序列平滑；(b) 1 步预测误差 $v_t$ 稳定、围绕 0，属样本外 1 步预测误差。

11.1.3 预测误差的性质（PDF p.584–586）

给定与 $y_t$ 独立的初值，$v_t$ 是 $\{y_1,\dots,y_t\}$ 的线性函数：$v_1=y_1-\mu_{1|0}$，$v_2=y_2-\mu_{1|0}-K_1(y_1-\mu_{1|0})$，$v_3=y_3-\mu_{1|0}-K_2(y_2-\mu_{1|0})-K_1(1-K_2)(y_1-\mu_{1|0})$……矩阵形式 $v=K(y-\mu_{1|0}1_T),\tag{11.15}$ $K$ 为单位下三角阵，$k_{i,i-1}=-K_{i-1}$，$k_{ij}=-(1-K_{i-1})(1-K_{i-2})\cdots(1-K_{j+1})K_j$（$j\le i-2$）。Kalman 增益 $K_t$ 不依赖 $\mu_{1|0}$ 和数据，只依赖 $\Sigma_{1|0},\sigma_e^2,\sigma_\eta^2$。
含义一：正态下 $\{v_t\}$ 相互独立。证明：变换雅可比为 1，$p(v)=p(y)=p(y_1)\prod_{j\ge2}p(y_j|F_{j-1})=\prod_jp(v_j)$。
含义二：Kalman 滤波给出 $\Omega=\mathrm{Cov}(y)$ 的 Cholesky 分解：$\mathrm{Cov}(v)=K\Omega K'=\mathrm{diag}\{V_1,\dots,V_T\}$（与第 10 章 Cholesky 正交化同理）。
状态误差递推：$x_t=\mu_t-\mu_{t|t-1}$，$\mathrm{Var}(x_t|F_{t-1})=\Sigma_{t|t-1}$， $v_t=x_t+e_t,\qquad x_{t+1}=L_tx_t+\eta_t-K_te_t,\tag{11.16}$ $L_t=1-K_t=\sigma_e^2/V_t$，$x_1=\mu_1-\mu_{1|0}$。这是以 $x_t$ 为状态、$v_t$ 为观测的时变状态空间模型。

11.1.4 状态平滑（State Smoothing）（PDF p.586–590）

目标：$\mu_t|F_T\sim N(\mu_{t|T},\Sigma_{t|T})$（平滑状态与平滑状态方差）。利用：$\{v_t,\dots,v_T\}$ 互相独立、与 $F_{t-1}$ 独立、均值 0、方差 $V_j$；$F_T$ 等价于 $\{F_{t-1},v_t,\dots,v_T\}$。由定理 11.1(3)： $\mu_{t|T}=\mu_{t|t-1}+\sum_{j=t}^T\mathrm{Cov}(\mu_t,v_j)V_j^{-1}v_j.\tag{11.17}$
协方差：$\mathrm{Cov}(\mu_t,v_j)=\mathrm{Cov}(x_t,v_j)$；$\mathrm{Cov}(x_t,v_t)=\Sigma_{t|t-1}$，$\mathrm{Cov}(x_t,v_{t+1})=\Sigma_{t|t-1}L_t$，……，$\mathrm{Cov}(x_t,v_T)=\Sigma_{t|t-1}\prod_{j=t}^{T-1}L_j$。故 $\mu_{t|T}=\mu_{t|t-1}+\Sigma_{t|t-1}q_{t-1}$， $q_{t-1}=\frac{v_t}{V_t}+L_t\frac{v_{t+1}}{V_{t+1}}+L_tL_{t+1}\frac{v_{t+2}}{V_{t+2}}+\cdots+\Big(\prod_{j=t}^{T-1}L_j\Big)\frac{v_T}{V_T},\tag{11.18}$ 满足 $q_{t-1}=v_t/V_t+L_tq_t$，$q_T=0$（11.19）。后向平滑递推： $q_{t-1}=V_t^{-1}v_t+L_tq_t,\qquad \mu_{t|T}=\mu_{t|t-1}+\Sigma_{t|t-1}q_{t-1},\quad t=T,\dots,1.\tag{11.20}$
平滑状态方差：由定理 11.1(4)， $\Sigma_{t|T}=\Sigma_{t|t-1}-\sum_{j=t}^T[\mathrm{Cov}(\mu_t,v_j)]^2V_j^{-1}=\Sigma_{t|t-1}-\Sigma_{t|t-1}^2M_{t-1},\tag{11.21, 11.22}$ $M_{t-1}=\frac1{V_t}+L_t^2\frac1{V_{t+1}}+\cdots+\big(\prod_{j=t}^{T-1}L_j^2\big)\frac1{V_T}=\mathrm{Var}(q_{t-1})$，$M_T=0$。后向递推 $M_{t-1}=V_t^{-1}+L_t^2M_t,\qquad \Sigma_{t|T}=\Sigma_{t|t-1}-\Sigma_{t|t-1}^2M_{t-1},\quad t=T,\dots,1.\tag{11.23}$
例 11.1 续：图 11.3 为滤波状态 $\mu_{t|t}$ 及 95% 逐点区间，图 11.4 为平滑状态 $\mu_{t|T}$ 及区间。平滑状态更平滑，区间更窄；$\mu_{1|1}$ 区间宽度取决于 $\Sigma_{1|0}$。
算法复杂度：前向滤波 + 后向平滑均为 $O(T)$（标量情形），不需对 $T\times T$ 协方差阵求逆。

11.1.5 缺失值（Missing Values）（PDF p.590–591）

设 $\{y_t\}_{t=\ell+1}^{\ell+h}$ 缺失。保持原时间刻度：$\mu_t=\mu_{\ell+1}+\sum_{j=\ell+1}^{t-1}\eta_j$，故 $E(\mu_t|F_{t-1})=\mu_{\ell+1|\ell}$，$\mathrm{Var}(\mu_t|F_{t-1})=\Sigma_{\ell+1|\ell}+(t-\ell-1)\sigma_\eta^2$，即 $\mu_{t|t-1}=\mu_{t-1|t-2},\qquad \Sigma_{t|t-1}=\Sigma_{t-1|t-2}+\sigma_\eta^2,\quad t=\ell+2,\dots,\ell+h.\tag{11.24}$ 等价于在缺失期取 $v_t=0$、$K_t=0$ 继续运行 (11.14)——没有新数据就没有新息与增益。（实用价值：处理停牌、节假日不对齐、不同频率数据。）

11.1.6 初始化的影响（Effect of Initialization）（PDF p.591–592）

$v_1=y_1-\mu_{1|0}$，$V_1=\Sigma_{1|0}+\sigma_e^2$， $\mu_{2|1}=\mu_{1|0}+\frac{\Sigma_{1|0}}{\Sigma_{1|0}+\sigma_e^2}(y_1-\mu_{1|0})$，$\Sigma_{2|1}=\frac{\Sigma_{1|0}}{\Sigma_{1|0}+\sigma_e^2}\sigma_e^2+\sigma_\eta^2$。令 $\Sigma_{1|0}\to\infty$：$\mu_{2|1}=y_1$，$\Sigma_{2|1}=\sigma_e^2+\sigma_\eta^2$，等价于视 $y_1$ 固定、$\mu_1\sim N(y_1,\sigma_e^2)$。称扩散初始化（diffuse initialization）——初始条件极不确定。
对平滑：$t=T,\dots,2$ 不受影响；对 $\mu_1$：$L_1=\sigma_e^2/V_1$， $\mu_{1|T}=\mu_{1|0}+\frac{\Sigma_{1|0}}{\Sigma_{1|0}+\sigma_e^2}(v_1+\sigma_e^2q_1)\to y_1+\sigma_e^2q_1$； $\Sigma_{1|T}=\frac{\Sigma_{1|0}}{\Sigma_{1|0}+\sigma_e^2}\sigma_e^2-\big(\frac{\Sigma_{1|0}}{\Sigma_{1|0}+\sigma_e^2}\big)^2\sigma_e^4M_1\to\sigma_e^2-\sigma_e^4M_1$。
建议：对 $\mu_1$ 所知甚少时用扩散初始化；若难以接受无穷方差的随机变量，可把 $\mu_1$ 当作额外参数与其他参数联合估计（与第 2、8 章精确 ML 密切相关）。

11.1.7 估计（Estimation）（PDF p.592）

由预测误差分解（prediction error decomposition），正态似然 $p(y_1,\dots,y_T|\sigma_e,\sigma_\eta)=p(y_1)\prod_{t=2}^Tp(v_t|F_{t-1}),\quad y_1\sim N(\mu_{1|0},V_1),\ v_t\sim N(0,V_t),$ $\ln L(\sigma_e,\sigma_\eta)=-\frac T2\ln(2\pi)-\frac12\sum_{t=1}^T\Big(\ln V_t+\frac{v_t^2}{V_t}\Big).\tag{11.25}$ 对数似然（含缺失值情形）可由 Kalman 滤波递推计算，再数值最优化。Matlab、RATS、S-Plus 均支持；本章用 Koopman, Shephard & Doornik (1999) 的 SsfPack（S-Plus 与 OX 中可用，免费）。

11.1.8 所用 S-Plus 命令（PDF p.592–596）

SsfPack 记号（表 11.1）：$\delta$→mDelta，$\Phi$→mPhi，$\Omega$→mOmega，$\Sigma$→mSigma。命令（表 11.2）：SsfFit（ML 估计）、CheckSsf（创建 ssf 对象）、KalmanFil（滤波）、KalmanSmo（平滑）、SsfMomentEst(task="STFIL")（滤波状态及方差）、SsfMomentEst(task="STSMO")（平滑状态及方差）、SsfCondDens(task="STSMO")（平滑状态无方差）。
流程：读 aa-rv-0304.txt，取对数；初值 ltm.start=c(3,1)；P1=-1（−1 表示扩散初始化，$\Sigma_{1|0}$ 很大），a1=0；定义函数返回 list(mPhi=c(1,1), mOmega=diag(c(sigma.eta^2, sigma.e^2)), mSigma=c(P1,a1))；SsfFit(..., lower=c(0,0), upper=c(100,100)) 得 $(\hat\sigma_\eta,\hat\sigma_e)=(0.07350827, 0.48026284)$，mOmega 对角为 0.0054035 与 0.2306524。随后 KalmanFil（输出 mOut、innov、std.innov、mGain、loglike 等）、KalmanSmo、SsfMomentEst 画滤波/平滑状态及 ±2 标准差区间。
模型检验：标准化预测误差 $\tilde v_t=v_t/\sqrt{V_t}$，$Q(25)=23.37(0.56)$，25 阶 ARCH LM 检验 18.48(0.82)，模型充分（archTest、autocorTest）。

11.2 线性状态空间模型（Linear State-Space Models）（PDF p.596–597）

许多经济金融模型可写成状态空间形式：ARIMA、含不可观测成分的动态线性模型、时变回归、随机波动率模型。一般高斯线性状态空间模型： $s_{t+1}=d_t+T_ts_t+R_t\eta_t,\tag{11.26}$ $y_t=c_t+Z_ts_t+e_t,\tag{11.27}$ $s_t$ 为 $m$ 维状态，$y_t$ 为 $k$ 维观测，$d_t,c_t$ 为确定性向量，$T_t$（$m\times m$）、$Z_t$（$k\times m$）为系数阵，$R_t$（$m\times n$）常由单位阵部分列组成；$\eta_t\sim N(0,Q_t)$，$e_t\sim N(0,H_t)$，相互独立（可放宽）；$s_1\sim N(\mu_{1|0},\Sigma_{1|0})$ 与噪声独立。
(11.27) 为测量/观测方程；(11.26) 为状态/转移方程（一阶马尔可夫）。$T_t,R_t,Q_t,Z_t,H_t$ 称系统矩阵（system matrices），常稀疏，可为参数 $\theta$ 的函数，用 ML 估计。
紧凑形式： $\begin{bmatrix}s_{t+1}\\y_t\end{bmatrix}=\delta_t+\Phi_ts_t+u_t,\quad \delta_t=\begin{bmatrix}d_t\\c_t\end{bmatrix},\ \Phi_t=\begin{bmatrix}T_t\\Z_t\end{bmatrix},\ u_t=\begin{bmatrix}R_t\eta_t\\e_t\end{bmatrix},\tag{11.28}$ $\Omega_t=\mathrm{Cov}(u_t)=\mathrm{diag}(R_tQ_tR_t',H_t)$。扩散初始化：$\Sigma_{1|0}=\Sigma_*+\lambda\Sigma_\infty$，$\lambda$ 很大（可趋于无穷）。SsfPack 中用 $\Sigma=\begin{bmatrix}\Sigma_{1|0}\\\mu_{1|0}'\end{bmatrix}_{(m+1)\times m}$。系统矩阵可时不变也可时变。

11.3 模型转换（Model Transformation）（PDF p.597–）

11.3.1 时变系数 CAPM（PDF p.597–599）

模型 $r_t=\alpha_t+\beta_tr_{M,t}+e_t,\ e_t\sim N(0,\sigma_e^2);\quad \alpha_{t+1}=\alpha_t+\eta_t,\ \eta_t\sim N(0,\sigma_\eta^2);\quad \beta_{t+1}=\beta_t+\epsilon_t,\ \epsilon_t\sim N(0,\sigma_\epsilon^2),\tag{11.29}$ $r_t$、$r_{M,t}$ 为资产与市场超额收益，新息互相独立；α、β 为随机游走。状态空间形式：$s_t=(\alpha_t,\beta_t)'$，$T_t=R_t=I_2$，$d_t=c_t=0$，$Z_t=(1,r_{M,t})$，$H_t=\sigma_e^2$，$Q_t=\mathrm{diag}\{\sigma_\eta^2,\sigma_\epsilon^2\}$；紧凑形式 $\delta_t=0$，$u_t=(\eta_t,\epsilon_t,e_t)'$，$\Phi_t=\begin{bmatrix}1&0\\0&1\\1&r_{M,t}\end{bmatrix}$，$\Omega=\mathrm{diag}(\sigma_\eta^2,\sigma_\epsilon^2,\sigma_e^2)$；扩散初始化 $\Sigma=\begin{bmatrix}-1&0\\0&-1\\0&0\end{bmatrix}$。
SsfPack 时变设定：需要数据矩阵 $X$（mX，存 $Z_t$ 的时变值）与索引矩阵（表 11.3：$J_\delta$→mJDelta，$J_\Phi$→mJPhi，$J_\Omega$→mJOmega）。$J_\Phi$ 与 $\Phi_t$ 同维，全设 −1，仅时变元素处填写 $X$ 中对应列号。
示例：GM 1990-01 至 2003-12 月度简单超额收益，S&P 500 超额收益为市场，设 $(\sigma_\eta,\sigma_\epsilon,\sigma_e)=(0.02,0.04,0.1)$：X.mtx=cbind(1,sp)；Phi.t=rbind(diag(2),rep(0,2))；Sigma=-Phi.t；Omega=diag(c(.02^2,.04^2,.1^2))；JPhi=matrix(-1,3,2); JPhi[3,1]=1; JPhi[3,2]=2；组成 list(mPhi, mOmega, mJPhi, mSigma, mX)（mX 为 168×2）。

11.3.2 ARMA 模型（PDF p.600–）

零均值 ARMA($p,q$)：$\phi(B)y_t=\theta(B)a_t$（11.30）。令 $m=\max(p,q+1)$，写成 $y_t=\sum_{i=1}^m\phi_iy_{t-i}+a_t-\sum_{j=1}^{m-1}\theta_ja_{t-j},\tag{11.31}$ $\phi_i=0$（$i>p$），$\theta_j=0$（$j>q$）。三种方法：
Akaike 方法（1975）：状态为产生预测所需的最小变量集合，$s_t=(y_{t|t},y_{t+1|t},\dots,y_{t+m-1|t})'$，$y_{t|t}=y_t$。观测方程 $y_t=Zs_t$，$Z=(1,0,\dots,0)$（11.32）。
- $s_{1,t+1}=y_{t+1}=y_{t+1|t}+a_{t+1}=s_{2t}+a_{t+1}$（11.33）。
- 用 MA 表示 $y_t=\sum\psi_ia_{t-i}$，$\psi_1=\phi_1-\theta_1$，$\psi_2=\phi_1\psi_1+\phi_2-\theta_2$，…，$\psi_{m-1}=\sum_{i=1}^{m-1}\phi_i\psi_{m-1-i}-\theta_{m-1}$（11.34）。
- 预测更新公式：$y_{t+j|t+1}=y_{t+j|t}+\psi_{j-1}a_{t+1}$（$j>0$）（11.35）——新信息体现在新息 $a_{t+1}$，以权重 $\psi_{j-1}$ 修正预测。
- $y_{t+m|t+1}=\sum_{i=1}^m\phi_iy_{t+m-i|t}+\psi_{m-1}a_{t+1}$（11.36）。
- 合并得转移方程 $\begin{bmatrix}y_{t+1}\\y_{t+2|t+1}\\\vdots\\y_{t+m|t+1}\end{bmatrix}=\begin{bmatrix}0&1&0&\cdots&0\\0&0&1&\cdots&0\\\vdots&&&\ddots&\vdots\\0&0&0&\cdots&1\\\phi_m&\phi_{m-1}&\phi_{m-2}&\cdots&\phi_1\end{bmatrix}\begin{bmatrix}y_t\\y_{t+1|t}\\\vdots\\y_{t+m-1|t}\end{bmatrix}+\begin{bmatrix}1\\\psi_1\\\vdots\\\psi_{m-1}\end{bmatrix}a_{t+1},\tag{11.37}$ 即 $s_{t+1}=Ts_t+R\eta_t$，$\eta_t=a_{t+1}\sim N(0,\sigma_a^2)$（11.38）。
Harvey 方法（1993, §4.4）：$s_{1t}=y_t$，其余递推定义：$y_{t+1}=\phi_1s_{1t}+s_{2t}+\eta_t$，$s_{2t}=\sum_{i=2}^m\phi_iy_{t+1-i}-\sum_{j=1}^{m-1}\theta_ja_{t+1-j}$；$s_{2,t+1}=\phi_2s_{1t}+s_{3t}-\theta_1\eta_t$；…；$s_{m,t+1}=\phi_ms_{1t}-\theta_{m-1}\eta_t$。于是 $s_{t+1}=Ts_t+R\eta_t,\quad y_t=Zs_t,\tag{11.39, 11.40}$ $T=\begin{bmatrix}\phi_1&1&0&\cdots&0\\\phi_2&0&1&\cdots&0\\\vdots&&&\ddots&\\\phi_{m-1}&0&0&\cdots&1\\\phi_m&0&0&\cdots&0\end{bmatrix},\quad R=\begin{bmatrix}1\\-\theta_1\\\vdots\\-\theta_{m-1}\end{bmatrix},$ 无测量误差；优点是 AR、MA 系数直接出现在系统矩阵中。
Aoki 方法（1987, 第 4 章）：先看 MA 模型 $y_t=\theta(B)a_t$，取 $s_t=(a_{t-q},\dots,a_{t-1})'$（原文下标印刷有误），转移矩阵为移位矩阵……（续）（接上，PDF p.603–606）MA 模型的 Aoki 形式：
\[s_{t+1}=\begin{bmatrix}0&1&\cdots&0\\\vdots&&\ddots&\vdots\\0&0&\cdots&1\\0&0&\cdots&0\end{bmatrix}s_t+\begin{bmatrix}0\\\vdots\\0\\1\end{bmatrix}a_t,\qquad y_t=(-\theta_q,-\theta_{q-1},\dots,-\theta_1)s_t+a_t,\tag{11.41}\]
$a_t$ 同时出现在状态方程和观测方程中。
- AR 模型 $\phi(B)z_t=a_t$ 的两种 Aoki 形式：(i) $s_t=(z_{t-p+1},\dots,z_t)'$，伴随矩阵（companion matrix）转移，末行为 $(\phi_p,\dots,\phi_1)$，冲击向量 $(0,\dots,0,1)'a_{t+1}$，$z_t=(0,\dots,0,1)s_t$（11.42）；(ii) 把最后一个元素换为 $z_t-a_t$，冲击向量 $(0,\dots,1,\phi_1)'a_t$，$z_t=(0,\dots,0,1)s_t+a_t$（11.43）。
- ARMA($p,q$)（设 $q<p$）：引入辅助变量 $z_t=a_t/\phi(B)$，则 $\phi(B)z_t=a_t$，$y_t=\theta(B)z_t$。用 (11.42) 作转移方程时观测方程 $y_t=(-\theta_{p-1},\dots,-\theta_1,1)s_t$（11.44）；用 (11.43) 时 $y_t=(-\theta_{p-1},\dots,-\theta_1,1)s_t+a_t$（11.45）。
小结：ARMA 有许多状态空间表示，各有优劣，估计和预测时任选其一即可。反之，对时不变状态空间模型，由 Cayley–Hamilton 定理可证观测 $y_t$ 服从 ARMA($m,m$)（$m$ 为状态维数）。
SsfPack：GetSsfArma 用 Harvey 方法。例：AR(1) $y_t=0.6y_{t-1}+a_t$，$a_t\sim N(0,0.4^2)$：mPhi=(0.6, 1)'，mOmega=diag(0.16, 0)，mSigma=(0.25, 0)'——平稳 AR(1) 取 $\Sigma_{1|0}=\mathrm{Var}(y_t)=0.16/(1-0.36)=0.25$，$\mu_{1|0}=0$。ARMA(2,1) $y_t=1.2y_{t-1}-0.35y_{t-2}+a_t-0.25a_{t-1}$，$a_t\sim N(0,1.1^2)$：$T=\begin{bmatrix}1.2&1\\-0.35&0\end{bmatrix}$，$Z=(1,0)$，mOmega 上块 $=\sigma_a^2RR'=\begin{bmatrix}1.21&-0.3025\\-0.3025&0.075625\end{bmatrix}$，mSigma 为 $(s_{1t},s_{2t})$ 的无条件协方差 $\begin{bmatrix}4.0607&-1.4874\\-1.4874&0.5731\end{bmatrix}$，$s_{1t}=y_t$，$s_{2t}=-0.35y_{t-1}-0.25a_t$（原文写作 $-0.25y_{t-2}$，应为 MA 项）。注意：SsfPack 中 MA 多项式为 $\theta(B)=1+\theta_1B+\cdots$，与文献常用的 $1-\theta_1B-\cdots$ 符号相反（故输入 ma=-0.25）。

11.3.3 线性回归模型（PDF p.606–608）

$y_t=x_t'\beta+e_t$，令 $s_t=\beta$（常数）， $\begin{bmatrix}s_{t+1}\\y_t\end{bmatrix}=\begin{bmatrix}I_p\\x_t'\end{bmatrix}s_t+\begin{bmatrix}0_p\\e_t\end{bmatrix},\tag{11.46}$ $T_t=I_p$，$Z_t=x_t'$，$d_t=c_t=0$，$Q_t=0$，$H_t=\sigma_e^2$。状态固定，应用扩散初始化。（此时 Kalman 滤波即递推最小二乘（RLS）。）
推广为随机系数：$\beta_{t+1}=\beta_t+R_t\eta_t$，$\eta_t\sim N(0,1)$，$R_t=(\sigma_1,\dots,\sigma_p)'$，$\sigma_i=0$ 则 $\beta_i$ 不变。
SsfPack：GetSsfReg(cbind(1,sp)) 生成市场模型的状态空间形式：mPhi=[I_2; 0 0]，mOmega=diag(0,0,1)，mSigma 扩散，mJPhi 第三行为 (1,2)，mX 为 168×2。

11.3.4 带 ARMA 误差的线性回归（PDF p.608–609）

$y_t=x_t'\beta+z_t$，$\phi(B)z_t=\theta(B)a_t$（11.47）；$x_t=1$ 时即非零均值 ARMA。设 $s_t$ 为 $z_t$ 的状态（如 Harvey 形式 11.39），定义 $s_t^*=(s_t',\beta_t')'$，$\beta_t=\beta$（11.48）， $s_{t+1}^*=T^*s_t^*+R^*\eta_t,\quad y_t=Z_t^*s_t^*,\tag{11.49, 11.50}$ $Z_t^*=(1,0,\dots,0,x_t')_{1\times(m+k)}$，$T^*=\mathrm{diag}(T,I_k)$，$R^*=(R',0')'$。
SsfPack GetSsfRegArma(X, ar=c(1.2,-0.35), ma=c(-0.25))：例 $y_t=\beta_0+\beta_1x_t+z_t$，$z_t=1.2z_{t-1}-0.35z_{t-2}+a_t-0.25a_{t-1}$，输出 mPhi（5×4，前两列为 ARMA 部分，后两列为单位阵，末行 (1,0,0,0) 并通过 mJPhi 第 5 行 (−1,−1,1,2) 指向 X 的两列）、mOmega（左上 $\begin{bmatrix}1&-0.25\\-0.25&0.0625\end{bmatrix}$）、mSigma（ARMA 部分无条件协方差 $\begin{bmatrix}3.356&-1.229\\-1.229&0.4736\end{bmatrix}$，回归系数部分 −1 扩散）。

11.3.5 标量不可观测成分模型（Scalar Unobserved Component Model）（PDF p.610–611）

结构时间序列模型（STSM）： $y_t=\mu_t+\gamma_t+\omega_t+e_t,\tag{11.51}$ 趋势 $\mu_t$、季节 $\gamma_t$、周期 $\omega_t$、不规则成分 $e_t$。
趋势（可能双单位根）： $\mu_{t+1}=\mu_t+\beta_t+\eta_t,\ \eta_t\sim N(0,\sigma_\eta^2);\quad \beta_t=\beta_{t-1}+\varsigma_t,\ \varsigma_t\sim N(0,\sigma_\varsigma^2),\tag{11.52}$ $\mu_1,\beta_1\sim N(0,\xi)$，$\xi$ 很大（如 $10^8$，Kitagawa & Gersch 1996）。$\sigma_\varsigma=0$：带漂移 $\beta_1$ 的随机游走；$\sigma_\varsigma=\sigma_\eta=0$：线性确定性趋势。
季节：$(1+B+\cdots+B^{s-1})\gamma_t=\omega_t$，$\omega_t\sim N(0,\sigma_\omega^2)$（11.53），$s$ 为季节周期；$\sigma_\omega=0$ 时季节确定性。
周期： $\begin{bmatrix}\omega_{t+1}\\\omega_{t+1}^*\end{bmatrix}=\delta\begin{bmatrix}\cos\lambda_c&\sin\lambda_c\\-\sin\lambda_c&\cos\lambda_c\end{bmatrix}\begin{bmatrix}\omega_t\\\omega_t^*\end{bmatrix}+\begin{bmatrix}\varepsilon_t\\\varepsilon_t^*\end{bmatrix},\tag{11.54}$ $(\varepsilon_t,\varepsilon_t^*)'\sim N(0,\sigma_\varepsilon^2(1-\delta^2)I_2)$，$\omega_0,\omega_0^*\sim N(0,\sigma_\varepsilon^2)$ 独立，$\delta\in(0,1]$ 为阻尼因子，$\lambda_c=2\pi/q$（$q$ 为周期）；$\delta=1$ 为确定性正余弦波。
SsfPack GetSsfStsm，最多 10 个周期成分；参数对应（表 11.4）：irregular→$\sigma_e$，level→$\sigma_\eta$，slope→$\sigma_\varsigma$，seasonalDummy/seasonalTrig/seasonalHS→$\sigma_\omega,s$，Cycle0…Cycle9→$\sigma_\varepsilon,\lambda_c,\delta$。例：局部趋势模型 $\sigma_e=0.4,\sigma_\eta=0.2$：GetSsfStsm(irregular=0.4, level=0.2)，mPhi=(1,1)'，mOmega=diag(0.04, 0.16)，mSigma=(−1, 0)'。

11.4 Kalman 滤波与平滑（Kalman Filter and Smoothing）（PDF p.611–620）

推导沿 11.1 节思路，参考 Durbin & Koopman (2001, 第 4 章)；应用型读者可跳过。

11.4.1 Kalman 滤波（PDF p.611–614）

记 $s_j|F_i\sim N(s_{j|i},\Sigma_{j|i})$。由 (11.26)： $s_{t+1|t}=d_t+T_ts_{t|t},\tag{11.55}\qquad \Sigma_{t+1|t}=T_t\Sigma_{t|t}T_t'+R_tQ_tR_t'.\tag{11.56}$
$y_{t|t-1}=c_t+Z_ts_{t|t-1}$，预测误差 $v_t=y_t-c_t-Z_ts_{t|t-1}=Z_t(s_t-s_{t|t-1})+e_t,\tag{11.57}$ 性质：$E(v_t|F_{t-1})=0$；与 $F_{t-1}$ 独立；$\{v_t\}$ 为独立正态向量序列。$V_t=Z_t\Sigma_{t|t-1}Z_t'+H_t$（11.58）。
由定理 11.1：$s_{t|t}=s_{t|t-1}+C_tV_t^{-1}v_t$（11.59），$C_t=\mathrm{Cov}(s_t,v_t|F_{t-1})=\Sigma_{t|t-1}Z_t'$（假设 $H_t$ 可逆从而 $V_t$ 可逆）。
合并：$s_{t+1|t}=d_t+T_ts_{t|t-1}+K_tv_t$（11.60），Kalman 增益 $K_t=T_t\Sigma_{t|t-1}Z_t'V_t^{-1}$（11.61）；$\Sigma_{t|t}=\Sigma_{t|t-1}-\Sigma_{t|t-1}Z_t'V_t^{-1}Z_t\Sigma_{t|t-1}$（11.62）；$\Sigma_{t+1|t}=T_t\Sigma_{t|t-1}L_t'+R_tQ_tR_t'$，$L_t=T_t-K_tZ_t$（11.63）。
一般 Kalman 滤波算法（给定 $s_{1|0},\Sigma_{1|0}$）： $\begin{aligned}v_t&=y_t-c_t-Z_ts_{t|t-1},\quad V_t=Z_t\Sigma_{t|t-1}Z_t'+H_t,\quad K_t=T_t\Sigma_{t|t-1}Z_t'V_t^{-1},\\ L_t&=T_t-K_tZ_t,\quad s_{t+1|t}=d_t+T_ts_{t|t-1}+K_tv_t,\quad \Sigma_{t+1|t}=T_t\Sigma_{t|t-1}L_t'+R_tQ_tR_t'.\end{aligned}\tag{11.64}$
含同期滤波量的版本：$v_t$；$C_t=\Sigma_{t|t-1}Z_t'$；$V_t=Z_tC_t+H_t$；$s_{t|t}=s_{t|t-1}+C_tV_t^{-1}v_t$；$\Sigma_{t|t}=\Sigma_{t|t-1}-C_tV_t^{-1}C_t'$；$s_{t+1|t}=d_t+T_ts_{t|t}$；$\Sigma_{t+1|t}=T_t\Sigma_{t|t}T_t'+R_tQ_tR_t'$。
稳态（steady state）：时不变模型中 $\Sigma_{t|t-1}\to\Sigma^*$，满足离散代数 Riccati 方程 $\Sigma^*=T\Sigma^*T'-T\Sigma^*Z'V^{-1}Z\Sigma^*T'+RQR',\quad V=Z\Sigma^*Z'+H.$ 达到稳态后 $V_t,K_t,\Sigma_{t+1|t}$ 均为常数，可大幅节省计算。
复杂度：每步主要为 $m\times m$ 矩阵乘法与 $k\times k$ 求逆，$O(T(m^3+k^3))$。

11.4.2 状态估计误差与预测误差（PDF p.614）

$x_t=s_t-s_{t|t-1}$，$\mathrm{Var}(x_t|F_{t-1})=\Sigma_{t|t-1}$；$v_t=Z_tx_t+e_t$；$x_{t+1}=T_tx_t+R_t\eta_t-K_tv_t=L_tx_t+R_t\eta_t-K_te_t$。即 $v_t=Z_tx_t+e_t,\qquad x_{t+1}=L_tx_t+R_t\eta_t-K_te_t,\tag{11.65}$ $x_1=s_1-s_{1|0}$。同样可证 $\{v_t\}$ 相互独立，且 $\{v_t,\dots,v_T\}$ 与 $F_{t-1}$ 独立。

11.4.3 状态平滑（PDF p.615–617）

$s_{t|T}=s_{t|t-1}+\sum_{j=t}^T\mathrm{Cov}(s_t,v_j)V_j^{-1}v_j$（11.66），$\mathrm{Cov}(s_t,v_j)=E(s_tx_j')Z_j'$（11.67），$E(s_tx_t')=\Sigma_{t|t-1}$，$E(s_tx_{t+1}')=\Sigma_{t|t-1}L_t'$，……，$E(s_tx_T')=\Sigma_{t|t-1}L_t'\cdots L_{T-1}'$（11.68）。
得 $s_{t|T}=s_{t|t-1}+\Sigma_{t|t-1}q_{t-1}$（11.69），$q_{t-1}=Z_t'V_t^{-1}v_t+L_t'Z_{t+1}'V_{t+1}^{-1}v_{t+1}+\cdots+L_t'\cdots L_{T-1}'Z_T'V_T^{-1}v_T$，后向递推 $q_{t-1}=Z_t'V_t^{-1}v_t+L_t'q_t$，$q_T=0$（11.70）。
固定区间平滑器（fixed interval smoother）（de Jong 1989）： $q_{t-1}=Z_t'V_t^{-1}v_t+L_t'q_t,\quad s_{t|T}=s_{t|t-1}+\Sigma_{t|t-1}q_{t-1},\quad t=T,\dots,1.\tag{11.71}$
平滑状态协方差：$\Sigma_{t|T}=\Sigma_{t|t-1}-\Sigma_{t|t-1}M_{t-1}\Sigma_{t|t-1}$，$M_{t-1}=Z_t'V_t^{-1}Z_t+L_t'M_tL_t$（11.72），$M_T=0$，$M_t=\mathrm{Var}(q_t)$（11.73）。
合并的后向递推： $\begin{aligned}q_{t-1}&=Z_t'V_t^{-1}v_t+L_t'q_t,& s_{t|T}&=s_{t|t-1}+\Sigma_{t|t-1}q_{t-1},\\ M_{t-1}&=Z_t'V_t^{-1}Z_t+L_t'M_tL_t,& \Sigma_{t|T}&=\Sigma_{t|t-1}-\Sigma_{t|t-1}M_{t-1}\Sigma_{t|t-1},\end{aligned}\quad t=T,\dots,1.\tag{11.74}$
实施两步：先前向 Kalman 滤波 (11.64) 并存储 $v_t,V_t,K_t,s_{t|t-1},\Sigma_{t|t-1}$；再后向用 (11.74)。

11.4.4 扰动平滑（Disturbance Smoothing）（PDF p.617–620）

平滑扰动 $e_{t|T}=E(e_t|F_T)$，$\eta_{t|T}=E(\eta_t|F_T)$，用于模型检验（异常值、结构突变诊断）等。
$e_{t|T}=\sum_{j\ge t}E(e_tv_j')V_j^{-1}v_j$（11.75）；$E(e_tv_t')=H_t$，$E(e_tv_j')=E(e_tx_j')Z_j'$（$j>t$）（11.76）；$E(e_tx_{t+1}')=-H_tK_t'$，$E(e_tx_{t+2}')=-H_tK_t'L_{t+1}'$，…（11.77）。于是 $e_{t|T}=H_t(V_t^{-1}v_t-K_t'q_t)=H_to_t,\tag{11.78}$ $o_t=V_t^{-1}v_t-K_t'q_t$ 称平滑测量误差（smoothing measurement error）。
$\eta_{t|T}=\sum_jE(\eta_tv_j')V_j^{-1}v_j$（11.79），$E(\eta_tv_{t+1}')=Q_tR_t'Z_{t+1}'$，$E(\eta_tx_{t+2}')=Q_tR_t'L_{t+1}'$，…，得 $\eta_{t|T}=Q_tR_t'q_t.\tag{11.80}$
Koopman (1993) 由此给出平滑状态的新前向递推：$s_{t+1|T}=d_t+T_ts_{t|T}+R_tQ_tR_t'q_t$（11.81），初值 $s_{1|T}=s_{1|0}+\Sigma_{1|0}q_0$（只需存 $q_t$，省内存）。
协方差：$\mathrm{Var}(e_t|F_T)=H_t-H_t(V_t^{-1}+K_t'M_tK_t)H_t=H_t-H_tN_tH_t$，$N_t=V_t^{-1}+K_t'M_tK_t$；$\mathrm{Var}(\eta_t|F_T)=Q_t-Q_tR_t'M_tR_tQ_t$。
扰动平滑算法汇总（$t=T,\dots,1$，$q_T=0,M_T=0$）： $\begin{aligned}e_{t|T}&=H_t(V_t^{-1}v_t-K_t'q_t),\quad \eta_{t|T}=Q_tR_t'q_t,\quad q_{t-1}=Z_t'V_t^{-1}v_t+L_t'q_t,\\ \mathrm{Var}(e_t|F_T)&=H_t-H_t(V_t^{-1}+K_t'M_tK_t)H_t,\quad \mathrm{Var}(\eta_t|F_T)=Q_t-Q_tR_t'M_tR_tQ_t,\\ M_{t-1}&=Z_t'V_t^{-1}Z_t+L_t'M_tL_t.\end{aligned}\tag{11.82}$

11.5 缺失值（Missing Values）（PDF p.620）

情形一：$y_t$ 在 $t=\ell+1,\dots,\ell+h$ 整体缺失：设 $v_t=0$，$K_t=0$，滤波照常：$s_{t+1|t}=d_t+T_ts_{t|t-1}$，$\Sigma_{t+1|t}=T_t\Sigma_{t|t-1}T_t'+R_tQ_tR_t'$；平滑中 $q_{t-1}=T_t'q_t$，$M_{t-1}=T_t'M_tT_t$。
情形二：$y_t$ 部分分量缺失：$y_t^*=Jy_t$，$J$ 的行是 $I_k$ 的部分行（选择矩阵），观测方程变为 $y_t^*=c_t^*+Z_t^*s_t+e_t^*$，$c_t^*=Jc_t$，$Z_t^*=JZ_t$，$H_t^*=JH_tJ'$，滤波与平滑在该时点用修改后的方程即可。易于处理缺失值是状态空间模型的一大优点。

11.6 预测（Forecasting）（PDF p.621–622）

预测原点 $t$，最小均方误差预测 $y_t(j)=E(y_{t+j}|F_t)$，$j=1,\dots,h$。可通过把 $y_{t+1},\dots,y_{t+h}$ 视为缺失值、用 Kalman 滤波得到。
1 步：$y_t(1)=c_{t+1}+Z_{t+1}s_{t+1|t}$，误差 $e_t(1)=Z_{t+1}(s_{t+1}-s_{t+1|t})+e_{t+1}$，$\mathrm{Var}[e_t(1)]=Z_{t+1}\Sigma_{t+1|t}Z_{t+1}'+H_{t+1}=V_{t+1}$。
$j$ 步： $y_t(j)=c_{t+j}+Z_{t+j}s_{t+j|t},\tag{11.83}\qquad \mathrm{Var}[e_t(j)]=Z_{t+j}\Sigma_{t+j|t}Z_{t+j}'+H_{t+j},\tag{11.84}$ $s_{t+j+1|t}=d_{t+j}+T_{t+j}s_{t+j|t},\qquad \Sigma_{t+j+1|t}=T_{t+j}\Sigma_{t+j|t}T_{t+j}'+R_{t+j}Q_{t+j}R_{t+j}'.\tag{11.85}$ 正是 $v_{t+j}=0$、$K_{t+j}=0$ 时的 (11.64)。
预测误差 $\{v_t\}$ 用于似然估计；标准化预测误差 $D_t^{-1/2}v_t$（$D_t=\mathrm{diag}\{V_t(1,1),\dots,V_t(k,k)\}$）用于模型检验。

11.7 应用（Application）（PDF p.622–）

例 11.2（GM 月度简单超额收益 %，1990-01 至 2003-12，S&P 500 超额收益为市场）：
- 固定系数市场模型 $r_t=\alpha+\beta r_{M,t}+e_t$（11.86），OLS：截距 0.1982（标准误 0.6302，t=0.31），斜率 1.0457（0.1453，t=7.20）；$R^2$=0.2378，调整 $R^2$=0.2332，DW=2.029，Jarque–Bera 2.53(0.28)，Ljung–Box 24.21(0.34)，残差标准误 8.13（166 自由度）。拟合 $r_t=0.20+1.0457r_{M,t}+e_t$，$\hat\sigma_e=8.13$，模型充分。
- SsfPack 估计同一模型：GetSsfReg(mX)，把 mOmega[3,3] 设为 $\exp(\text{parm})$（以对数参数化保证为正），SsfFit(c.start=10, gm, "reg.m", mX=X.mtx)，$\sqrt{\exp(\hat\theta)}=8.130114$；SsfMomentEst(task="STSMO") 得平滑状态 (0.1982025, 1.045702)，标准差 (0.6302091, 0.1453139)——与 OLS 完全一致。
- 时变 CAPM（11.3.1 节）：函数 tv.capm 以对数参数化三个方差，SsfFit(tv.start=c(0,0,0), ...) 得 $(\hat\sigma_\eta,\hat\sigma_\epsilon,\hat\sigma_e)=(4.91\times10^{-5},\ 0.0122,\ 8.125)$；用 SsfCondDens 计算平滑状态与平滑响应（不含方差）……（续）（PDF p.624–625）$\hat\sigma_\eta=4.91\times10^{-5}$ 与 $\hat\sigma_\epsilon=1.22\times10^{-2}$ 都接近 0，说明 GM 的 $\alpha_t,\beta_t$ 基本恒定，与固定系数模型拟合良好一致。图 11.5：(a) 超额收益；(b) 期望收益 $r_{t|T}$；(c) $\alpha_t$ 估计（约 0.206，几乎不变）；(d) $\beta_t$ 估计（在 1.02–1.06 间小幅变动）。纵轴刻度很窄，确认固定系数模型足够。
例 11.3（Johnson & Johnson 1960–1980 季度 EPS 对数，第 2 章数据）：不可观测成分模型
\[y_t=\mu_t+\gamma_t+e_t,\ e_t\sim N(0,\sigma_e^2);\quad \mu_{t+1}=\mu_t+\eta_t;\quad (1+B+B^2+B^3)\gamma_t=\omega_t,\tag{11.87}\]
即 $\gamma_t=-\sum_{j=1}^3\gamma_{t-j}+\omega_t$，三个参数 $\sigma_e,\sigma_\eta,\sigma_\omega$。状态空间形式：
\[\begin{bmatrix}\mu_{t+1}\\\gamma_{t+1}\\\gamma_t\\\gamma_{t-1}\end{bmatrix}=\begin{bmatrix}1&0&0&0\\0&-1&-1&-1\\0&1&0&0\\0&0&1&0\end{bmatrix}\begin{bmatrix}\mu_t\\\gamma_t\\\gamma_{t-1}\\\gamma_{t-2}\end{bmatrix}+\begin{bmatrix}1&0\\0&1\\0&0\\0&0\end{bmatrix}\begin{bmatrix}\eta_t\\\omega_t\end{bmatrix},\quad y_t=[1,1,0,0]s_t+e_t,\]
$\mathrm{Cov}(\eta_t,\omega_t)=\mathrm{diag}(\sigma_\eta^2,\sigma_\omega^2)$。用 GetSsfStsm(irregular, level, seasonalDummy=c(σω,4))，对数参数化 SsfFit，得 $(\hat\sigma_e,\hat\sigma_\eta,\hat\sigma_\omega)=(2.04\times10^{-6},\ 7.27\times10^{-2},\ 2.93\times10^{-2})$；mOmega 对角 0.005285、0.000859、…、$4.18\times10^{-12}$。
- SsfMomentEst(task="STSMO") 得平滑趋势 $\mu_{t|T}$ 与季节 $\gamma_{t|T}$（$T=84$）及 95% 区间（图 11.6）：季节模式随时间演变。图 11.7：(a) Kalman 1 步预测误差，(b) 平滑响应残差（KalmanSmo(...)$response.residuals）。
- 注意：成分分解不唯一，依赖模型设定和约束（如改用 seasonalTrig 会得到另一分解），解释需谨慎；但只要是有效分解，对预测无影响。

第 11 章习题（PDF p.629–631）

11.1：ARMA(1,1) $y_t-0.8y_{t-1}=a_t+0.4a_{t-1}$，$a_t\sim N(0,0.49)$，分别用 Akaike、Harvey、Aoki 方法写成状态空间形式。
11.2：Alcoa 20 分钟收益构造的已实现波动率：拟合 ARIMA(0,1,1)；估计局部趋势模型的 $\sigma_e,\sigma_\eta$，画滤波与平滑状态及 95% 区间。
11.3：Pfizer 与 S&P 500 1990–2003 月度超额收益：固定系数市场模型；时变 CAPM，估计 $\alpha_t,\beta_t$ 新息标准差并画平滑估计。
11.4：AR(3) 加测量误差 $y_t=x_t+e_t$ 写成状态空间；若 $E(e_t)=c\ne0$ 如何修改。
11.5：美国 PPI（1947-01 至 2009-11）对数差分去均值后 AR(3)；假定有独立测量误差，用状态空间估计状态新息方差与 $\sigma_e^2$，画平滑 $x_t$ 与滤波响应残差。

参考文献（PDF p.631）：Akaike (1975)、Aoki (1987)、Anderson & Moore (1979)、de Jong (1989)、Durbin & Koopman (2001)、Hamilton (1994)、Harvey (1993)、Kalman (1960)、Kim & Nelson (1999)、Koopman (1993)、Koopman, Shephard & Doornik (1999)、Shumway & Stoffer (2000)、West & Harrison (1997)。PDF p.632 空白。

第 11 章本章要点

线性高斯状态空间模型 = 状态方程 $s_{t+1}=d_t+T_ts_t+R_t\eta_t$ + 观测方程 $y_t=c_t+Z_ts_t+e_t$；ARIMA、时变系数回归、不可观测成分（趋势/季节/周期）、带 ARMA 误差的回归都可纳入。
局部趋势模型等价于 ARIMA(0,1,1)（简单指数平滑），$\theta$ 与 $\sigma_e^2,\sigma_\eta^2$ 一一对应（$\theta>0$ 时）。
Kalman 滤波来自多元正态条件分布公式：预测–更新两步，增益 $K_t$ 是状态对新息的回归系数；新息 $v_t$ 相互独立，给出协方差阵的 Cholesky 分解，从而似然可写为 $-\frac12\sum(\ln|V_t|+v_t'V_t^{-1}v_t)$。
平滑（固定区间平滑器、扰动平滑）用后向递推 $q_{t-1},M_{t-1}$；平滑估计比滤波更平滑、区间更窄。
缺失值：令 $v_t=0,K_t=0$；部分缺失用选择矩阵 $J$；预测 = 把未来视为缺失值。
扩散初始化处理未知初值；时不变模型存在稳态解（Riccati 方程）。
ARMA 的状态空间表示不唯一（Akaike、Harvey、Aoki）；SsfPack 的 MA 符号约定与文献相反。

与量化交易的关联

时变 β 与动态对冲：时变 CAPM（随机游走 α、β）+ Kalman 滤波是估计时变市场 β、因子暴露的标准工具；配对交易中用 Kalman 滤波在线估计动态对冲比 $\gamma_t$（把 8.8 节的静态协整回归改为随机游走系数回归）是业界常见做法。滤波估计（只用历史信息）可用于回测，平滑估计会引入未来信息，回测中不能使用——这是常见的前视偏差陷阱。
波动率估计：已实现波动率 + 局部趋势模型分离微观结构噪声与真实波动，用于日内波动预测和期权定价输入。
信号去噪与趋势提取：局部水平/局部线性趋势模型等价于自适应指数平滑，可用于价格趋势信号、均线类策略的理论基础；信噪比 $\sigma_\eta/\sigma_e$ 决定平滑程度。
数据工程：缺失值处理（停牌、节假日、不同市场交易日历不一致）、混频数据、不规则采样都可在状态空间框架下统一处理。
宏观与基本面：不可观测成分模型分解季度盈利的趋势与季节，用于盈利预测与季节调整；动态因子模型、利率期限结构（动态 Nelson–Siegel）均是状态空间模型。
估计：任何带潜变量的线性高斯模型都可用预测误差分解构造精确似然。

第 12 章马尔可夫链蒙特卡罗方法及其应用（Markov Chain Monte Carlo Methods with Applications）（PDF p.633 起）

本章引言（PDF p.633–634）

计算能力与方法进步使 MCMC 与**数据增广（data augmentation）**成为可能，本章介绍其思想，讨论基于 Gibbs 抽样的贝叶斯推断及在金融计量中的应用（参考 Carlin & Louis 2000；Gelman, Carlin, Stern & Rubin 2003）。作者认为贝叶斯推断与 MCMC 几乎适用于金融计量的所有问题。
马尔可夫过程：$\{X_t\}$ 取值于空间 $\Omega$，给定 $X_t$ 时未来 $X_h$（$h>t$）与过去 $X_s$（$s<t$）无关：$P(X_h|X_s,s\le t)=P(X_h|X_t)$；离散时间 $P(X_h|X_t,X_{t-1},\dots)=P(X_h|X_t)$。转移概率函数 $P_t(\theta,h,A)=P(X_h\in A|X_t=\theta)$；若只依赖 $h-t$ 而不依赖 $t$，称有平稳转移分布。

12.1 马尔可夫链模拟（Markov Chain Simulation）（PDF p.634–635）

推断需要后验 $P(\theta|X)$。思想：在参数空间 $\Theta$ 上构造一个马尔可夫过程，使其平稳分布为 $P(\theta|X)$，运行足够长时间使当前值分布足够接近平稳分布。给定目标分布，可构造许多满足条件的链；用此法得到 $P(\theta|X)$ 的方法统称 MCMC。
历史脉络：EM 算法（Dempster, Laird & Rubin 1977）处理缺失值——M 步：若缺失值已知，用完全数据方法做 ML 估计；E 步：给定数据与拟合模型，求缺失值的条件期望并填补；从任意值开始反复迭代至收敛。Tanner & Wong (1987) 两方面推广：(1) 迭代模拟——用从条件分布的随机抽取代替条件期望；(2) 数据增广——加入辅助变量，常可简化或加速模拟。

12.2 Gibbs 抽样（Gibbs Sampling）（PDF p.635–637）

Geman & Geman (1984)、Gelfand & Smith (1990)，最流行的 MCMC 方法。"参数"含义广：缺失数据点可视为参数；不可观测量（如 $N$ 笔成交对应的"真实"价格）可视为 $N$ 个参数——与数据增广相关。
设三个参数 $\theta_1,\theta_2,\theta_3$，数据 $X$，模型 $M$。似然难求，但三个完全条件分布（full conditional distributions） $f_1(\theta_1|\theta_2,\theta_3,X,M),\quad f_2(\theta_2|\theta_3,\theta_1,X,M),\quad f_3(\theta_3|\theta_1,\theta_2,X,M)\tag{12.1}$ 可抽样（不必知道确切形式，只需能抽随机数）。
算法：给定任意初值 $\theta_{2,0},\theta_{3,0}$：(1) 从 $f_1(\theta_1|\theta_{2,0},\theta_{3,0})$ 抽 $\theta_{1,1}$；(2) 从 $f_2(\theta_2|\theta_{3,0},\theta_{1,1})$ 抽 $\theta_{2,1}$；(3) 从 $f_3(\theta_3|\theta_{1,1},\theta_{2,1})$ 抽 $\theta_{3,1}$。完成一次 Gibbs 迭代；重复 $m$ 次。正则条件下（实质要求从任意初值出发能访问整个参数空间；理论见 Tierney 1994），$m$ 足够大时 $(\theta_{1,m},\theta_{2,m},\theta_{3,m})$ 近似为联合后验 $f(\theta_1,\theta_2,\theta_3|X,M)$ 的一次抽样。
实践：运行 $n$ 次，丢弃前 $m$ 次，得 Gibbs 样本 $\{(\theta_{1,j},\theta_{2,j},\theta_{3,j})\}_{j=m+1}^n$（12.2）。点估计与方差 $\hat\theta_i=\frac{1}{n-m}\sum_{j=m+1}^n\theta_{i,j},\qquad \hat\sigma_i^2=\frac{1}{n-m-1}\sum_{j=m+1}^n(\theta_{i,j}-\hat\theta_i)^2.\tag{12.3}$ 检验 $H_0:\theta_1=\theta_2$：对 $\theta_{1,j}-\theta_{2,j}$ 求均值与方差，用 $t=\hat\theta/\hat\sigma$。
备注：被丢弃的前 $m$ 个称 burn-in 样本，保证剩余样本接近联合分布的随机样本。另一做法：用不同初值运行许多较短的链，取每条链最后一次抽样组成 Gibbs 样本。
优点：通过完全条件分布把高维估计分解为低维问题（极端情形 $N$ 个一元条件分布）。但参数高度相关时应联合抽样（如 $\theta_1,\theta_2$ 高度相关，用 $f(\theta_1,\theta_2|\theta_3)$ 和 $f_3(\theta_3|\theta_1,\theta_2)$），否则收敛慢（Liu, Wong & Kong 1994，即"分块 Gibbs"）。
收敛诊断：理论只说 $m$ 足够大时收敛，无具体指导；诊断方法很多但无共识，没有方法能 100% 保证收敛（Carlin & Louis 2000；Gelman et al. 2003）。实践中应用不同初值重复多次确认收敛（可补充：Gelman–Rubin $\hat R$、迹图、自相关与有效样本量）。

12.3 贝叶斯推断（Bayesian Inference）（PDF p.637–641）

完全条件分布在文献中称条件后验分布。

12.3.1 后验分布（Posterior Distributions）（PDF p.637–638）

两种推断范式：经典（极大似然）与贝叶斯（先验 + 数据 → 后验）。本书前面都是经典方法，但所有问题都有贝叶斯解法，借助 MCMC 已可行，多数情况下结果相近；有时贝叶斯更优，例如 VaR 计算中可自然纳入参数不确定性，代价是计算量大。
先验 $P(\theta)$，似然 $f(X|\theta)$： $f(\theta|X)=\frac{f(X|\theta)P(\theta)}{f(X)},\quad f(X)=\int f(X|\theta)P(\theta)d\theta,\tag{12.4}$ $f(\theta|X)\propto f(X|\theta)P(\theta).\tag{12.5}$ 基于似然的推断等价于取常数先验的贝叶斯方法。

12.3.2 共轭先验（Conjugate Prior Distributions）（PDF p.638–641）

先验与后验同属一族即共轭先验；在 MCMC 中意味着条件后验有闭式，可直接用常规随机数生成器抽样（DeGroot 1970 第 9 章）。

Result 12.1（正态均值，方差已知）：$x_i\sim N(\mu,\sigma^2)$，$\mu\sim N(\mu_o,\sigma_o^2)$，则后验 $N(\mu_*,\sigma_*^2)$： $\mu_*=\frac{\sigma^2\mu_o+n\sigma_o^2\bar x}{\sigma^2+n\sigma_o^2},\qquad \sigma_*^2=\frac{\sigma^2\sigma_o^2}{\sigma^2+n\sigma_o^2}.$ 用精度 $\eta=1/\sigma^2$：$\eta_*=\eta_o+n\eta$，$\mu_*=\frac{\eta_o}{\eta_*}\mu_o+\frac{n\eta}{\eta_*}\bar x$。即后验精度 = 先验精度 + 数据精度（$\bar x$ 为充分统计量，精度 $n\eta$）；后验均值是按精度加权的平均；$n$ 增大时先验影响递减。
Result 12.1a（多元）：$x_i\sim N(\mu,\Sigma)$（$\Sigma$ 已知），$\mu\sim N(\mu_o,\Sigma_o)$，后验 $N(\mu_*,\Sigma_*)$： $\Sigma_*^{-1}=\Sigma_o^{-1}+n\Sigma^{-1},\qquad \mu_*=\Sigma_*(\Sigma_o^{-1}\mu_o+n\Sigma^{-1}\bar x).$ （回归模型 MCMC 中很有用，Box & Tiao 1973。）
Gamma 分布 $f(\eta|\alpha,\beta)=\frac{\beta^\alpha}{\Gamma(\alpha)}\eta^{\alpha-1}e^{-\beta\eta}$，$E=\alpha/\beta$，$\mathrm{Var}=\alpha/\beta^2$。Result 12.2（正态精度，均值已知）：$\eta\sim\Gamma(\alpha,\beta)$，后验 $\Gamma\big(\alpha+n/2,\ \beta+\sum(x_i-\mu)^2/2\big)$。
Beta 分布 $f(\theta|\alpha,\beta)=\frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\theta^{\alpha-1}(1-\theta)^{\beta-1}$，$E=\alpha/(\alpha+\beta)$，$\mathrm{Var}=\frac{\alpha\beta}{(\alpha+\beta)^2(\alpha+\beta+1)}$。Result 12.3（Bernoulli）：Beta($\alpha,\beta$) 先验，后验 Beta($\alpha+\sum x_i$, $\beta+n-\sum x_i$)。
Result 12.4（Poisson）：$\lambda\sim\Gamma(\alpha,\beta)$，后验 $\Gamma(\alpha+\sum x_i,\ \beta+n)$。
Result 12.5（指数分布）：$\lambda\sim\Gamma(\alpha,\beta)$，后验 $\Gamma(\alpha+n,\ \beta+\sum x_i)$。
负二项分布 $p(n|m,\lambda)=\binom{m+n-1}{n}\lambda^m(1-\lambda)^n$；金融例子：公司有 $m$ 个职位，每个 MBA 应聘者独立以概率 $\lambda$ 合适，面试总数 $Y$，则 $X=Y-m$ 为负二项。Result 12.6：$m$ 固定，$\lambda\sim$ Beta($\alpha,\beta$)，后验 Beta($\alpha+mn$, $\beta+\sum x_i$)。
Result 12.7（正态，均值与精度都未知，正态–Gamma 先验）：$\mu|\eta\sim N(\mu_o,(\tau_o\eta)^{-1})$，$\eta\sim\Gamma(\alpha,\beta)$。则 $\mu|\eta$ 后验为正态，均值 $\mu_*=\frac{\tau_o\mu_o+n\bar x}{\tau_o+n}$，精度 $(\tau_o+n)\eta$；$\eta$ 边际后验 $\Gamma(\alpha+n/2,\beta_*)$， $\beta_*=\beta+\frac12\sum_{i=1}^n(x_i-\bar x)^2+\frac{\tau_on(\bar x-\mu_o)^2}{2(\tau_o+n)}.$
逆卡方分布：$1/Y\sim\chi^2_v$，密度 $f(y|v)=\frac{2^{-v/2}}{\Gamma(v/2)}y^{-(v/2+1)}e^{-1/(2y)}$，$E(Y)=1/(v-2)$（$v>2$），$\mathrm{Var}(Y)=\frac{2}{(v-2)^2(v-4)}$（$v>4$）。Result 12.8（零均值正态的方差）：先验 $v\lambda/\sigma^2\sim\chi^2_v$，后验 $(v\lambda+\sum a_i^2)/\sigma^2\sim\chi^2_{v+n}$。

12.4 替代算法（Alternative Algorithms）（PDF p.642–644）

条件后验无闭式时的抽样方法。

12.4.1 Metropolis 算法（PDF p.642–643）

适用于条件后验已知到一个归一化常数的情形（Metropolis & Ulam 1949；Metropolis et al. 1953）。步骤：
1. 取初值 $\theta_0$ 使 $f(\theta_0|X)>0$；
2. 对 $t=1,2,\dots$：(a) 从跳跃分布/建议分布（jumping / proposal distribution） $J_t(\theta^*|\theta_{t-1})$ 抽候选 $\theta^*$，$J_t$ 必须对称：$J_t(\theta_i|\theta_j)=J_t(\theta_j|\theta_i)$；(b) 计算 $r=f(\theta^*|X)/f(\theta_{t-1}|X)$；(c) 以概率 $\min(r,1)$ 令 $\theta_t=\theta^*$，否则 $\theta_t=\theta_{t-1}$。
正则条件下 $\{\theta_t\}$ 依分布收敛到 $f(\theta|X)$。只用比值，不需归一化常数；需能计算 $r$、从跳跃分布抽样、抽均匀随机数决定接受与否。直观：跳到后验密度更高处必接受；更低处以概率 $r$ 接受。
对称跳跃分布例：以当前值为中心的正态或 Student-t（随机游走 Metropolis）。

12.4.2 Metropolis–Hastings 算法（PDF p.643）

Hastings (1970) 推广：跳跃分布不必对称，接受比改为 $r=\frac{f(\theta^*|X)/J_t(\theta^*|\theta_{t-1})}{f(\theta_{t-1}|X)/J_t(\theta_{t-1}|\theta^*)}=\frac{f(\theta^*|X)J_t(\theta_{t-1}|\theta^*)}{f(\theta_{t-1}|X)J_t(\theta^*|\theta_{t-1})}.$ 提高效率的方法见 Tierney (1994)。

12.4.3 Griddy Gibbs（PDF p.643–644）

金融模型常含非线性参数（ARMA 的 MA 参数、GARCH 参数），条件后验无闭式。Tanner (1996) 的 Griddy Gibbs 用于一元条件后验，适用面广但可能低效。对标量 $\theta_i$（$\theta_{-i}$ 为其余参数）：
1. 在适当区间取格点 $\theta_{i1}\le\cdots\le\theta_{im}$，计算 $w_j=f(\theta_{ij}|X,\theta_{-i})$；
2. 用 $\{w_j\}$ 近似逆 CDF；
3. 抽 U(0,1) 并经近似逆 CDF 变换得 $\theta_i$ 的抽样。
说明：不需归一化常数；最简近似为离散分布 $p(\theta_{ij})=w_j/\sum_vw_v$；区间选择需检查——若 Gibbs 抽样直方图在端点处有明显概率质量则扩大区间，若概率集中在内部则缩短（过宽导致多数 $w_j\approx0$，低效）。Griddy Gibbs 或 MH 可嵌入 Gibbs 抽样中抽部分参数（Metropolis-within-Gibbs）。

12.5 带时间序列误差的线性回归（Linear Regression with Time Series Errors）（PDF p.644–648）

模型（第 2 章曾用 SCA 估计）： $y_t=x_t'\beta+z_t,\quad z_t=\phi z_{t-1}+a_t,\quad a_t\sim\text{iid }N(0,\sigma^2),\quad t=1,\dots,n,\tag{12.6}$ $x_t=(1,x_{1t},\dots,x_{kt})'$（可含 $y_t$ 滞后），$\theta=(\beta',\phi,\sigma^2)'$。
思路：在回归估计与时间序列估计间迭代——已知时间序列模型则回归易估；已知回归则 $z_t=y_t-x_t'\beta$ 可估 AR(1)。需要 $f(\beta|Y,X,\phi,\sigma^2)$、$f(\phi|Y,X,\beta,\sigma^2)$、$f(\sigma^2|Y,X,\beta,\phi)$。
独立共轭先验： $\beta\sim N(\beta_o,\Sigma_o),\quad \phi\sim N(\phi_o,\sigma_o^2),\quad \frac{v\lambda}{\sigma^2}\sim\chi^2_v,\tag{12.7}$ **超参数（hyperparameters）**通常取 $\beta_o=0$，$\phi_o=0$，$\Sigma_o$ 为大对角阵（弱信息先验）。
$\beta$ 的条件后验：给定 $\phi$，做准差分 $y_{o,t}=y_t-\phi y_{t-1}$，$x_{o,t}=x_t-\phi x_{t-1}$，得 $y_{o,t}=\beta'x_{o,t}+a_t$（$t=2,\dots,n$）（12.8）。LS 估计 $\hat\beta=(\sum x_{o,t}x_{o,t}')^{-1}\sum x_{o,t}y_{o,t}\sim N(\beta,\sigma^2(\sum x_{o,t}x_{o,t}')^{-1})$。由 Result 12.1a： $(\beta|Y,X,\phi,\sigma^2)\sim N(\beta_*,\Sigma_*),\quad \Sigma_*^{-1}=\frac{\sum_{t=2}^nx_{o,t}x_{o,t}'}{\sigma^2}+\Sigma_o^{-1},\quad \beta_*=\Sigma_*\Big(\frac{\sum x_{o,t}x_{o,t}'}{\sigma^2}\hat\beta+\Sigma_o^{-1}\beta_o\Big).\tag{12.9}$
$\phi$ 的条件后验：给定 $\beta$，$z_t=y_t-\beta'x_t$，$\hat\phi=\sum z_{t-1}z_t/\sum z_{t-1}^2$，方差 $\sigma^2/\sum z_{t-1}^2$。由 Result 12.1： $\sigma_*^{-2}=\frac{\sum_{t=2}^nz_{t-1}^2}{\sigma^2}+\sigma_o^{-2},\qquad \phi_*=\sigma_*^2\Big(\frac{\sum z_{t-1}^2}{\sigma^2}\hat\phi+\sigma_o^{-2}\phi_o\Big).\tag{12.10}$
$\sigma^2$ 的条件后验：$a_t=z_t-\phi z_{t-1}$，由 Result 12.8： $\frac{v\lambda+\sum_{t=2}^na_t^2}{\sigma^2}\sim\chi^2_{v+(n-1)}.\tag{12.11}$
Gibbs 步骤：(1) 设定先验超参数；(2) 任意初值（如忽略时序误差的 OLS）；(3) 从 (12.9) 抽 $\beta$；(4) 从 (12.10) 抽 $\phi$；(5) 从 (12.11) 抽 $\sigma^2$；重复 3–5 多次，样本均值作点估计。
例 12.1（美国 1 年与 3 年期国债恒定期限周利率，1962-01-05 至 2009-04-10，圣路易斯联储）：因单位根，用周变化 $c_{3t}=r_{3t}-r_{3,t-1}$、$c_{1t}=r_{1t}-r_{1,t-1}$（%）。第 2 章用 MA(1) 误差，这里用 AR(2) 误差。传统方法（R）：
\[c_{3t}=0.782c_{1t}+z_t,\quad z_t=0.183z_{t-1}-0.036z_{t-2}+a_t,\quad \sigma_a=0.068,\tag{12.12}\]
标准误 0.0075、0.0201、0.0201；除滞后 4、6 处残差 ACF 边际显著外模型充分。写成 $c_{3t}=\beta c_{1t}+z_t$，$z_t=\phi_1z_{t-1}+\phi_2z_{t-2}+a_t$（12.13），先验 $\beta\sim N(0,4)$，$\phi\sim N[0,\mathrm{diag}(0.25,0.16)]$，$(v\lambda)/\sigma^2=(10\times0.05)/\sigma^2\sim\chi^2_{10}$。OLS 两步法初值；样本量 2466；迭代 2100 次，丢弃前 100 次。
- 表 12.1 后验均值（标准误）：$\beta$ 0.793 (0.008)，$\phi_1$ 0.184 (0.019)，$\phi_2$ −0.036 (0.021)，$\sigma^2$ 0.00479 (0.00013)，即 $\sigma\approx0.069$。图 12.1 抽样迹图稳定；图 12.2 为边际后验直方图。换初值结果相似，判断已收敛；后验均值接近 (12.12)，因样本大、模型简单。

12.6 缺失值与异常值（Missing Values and Outliers）（PDF p.648–）

加性异常值（additive outlier, AO）： $y_t=\begin{cases}x_h+\omega,&t=h\\x_t,&\text{otherwise}\end{cases}\tag{12.14}$ $\omega$ 为异常幅度，$x_t$ 为无异常序列；如记录错误（录入、测量误差）。异常值会导致参数估计严重偏倚和模型误设。
思想：把 $x_h$ 当作缺失值，求其在其余数据下的条件分布；若观测 $y_h$ 在该分布下很可能出现则非异常，概率很小则判为 AO。故缺失值处理与 AO 检测基于同一思想。缺失值可用 Kalman 滤波（第 11 章，Jones 1980）或 MCMC（McCulloch & Tsay 1994a）处理；异常值文献：Chang, Tiao & Chen (1988)、Tsay (1988)、Tsay, Peña & Pankratz (2000)，异常值按影响分四类，此处只讨论 AO。

12.6.1 缺失值（PDF p.649–）

AR($p$)：$x_t=\phi_1x_{t-1}+\cdots+\phi_px_{t-p}+a_t$（12.15），$x_h$ 缺失（$1<h<n$）。参数 $\theta=(\phi',x_h,\sigma^2)'$，把 $x_h$ 当未知参数。先验 $\phi\sim N(\phi_o,\Sigma_o)$，$x_h\sim N(\mu_o,\sigma_o^2)$，$v\lambda/\sigma^2\sim\chi^2_v$。$\phi$ 与 $\sigma^2$ 的条件后验同 12.5 节；$x_h$ 的条件后验为一元正态 $N(\mu_*,\sigma_h^2)$，可通过线性回归得到。
给定模型与数据，$x_h$ 只与 $\{x_{h-p},\dots,x_{h-1},x_{h+1},\dots,x_{h+p}\}$ 相关：
1. $t=h$：令 $y_h=\phi_1x_{h-1}+\cdots+\phi_px_{h-p}$，$b_h=-a_h$，则 $y_h=x_h+b_h=\phi_0x_h+b_h$，$\phi_0=1$；
2. $t=h+j$（$j=1,\dots,p$）：令 $y_{h+j}=x_{h+j}-\phi_1x_{h+j-1}-\cdots-\phi_{j-1}x_{h+1}-\phi_{j+1}x_{h-1}-\cdots-\phi_px_{h+j-p}$，$b_{h+j}=a_{h+j}$，则 $y_{h+j}=\phi_jx_h+b_{h+j}$。于是得到 $p+1$ 个方程 $y_{h+j}=\phi_jx_h+b_{h+j},\quad j=0,\dots,p,\tag{12.16}$ （正态对称，$a_h$ 与 $-a_h$ 同分布）这是含 $p+1$ 个数据点的无截距简单线性回归。（续）（PDF p.651）$x_h$ 的 LS 估计与方差： $\hat x_h=\frac{\sum_{j=0}^p\phi_jy_{h+j}}{\sum_{j=0}^p\phi_j^2},\qquad \mathrm{Var}(\hat x_h)=\frac{\sigma^2}{\sum_{j=0}^p\phi_j^2}.$ $p=1$ 时 $\hat x_h=\frac{\phi_1}{1+\phi_1^2}(x_{h-1}+x_{h+1})$，称 $x_h$ 的滤波值（filtered value）；高斯 AR(1) 时间可逆，故两侧邻居等权。
由 Result 12.1，$x_h$ 后验为正态： $\mu_*=\frac{\sigma^2\mu_o+\sigma_o^2(\sum_{j=0}^p\phi_j^2)\hat x_h}{\sigma^2+\sigma_o^2\sum_j\phi_j^2},\qquad \sigma_*^2=\frac{\sigma^2\sigma_o^2}{\sigma^2+\sigma_o^2\sum_j\phi_j^2}.\tag{12.17}$
连续缺失（patch）：(1) 直接推广——如 $x_h,x_{h+1}$ 都缺失，它们与 $\{x_{h-p},\dots,x_{h-1};x_{h+2},\dots,x_{h+p+1}\}$ 相关，构造多元线性回归，结合先验得 $(x_h,x_{h+1})$ 的二元正态后验，Gibbs 中联合抽取；(2) 在一次 Gibbs 迭代中多次使用单个缺失值公式，逐个抽取。因相邻缺失值相关，联合抽取更好，尤其缺失段较长时；缺失少时逐个抽取亦可。
备注：上述假设 $h-p\ge1$、$h+p\le n$；靠近样本端点时需调整回归的数据点数。

12.6.2 异常值检测（Outlier Detection）（PDF p.652–656）

MCMC 框架下 AO 检测很直接。除幅度相近的成片 AO 外，McCulloch & Tsay (1994a) 的简单 Gibbs 抽样效果良好（Justel, Peña & Tsay 2001）。对其他模型可借助 MH 或 Griddy Gibbs 抽非线性参数。
模型： $y_t=\delta_t\beta_t+x_t,\quad t=1,\dots,n,\tag{12.18}$ $\delta_t$ 为独立 Bernoulli，$P(\delta_t=1)=\epsilon$；$\beta_t$ 为独立的异常幅度；$x_t=\phi_0+\phi_1x_{t-1}+\cdots+\phi_px_{t-p}+a_t$ 为无异常 AR($p$)。每个时点都可能出现 AO，概率 $\epsilon$。$n$ 个数据却有 $2n+p+3$ 个参数：$\phi$、$\delta$、$\beta$、$\sigma^2$、$\epsilon$——$\delta,\beta$ 由数据增广引入。
共轭先验：$\phi\sim N(\phi_o,\Sigma_o)$，$v\lambda/\sigma^2\sim\chi^2_v$，$\epsilon\sim\mathrm{Beta}(\gamma_1,\gamma_2)$，$\beta_t\sim N(0,\xi^2)$。需要 $f(\phi|\cdot)$、$f(\delta_h|\cdot)$、$f(\beta_h|\cdot)$、$f(\epsilon|Y,\delta)$、$f(\sigma^2|\cdot)$。
$\phi$ 与 $\sigma^2$：给定 $\delta,\beta$，$x_t=y_t-\delta_t\beta_t$，$\hat\phi=(\sum_{t=p+1}^n\mathbf{x}_{t-1}\mathbf{x}_{t-1}')^{-1}\sum\mathbf{x}_{t-1}x_t$，$\mathbf{x}_{t-1}=(1,x_{t-1},\dots,x_{t-p})'$，后验同 (12.9) 形式；$\sigma^2$：$(v\lambda+\sum_{t=p+1}^na_t^2)/\sigma^2\sim\chi^2_{v+(n-p)}$，$a_t=x_t-\phi'\mathbf{x}_{t-1}$。（$\epsilon$ 的后验为 Beta($\gamma_1+\sum\delta_t$, $\gamma_2+n-\sum\delta_t$)，由 Result 12.3，原文未单列。）
$\delta_h$：只与 $j=h-p,\dots,h+p$ 的 $y_j,\beta_j,\delta_j$（$j\ne h$）、$\phi$、$\sigma^2$ 有关。令 $x_j^*=x_j$（$j\ne h$），$x_h^*=y_h$， $w_j=x_j^*-\phi_0-\phi_1x_{j-1}^*-\cdots-\phi_px_{j-p}^*$，$j=h,\dots,h+p$。
- Case I（$\delta_h=0$，概率 $1-\epsilon$）：$w_j=a_j\sim N(0,\sigma^2)$。
- Case II（$\delta_h=1$，概率 $\epsilon$）：$x_h^*=x_h+\beta_h$，$w_h\sim N(\beta_h,\sigma^2)$，$w_j\sim N(-\phi_{j-h}\beta_h,\sigma^2)$（$j>h$）；令 $\psi_0=-1$，$\psi_i=\phi_i$，统一为 $w_j\sim N(-\psi_{j-h}\beta_h,\sigma^2)$。
- 令 $m=\min(n,h+p)$， $P(\delta_h=1|\cdot)=\frac{\epsilon\exp[-\sum_{j=h}^m(w_j+\psi_{j-h}\beta_h)^2/(2\sigma^2)]}{\epsilon\exp[-\sum_{j=h}^m(w_j+\psi_{j-h}\beta_h)^2/(2\sigma^2)]+(1-\epsilon)\exp[-\sum_{j=h}^mw_j^2/(2\sigma^2)]},\tag{12.19}$ 即按先验概率加权比较两种情形下的似然。
$\beta_h$：$\delta_h=0$ 时 $\beta_h\sim N(0,\xi^2)$（先验）；$\delta_h=1$ 时由回归 $w_j=-\psi_{j-h}\beta_h+a_j$ 得 $\hat\beta_h=\frac{\sum_{j=h}^m-\psi_{j-h}w_j}{\sum_{j=h}^m\psi_{j-h}^2}$，方差 $\sigma^2/\sum\psi_{j-h}^2$，后验 $\beta_h^*=\frac{-(\sum_{j=h}^m\psi_{j-h}w_j)\xi^2}{\sigma^2+(\sum_{j=h}^m\psi_{j-h}^2)\xi^2},\qquad \sigma_{h*}^2=\frac{\sigma^2\xi^2}{\sigma^2+(\sum\psi_{j-h}^2)\xi^2}.$
例 12.2：美国 3 年期国债周变化 $c_{3t}$，1988-03-18 至 1999-09-10，600 个观测（例 12.1 子序列，图 12.3(a)）。PACF 建议 AR(3)：$c_{3t}=0.227c_{3,t-1}+0.006c_{3,t-2}+0.114c_{3,t-3}+a_t$，$\sigma^2=0.0128$（原文第三项下标误为 $t-2$），标准误 0.041、0.042、0.041，$Q(12)=11.4$ 不显著。
- Gibbs 联合估计 + AO 检测，先验 $\phi\sim N(0,0.25I_3)$，$v\lambda/\sigma^2=5\times0.00256/\sigma^2\sim\chi^2_5$，$\gamma_1=5,\gamma_2=95$（预期 5% 为 AO），$\xi^2=0.1$（$0.00256\approx\sigma^2/5$，$\xi^2\approx9\sigma^2$）。初值 $\epsilon=0.05$，$\sigma^2=0.012$，$\phi=(0.2,0.02,0.1)$；1050 次迭代丢弃前 50。
- 结果：$c_{3t}=0.252c_{3,t-1}+0.003c_{3,t-2}+0.110c_{3,t-3}+a_t$，$\sigma^2=0.0118$，后验标准差 0.046、0.045、0.046、0.0008，与 ML 相近。图 12.3(b) 各点为 AO 的后验概率，(c) 后验异常幅度。$t=323$（1994-05-20）AO 概率 0.83，幅度 −0.304（$c_{3t}$ 从 0.24 变为 −0.34，两周内约下降 0.6%）；$t=201$（1992-01-17）概率 0.58，幅度 0.176（从 −0.02 变为 0.33）。
备注：Gibbs 检测计算量大，但联合估计参数与异常值；传统方法把估计与检测分开，快但多个异常值时可能虚假检测。SCA 也识别出 $t=323$、$t=201$ 为最显著的两个 AO，幅度 −0.39、0.36。

12.7 随机波动率模型（Stochastic Volatility Models）（PDF p.656–668）

MCMC 在金融中的重要应用是估计 SV 模型（Jacquier, Polson & Rossi 1994，JPR）。

单变量 SV： $r_t=\beta_0+\beta_1x_{1t}+\cdots+\beta_px_{pt}+a_t,\quad a_t=\sqrt{h_t}\epsilon_t,\tag{12.20}$ $\ln h_t=\alpha_0+\alpha_1\ln h_{t-1}+v_t,\tag{12.21}$ $x_{it}$ 在 $t-1$ 时可得（可含滞后收益），$\epsilon_t\sim N(0,1)$，$v_t\sim N(0,\sigma_v^2)$，两者独立；取对数保证 $h_t>0$；$|\alpha_1|<1$ 保证对数波动平稳；可用更高阶 AR。
参数 $\beta$、$\omega=(\alpha_0,\alpha_1,\sigma_v^2)'$，不可观测波动 $H=(h_1,\dots,h_n)'$ 为辅助变量。ML 困难：似然是对 $n$ 维 $H$ 的混合积分 $f(R|X,\beta,\omega)=\int f(R|X,\beta,H)f(H|\omega)dH$。贝叶斯框架下 $H$ 是增广参数，先验 $p(\beta,\omega)=p(\beta)p(\omega)$，Gibbs 抽取 $f(\beta|R,X,H,\omega)$、$f(H|R,X,\beta,\omega)$、$f(\omega|R,X,\beta,H)$。

12.7.1 单变量模型估计（PDF p.657–663）

$\beta$：给定 $H$，均值方程为异方差回归，除以 $\sqrt{h_t}$：$r_{o,t}=x_{o,t}'\beta+\epsilon_t$（12.22），$r_{o,t}=r_t/\sqrt{h_t}$，$x_{o,t}=x_t/\sqrt{h_t}$。先验 $N(\beta_o,A_o)$，后验 $N(\beta_*,A_*)$：$A_*^{-1}=\sum x_{o,t}x_{o,t}'+A_o^{-1}$，$\beta_*=A_*(\sum x_{o,t}r_{o,t}+A_o^{-1}\beta_o)$。
$h_t$ 逐个抽取： $f(h_t|\cdot)\propto f(a_t|h_t)f(h_t|h_{t-1})f(h_{t+1}|h_t)\propto h_t^{-1.5}\exp\Big[-\frac{(r_t-x_t'\beta)^2}{2h_t}-\frac{(\ln h_t-\mu_t)^2}{2\sigma^2}\Big],\tag{12.23}$ $\mu_t=\frac{\alpha_0(1-\alpha_1)+\alpha_1(\ln h_{t+1}+\ln h_{t-1})}{1+\alpha_1^2}$，$\sigma^2=\frac{\sigma_v^2}{1+\alpha_1^2}$。推导用到：$a_t|h_t\sim N(0,h_t)$；$\ln h_t|\ln h_{t-1}$ 与 $\ln h_{t+1}|\ln h_t$ 的正态性；$d\ln h_t=h_t^{-1}dh_t$；配方恒等式 $(x-a)^2A+(x-b)^2C=(x-c)^2(A+C)+(a-b)^2AC/(A+C)$，$c=(Aa+Cb)/(A+C)$（Box & Tiao 1973, p.418 引理 1 的标量版），取 $A=1$，$a=\alpha_0+\alpha_1\ln h_{t-1}$（原文漏 $\alpha_1$），$C=\alpha_1^2$，$b=(\ln h_{t+1}-\alpha_0)/\alpha_1$。JPR 用 Metropolis 抽 $h_t$；本节用 Griddy Gibbs，$h_t$ 的范围取样本无条件方差的倍数。
$\omega$：分块 $\alpha=(\alpha_0,\alpha_1)'$ 与 $\sigma_v^2$。给定 $H$，$\ln h_t$ 为 AR(1)：先验 $N(\alpha_o,C_o)$，后验 $C_*^{-1}=\sum_{t=2}^nz_tz_t'/\sigma_v^2+C_o^{-1}$，$\alpha_*=C_*(\sum z_t\ln h_t/\sigma_v^2+C_o^{-1}\alpha_o)$，$z_t=(1,\ln h_{t-1})'$；$v_t=\ln h_t-\alpha_0-\alpha_1\ln h_{t-1}$，先验 $m\lambda/\sigma_v^2\sim\chi^2_m$，后验 $(m\lambda+\sum_{t=2}^nv_t^2)/\sigma_v^2\sim\chi^2_{m+n-1}$。
端点备注：(12.23) 适用于 $1<t<n$。简单处理：$h_1$ 固定，$t=n$ 用 $\ln h_n\sim N(\alpha_0+\alpha_1\ln h_{n-1},\sigma_v^2)$；或用 $h_{n+1}$ 的预测（原点 $n-1$ 的 2 步预测 $\alpha_0+\alpha_1(\alpha_0+\alpha_1\ln h_{n-1})$）和 $h_0$ 的后向预测（时间可逆：$\ln h_t-\eta=\alpha_1(\ln h_{t+1}-\eta)+v_t^*$，$\eta=\alpha_0/(1-\alpha_1)$，2 步后向预测 $\alpha_1^2(\ln h_2-\eta)$）。
备注：(12.23) 也可由 AR(1) 缺失值结果得到：把 $\ln h_t$ 当缺失值，构造两观测的回归 $y_t=\ln h_t+b_t$（$y_t=\alpha_0+\alpha_1\ln h_{t-1}$）与 $y_{t+1}=\alpha_1\ln h_t+b_{t+1}$（$y_{t+1}=\ln h_{t+1}-\alpha_0$），LS 估计正是 $\mu_t$，方差 $\sigma_v^2/(1+\alpha_1^2)$；(12.23) 即 $a_t\sim N(0,h_t)$ 与 $\widehat{\ln h_t}\sim N(\ln h_t,\sigma_v^2/(1+\alpha_1^2))$ 之积再做变量变换。易推广到 AR($p$)（假设前 $p$ 个 $h_t$ 固定）。$h_t$ 初值可由第 3 章波动率模型拟合得到。
例 12.3（S&P 500 月度对数收益 %，1962-01 至 2009-12，575 个观测，用每月首个交易日收盘指数；图 12.4）：
- 高斯 GARCH(1,1)：$r_t=0.552+a_t$，$h_t=0.878+0.125a_{t-1}^2+0.837h_{t-1}$（12.26），t 值均 >2.56；$Q(12)=10.04(0.61)$，平方 6.14(0.91)。
- SV 模型 $r_t=\mu+a_t$，$a_t=\sqrt{h_t}\epsilon_t$，$\ln h_t=\alpha_0+\alpha_1\ln h_{t-1}+v_t$（12.27）。先验 $\mu\sim N(0,4)$，$\alpha\sim N[(0,0.6)',\mathrm{diag}(0.25,0.04)]$，$10\times0.1/\sigma_v^2\sim\chi^2_{10}$。初值：$h_{0t}$ 取 GARCH 拟合值，$\alpha,\sigma_v^2$ 取 $\ln h_{0t}$ 的 LS 估计，$\mu$ 取样本均值。Griddy Gibbs 400 格点，第 $j$ 次迭代范围 $[\eta_{1t},\eta_{2t}]$，$\eta_{1t}=0.6\max(h_{j-1,t},h_{0t})$，$\eta_{2t}=1.4\min(h_{j-1,t},h_{0t})$（原文如此，疑应为 min/max 互换）。2500 次迭代丢弃前 500。
- 后验均值（标准误）：$\mu$ 0.409 (0.157)，$\alpha_0$ 0.454 (0.068)，$\alpha_1$ 0.837 (0.025)，$\sigma_v^2$ 0.086 (0.007)。$\alpha_1=0.837$ 说明波动强持续，小于 JPR 用日数据所得值。图 12.5 先验（虚线，较无信息）与后验（实线）密度，$\mu$、$\sigma_v^2$ 后验很集中；图 12.6 SV 后验均值波动率与 GARCH 拟合形态相似。不同初值、先验、迭代次数结果稳定；Griddy Gibbs 的结果与效率依赖 $h_t$ 范围设定。

12.7.2 多元随机波动率模型（PDF p.663–668）

用第 10 章 Cholesky 分解，二元情形：$b_{1t}=a_{1t}$，$b_{2t}=a_{2t}-q_{21,t}b_{1t}$（回归 $a_{2t}=q_{21,t}a_{1t}+b_{2t}$）， $\Sigma_t=\begin{bmatrix}1&0\\q_{21,t}&1\end{bmatrix}\begin{bmatrix}g_{11,t}&0\\0&g_{22,t}\end{bmatrix}\begin{bmatrix}1&q_{21,t}\\0&1\end{bmatrix},\tag{12.28}$ $g_{ii,t}=\mathrm{Var}(b_{it}|F_{t-1})$，$b_{1t}\perp b_{2t}$。
模型： $r_t=\beta_0+\beta_1x_t+a_t,\tag{12.29}\quad \ln g_{ii,t}=\alpha_{i0}+\alpha_{i1}\ln g_{ii,t-1}+v_{it},\tag{12.30}\quad q_{21,t}=\gamma_0+\gamma_1q_{21,t-1}+u_t,\tag{12.31}$ $v_{1t},v_{2t},u_t$ 为独立高斯白噪声，方差 $\sigma_{iv}^2,\sigma_u^2$。传统参数 $\beta,\alpha_i,\sigma_{iv}^2,\gamma,\sigma_u^2$；增广参数 $Q,G_1,G_2$。
Gibbs：(1) 用 (12.22) 逐行抽 $\beta_0,\beta_1$；(2) 用 (12.23)（$a_t$ 换成 $a_{1t}$）抽 $g_{11,t}$；(3) 同单变量法抽 $\alpha_1,\sigma_{1v}^2$；对 $\alpha_2,\sigma_{2v}^2,g_{22,t}$ 先算 $b_{2t}=a_{2t}-q_{21,t}a_{1t}\sim N(0,g_{22,t})$ 再同法处理。另需：
- $\gamma$：AR(1)，先验 $N(\gamma_o,D_o)$，$D_*^{-1}=\sum z_tz_t'/\sigma_u^2+D_o^{-1}$，$\gamma_*=D_*(\sum z_tq_{21,t}/\sigma_u^2+D_o^{-1}\gamma_o)$，$z_t=(1,q_{21,t-1})'$；
- $\sigma_u^2$：$(m\lambda+\sum u_t^2)/\sigma_u^2\sim\chi^2_{m+n-1}$，$u_t=q_{21,t}-\gamma_0-\gamma_1q_{21,t-1}$；
- $q_{21,t}$： $f(q_{21,t}|\cdot)\propto g_{22,t}^{-0.5}\exp\Big[-\frac{(a_{2t}-q_{21,t}a_{1t})^2}{2g_{22,t}}\Big]\exp\Big[-\frac{(q_{21,t}-\mu_t)^2}{2\sigma^2}\Big],\tag{12.32}$ $\mu_t=\frac{\gamma_0(1-\gamma_1)+\gamma_1(q_{21,t-1}+q_{21,t+1})}{1+\gamma_1^2}$，$\sigma^2=\sigma_u^2/(1+\gamma_1^2)$。第一项是均值 $a_{2t}/a_{1t}$、方差 $g_{22,t}/a_{1t}^2$ 的正态，第二项是 $N(\mu_t,\sigma^2)$，由 Result 12.1 得闭式正态后验： $\frac1{\sigma_*^2}=\frac{a_{1t}^2}{g_{22,t}}+\frac{1+\gamma_1^2}{\sigma_u^2},\qquad \mu_*=\sigma_*^2\Big(\frac{1+\gamma_1^2}{\sigma_u^2}\mu_t+\frac{a_{1t}^2}{g_{22,t}}\cdot\frac{a_{2t}}{a_{1t}}\Big).$
例 12.4（IBM 与 S&P 500 月度对数收益，1962-01 至 2009-12，$r_t=(IBM_t,SP_t)'$，图 12.7）：
- Cholesky 时变相关 GARCH：$r_t=\beta_0+a_t$（12.33），$g_{11,t}=\alpha_{10}+\alpha_{11}g_{11,t-1}+\alpha_{12}a_{1,t-1}^2$（12.34），$g_{22,t}=\alpha_{20}+\alpha_{22}b_{2,t-1}^2$（12.35），$q_{21,t}=\gamma_0$（12.36）。表 12.2(a)：$\beta_{01}=0.69(0.30)$，$\beta_{02}=0.49(0.18)$，$\alpha_{10}=3.98(1.22)$，$\alpha_{11}=0.80(0.04)$，$\alpha_{12}=0.12(0.03)$，$\alpha_{20}=10.67(0.53)$，$\alpha_{22}=0.12(0.04)$，$\gamma_0=0.37(0.01)$。
- BEKK(1,1)：$\beta_0=(0.70,0.54)'$，$A=\begin{bmatrix}0.80&0\\0.83&0.01\end{bmatrix}$，$A_1=\begin{bmatrix}0.07&0.33\\-0.06&0.43\end{bmatrix}$，$B_1=\begin{bmatrix}1.00&-0.12\\0.01&0.90\end{bmatrix}$（Matlab 估计）。
- SV 模型：同均值方程，$\ln g_{ii,t}=\alpha_{i0}+\alpha_{i1}\ln g_{ii,t-1}+v_{it}$（12.37–12.38），$q_{21,t}=\gamma_0+u_t$（12.39）。先验 $\beta_{i0}\sim N(0,4)$，$\alpha_i\sim N[(0,0.7)',\mathrm{diag}(0.25,0.04)]$，$\gamma_0\sim N(0,1)$，$10\times0.1/\sigma_{iv}^2\sim\chi^2_{10}$，$5\times0.2/\sigma_u^2\sim\chi^2_5$。初值取自 BEKK；$t=1$ 值固定；2500 次迭代丢弃 500；$g_{ii,t}$ 用 500 格点 Griddy Gibbs。表 12.2(b) 后验均值（标准误）：$\beta_{01}$ 0.53(0.26)，$\beta_{02}$ 0.51(0.17)，$\alpha_{10}$ 0.75(0.11)，$\alpha_{11}$ 0.80(0.03)，$\sigma_{1v}^2$ 0.07(0.01)，$\alpha_{20}$ 0.43(0.06)，$\alpha_{21}$ 0.81(0.03)，$\sigma_{2v}^2$ 0.07(0.01)，$\gamma_0$ 0.38(0.03)，$\sigma_u^2$ 0.07(0.01)。
- 收敛检查：不同初值与迭代次数，500+2000 与 500+1000 两个 Gibbs 样本的 $g_{11,t},g_{22,t},q_{21,t},\sigma_{22,t},\sigma_{21,t},\rho_{21,t}$ 后验均值散点图都贴近 $y=x$（图 12.8），结果稳定。
- 比较三模型：(1) 均值方程基本相同；(2) IBM 条件方差（图 12.9）三者都显示波动聚集和上升趋势，但 GARCH 峰值更高且 1993 年多一个峰；(3) S&P 500 条件方差（图 12.10）GARCH 在 1993 年附近出现额外峰，单变量分析（图 12.6）没有——是二元 GARCH 依赖 IBM 收益而产生的虚假波动峰，SV 与 BEKK 都没有；SV 的 S&P 波动与单变量结果相似；(4) 条件相关（图 12.11）差异很大：Cholesky-GARCH 相关平滑且恒正，均值 0.59、标准差 0.07、范围 (0.411, 0.849)；BEKK 在 1993 年附近出现小负值，均值 0.59、标准差 0.13、范围 (−0.020, 0.877)；SV 逐月变化剧烈，均值 0.60、标准差 0.14、范围 (−0.161, 0.839)，在若干孤立时期出现负相关——因 SV 中 $q_{21,t}$ 含随机冲击 $u_t$。
备注：Gibbs 估计可推广到其他二元 SV 模型，所需条件后验需扩展但思路相同。

12.8 SV 估计的新方法（New Approach to SV Estimation）（PDF p.669–）

在 Kalman 滤波框架内用前向滤波后向抽样（FFBS, forward filtering and backward sampling）提高 Gibbs 效率：借助正态混合联合抽取整个波动过程，大幅减少计算时间；可推广到带杠杆效应与跳跃的随机扩散模型。
重参数化： $r_t=x_t'\beta+\sigma_0\exp(z_t/2)\epsilon_t,\tag{12.40}\qquad z_{t+1}=\alpha z_t+\eta_t,\tag{12.41}$ $z_t$ 为零均值对数波动，$(\epsilon_t,\eta_t)'$ 二元正态，协方差 $\Sigma=\begin{bmatrix}1&\rho\sigma_\eta\\\rho\sigma_\eta&\sigma_\eta^2\end{bmatrix}$；$\rho$ 刻画杠杆效应，通常为负（负收益推高波动）。与前面模型关系：$z_t=\ln h_t-\ln\sigma_0^2$，$\sigma_0^2=\exp\{E[\ln h_t]\}$。波动率 $\sigma_0\exp(z_t/2)$ 恒正；关键是 $\eta_t$ 是 $z_{t+1}$ 的新息、与 $z_t$ 独立——这一时间错位使杠杆效应可处理；若写成 $z_t=\alpha z_{t-1}+\eta_t$，则 $\eta_t$ 与 $\epsilon_t$ 相关会导致 $z_t$ 与 $\epsilon_t$ 相关，产生识别问题。
备注：等价写法 $r_t=x_t'\beta+\sigma_0\exp(z_{t-1}/2)\epsilon_t$，$z_t=\alpha z_{t-1}+\eta_t$；或 $r_t=x_t'\beta+\exp(z_{t-1}^*/2)\epsilon_t$，$z_t^*=\alpha_0+\alpha z_{t-1}^*+\eta_t$，$E(z_t^*)=\alpha_0/(1-\alpha)\ne0$。
参数 $\beta,\sigma_0,\alpha,\rho,\sigma_\eta$ 与 $z=(z_1,\dots,z_n)'$（设 $z_1$ 已知）。条件后验：
1. $\beta$：同 12.7.1，把 $\sqrt{h_t}$ 换为 $\sigma_0\exp(z_t/2)$。
2. $\alpha$：给定 $z,\sigma_\eta^2$ 为 AR(1) 系数，正态先验下后验易得。
3. $\sigma_0^2$：$v_t=(r_t-x_t'\beta)\exp(-z_t/2)=\sigma_0\epsilon_t$ iid $N(0,\sigma_0^2)$；先验 $m\lambda/\sigma_0^2\sim\chi^2_m$，后验 $(m\lambda+\sum_{t=1}^nv_t^2)/\sigma_0^2\sim\chi^2_{m+n}$。
4. $(\rho,\sigma_\eta^2)$：给定其余，可得 $b_t=(\epsilon_t,\eta_t)'$，似然 $\propto|\Sigma|^{-(n-1)/2}\exp[-\frac12\mathrm{tr}(\Sigma^{-1}\sum b_tb_t')]$，但 $\rho,\sigma_\eta^2$ 无法分离。采用 Jacquier, Polson & Rossi (2004) 重参数化 $\Sigma=\begin{bmatrix}1&\phi\\\phi&\omega+\phi^2\end{bmatrix}$，$\omega=\sigma_\eta^2(1-\rho^2)$，$|\Sigma|=\omega$，$\Sigma^{-1}=\frac1\omega S+\begin{bmatrix}1&0\\0&0\end{bmatrix}$，$S=\begin{bmatrix}\phi^2&-\phi\\-\phi&1\end{bmatrix}$。令 $e=(\epsilon_2,\dots,\epsilon_n)'$，$\eta=(\eta_2,\dots,\eta_n)'$，$R=(e,\eta)'(e,\eta)$，似然 $\ell(\phi,\omega)\propto\omega^{-(n-1)/2}\exp[-\frac1{2\omega}\mathrm{tr}(SR)]$。共轭先验 $\omega\sim IG(\gamma_0/2,\gamma_1/2)$，$\phi|\omega\sim N(0,\omega/2)$，则 $\phi\sim N\Big(\tilde\phi,\frac{\omega}{2+e'e}\Big),\ \tilde\phi=\frac{e'\eta}{2+e'e};\qquad \omega\sim IG\Big(\frac{n+1+\gamma_0}{2},\ \frac12\Big[\gamma_1+\eta'\eta-\frac{(e'\eta)^2}{2+e'e}\Big]\Big).$ 再得 $\sigma_\eta^2=\omega+\phi^2$，$\rho=\phi/\sigma_\eta$。IG($\alpha,\beta$) 密度 $f(\omega)=\frac{\beta^\alpha}{\Gamma(\alpha)}\omega^{-(\alpha+1)}e^{-\beta/\omega}$。
5. $z$ 的联合抽取：$(r_t-x_t'\beta)^2/\sigma_0^2=\exp(z_t)\epsilon_t^2$，令 $y_t=\ln[(r_t-x_t'\beta)^2/\sigma_0^2]$， $y_t=z_t+\epsilon_t^*,\quad \epsilon_t^*=\ln(\epsilon_t^2),\tag{12.42}$ 以 (12.42) 为观测方程、(12.41) 为状态方程（原文误写为 12.40），构成线性状态空间模型，但 $\epsilon_t^*\sim\ln\chi^2_1$ 非正态。Kim, Shephard & Chib (1998)（KSC）用 7 个正态的混合近似：$f(\epsilon_t^*)\approx\sum_{i=1}^7p_iN(\mu_i,\varpi_i^2)$（另见 Chib, Nardari & Shephard 2002）。表 12.3： | $i$ | $p_i$ | $\mu_i$ | $\varpi_i^2$ | |---|---|---|---| | 1 | 0.00730 | −11.4004 | 5.7960 | | 2 | 0.10556 | −5.2432 | 2.6137 | | 3 | 0.00002 | −9.8373 | 5.1795 | | 4 | 0.04395 | 1.5075 | 0.1674 | | 5 | 0.34001 | −0.6510 | 0.6401 | | 6 | 0.24566 | 0.5248 | 0.3402 | | 7 | 0.25750 | −2.3586 | 1.2626 | 图 12.12：基于 100,000 个观测，$\ln\chi^2_1$ 密度（实线）与 7 正态混合（虚线）几乎重合。（续）
为什么需要高斯状态空间模型：可以联合、高效地抽取整个对数波动序列 $z$。考虑无杠杆（$\eta_t$ 与 $e_t$ 不相关）的特殊模型 $z_{t+1}=\alpha z_t+\eta_t,\ \eta_t\sim\text{iid }N(0,\sigma_\eta^2),\tag{12.43}\qquad y_t=c_t+z_t+e_t,\ e_t\sim\text{ind. }N(0,H_t),\tag{12.44}$ $(c_t,H_t)$ 取表 12.3 中某个 $(\mu_i,\varpi_i^2)$。Kalman 滤波： $v_t=y_t-c_t-z_{t|t-1},\ V_t=\Sigma_{t|t-1}+H_t,\ z_{t|t}=z_{t|t-1}+\Sigma_{t|t-1}V_t^{-1}v_t,\ \Sigma_{t|t}=\Sigma_{t|t-1}-\Sigma_{t|t-1}^2V_t^{-1},\ z_{t+1|t}=\alpha z_{t|t},\ \Sigma_{t+1|t}=\alpha^2\Sigma_{t|t}+\sigma_\eta^2.\tag{12.45}$

前向滤波后向抽样（FFBS）（PDF p.675–677）：

分解联合后验（利用马尔可夫性：给定 $z_{t+1}$，$z_t$ 与 $z_{t+j}$（$j>1$）独立）： $p(z|F_n)=p(z_n|F_n)\,p(z_{n-1}|z_n,F_n)\,p(z_{n-2}|z_{n-1},F_n)\cdots p(z_2|z_3,F_n).\tag{12.46}$
$p(z_n|F_n)=N(z_{n|n},\Sigma_{n|n})$。关键性质：$p(z_{n-1}|z_n,F_n)=p(z_{n-1}|z_n,F_{n-1},v_n)=p(z_{n-1}|z_n,F_{n-1})$（12.47–12.48），因 $z_{n-1}$ 与 $v_n$（给定 $z_n$）独立。由 Kalman 滤波 $\begin{bmatrix}z_{n-1}\\z_n\end{bmatrix}\Big|F_{n-1}\sim N\left(\begin{bmatrix}z_{n-1|n-1}\\z_{n|n-1}\end{bmatrix},\begin{bmatrix}\Sigma_{n-1|n-1}&\alpha\Sigma_{n-1|n-1}\\\alpha\Sigma_{n-1|n-1}&\Sigma_{n|n-1}\end{bmatrix}\right),\tag{12.49}$ 由定理 11.1：$p(z_{n-1}|z_n,F_n)=N(\mu_{n-1}^*,\Sigma_{n-1}^*)$（12.50），$\mu_{n-1}^*=z_{n-1|n-1}+\alpha\Sigma_{n-1|n-1}\Sigma_{n|n-1}^{-1}(z_n-z_{n|n-1})$，$\Sigma_{n-1}^*=\Sigma_{n-1|n-1}-\alpha^2\Sigma_{n-1|n-1}^2\Sigma_{n|n-1}^{-1}$。
一般地 $p(z_t|z_{t+1},F_n)=p(z_t|z_{t+1},F_t)$，由 (12.51) 的二元正态得 $N(\mu_t^*,\Sigma_t^*)$：$\mu_t^*=z_{t|t}+\alpha\Sigma_{t|t}\Sigma_{t+1|t}^{-1}(z_{t+1}-z_{t+1|t})$，$\Sigma_t^*=\Sigma_{t|t}-\alpha^2\Sigma_{t|t}^2\Sigma_{t+1|t}^{-1}$。
流程：给定 $z_{1|0},\Sigma_{1|0}$，前向运行 Kalman 滤波，然后从 $z_n$ 开始后向递归抽样，得到 $z$ 的一个联合实现。称 FFBS（Carter & Kohn 1994；Frühwirth-Schnatter 1994）。因 $z_t$ 序列相关，联合抽取比逐个抽取高效得多。
备注：FFBS 适用于一般线性高斯状态空间模型，核心恒等式 $p(S_t|S_{t+1},F_n)=p(S_t|S_{t+1},F_t,v_{t+1},\dots,v_n)=p(S_t|S_{t+1},F_t)$。
混合指示变量：为确定 $(c_t,H_t)$，引入独立指示变量 $I_t\in\{1,\dots,7\}$。给定 $z_t$，令 $q_{it}=\phi[(y_t-z_t-\mu_i)/\varpi_i]$（原文记为标准正态 CDF，按似然含义应为密度 $\varpi_i^{-1}\phi(\cdot)$），作为 $I_t$ 的似然；以 $p_i$ 为先验，后验 $p_{it}=p_iq_{it}/\sum_jp_jq_{jt}$；抽到 $I_t=j$ 则 $c_t=\mu_j$，$H_t=\varpi_j^2$。由此用近似线性高斯状态空间模型联合抽 $z$，Gibbs 抽样估计单变量 SV 很高效。
杠杆效应的处理（Artigas & Tsay 2004）：(12.42) 的平方变换丢失 $\eta_t$ 与 $\epsilon_t$ 的相关性，无法估计杠杆。当 $\rho\ne0$：$\eta_t=\rho\sigma_\eta\epsilon_t+\eta_t^*$，$\eta_t^*$ 与 $\epsilon_t$ 独立，$\mathrm{Var}(\eta_t^*)=\sigma_\eta^2(1-\rho^2)$；代入 $\epsilon_t=(r_t-x_t'\beta)\exp(-z_t/2)/\sigma_0$： $z_{t+1}=\alpha z_t+\frac{\rho\sigma_\eta(r_t-x_t'\beta)}{\sigma_0}\exp\Big(-\frac{z_t}2\Big)+\eta_t^*=G(z_t)+\eta_t^*,\tag{12.52}$ 非线性转移方程，Kalman 滤波不再直接适用。用时变线性化（类似扩展 Kalman 滤波）修改 (12.45) 后两式： $z_{t+1|t}=G(z_{t|t}),\qquad \Sigma_{t+1|t}=g(z_{t|t})^2\Sigma_{t|t}+\sigma_\eta^2(1-\rho^2),\tag{12.53}$ $g(z_{t|t})=\partial G(x)/\partial x|_{x=z_{t|t}}$（原文称"平滑状态"，实为滤波状态）。
例 12.5（S&P 500 月度对数收益，1962-01 至 2004-11，515 个观测，例 12.3 的子序列）：模型 $r_t=\mu+\sigma_o\exp(z_t/2)\epsilon_t$，$z_{t+1}=\alpha z_t+\eta_t$（12.54）；模型 1 无杠杆，模型 2 有 $\mathrm{corr}(\epsilon_t,\eta_t)=\rho$。Matlab 程序，FFBS–Gibbs 2000+8000 次迭代（前 2000 为 burn-in）。表 12.4：
- 有杠杆：$\mu$ 0.0081 (0.0274)，$\sigma_o$ 0.0764 (0.0255)，$\alpha$ −0.0616 (0.1186)，$\sigma_\eta$ 2.5639 (0.3924)，$\rho$ −0.3892 (0.0292)；
- 无杠杆：$\mu$ 0.0080 (0.0279)，$\sigma_o$ 0.0775 (0.0266)，$\alpha$ −0.0613 (0.1164)，$\sigma_\eta$ 2.5827 (0.3783)。
- （注：$\alpha\approx-0.06$、$\sigma_\eta\approx2.56$ 与例 12.3 的高持续性差异很大，疑原表列错位或参数含义不同，编写时宜核对；按原文记录。）$\hat\rho=-0.39$ 与文献常见值接近。图 12.13 两模型波动率后验均值非常接近，与例 12.3 形态和量级相似（注意图 12.6 是百分比收益的条件方差，图 12.13 是对数收益的条件标准差）。

12.9 马尔可夫转换模型（Markov Switching Models）（PDF p.680–686）

MCMC 相对传统似然法的另一优势场景。McCulloch & Tsay (1994b) 用 Gibbs 估计各状态波动率不变的 Markov 转换模型，应用于美国季调实际 GNP 季度增长率，发现"收缩"与"扩张"期动态显著不同。本节关注波动率转换。
两状态、不同风险溢价与 GARCH 动态的模型： $r_t=\begin{cases}\beta_1\sqrt{h_t}+\sqrt{h_t}\epsilon_t,\ h_t=\alpha_{10}+\alpha_{11}h_{t-1}+\alpha_{12}a_{t-1}^2,&s_t=1,\\ \beta_2\sqrt{h_t}+\sqrt{h_t}\epsilon_t,\ h_t=\alpha_{20}+\alpha_{21}h_{t-1}+\alpha_{22}a_{t-1}^2,&s_t=2,\end{cases}\tag{12.55}$ $a_t=\sqrt{h_t}\epsilon_t$，$\epsilon_t\sim N(0,1)$，$\alpha_{ij}$ 满足无条件方差存在的条件。转移概率 $P(s_t=2|s_{t-1}=1)=e_1,\quad P(s_t=1|s_{t-1}=2)=e_2,\quad 0<e_i<1,\tag{12.56}$ $e_i$ 小表示倾向停留在状态 $i$，期望持续期 $1/e_i$。识别约束 $\beta_2>\beta_1$（状态 2 风险溢价更高，仅为唯一标记状态）。$\alpha_{1j}=\alpha_{2j}$ 时退化为所有状态同一 GARCH；若 $\beta_i\sqrt{h_t}$ 换成 $\beta_i$ 则为简单 Markov 转换 GARCH。(12.55) 是 Markov 转换 GARCH-M 模型。
设 $h_1$ 固定为样本方差（也可作为参数估计，大样本下影响可忽略）。传统参数 $\beta=(\beta_1,\beta_2)'$、$\alpha_i=(\alpha_{i0},\alpha_{i1},\alpha_{i2})'$、$e=(e_1,e_2)'$；状态向量 $S=(s_1,\dots,s_n)'$ 为增广参数；给定 $h_1,\alpha_i,S$，$H$ 可递推计算。
收益依赖波动率意味着收益有序列相关、有一定可预测性；但未来状态未知，预测是各状态配置的混合，点预测不确定性高。
似然是对所有状态配置的混合，复杂；Gibbs 只需 $f(\beta|R,S,H,\alpha_1,\alpha_2)$、$f(\alpha_i|R,S,H,\alpha_{j\ne i})$、$P(S|R,h_1,\alpha_1,\alpha_2)$、$f(e_i|S)$。先验：$\beta_i\sim N(\beta_{io},\sigma_{io}^2)$，$e_i\sim\mathrm{Beta}(\gamma_{i1},\gamma_{i2})$，$\alpha_{ij}$ 在适当区间上均匀（非线性参数，用 Griddy Gibbs）。
1. $\beta_i$：只依赖状态 $i$ 的数据。令 $r_{it}=r_t/\sqrt{h_t}$（$s_t=i$），则 $r_{it}=\beta_i+\epsilon_t$；$\bar r_i=\sum_{s_t=i}r_{it}/n_i$，后验正态：$1/\sigma_{i*}^2=n_i+1/\sigma_{io}^2$，$\beta_i^*=\sigma_{i*}^2(n_i\bar r_i+\beta_{io}/\sigma_{io}^2)$。
2. $\alpha_{ij}$：逐个用 Griddy Gibbs，对数条件后验 $\propto-\frac12\sum_{s_t=i}[\ln h_t+(r_t-\beta_i\sqrt{h_t})^2/h_t]$（$h_t$ 含 $\alpha_{ij}$），在区间上取格点（如 $0\le\alpha_{11}<1-\alpha_{12}$）。
3. $e_i$：$\ell_1$ 为 1→2 转换次数，$\ell_2$ 为 2→1 次数，$n_i$ 为状态 $i$ 数据点数，由 Result 12.3，后验 Beta($\gamma_{i1}+\ell_i$, $\gamma_{i2}+n_i-\ell_i$)。
4. $s_j$：逐个抽取。$P(s_j|\cdot)\propto\prod_{t=j}^nf(a_t|H)P(s_j|S_{-j})$，$P(s_j=i|S_{-j})=P(s_j=i|s_{j-1},s_{j+1})$ 由转移概率计算；设 $s_j=i$ 后递推 $h_t$（$t\ge j$），似然 $L(s_j=i)\propto\exp(f_{ji})$，$f_{ji}=\sum_{t=j}^n-\frac12[\ln h_t+a_t^2/h_t]$，$a_t=r_t-\beta_{s_t}\sqrt{h_t}$。 $P(s_j=1|\cdot)=\frac{P(s_j=1|s_{j-1},s_{j+1})L(s_j=1)}{P(s_j=1|s_{j-1},s_{j+1})L(s_j=1)+P(s_j=2|s_{j-1},s_{j+1})L(s_j=2)},$ 再用 U[0,1] 抽取。
备注：$e_1,e_2$ 小时 $s_j$ 与 $s_{j+1}$ 高度相关，联合抽取若干 $s_j$ 更有效，但需枚举的状态配置数随联合个数快速增长。（注：GARCH 的路径依赖使改变 $s_j$ 影响其后全部 $h_t$，因此似然要从 $j$ 累乘到 $n$，计算量 $O(n^2)$。）
例 12.6（GE 月度对数收益 %，1926-01 至 1999-12，888 个观测，图 12.14(a)）：
- GARCH-M 基准：$r_t=0.182\sqrt{h_t}+a_t$，$h_t=0.546+1.740h_{t-1}-0.775h_{t-2}+0.025a_{t-1}^2$（12.57），p 值均 <0.0006，Ljung–Box 无问题；风险溢价为正且显著。写成 $(1-1.765B+0.775B^2)a_t^2=0.546+(1-0.025B)\eta_t$（$\eta_t=a_t^2-h_t$），AR 多项式分解 $(1-0.945B)(1-0.820B)$，两实根模小于 1，无条件方差 $0.546/(1-1.765+0.775)\approx49.64$。
- Markov 转换：先验 $\beta_1\sim N(0.3,0.09)$，$\beta_2\sim N(1.3,0.09)$，$e_i\sim\mathrm{Beta}(5,95)$；初值 $e_i=0.1$，$s_1$ 等概率 Bernoulli、其余按初始转移概率依次生成，$\alpha_1=(1.0,0.6,0.2)'$，$\alpha_2=(2,0.7,0.1)'$。Griddy Gibbs 400 格点，范围 $\alpha_{i0}\in[0,6]$，$\alpha_{i1}\in[0,1]$，$\alpha_{i2}\in[0,0.5]$，约束 $\alpha_{i1}+\alpha_{i2}<1$。5000+2000 次迭代，用后 2000 次。
- 表 12.5 后验均值（标准差）：状态 1：$\beta_1$ 0.111 (0.043)，$e_1$ 0.089 (0.012)，$\alpha_{10}$ 2.070 (1.001)，$\alpha_{11}$ 0.844 (0.038)，$\alpha_{12}$ 0.033 (0.033)；状态 2：$\beta_2$ 0.247 (0.050)，$e_2$ 0.112 (0.014)，$\alpha_{20}$ 2.740 (1.073)，$\alpha_{21}$ 0.869 (0.031)，$\alpha_{22}$ 0.068 (0.024)；差值：$\beta_2-\beta_1$ 0.135 (0.063)（5% 显著），$e_2-e_1$ 0.023 (0.019)，$\alpha_{20}-\alpha_{10}$ 0.670 (1.608)，$\alpha_{21}-\alpha_{11}$ 0.026 (0.050)，$\alpha_{22}-\alpha_{12}$ −0.064 (0.043)。
- 结论：风险溢价差异显著；波动参数后验均值差异不显著，但后验分布形状不同（图 12.15、12.16 直方图）；图 12.17：状态 1 的持续性 $\alpha_{11}+\alpha_{12}$ 常触及边界 1.0，状态 2 不会；两状态期望持续期约 11 个月与 9 个月（$1/0.089$、$1/0.112$）；图 12.14(b) 为各观测处于状态 2 的后验概率。图 12.18：两模型拟合波动率形态相似且与平方收益一致，简单 GARCH-M 更平滑、估计波动更低。

12.10 预测（Forecasting）（PDF p.686–689）

MCMC 框架下预测很简单：在每次 Gibbs 迭代中用该次参数模拟预测期的实现。以单变量 SV（12.20–12.21）为例，第 $j$ 次迭代参数 $\beta_j,\alpha_j,\sigma_{v,j}^2$： $r_t=\beta_{0,j}+\beta_{1,j}x_{1t}+\cdots+\beta_{p,j}x_{pt}+a_t,\tag{12.58}\qquad \ln h_t=\alpha_{0,j}+\alpha_{1,j}\ln h_{t-1}+v_t,\ \mathrm{Var}(v_t)=\sigma_{v,j}^2.\tag{12.59}$ 步骤：抽 $v_{n+1}\sim N(0,\sigma_{v,j}^2)$ 由 (12.59) 算 $h_{n+1,j}$；抽 $\epsilon_{n+1}\sim N(0,1)$ 得 $a_{n+1,j}=\sqrt{h_{n+1,j}}\epsilon_{n+1}$（原文漏开方），由 (12.58) 算 $r_{n+1,j}$；对 $i=2,\dots,\ell$ 依次重复。（解释变量需已知或可依次预测。）
运行 $M+N$ 次迭代时只需对后 $N$ 次计算预测，得随机样本 $\{r_{n+1,j},\dots,r_{n+\ell,j}\}_{j=1}^N$ 与 $\{h_{n+1,j},\dots,h_{n+\ell,j}\}_{j=1}^N$；点预测为样本均值，样本标准差衡量预测误差。可用重要性抽样提高波动率预测效率（Gelman et al. 2003）。
例 12.7（例 12.3 续，S&P 500 月度对数收益 1962–1999，原点 1999-12）：表 12.6，1–5 步：
- 收益：GARCH 均为 0.66；SVM 为 0.53、0.78、0.92、0.88、0.84。
- 波动率：GARCH 17.98、18.12、18.24、18.34、18.42（逐渐升向无条件方差 $3.349/(1-0.086-0.735)=18.78$，此处 GARCH 参数为该样本期的估计，与式 12.26 数值不同）；SVM 19.31、19.36、19.35、19.65、20.13（2000+2000 次迭代）。
- SV 预测高于 GARCH：SV 的 MCMC 预测纳入了参数不确定性，GARCH 把参数当固定已知——这是 GARCH 相对期权隐含波动率倾向低估波动的原因之一。
备注：MCMC 实际给出波动率的预测分布，比点预测信息更多，可直接得到 VaR 所需分位数。

12.11 其他应用（Other Applications）（PDF p.689）

Zhang, Russell & Tsay (2008) 分析买卖报价的信息决定因素；McCulloch & Tsay (2001) 估计 IBM 交易数据的分层模型；Eraker (2001)、Elerian, Chib & Shephard (2001) 估计扩散方程；VaR 计算中自然地评估预测分布。关键问题不是能否用，而是效率。

第 12 章习题（PDF p.690–691）

12.1：$x\sim N(\mu,4)$，先验 $\mu\sim N(0,25)$，求单个观测下 $\mu$ 的后验（Result 12.1 直接应用：后验均值 $25x/29$，方差 $100/29$）。
12.2：回归误差为 AR($p$) 时，在共轭先验 $\beta\sim N(\beta_o,\Sigma_o)$，$\phi\sim N(\phi_o,A_o)$，$v\lambda/\sigma^2\sim\chi^2_v$ 下推导三个条件后验。
12.3：AR($p$) 中 $x_h,x_{h+1}$ 两个缺失值、联合正态先验，求其条件后验（二元正态）。
12.4：Ford 1965–2008 月度对数收益：GARCH 与 SV 模型比较（m-fsp6508.txt）。
12.5：Cisco 2001–2008 日对数收益 SV 模型，得 1 步波动率预测分布，并计算 100 万美元多头、概率 0.01 的次日 VaR。
12.6：Ford 与 S&P 指数二元 SV 模型，讨论波动关系并计算 Ford 的时变 β。
12.7：P&G 与价值加权指数 1965–2008：二元 SV、BEKK(1,1) 并比较。
12.8：30 年期抵押贷款利率与 3 月期国库券利率（1971-04 至 2009-09）：时间序列误差回归，再用 MCMC 重估并比较。

参考文献（PDF p.691–692）：Artigas & Tsay (2004)、Box & Tiao (1973)、Carlin & Louis (2000)、Carter & Kohn (1994)、Chang, Tiao & Chen (1988)、Chib, Nardari & Shephard (2002)、DeGroot (1970)、Dempster, Laird & Rubin (1977)、Elerian, Chib & Shephard (2001)、Eraker (2001)、Frühwirth-Schnatter (1994)、Gelfand & Smith (1990)、Gelman et al. (2003)、Geman & Geman (1984)、Hastings (1970)、Jacquier, Polson & Rossi (1994, 2004)、Jones (1980)、Justel, Peña & Tsay (2001)、Kim, Shephard & Chib (1998)、Liu, Wong & Kong (1994)、McCulloch & Tsay (1994a, 1994b, 2001)、Metropolis & Ulam (1949)、Metropolis et al. (1953)、Tanner (1996)、Tanner & Wong (1987)、Tierney (1994)、Tsay (1988)、Tsay, Peña & Pankratz (2000)、Zhang, Russell & Tsay (2008)。

第 12 章本章要点

MCMC：构造平稳分布为后验 $P(\theta|X)$ 的马尔可夫链；源于 EM 算法 → 迭代模拟 + 数据增广（缺失值、潜在波动、状态、异常指示都当作参数）。
Gibbs 抽样：依次从完全条件分布抽样；burn-in；高相关参数应分块联合抽取；收敛诊断需多初值重复。
共轭先验（正态–正态、正态–Gamma、Beta–Bernoulli、Gamma–Poisson/指数、Beta–负二项、逆卡方–方差）给出闭式条件后验；后验精度 = 先验精度 + 数据精度。
无闭式时：Metropolis（对称建议）、Metropolis–Hastings（一般建议，接受比含建议密度修正）、Griddy Gibbs（一元格点逆 CDF）。
应用模板：时序误差回归（准差分 + 共轭）、缺失值（$p+1$ 点回归求滤波值）、加性异常值（Bernoulli 指示 + 幅度的数据增广，后验 AO 概率）、单/多元 SV（逐点抽 $h_t$ 或 FFBS 联合抽取、KSC 七正态混合、杠杆效应的线性化）、Markov 转换 GARCH-M（逐个抽状态）。
MCMC 预测：逐迭代模拟未来路径，得到纳入参数不确定性的预测分布，可直接用于 VaR。

与量化交易的关联

风险管理：MCMC 预测分布天然包含参数不确定性，VaR/ES 更保守、更贴近现实；例 12.7 显示 GARCH 点估计会低估波动。
波动率建模与期权：SV 模型（含杠杆 $\rho\approx-0.4$）比 GARCH 更贴近连续时间期权定价模型（Heston 等的离散版本）；FFBS + KSC 混合是 SV 估计的标准工业算法，也可用于跳跃扩散。多元 SV 给出时变相关与时变 β（例 12.4、习题 12.6），可用于动态对冲。
数据清洗：Gibbs 异常值检测与缺失值插补可用于行情数据清洗（录入错误、坏点、停牌缺失），比"3σ 剔除"更有统计依据，且能同时估计模型。
区制识别：Markov 转换 GARCH-M 可识别高/低风险溢价区制，用于择时、风险预算切换与区制条件下的资产配置；状态后验概率可作为策略信号或风控开关（但需注意后验概率使用全样本平滑会引入前视偏差，实盘应使用滤波概率）。
贝叶斯回归：时序误差回归的 Gibbs 框架可推广到因子模型的贝叶斯估计、收缩先验（与 Black–Litterman 的贝叶斯思想一脉相承），在小样本因子检验中有用。
局限：计算量大，收敛诊断无保证；对高频、大截面问题需权衡效率（作者自己也指出"关键是效率"）。

索引及尾页（PDF p.693–707）

PDF p.693–697 为全书索引（Index），按字母顺序列出术语及页码（原书页码），例如 ACD 模型 255、ARCH 模型 115、Cholesky decomposition 400/458/517、Cointegration 91/428、Dynamic conditional correlation model 531、Forward filtering and backward sampling 655、Gibbs sampling 615、Kalman filter 563/593、Pairs trading 446、Stochastic volatility model 153/637、Threshold cointegration 444、Value at risk 326/546 等；还按数据集列出各示例数据的出处页码（如 IBM、S&P 500、GE、Cisco、Intel、BHP、Vale、恒生、日经等），以及 R、S-Plus、OX 命令索引（如 ca.jo 未列入、factanal 495、princomp 488、mfactor 500、mgarch 507、coint 439、VECM 439）。
PDF p.698–707 为空白页（文本抽取为空）。

第 8 章 多元时间序列分析及其应用（Multivariate Time Series Analysis and Its Applications）（接上一块）

8.6.3 协整检验：Johansen 检验（PDF p.457）（接上一块）

8.6.4 协整 VAR 模型的预测（PDF p.457）

8.6.5 实例：美国 3 个月与 6 个月国库券利率（PDF p.458–462）

8.7 门限协整与套利（Threshold Cointegration and Arbitrage）（PDF p.462–466）

8.7.1 多元门限模型（Multivariate Threshold Model）（PDF p.464）

8.7.2 数据（PDF p.465）

8.7.3 估计（PDF p.465–467）

8.8 配对交易（Pairs Trading）（PDF p.466–476）

8.8.1 理论框架（PDF p.466–468）

8.8.2 交易策略（PDF p.468–469）

8.8.3 简单示例：BHP 与 VALE（PDF p.469–476）

附录 A 向量与矩阵复习（Review of Vectors and Matrices）（PDF p.476–483）

附录 B 多元正态分布（Multivariate Normal Distributions）（PDF p.480–481）

附录 C 部分 SCA 命令（PDF p.481–482）

第 8 章习题（PDF p.482–484）

第 8 章（后半）本章要点

与量化交易的关联

推荐习题

第 9 章 主成分分析与因子模型（Principal Component Analysis and Factor Models）（PDF p.487 起）

本章引言（PDF p.487–488）

9.1 因子模型（A Factor Model）（PDF p.488–489）

9.2 宏观经济因子模型（Macroeconometric Factor Models）（PDF p.490–496）

9.2.1 单因子模型：市场模型（PDF p.490–494）

9.2.2 多因子模型（PDF p.494–496）

9.3 基本面因子模型（Fundamental Factor Models）（PDF p.496 起）

9.3.1 BARRA 因子模型（PDF p.497–）

9.3.2 Fama–French 方法（PDF p.502–503）

9.4 主成分分析（Principal Component Analysis）（PDF p.503–508）

9.4.1 PCA 理论（PDF p.503–505）

9.4.2 经验 PCA（Empirical PCA）（PDF p.505–508）

9.5 统计因子分析（Statistical Factor Analysis）（PDF p.509–518）

9.5.1 估计（PDF p.510–512）

9.5.2 因子旋转（Factor Rotation）（PDF p.512）

9.5.3 应用（PDF p.512–518）

9.6 渐近主成分分析（Asymptotic Principal Component Analysis）（PDF p.518–）

9.6.1 因子个数的选择（PDF p.519）

9.6.2 实例（PDF p.520–523）

第 9 章习题（PDF p.521–523）

第 9 章 本章要点

与量化交易的关联

推荐习题

第 10 章 多元波动率模型及其应用（Multivariate Volatility Models and Their Applications）（PDF p.525 起）

本章引言（PDF p.525–526）

10.1 指数加权估计（Exponentially Weighted Estimate）（PDF p.526–530）

10.2 若干多元 GARCH 模型（Some Multivariate GARCH Models）（PDF p.530–536）

10.2.1 对角向量化（DVEC）模型（PDF p.530–533）

10.2.2 BEKK 模型（PDF p.533–536）

10.3 重参数化（Reparameterization）（PDF p.536–）

10.3.1 使用相关系数（Use of Correlations）（PDF p.536–537）

10.3.2 Cholesky 分解（PDF p.537–）

10.4 二元收益的 GARCH 模型（GARCH Models for Bivariate Returns）（PDF p.541–551）

10.4.1 常相关模型（Constant-Correlation Models）（PDF p.541–545）

10.4.2 时变相关模型（Time-Varying Correlation Models）（PDF p.545–551）

10.4.3 动态相关模型（Dynamic Correlation Models, DCC）（PDF p.551–557）

10.5 高维波动率模型（Higher Dimensional Volatility Models）（PDF p.557–）

10.6 因子–波动率模型（Factor–Volatility Models）（PDF p.563–566）

10.7 应用：多资产 VaR（Application）（PDF p.566–568）

10.8 多元 t 分布（Multivariate t Distribution）（PDF p.568–569）

附录：估计说明（Appendix: Some Remarks on Estimation）（PDF p.569–574）

第 10 章习题（PDF p.574–575）

第 10 章 本章要点

与量化交易的关联

推荐习题

第 11 章 状态空间模型与卡尔曼滤波（State-Space Models and Kalman Filter）（PDF p.577 起）

本章引言（PDF p.577）

11.1 局部趋势模型（Local Trend Model）（PDF p.578–）

11.1.1 统计推断（Statistical Inference）（PDF p.581–582）

11.1.2 Kalman 滤波（Kalman Filter）（PDF p.582–）

11.1.3 预测误差的性质（PDF p.584–586）

11.1.4 状态平滑（State Smoothing）（PDF p.586–590）

11.1.5 缺失值（Missing Values）（PDF p.590–591）

11.1.6 初始化的影响（Effect of Initialization）（PDF p.591–592）

11.1.7 估计（Estimation）（PDF p.592）

11.1.8 所用 S-Plus 命令（PDF p.592–596）

11.2 线性状态空间模型（Linear State-Space Models）（PDF p.596–597）

11.3 模型转换（Model Transformation）（PDF p.597–）

11.3.1 时变系数 CAPM（PDF p.597–599）

11.3.2 ARMA 模型（PDF p.600–）

11.3.3 线性回归模型（PDF p.606–608）

第 8 章多元时间序列分析及其应用（Multivariate Time Series Analysis and Its Applications）（接上一块）

第 9 章主成分分析与因子模型（Principal Component Analysis and Factor Models）（PDF p.487 起）

第 9 章本章要点

第 10 章多元波动率模型及其应用（Multivariate Volatility Models and Their Applications）（PDF p.525 起）

第 10 章本章要点

第 11 章状态空间模型与卡尔曼滤波（State-Space Models and Kalman Filter）（PDF p.577 起）

第 11 章本章要点

第 12 章马尔可夫链蒙特卡罗方法及其应用（Markov Chain Monte Carlo Methods with Applications）（PDF p.633 起）

第 12 章本章要点