时间序列分析方法之卡尔曼滤波

第十三章 卡尔曼滤波
在本章中，我们介绍一种被称为卡尔曼滤波的十分有用的工具。

卡尔曼滤波的根本思想是将动态系统表示成为一种称为状态空间表示的特殊情形。

卡尔曼滤波是对系统线性投影进行序列更新的算法。

除了一般的优点以外，这种算法对计算确切的有限样本预测、计算Gauss ARMA 模型确实切似然函数、估计具有时变参数的自回归模型等，都提供了重要方法。

§13.1 动态系统的状态空间表示
我们已经介绍过一些随机过程的动态表示方法，下面我们在以前的假设根底上，继续分析动态系统的表示方法。

13 继续使用的假设
假设表示时刻观测到的n 维随机向量，一类非常丰富的描述动态性的模型可以利用一些可能无法观测的被称为状态向量(state vector)的r 维向量表示，因此表示动态性的状态空间表示(state-space representation)由以下方程系统给出：
状态方程(state model) (13.1) 量测方程(observation model) (13.2)
这里，和分别是阶数为，和的参数矩阵，是的外生或者前定变量。

方程(13.1)被称为状态方程(state model)，方程(13.2)被称为量测方程(observation model)，维向量和维向量都是向量白噪声，满足：
(13.3) (13.4)
这里和是和阶矩阵。

假设扰动项和对于所有阶滞后都是不相关的，即对所有和，有： (13.5)
t x 是外生或者前定变量的假定意味着，在除了包含在121,,,y y y --t t 内的信息以外，t
x 没有为s t +ξ和s t +w ( ,2,1,0=s )提供任何新的信息。

例如，t x 可以包括t y 的滞后值，也可以包括与τξ和τw (任意τ)不相关的变量。

方程系统中方程(13.1)至方程(13.5)可以表示有限观测值的序列},,,{21T y y y ，这时需要状态向量初始值1ξ。

假设1ξ与t v 和t w 的任何实现都不相关：
0ξv =')(1
t E ，对任意T t ,,2,1 = (13.6) 0ξw =')(1
t E ，对任意T t ,,2,1 = (13.7) 状态方程(13.1)说明，t ξ可以表示成为},,,,{321t v v v ξ 的线性函数：
1122221ξF v F v F v F v ξ----+++++=t t t t t t ，T t ,,3,2 = (13.8)
因此，方程(13.6)和方程(13.3)意味着t v 与所有ξ的滞后值都是不相关的：
0ξv =')(τ
t E ，1,,2,1 --=t t τ (13.9) 类似地，可以得到：
0ξw =')(τ
t E ，T ,,2,1 =τ (13.10) 0w ξH x A w y w t ='+'+'='])([)(t t t t E E τ
，1,,2,1 --=t t τ (13.11)
0y v =')(τt E ，1,,2,1 --=t t τ (13.12)
上述系统是相当灵活的，它的一些结论也可以推广到t v 与t w 相关的系统中，而且系数矩阵),,,,(R H A Q F 也可以是时间的函数。

如果我们仅仅关注到上述系统的根本形式，那么下面的论述将是十分清晰的。

13.1.2 状态空间表示的例子 考虑一元)(p AR 过程：
111211)()()(++--++-++-+-=-t p t p t t t y y y y εμφμφμφμ
⎩
⎨
⎧≠==ττ
σεετt t E t ,0,)(2 这个)(p AR 过程可以表示成为下面的状态空间模型形式：
状态方程(p r =)
⎥
⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡+⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡---⎥⎥⎥
⎥⎥⎥⎦
⎤
⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎢⎣⎡---++---+-+00010
00010000111112
1
21
t p t t t p p p t t t y y y y y y εμμμφφφφμμμ (13.13) 量测方程：
⎥⎥⎥⎥⎥⎦
⎤
⎢⎢⎢⎢⎢⎣⎡---+=+--μμμμ11]001[p t t t t y y y y (13.14) 对应地，我们指定：
⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡---=+--μμμ11p t t t t y y y ξ，⎥⎥⎥
⎥⎥⎥⎦
⎤
⎢⎢⎢
⎢⎢⎢⎣⎡=-01000010000112
1 p p φφφφF ，⎥
⎥⎥⎥
⎦⎤⎢⎢⎢⎢⎣⎡=00 t t εv ，⎥
⎥
⎥⎥
⎥
⎦⎤
⎢⎢⎢⎢
⎢⎣⎡=000000
002 σQ
这里变量和参数矩阵对应为：
t t y =y ，μ='A ，1=t x ，[]001 ='H ，0=t w ，0=R
注意到这里的状态方程只是一个一阶向量自回归方程，量测方程只是一个简单的等式。

因此，我们已经看到，状态空间表示只是总结)(p AR 过程的另外一种方式。

将)(p AR 过程表示成为这种方式的原因在于，这样可以获得归纳)(p AR 过程动态性的适宜方式，这是我们对任何系统状态空间表示感兴趣的根本原因。

另外一个例子是，我们考虑一元)1(MA 过程：
1-++=t t t y εθεμ
对应地，它可以表示成为状态空间模型形式为： 状态方程(2=r )：
⎥⎦⎤⎢⎣⎡+⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣
⎡=⎥⎦⎤⎢⎣⎡+-+00100111t t t t t εεεεε
量测方程(1=n )： []⎥⎦⎤
⎢⎣⎡+=+t t t y εεθμ11
这里：
⎥⎦⎤⎢⎣⎡=-1t t t εεξ，⎥⎦⎤⎢⎣⎡=0100F ，⎥⎦⎤
⎢⎣⎡=+01t t εv ，⎥⎦
⎤
⎢⎣⎡=0002σQ t t y =y ，μ='A ，1=t x ，[]θ1='H ，0=t w ，0=R 将给定系统表示成为状态方程的方式有多种。

例如，可以将)1(MA 过程表示成为下面类
型的状态空间模型：
状态方程(2=r )：
⎥⎦⎤⎢⎣⎡+⎥⎦⎤⎢⎣⎡+⎥⎦⎤⎢⎣⎡=⎥⎦⎤⎢⎣⎡+++-++111110010t t t t t t t t εθεεθεθεεθεθε
量测方程(1=n )：
[]⎥⎦
⎤
⎢⎣⎡++=-t t t t y εθεθεμ101 显然上面的)1(MA 过程、两种状态空间模型表示都是具有相同特征的过程表示，这三种表示都具有相同的预测和相同的似然函数值，也就无须讨论哪一种方式更为适宜。

更一般地，一元),(q p ARMA 模型可以通过定义}1,max {+≡q p r 进行状态空间模型表示：
1
122112211)
()()(+-------+++++-++-+-=-r t r t t t r t r t t t y y y y εθεθεθεμφμφμφμ 5)
这里的参数约束是：当p j >时，0=j φ；当q j >时，0=j θ。

考虑以下状态空间模型表示为： 状态方程(}1,max {+≡q p r )：
⎥
⎥⎥⎥⎦
⎤
⎢⎢⎢⎢⎣⎡+⎥⎥⎥
⎥⎥
⎥⎦
⎤
⎢⎢⎢
⎢⎢
⎢⎣⎡=+-+0001000010000111211
t t r r t εφφφφξξ6) 量测方程(1=n )：
[]t r t y ξ1211-+=θθθμ (13.15)
为了验证方程67)表示了系统与方程(13.15)一致，假设t j ξ表示向量t ξ的第j 个元素，因
此状态方程的第2行表示： t t 112ξξ=+
第3行说明： 11213-+==t t t ξξξ 更一般地，第j 行表示：
1111+-+=t j t j L ξξ
因此状态方程的第1行意味着：
1,1123211,1)(+-++++++=t t r r t L L L εξφφφφξ
或者：
11,1221)1(++=----t t r r L L L εξφφφ 8)
量测方程说明：
t r r t L L L y ,111221)1(ξ--+++++=θθθμ 9)
在方程9)两端乘以算子多项式)1(221r r L L L φφφ---- ，并利用方程8)，可以得到：
t
r r t r r L L L y L L L εθθθμφφφ)1()()1(11221221--++++=
-----
这就是原来的),(q p ARMA 模型，即方程5)。

状态空间形式是描述随机过程的和，或者测量误差结果的模型的非常适宜的方式。

例如，Fama 和Gibbons (1982)开始着手研究事前实际利率(ex ante real interest rate )行为 (事前实际利率是名义利率t i 减去预期通货膨胀率e t π)。

由于经济计量学家通过证券市场推断的预期通货膨胀率的数据，因此这个变量不是可以观测的。

因此在这种应用中状态变量是一个标量，即：μπξ--=e t t t i ，这里μ表示平均事前实际利率。

Fama 和Gibbons (1982)假设事前实际利率服从)1(AR 过程：
11+++=t t t v ξφξ (13.20)
经济计量学家可以观测到事后实际利率(名义利率t i 减去真实通货膨胀率t π)，这可以表示为：
t t t t e t e t t t t w i i ++=-+-=-ξμππππ)()( (13.21)
这里)(t e t t w ππ-≡是人们预测通货膨胀率时的误差。

如果人们以最优的方式形成通货膨胀率预测，那么t w 与自身的滞后值和事前实际利率是无关的。

因此方程(13.20)和方程
(13.21)是状态空间模型，这里1==n r ，φ=F ，t t t i y π-=，
μ='t x A ，1=H ，t e t t ππ-=w 。

状态空间模型框架的另外一个有趣例子是Stock 和Waston (1991)的研究，他们假设存在
表示经济周期状态的不可观测变量t C 。

假设),,,(21nt t t y y y 是n 个可以观测的宏观经济变量，每个都受到经济周期的影响，并且具有与i j y jt ≠,中移动不相关的奇异成分(表示为
t i χ)。

如果经济周期和每个奇异成分可以利用一元)1(AR 过程描述，那么]1)1[(⨯+n 维状态
向量是：
⎥⎥⎥⎥⎥
⎥⎦
⎤⎢⎢⎢⎢⎢⎢⎣⎡=t n t t t t C χχχ 21ξ 2)
该状态变量具有的状态方程为： ⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡+⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣
⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡++++++++112111,21112111000000000000t n t t t C t n t t t C C C C t n t t t v v v v C C χχχφφφφχχχ3) 量测方程为： ⎥⎥
⎥
⎥⎥
⎥
⎦
⎤
⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡+⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡t n t t t n
n t n t t t C y y y y χχχγγγγμμμμ 21321321321100000010001 4)
因此，参数i γ描述第i 个序列对经济周期反响的敏感性。

为了出现和描述p 2)中的t C 和t i χ替换为)1(⨯p 阶向量),,,(11'+--p t t t C C C 和),,,(11'+--p t i t i t i χχχ ，这时t ξ是
]1)1[(⨯+p n 3)中的标量i φ需要利用)(p p ⨯阶矩i F 阵替换，该矩阵结构与方程(13.13)类似。

还需要在量测方程(13.24)中H '的列中参加阶数为)]1(([-⨯p n 的零子块。

§13.2 卡尔曼滤波的推导
卡尔曼滤波是估计状态空间模型的重要方法，也是应用广泛的参数估计方法。

下面我们介绍卡尔曼滤波的有关公式。

13.2.1 卡尔曼滤波的回忆 Overview of the Kalman Filter
考虑上述讨论的状态空间模型的一般形式，为了方便，我们将使用的一些关键方程在这里重复表示如下： )
1()
1()()1(1
1⨯⨯⨯⨯+=++r r r r r t t t v ξF
ξ
(1)()(1)()(1)(1)
t t t
n n k k n r r n ''=++⨯⨯⨯⨯⨯⨯t y A x H ξw
⎪⎩⎪⎨⎧≠=='⨯τ
ττ
t t E r r t ,,)(0v v Q ，
⎪⎩⎪⎨⎧≠=='⨯τ
ττt t E n n t ,,)(0w w R 假设我们已经得到了观测值T y y y ,,,21 ，T x x x ,,,21 ；一个最终目标是基于这些观
测值估计系统的所有未知参数。

但是，目前我们暂时假设参数矩阵),,,,(R H A Q F 的特定数值都是确定性的。

如何估计这些参数在后面的内容中讨论。

卡尔曼滤波具有多种应用。

它的根本动因是作为一种计算状态向量基于时刻t 观测到的数据进行最小二乘预测的算法。

)|(ˆˆ1|1t t t t E Υξξ++≡ (13.25)
这里：
1111
{,,,,,,,}t t t t t --'''''''Y =y y y x x x 这里)|(ˆ1t t E Υξ+表示1+t ξ基于t
Υ和常数的线性投影。

卡尔曼滤波是采用叠代算法计算这些预测的，按顺序分别产生0|1ˆξ，1|2ˆξ，……，1|ˆ-T T ξ。

与这些预测有关的是均方误差矩
阵，可以由一个r r ⨯阶矩阵表示：
])ˆ)(ˆ([|11|11|1'--≡+++++t t t t t t t t E ξξξξP (13.26)
13.2.2 叠代的开始
叠代首先从0|1ˆξ开始，0
|1ˆξ表示在没有y 和x 观测值的根底上对1ξ的预测。

这就是1ξ的无条件均值：
)(ˆ1
0|1ξξE = 与此相关的MSE 为：
])()((([1111|1'--=ξξξξP E E E t
例如，对)1(MA 系统的状态空间表示，状态向量为：
⎥⎦
⎤⎢⎣⎡=-1t t t εεξ
这时有：
⎥
⎦
⎤⎢⎣⎡=⎥⎦⎤⎢⎣⎡=-00ˆ10|1t t E εεξ []⎥⎦
⎤⎢⎣⎡=⎪
⎪⎭⎫ ⎝⎛⎥⎦⎤⎢⎣⎡=220101|100σσεεεεE t
P ，)(2
2t E εσ= 更一般地，如果矩阵F 的特征根都落在单位圆内，那么状态方程表示的过程}{t ξ是协方差平稳的。

因此，对状态方程两端取无条件数学期望，可以得到：
)()(1t t E E ξF ξ=+
由于过程}{t ξ是协方差平稳的，那么有：
0ξF I =-)()(t r E
由于矩阵F 没有单位根，因此矩阵)(F I -r 是非奇异的，因此这个方程存在唯一零解，也就是有：0ξ=)(t E 。

那么}{t ξ的无条件方差也可以类似地得到，取矩阵的转置并取数学期望，可以得到(由于存在正交性，下面的交叉项的数学期望为零)：
)()()])([()(111111++++++'+''='+''+='t t t t t t t t t t E E E E v v F ξξF v F ξv ξF ξξ
假设矩阵Σ表示}{t ξ的协方差矩阵，那么有：
Q F ΣF Σ+'=
这个方程的解可以表示为：
)vec()]([)vec(12Q F F I Σ-⊗-=r
因此，一般情况下，如果矩阵F 的特征根都落在单位圆内，因此卡尔曼滤波的叠代可以
从0ξ=0
|1ˆ和0|1P 开始，这里的0|1P 表示成为列向量可以从下式得到： )vec()]([)vec(10|12Q F F I P -⊗-=r
如果矩阵F 的局部特征根落在单位圆上或者单位圆外，或者1ξ无法从状态方程中获得，
这时0
|1ˆξ可以利用分析者对初始1ξ的最优猜想来替代，而0|1P 是归纳这种预测置信区间的正定矩阵。

0|1P 中对角线上比拟大的数值对应着对1ξ真实取值较高的非确定性。

13.2.3 预测t y
给定开始的初值0|1ˆξ和0|1P ，下一步是计算下一个时期类似的数量1|2ˆξ和1|2P 。

由于计算对T t ,,3,2 =都具有相同的形式，因此我们讨论在时刻t 的一般形式。

给定1|ˆ-t t ξ和1|-t t P ，目的是计算t t |1ˆ+ξ和t t |1+P 。

首先，我们需要注意到，我们假设除了包含在1-Y t 内的信息以外，t x 不再包含关于t ξ的信息，因此有：
1|11ˆ)|(ˆ),|(ˆ---==t t t t t t t E E ξΥξΥx ξ
下面我们考虑对t y 的预测：
),|(ˆˆ1
1|--≡t t t t t E Υx y y 注意到根据状态方程，可以得到：
t
t t t t E ξH x A ξx y '+'=),|(ˆ 因此，根据投影的叠代定律，有：
1|11|ˆ),|(ˆˆ---'+'='+'=t t t t t t t t t E ξH x A Υx ξH x A y
这个预测的误差为：
t t t t t t t t t t t t t w ξξH ξH x A w ξH x A y y +-'='-'-+'+'=----)ˆ(ˆˆ1|1|1|
因此预测的MSE 为：
][])ˆ()ˆ([])ˆ)(ˆ[(1|1|1|1|t
t t t t t t t t t t t t t E E E w w H ξξξξH y y y y '+'--'='------ 由于0ξξw ='--])ˆ([1
|t t t t E ，因此上式中交叉项为零。

这个正交条件需要根据假设和投影性质加以验证。

这时可以将MSE 表示为：
R H P H y y y
y +'='-----]])ˆ)(ˆ[(1|1|1|t t t t t t t t E 13.2.4 关于t ξ推断的更新
给定开始的初值0|1ˆξ和0|1P ，下一步是计算下一个时期类似的数量1
|2ˆξ和1|2P 。

由于计算对
13.2.5 产生1+t ξ的预测
给定开始的初值0|1ˆξ和0|1P ，下一步是计算下一个时期类似的数量1
|2ˆξ和1|2P 。

由于计算对
13.2.6 归纳和注释
给定开始的初值0|1ˆξ和0|1P ，下一步是计算下一个时期类似的数量1
|2ˆξ和1|2P 。

由于计算对
§13.3 基于状态空间表示的预测 Forecasts Based on the State-Space Representation §13.4 参数的极大似然估计 Maximum Likelihood Estimation of Parameters §13.5 稳态卡尔曼滤波 The Steady-State Kalman Filter 13.5.1 卡尔曼滤波的收敛性 §13.6 光滑 Smoothing
§13.7 利用卡尔曼滤波的统计推断 Statistical Inference with Kalman Filter 我们已经介绍过一些随机过程的动态表示方法，下面我们在以前的假设根底上，继续分析动
§13.8 时变参数 Time Varying Parameter
我们已经介绍过一些随机过程的动态表示方法，下面我们在以前的假设根底上，继续分析动。