1、时间序列分析方法讲义 第 5 章 最大似然估计1第五章 最大似然估计在本章中我们开始讨论时间序列模型的参数估计方法，其中极大似然估计是一种最为常用的参数估计方法。我们仅仅讨论极大似然估计的原理和似然函数的推导，而对获取极大似然估计的算法不加以详述。5.1 引 言5.1.1 ARMA 模型的极大似然估计假设数据的真实生成过程是一个 过程，则该过程的数据生成机制为：),(qpARMqttttpttt YYc 121其中 是白噪声序列，满足：tsEts,0)(我们将要讨论如何利用 的观测值来估计母体参数：tY),( 22121qpc我们将要采用的方法是极大似然估计方法，因此需要获得似然函数的表达式。
2、假设获得了 个样本 ，如果能够计算出相应的联合概率密度函数：T),Ty ;(21),(1Yf上述函数可以视为在给定参数下样本发生的概率，因此合理的参数取值是使得上述概率最大，如此参数便称为极大似然估计。这时我们需要极大化上述联合概率密度。为此，我们假设噪声序列是高斯白噪声序列，即 ),0(.2Ndit虽然这个假设非常强，但是在这样假设下得到的参数估计 ，对于非 Gauss 过程来说也是很有意义的。具体求解极大似然估计的步骤是：一是先求出并计算似然函数，二是求似然函数的最大值。这里涉及到一些代表性的非线性数值优化问题。5.2 高斯 过程的似然函数)1(AR假设数据生成过程是一个具有高斯白噪声序列
3、的 过程：)1(ARtttYc1这时对应的参数向量为： 。我们首先寻求联合概率分布函数，也就是这),(2c些参数对应的似然函数。(1) 求上述过程似然函数的代表性过程是利用条件概率密度进行传递，所以需要先求出 的概率密度。它的均值和方差为：1Y，cE221)(Y由于它具有正态分析，因此对应的密度函数为： )1/(2exp)/(1),;();( 22211 cycyfyfYY(2) 在给定 的条件下， 的条件概率分布可以得到：Y时间序列分析方法讲义 第 5 章 最大似然估计2),(| 2112ycNyY对应的概率密度函数为： 21212| )(exp);|(1 ycfY(3) 类似地，在给定前两
4、个观测值的条件， 的条件概率密度函数为：3Y22123,| )();,|(213 yfY注意到上述条件概率分布中只依赖一阶滞后的条件观测值。(4) 最后一个样本的条件概率分布为： 212121,| )(exp);,|(1213 TTTY ycyyfT注意到上述条件概率分布中也只依赖一阶滞后的条件观测值。(5) 根据无条件密度函数与条件密度函数之间的关系，可以得到： Tt tYYTY yfyfyyf tT 21|1121, );|();();,( 11 经常对上述函数取对数，得到对数似然函数： );|(log);(log)( 1|211tYtYfftL(6) 将具体的密度函数代入上式，可以得到
5、过程的似然函数为：(ARTt ttycTTy2212212 )(log/)()log2/)1( )/ll) 可以将上述似然函数表示为更为紧凑的向量和矩阵形式。令均值向量和自协方差为和 ，注意到过程之间具有的自协方差函数表达形式，则有：，V2 11132132122 TTT这样一来，所观测到的样本可以当作多元正态母体 的一个简单抽样，具有的)(,N联合概率密度函数为： )(y)();(y 1Y2exp|22/1/Tf理论上可以对上述极大似然函数求导数，然后获得参数估计。但是，一般情况下的导数方程是非线性方程，难以获得精确的最大值估计。一种近似的方法是假设第一个观测值是确定性的，然后求解给定 时的
6、条件似然函数值，这时的目标函数是：1YTt ttTY ycTyfT 22122|, )( log/)(log/);|,(log12 上式最大值相当于求下式的最小值：时间序列分析方法讲义 第 5 章 最大似然估计3Tt ttyc221)(上式的最小值就是线性回归的最小二乘估计，满足方程： TtttTttt yyc212121类似地，噪声的方差为： Tt ttc21)(当样本容量足够大时，可以证明上述近似或者条件极大似然估计具有与精确极大似然估计一致的极限分布。5.3 高斯 过程的似然函数)(pAR对于一般的高阶自回归过程：，tptttt YYc21 ),0(.2Ndit此时所要估计的总体参数向量
7、是： 。,(21pc(1) 似然函数的估值 Evaluating the Likelihood Function假设我们获得了 个来自 过程的样本，假设前 个样本表示为T)(pAR),(21ppyy可以将这个向量当作 维 Gauss 变量的一个样本。这个向量的均值表示为 ，它的p每个分量都是： )/(1pc假设 是 的协方差矩阵，则有：pV2,Y 221 222 11212 )()()( )()()( ppp pp YEYEYE 对于一阶自回归过程而言( )，上述矩阵是一个标量， ；对于 阶1/Vp自回归而言： 032130211 1202 pppV这里 是 过程的第 个自协方差，可以按照以前
8、的介绍公式计算。j)(ARj由于自回归过程的条件相依性具有截断性质，因此我们将样本分为 p 个一组，样本中前 个观测值的联合概率分布为 ，密度为：),(2pNV时间序列分析方法讲义 第 5 章 最大似然估计4 )(yV)(yV);( ppppppYfp 122/12/2/1, ex|)()(|,1 对于样本中剩余的观测值 ，我们可以使用推断误差分解(prediction ),Tyerror decomposition)，将前 个观测值作为条件，则第 个观测值的条件分布为 Gauss 分t t布，且均值和方差分别为：，pttt yyc21 2只有 个最近的观测值与这个分布有关，因此，对于 ，则有
9、：p pt 22122 1,|11,| ex ,|,| 1121 )( );();( ptttt tttYYttYY yycyff tpttt 因此，整个样本的似然函数为： Tpt ptttYYppYTTY yyyyfyf ptptpT 1 21,|121, ,|, 1211 );();( 则对数似然函数形式为： Tpt ptttt pppTpt ptttt pppttYY yycy yycyTpfLpt1221 112 22112121,| )()(|log)l()2log(ll )()(|og)()2og(,|l)(11 )VV); 为了获得上述似然函数值，我们需要获得逆矩阵 ，为此我们有
10、下述命题：1命题 5.1 利用 表示矩阵 的第 位置的元素，则对任意 ，有：)(vij 1p),(ji pji1jipjkijkikijijpv110)( 这里 。证明：略。 End因为 是对角矩阵，因此也可以得到 时的元素 。1pVji)(pvij例如，对 过程而言， 是一个标量，取 ，得到：)(AR1pV1)()(2210101 kkv时间序列分析方法讲义 第 5 章 最大似然估计5因此有： 212V因此命题 5.1 确实可以重新得到 过程的方差表达式。)1(AR对于 的情形，利用命题 5.1 可以得到：p221112)(可以计算行列式值为： )()()(| 21221212 V并且有：
11、)(1)()2)(1)( 1,)()( 222122222 yyyy因此，对于 Gauss 条件下的 过程，确切的似然函数为：AR Tt ttt ycy yyL3221 22212 12 )()()(1 log(l)log()( (这里： )/(2) 条件极大似然估计 Conditional Maximum Likelihood Estimates由于目标函数形式比较复杂，因此对 过程的确切极大化必须使用数值算法。与)(pAR此对应，以前 个样本为条件的对数似然函数具有下述简单形式：p Tpt ptttttTYYY yycyTyfpt 12212211,|, )log()2log( ,|,lo
12、g111 )();注意到，极大化上式的参数 与极小化下式的参数选择是一致的：,1c Tpt ptttt yycy1 221)(因此这些参数 的条件极大似然估计是 基于常数和自身滞后值的普通),(p ty最小二乘回归估计， 的条件极大似然估计是这个回归方程平方残差的平均值：2 Tpt ptttt ycy1 2212 )( 显然上述条件似然函数与确切似然函数相比，缺少了初始样本的母体分布，这样就降低了样本发生的似然性，这就是条件似然函数与确切似然函数的差异。类似地，确切的极大似然估计和条件极大似然估计能够得到相同的大样本分布。时间序列分析方法讲义 第 5 章 最大似然估计6(3) 非 Gauss
13、分布时间序列的极大似然估计 Maximum Likelihood Estimation for Non-Gaussian Time Series根据线性回归模型的性质，如果假设随机过程关于二阶矩是遍历的，则我们知道普通最小二乘估计也是下面线性投影系数的一致估计(consistent estimate)： ),|(21ptttYYE同时这个 OLS 估计也使得 Gauss 条件似然函数达到最大。因此，即使一个过程不是Gauss 过程，但是我们错误地将它当作 Gauss 过程，并且极大化它的似然函数，则得到参数估计 仍然是母体参数的一致估计。),(21pc一个从错误指定的似然函数中获得的极大似然估
14、计(例如，当真实过程不是 Gauss 过程，而从 Gauss 过程假设下获得的 MLE)被称为“伪极大似然估计”(quasi-maximum likelihood estimate)。但是，有些时候，例如上述提到的情形， “伪极大似然估计”能够为人们感兴趣的参数提供一致估计。但是，如果真实数据不是服从 Gauss 过程的，在 Gauss 假设下获得参数估计的标准差则是不正确的，它不具有估计的一致性。另外，如果原始数据不是服从 Gauss 过程的，有时通过简单的变化，例如取对数等，可以获得 Gauss 分布过程。对于一个正的随机变量 ，Box and Cox (1964)提出了一类比较tY一般的
15、函数变换： 0,log1)(tttY一种方法是，选取特殊的 ，在 是高斯 ARMA 过程的假设下极大化 的似然函)(tY)(tY数，这时能够使得似然函数取最大值的 被认为是最合适的变换选择。但是，也有一些研究者对这种方法在实用中的效果提出了质疑。5.4 高斯 过程的似然函数)1(MA利用类似自回归过程似然函数的推导方法，我们开始类似地讨论移动平均过程的极大似然估计问题。5.4.1 条件似然函数 Conditional Likelihood Function如果基于 Y 的初始值，则自回归过程的条件似然函数计算是相当简单的。类似地，如果利用初始的 作为条件，移动平均过程的似然函数计算也是比较简单
16、的。考虑一个高斯 过程：)1(MA，ttt),0(.2Ndit利用参数向量 表示母体参数，如果 的数值是确定可知的，则：,21t)(| 11ttNY表示成为密度函数情形为： 21221| 1exp,|1 )()( tttY yyft 假设我们确定性地已知 ，并且获得观测值 ，则下一个阶段的误差也是确定性0已知的，可以确定为： 1y时间序列分析方法讲义 第 5 章 最大似然估计7利用上述条件概率分布公式得到： 21220120,| 1exp1,|12 )()( yyfY由于 是确定性可知的， 可以按照下式计算：12因此，已知初值 ，可以利用叠代方式从样本 中获得误差序列0 ,21Ty：,1T，1
17、tttyT,2则条件概率分布为： 22 10|010,| exp ;|;,| 1121 t tYttYY yfyyf ttt )()(则联合概率分布可以通过单独的条件概率密度的乘积得到： Tt ttYYY yyfyf ytttt 2 01210,|01| 121|, ;,|;| ;|,121012 )()( ) 则条件对数似然函数为： TttTyfLtt 122010|, )log()l(2 ;|,og)(121 )虽然上述函数比较简单，但由于是参数的非线性函数，无法直接利用解析方法获得参数的极大似然估计。因此，即使是 1 阶移动平均过程，条件极大似然估计也需要利用数值优化方法求解。从任意初始
18、误差 开始叠代，可以得到：00121 )()()()()( tttttt yyyy 如果 显著地比 1 小，则附加初始条件 带来的影响将快速地消失，则当样本|0数量比较大的时候，条件似然函数是无条件似然函数比较好的近似。相反地，如果 ，1|则附加 带来的结果将随着时间而累积，在这种情形下条件概率方式不是很好的近似0方法。这时如果数值方法模拟的结果是 ，则需要放弃这样的估计结果。这是需要利1|用 的倒数作为初值重新进行条件似然估计的数值计算。5.4.2 确切似然函数 Exact Likelihood Function对于高斯 过程，有两种比较合适的算法来计算确切似然函数。一种是使用后面)1(MA
19、介绍的卡尔曼滤波方法，另一种是方差或协方差矩阵的三角因子分解。下面我们开始介绍矩阵的三角因子分解方法。假设样本观测值可以利用 矩阵表示： ，这个向量具有均值为1T),(21Tyy和 的方差协方差矩阵 ：),( )(YE从高斯 过程中获取样本的方差协方差矩阵为：1MA )1(000)1()(02222 时间序列分析方法讲义 第 5 章 最大似然估计8则联合分布表示的似然函数是： )(y)();(y 1Y2exp|22/1/Tf似然性的推断误差分解可以从下面矩阵 的三角因子化中获得：AD这里的矩阵 是下面形式的下三角矩阵： 11)(00 001)(0(242422 n 矩阵 是下述形式的对角矩阵D
20、 )1(242426242 100 0110 n 将一个对角矩阵表示成为下三角矩阵的对角化，称为一个矩阵的三角因子化。将这样的矩阵分解带入到似然函数中，得到： )(yAD)(yAD);(yY 12/12/ exp|Tf注意到矩阵 是主对角线都是 1 的下三角矩阵，因此它的行列式为 ，这时有：|AD进一步定义： )(y1则似然函数可以表示成为： yD);(Y 21exp|22/1/Tf注意到变换形式也可以表示为： yA这个系统的第一行表明 ，而第 t 行意味着：1y12421 tttt y）（ ）（ 因此向量 可以从 开始，利用上式进行叠代得到， 。变量y1 Tt,32可以解释为 基于常数和 的
21、线性投影，而矩阵 的第 个对角元素是这tyt 12,ytt D个线性投影的 MSE： ）（ 124221)( tttYEd 时间序列分析方法讲义 第 5 章 最大似然估计9由于矩阵 是对角矩阵，其行列式是对角元素的乘积：DTttd1|而矩阵 可以利用矩阵 的主对角线元素取倒数得到。因此：Ttty12这样一来，似然函数可以简化为： TttTttdydf 1212/ exp)(;yY因此，高斯 过程的确切对数似然函数为：MATt TtttdyfL112)log(2)log()(log)( ;如果给定参数的初始数值，可以利用叠代方法计算上述似然函数。这样的计算可以对尝试参数的优劣进行比较。与条件似然
22、函数不同的是，上述似然函数的表达式对任何参数 都是有效的，而不取决于是否该过程具有可逆 表示。)1(MA5.5 高斯 过程的似然函数(q在讨论了一阶移动平均过程似然函数的基础上，我们继续讨论 过程的极大似然)(qMA估计问题。5.5.1 高斯 过程的条件似然函数 Conditional Likelihood Function)(qA考虑一个高斯 过程：M，qtttttY21 ),0(.2Ndit此时一个简单方法是考虑给定前 个 初值的条件似然函数，假设初值为：0210q利用这些初值，我们对时刻 ，利用下式进行叠代：Tt, qtttty 21利用向量 表示 向量 ，则条件对数似然函数为：0),(
23、110 TttYYTyyfLtt 1220|, )log()2l( ;|,og)(11 )这里的参数向量为： 。与 1 阶情形类似，这里要求下述特征方程),q的根落在单位圆外面： 012qzz上述条件是要求移动平均过程是可逆的，这时条件似然函数的参数估计对误差初值的选取没有严重的依赖性。5.5.2 高斯 过程的确切似然函数 Exact Likelihood Function)(qMA由于高斯过程的联合分布是多元正态概率分布，因此高斯 过程的确切似然函数)(qMA时间序列分析方法讲义 第 5 章 最大似然估计10可以表示为： )(y)();(y 1Y2exp|22/1/Tf这里的样本为 ，均值为
24、 ， 的方差协方差矩),(Ty , T阵为： qq0上述矩阵中的元素为： qjijiji|,0),(|这里 是 过程的 k 阶自协方差函数：k)(MA， qkqk, ,21,0),212 10与一阶情形是一样的，确切似然函数可以通过卡尔曼滤波或者矩阵 的三角因子化获得。仍然考虑三角因子化，这时矩阵 可以分解为：AD这里的矩阵 是下三角矩阵，矩阵 是对角矩阵，它们具有前面介绍预测时的表达式D类似。进一步分析可以发现，由于矩阵 具有的特殊谱结构，矩阵 具有性质：对于A， 。对数给定的数值参数，计算机程序可以十分容易地获得这些矩阵的数jqi0ija值。将这样的矩阵分解带入到 过程的似然函数中，得到：
25、)(qMAyD);(yY 21exp|22/1/Tf这里： ， 的元素可以叠代计算如下：Ay，1，122)(ay，33 qtttttttt yay)( ,2,1, 因此，高斯 过程的确切对数似然函数为：(qMATt TtttdTfL112)log()log()log)( y;Y给定参数数值，可以计算上述似然函数的函数值。但是，给定样本，计算上述函数的极大值的非线性优化过程却是比较复杂的。时间序列分析方法讲义 第 5 章 最大似然估计115.6 高斯 过程的似然函数),(qpARM将自回归过程和移动平均过程的方法结合起来，就可以讨论移动平均过程的条件似然函数和确切似然函数问题。5.6.1 条件似
26、然函数 Conditional Likelihood Function考虑一个高斯 过程：),(qpAR，qttttptttt YYc 2121,0(.Ndi我们的目标是估计母体参数向量 。),( 221 qpc此时，对自回归部分条件似然函数的近似依赖于初始的 ，对移动平均部分条件似然y函数的近似依赖于初始的 。对移动平均过程 的似然函数共同的近似即依赖),(ARM初始的 ，也依赖初始的 。y将初始的 和 当作给定的，则误差序列),(110py ,110q可以通过 ，利用下式叠代获得：,21T 2T，qtttpttttt yYc 2,因此条件似然函数是： TttYYTyyfLtt 122 01
27、,|, )log()2l( ;,|,og)(011 )Y一种选择是将 和 的初值设定为它们的预期值，即：y，pscy21 1,0ps，01,q然后对 ，进行上述的叠代计算。Tt,另外，Box and Jenkins (1976) 建议设定 为零，但将 设定为它们的真实数值。因此，y对上式的叠代可以从时点 开始，将 选取为真实的样本观测值，选取pt ,21py初值： 011qpp则对数条件似然函数可以计算为： Tptt qpppTTyyyfL 122 1112)log()log(2 0,0,;,|,l)( ) 如同移动平均过程的情形，这种近似只有当下述特征方程： 012qzz的特征根全部落在单位圆外面才能应用。
