《2022年时间序列分析讲义Kalman滤波定义 .pdf》由会员分享,可在线阅读,更多相关《2022年时间序列分析讲义Kalman滤波定义 .pdf(11页珍藏版)》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、时间序列分析方法讲义第 11章Kalman 滤波1 第十一章 Kalman 滤波本章我们介绍由R. E. Kalman 命名的一种十分有用的工具。其核心思想是将动态系统表示成为状态空间表示(state-space representation) 。Kalman 滤波是按顺序更新系统线性投影的一种算法。该算法对于计算例如时变系数回归等具有重要应用。11.1 动态系统的状态空间表示1. 状态空间模型和状态空间表示假设ty是一个1n维的在时刻t 可以观测的向量。很多情形下,ty可以利用一个不可观测的1r随机向量t表示,该向量称为状态向量(state vector)。则变量ty服从的动态过程可以利用状
2、态空间表示为:11tttFv(11.1) ttttyA xH w(11.2) 这里矩阵F、A和H分别是()rr、()nk和()nr维的参数矩阵。tx是(1)k维的外生或前定变量。方程 (11.1)被称为状态方程(state equation);方程 (11.2)被称为观测方程(observation equation)。这里(1)r维误差向量tv和(1)n维误差向量tw均为向量白噪声过程:,(),ttEtQv v0(11.3) ,(),ttEtRw w0(11.4) 其中Q和R分别是rr和nn维的正定矩阵。我们继续假设状态方程误差tv和观测方程误差tw在所有时间间隔上是不相关的,即对所有的t和
3、:()tE v w0(11.5) 这里所假设的tx是(1)k维的外生或前定变量的含义是,除了包含在111,ttyyy中的信息以外,tx没有提高任何关于ts和t sw(1, 2,s)的信息。例如,tx中可以包含变量ty的滞后变量,或者包含与ts和tsw(1, 2,s)不相关的解释变量。如果假设初始状态向量1是可以得到, 则上述方程 (11.1)(11.5)可以典型地用来描述有限个样本观测值过程:12,Tyyy。(1) 假设1与tv和tw的任何实现都无关,则有:对于1, 2,tT1()tE v 0,1()tE w 0(11.6) 状态方程 (11.1)意味着t可以表示为:2211221tttttt
4、vFvF vFvF,2, 3,tT因此方程 (11.6)意味着tv与滞后的t无关:()tE v 0,1,2,1tt(11.7) ()()ttEEw yw A xH w0,1,2,1tt(11.8) ()tE v y0,1,2,1tt(11.9) 因此,上述方程 (11.1)(11.5)表示的模型系统是相当灵活的。我们还可以将其推广为tv与tw相关的情形。另外,各种参数矩阵(F、Q、A、H或R)可以是时间的函数。但是,为了更为清楚和简洁,我们集中考虑状态空间模型的基本形式,即方程(11.1)(11.5)表示的名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - -
5、 - - - - - - 名师精心整理 - - - - - - - 第 1 页,共 11 页 - - - - - - - - - 时间序列分析方法讲义第 11章Kalman 滤波2 模型系统。2. 状态空间表示的例子例 11.1 AR(p)过程的状态空间表示考虑一个单变量的AR(p)过程:112111()()()tttptptyyyy(11.10) 2,()0,ttEt这个模型可以表示成为如下状态空间形式:状态方程 (rp) 1211112110000010000010ppttttttptpyyyyyy(11.11) 观测方程 (1n):11100ttttpyyyy(11.12) 与状态空间表
6、示的标准形式对照,相应的参数矩阵为:11ttttpyyy,121100001000010ppF1100ttv,200000000Qttyy,A,1tx,100H,0tw,0R上述状态方程是一个一阶向量差分方程,而观测方程是一个简单的恒等式。将自回归过程表示成为状态空间表示的主要原因是说明状态空间表示是归纳模型动态性质的一种重要表示方式和归纳方式。例 11.2 MA(1)过程的状态空间表示我们考虑一个一元MA(1)过程:1ttty(11.13) 这个模型可以表示成为如下状态空间形式。状态方程 (2r):名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - -
7、- - - - - 名师精心整理 - - - - - - - 第 2 页,共 11 页 - - - - - - - - - 时间序列分析方法讲义第 11章Kalman 滤波3 11100100ttttt(11.14) 观测方程 (1n):11ttty(11.15) 与状态空间模型表示的符号对应有:1ttt,0010F,110ttv,2000Q,ttyy,A,1tx,1H,0tw,0R应该注意到, 一个动态模型可以具有多种状态空间表示,例如 MA(1)过程还可以表示成为:状态方程 (2r):111110100tttttttt(11.16) 观测方程 (1n):110tttty(11.17) 显然
8、,MA(1)过程和相应的两种状态空间表示都刻画的是同一个过程,无论是哪个表示,我们都可以获得相同的预测和似然函数值,因此我们可以选择最方便处理的形式。更一般地,一个单变量ARMA(p, q)模型可以通过定义max( ,1)rp q来表示成为状态空间模型形式:112211221()()()tttrtrtttrtryyyy(11.18) 此处参数限制为:0j,jp;0j,jq;考虑下述状态空间表示:状态方程 (max(,1)rp q): 1211110000010000010rrttt(11.19) 观测方程 (1n):111trty(11.20) 为了验证方程(11.19)和方程 (11.20)
9、与方程 (11.18)表示同一个随机过程,假设jt表示向量t的第 j 个元素,因此状态方程的第二行表明:2,11tt名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 3 页,共 11 页 - - - - - - - - - 时间序列分析方法讲义第 11章Kalman 滤波4 状态方程的第三行表明:3,12,1,1ttt更一般地有:1,111jj ttL因此状态方程的第一行意味着:11,11211()rtrttLL或者:2121,11(1)rrttLLL利用观测方程,可以得到:21121
10、1(1)rtrtyLLL对上式两端乘以滞后算子多项式212(1)rrLLL,得到:22112121(1)()(1)rrrtrtLLLyLLL显然上式就是方程(11.18),因此可见状态空间表示与自回归移动平均过程是一致的。状态空间形式也非常适合模型化处理随机过程的和,或者描述测量误差序列。下面我们给出两个具体的例子。例 13.1 事前实际利率方程Fama and Gibbons (1982) 想研究事前 (ex ante real interest rate) 实际利率的行为,事前实际利率就是名义利率ti减去预期通货膨胀率et。因为经济计量学家没有从证券市场推断出的通货膨胀率数据,因此这个变量
11、是不可观测的。因此, 在这个例子中的状态变量是一个标量:ettti,这里表示平均事前实际利率。Fama and Gibbons 假设事前实际利率服从 AR(1) 过程:11tttv注意到,经济计量学家可以观测到事后实际利率(ex post real interest rate)( 名义利率ti减去实际通货膨胀率t),它可以表示成为:()()eettttttttiiw这里()etttw是人们预测通货膨胀率时的误差。如果经济个体以最优方式形成预测,则这个预测误差()etttw应该与自身滞后值和事前实际通货膨胀率不相关。因此,上述两个方程构成了状态空间表示,其中对应的参数和变量为:1rn,F,ttt
12、yi,tA x,1H,etttw例 13.2 经济周期模型Stock and Watson (1991) 给出了状态空间模型的另一个有趣的应用。该研究假设存在一个代表经济周期状态的不可观测变量tC,一组由n 个可观测的宏观经济变量构成的向量为12(,)ttntyyy,其中每个分量都受到经济周期阶段的影响。这些分量当中都包含着自己的奇异成分,利用it表示,该奇异成分与其他jty,ij中的变化不相关。如果经济周期和每个奇异成分都可以利用单变量AR(1)过程描述,则此时的(1) 1n维状态向量为:名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - -
13、 - - 名师精心整理 - - - - - - - 第 4 页,共 11 页 - - - - - - - - - 时间序列分析方法讲义第 11章Kalman 滤波5 12ttttntC状态方程为:,111,111112,12122,11000000000000C ttCtttttttn tntnntvCCvvv观测方程为:11122213332100010000001ttttttntnnntyCyyy因此, 系数i描述第 i 个时间序列对经济周期的灵敏性。为了出现 p 阶自回归,Stock and Watson (1991) 将上述方程中的周期经济周期状态和起义成分分别换为1p阶向量:11(,
14、)tttpC CC,,1,1(,)iti ti tp其余的系数矩阵也可以做出相应调整。11.2 卡尔曼滤波的推导1. 卡尔曼滤波的概述我们考虑一般形式的状态空间模型,为了方便,我们在此重新表述一下:(1)()(1)(1)11rrrrrtttFv(11.21) (1)()(1)()(1)(1)nnkknrrnttttyA xH w(11.22) (),(),rrttEtQv v0(11.23) (),(),tnntEtRw w0(11.24) 假设我们已经观测到12,Tyyy、12,Tx xx,目的是利用这些观测值估计系统中的任意未知参数。对目前来说,我们一直假设我们确定性地知道矩阵F、Q、A、
15、H和R的特定数值,后面我们将给出估计这些参数的具体方法和过程。卡尔曼滤波具有多种用途。这里的使用动机是将其作为一种基于时期t 所获得数据来计名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 5 页,共 11 页 - - - - - - - - - 时间序列分析方法讲义第 11章Kalman 滤波6 算状态向量的线性最小二乘预测的算法,我们所要预测的线性函数为:1|1?(|)ttttE(11.25) 这里:1111(,;,)tttttyyyx xx其中1?(|)ttE 是1t基于常数和t
16、的线性投影。卡尔曼滤波就是叠带预测这些投影,按照顺序生成序列1|0?、2|1?、|1?T T与这些预测相关的是均方误差矩阵,可以利用下面的rr阶矩阵表示:1|11|11|?()() ttttttttEP(11.26) 2. Starting the Recursion 迭代的开始开始迭代时我们需要计算1|0?,这是在没有关于y和x观测值的情形下对1的预测。此时最好的预测就是1的无条件均值。1|01?()E与此相关的MSE 为:1|01111()() EEEP例如,对 MA(1)过程的状态空间表示而言,状态向量为:1ttt对此有:1|010?0ttE211|0102000EP,22()tE更一般
17、地, 如果矩阵F的特征值都在单位圆内,则方程 (11.21)描述的自回归过程是协方差平稳过程。我们可以对这个方程两端求数学期望得到无条件均值:1ttE F 上述方程具有唯一解:tE 0类似地,我们可以求出t的无条件方差:111111()() ()ttttttttttEEEE FvFvF Fvv假设是t的无条件方差矩阵,则有:F FQ这个方程的解为:21vec( )()vec( )rIFFQ因此在一般情形下,如果假设矩阵F的特征值都在单位圆内,则卡尔曼滤波的迭代可以从零向量1|0?0和rr矩阵1|0P开始,这个矩阵的所有列拉直表示成为列向量为:21vec( )()vec( )rIFFQ。名师资料
18、总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 6 页,共 11 页 - - - - - - - - - 时间序列分析方法讲义第 11章Kalman 滤波7 如果矩阵F的特征值都在单位圆外或者单位圆上,如果初始值1不被认为是从该过程中提取的任意值,则1|0?可以取为人们认为最可能的猜测值,而矩阵1|0P则是一个归纳这种猜测置信程度的正定矩阵,其对角线元素越大,则越意味着关于这些数值猜测的不确定性越大。3. 预测ty给定初始值1|0?和1|0P,下一步是计算下一个阶段的这些类似数值2|1?和2
19、|1P。对2,3,tT的计算具有完全一致的形式,因此我们以一般的形式进行讨论:给定|1?t t和|1t tP,我们来计算1|?tt和1|ttP。首先,我们需要注意到, 我们假设tx中除了包含在1t中信息以外, 没有再提供关于t的信息。此时有:11|1?(|,)(|)tttttt tEE x其次,我们考虑预测ty的值:|11?(|,)t ttttEyy x根据方程 (11.22),我们得到:?(| ,)tttttE yx A x +H 因此,根据投影的迭代定律,有:|11|1?(| ,)t ttttttt tEyA x +H xA x + H (11.27) 预测误差为:|1|1|1?()tt
20、tttttt ttt ttyyA x +H wA xH H w该预测误差的均方期望(MSE):|1|1|1|1?()() ()() tt ttt ttt ttt tttEEEyyyyH Hw w其中交叉项为零(对此我们可以进一步验证):|1?() ttt tE w 0利用上述结算结果,可以得到均方期望(MSE):|1|1|1?()() tt ttt tt tE yyyyH PHR(11.28) 4. 更新关于t的推断下面关于t当前值的推断可以在观测值ty的基础上获得的:|1?(|,)(|)t tttttttEE y x(11.29) 这可以利用线性投影的更新获得:1|1|1|1|1|1|1?(
21、)() ()() ()t tt ttt ttt ttt ttt ttt tEE yyyyyyyy但是:|1|1|1|1|1|1|1?()() ()() ?()()tt ttt ttt ttt tttt ttt tt tEEEyyHwHPH(11.30) 于是可以得到:名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 7 页,共 11 页 - - - - - - - - - 时间序列分析方法讲义第 11章Kalman 滤波8 1|1|1|1|1?() ()t tt tt tt tttt
22、tPH H PHRyA xH (11.31) 这个更新投影的MSE 记为| t tP,可以表示为:|1|1|1|11|1|1|1|11|1|1|1|1?()() ?()() ()() ?()() ()() ()t ttt ttt ttt ttt ttt ttt ttt ttt ttt ttt tt tt tt tt tEEEEEP yyyyyyyy PPH H PHRH P(11.32) 5. 形成1t的预测下一步,我们使用状态空间方程获得1t的预测:1|11|?(|)?(|)(|)?ttttttttt tEEEFvF0将方程 (11.31)代入上式,得到:11|1|1|1|1?() ()tt
23、t tt tt tttt tFFPH H PHRyA xH (11.33) 上述系数被称为收益矩阵(gain matrix ),表示为:1|1|1()tt tt tKFPH H PHR可以将方程 (11.33)简化为:1|1|1?()ttt ttttt tFKyA xH (11.34) 进一步,我们可以获得这个预测的MSE 为:1|1|1|1|1|1|1|11|?()() ?()() ?()() ttttttttttt tttt ttt ttt tttt tEEEEPFvFFvFFFvvFP FQ(11.35) 将方程 (11.32)代入上式,得到:11|1|1|1|1()ttt tt tt
24、tt tPF PPH H PHRH PFQ(11.36) 6. 归纳和评注下面我们将获得的结论归纳和总结一下。卡尔曼滤波的迭代是从下面1的无条件均值和无条件方差开始的:1|01?()E,1|01111()() EEEP典型的计算初始值为:1|0?0,211|0vec()()vec( )rP IFFQ然后我们按照下面两个公式进行迭代计算:11|1|1|1|1?() ()ttt tt tt tttt tFFPH H PHRyA xH 11|1|1|1|1()ttt tt tt tt tPF PPH H PHRH PFQ直至迭代到1, 2,tT。名师资料总结 - - -精品资料欢迎下载 - - -
25、- - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 8 页,共 11 页 - - - - - - - - - 时间序列分析方法讲义第 11章Kalman 滤波9 显然,迭代得到的序列1|?tt表示基于常数和变量11,ttyyy、11,ttx xx的最优线性预测,而1|ttP给出了这个预测的MSE。类似地,对1ty的预测和预测的MSE 为:1|1111|?(|,)ttttttttEyyxA x+H 11|11|1|?()() ttttttttE yyyyH PHR值得注意的是,在关于1|ttP的迭代计算过程中,甚至不需要对序列1|?tt进行计
26、算,它与数据过程没有关系,只与模型的系数矩阵有关。另外一个描述1|ttP的迭代过程的方式有时是有用的。我们从状态方程中减去卡尔曼更新方程 (11.34),得到:11|1|1?()()ttttt ttttt ttF KyA xH v再利用观测方程得到:11|11?()()tttttt ttttFK HK wv乘以自身的转置矩阵,再求数学期望得到:11|11|1|1?()() ()()() ()ttttttttt ttt ttttEEFK HFHKK RKQ利用1|ttP的定义,得到:1|1()()tttt ttttPFK H PFHKK RKQ(11.37) 上述方程与卡尔曼更新方程产生相同的序
27、列。11.3 基于状态空间表示的预测使用已知的数值矩阵F、Q、A、H和R,已经样本观测值的实际数据,上面介绍的卡尔曼滤波的计算可以通过计算机实现。为了使得这些想法更为清晰,我们利用一个简单的例子解析地说明这些计算的结果。例 11.3 利用卡尔曼滤波获得MA(1) 过程的确切有限样本预测考虑 MA(1) 过程的状态空间表示:状态方程 (2r):11100100ttttt观测方程 (1n):11ttty状态空间模型的变量符号为:1ttt,0010F,110ttv,2000Q,ttyy,A,1tx,1H,0tw,0R我们按照如下方式获得卡尔曼滤波:(1) 卡尔曼滤波的初始值为:名师资料总结 - -
28、-精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 9 页,共 11 页 - - - - - - - - - 时间序列分析方法讲义第 11章Kalman 滤波10 1|00?0,21|0200P(2) 对第一阶段的观测变量预测为:1|01|0? yH 对应的 MSE 为:222211|01|0210?()10(1)0E yyH P HR当然,这些数值恰好是MA(1) 过程的无条件均值和无条件方差。为了了解当2,3,tT时迭代的过程,我们考虑更新过程的基本形式。注意到F的第一行元素均为零,所以对所有的t,1|?
29、tt的第一行元素也都为零。注意到:1|1|?ttttt t十分自然地,1|?0tt。我们可以得到1ty的预测为:1|1|?1?ttttt tt ty对于这个例子而言,卡尔曼滤波更新方程的MSE 为:21|10001010000ttt tt ttpPFP FQP经过对比可以知道,1|ttP的第 2 行和第 2 列元素1tp与| t tP的第 1 行和第 1列元素相等,它是|?t t的 MSE。21|?()ttt tpE1|ttP的第 1 行和第 1 列元素可以解释为1|?tt的 MSE。 我们已经说明这个预测总是零值,而且对所有 t, 它的 MSE是2。 由于1|ttP是对角矩阵,因此预测误差1
30、1|?()ttt与|?()tt t不相关。根据1ty预测的 MSE 公式,得到:211|1|21221?()10100tttttttEyyppH PHR我们仍然可以从预测1|?ttt ty的本质性质发现这个公式的启示。2211|1|2221|?()()()?()()ttttt tttt tEyyEEE显然上述公式与前面的推导相符。根据预测迭代公式,|?t t可以根据下面公式迭代生成:21| 122|1|1000000101?10100tttt tttttypp名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - -
31、- - - - - 第 10 页,共 11 页 - - - - - - - - - 时间序列分析方法讲义第 11章Kalman 滤波11 或者:2|1| 122?t tttttyp迭代的初始值为0|0?0。需要注意的是,这里的|?t t与下面的近似是不同的:1?ttty,0?0收益矩阵为tK:222222000101/()100ttttpppK最后,我们得到:222|22100011000t tttttppppP名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 11 页,共 11 页 - - - - - - - - -