高斯向量似然函数的一个显式公式

使用公式∈RNexp(-[u\'Au*+ 2h′Bu*+ h′Ch]）du=（2π）N/2（det（A））-1/2exp(-[h′（C）]- 文学士-1B）h]），其中N=kq是u=的尺寸*anddet（A）=det（λ′λ+Iq）kdet(Ohm)-qwe推断：pdf（Z）=exp(-[v（Z）′[C]- 文学士-1B′[v（Z）]（（2π）kdet(Ohm))T/2det（λ′λ+Iq）k/2（29）指数在Z中是二次的- 文学士-1B′=Θ-1′T[IT]- λ[λ′λ+Iq]-1λ′]Θ-1T Ohm-引理1- λ[λ′λ+Iq]-1λ′为正定义。以下引理在向量化中是众所周知的：引理2对于任意两个大小相同的矩阵M和N，v（M）′v（N）=Tr（M′N）（30），特别是当H和Ohm 分别是大小为T×T和k×k的对称矩阵，如果X是大小为T×k的矩阵，我们有v（X）′（H Ohm)v（X）=v（HX）′v（XOhm) = Tr（X′HX）Ohm)这在vec和Tr之间提供了一种联系。集合‘K=’K（θ）=λ′λ+IQ那么‘K是一个q×q矩阵。K=K（θ，T）=IT- λ′K-1λ′=IT- λ[λ′λ+Iq]λ′K是一个T×T矩阵。然后是Pdf（Z | X，θ，Φ，Ohm) =经验(-Tr（Z′Θ-1′TK（θ，T）Θ-1TZOhm-1） ）（（2π）kdet(Ohm))T/2det（λ′λ+Iq）k/2（31），由此我们证明了（2）。现在让我们考虑Φ和中的部分优化问题Ohm 给定θ。替代品Z=X- u - LXΦ- ... - LpXΦpin方程（31），问题是找到Φiminimizing:Tr（X- u - LXΦ- ... - LpXΦp′Θ-1′K（θ，T）Θ-1T（X- u - LXΦ- ... - LpXΦp）Ohm-1） =Tr（Θ）-1X- Θ-1Tu- Θ-1LXΦ- ... - Θ-1LpXΦp）′K（θ，T）（Θ-1X- Θ-1u - Θ-1LXΦ- ... - Θ-1LpXΦp）Ohm-1） （32）由于表达式是以u和Φi为单位的二次表达式，因此可以通过线性回归和修改的内积对其进行优化。我们将矩阵Xθ，lagof大小T×（k×（1+p））表示为（Θ-1T1 |Θ-1TX |Θ-1TLX |····|-1TLpX）。如果我们不包括常数项u，我们可以排除块Θ-1T1。S etXθ=Θ-1TX（33）那么uΦΦ...Φpopt=（X′θ，lagKXθ，lag）-1X′θ，lagKXθ是最佳选择。引理A.1证明了这一点。

我们注意到，这并不取决于Ohm.通过选择Φ，isTr（（Xθ）′KXθ上方的二次型（32）的最小值Ohm-1.-Tr（（Xθ）′KXθ，滞后（X′θ，滞后KXθ，滞后）-1X′θ，lagKXθOhm-1） （34）与VAR模型的最大似然参数f类似，使用与det和Tr有关的Jacobian公式的论证表明Ohm 使对数可能性最小化的方法是：Ohmopt（θ）=T（Xθ）′KXθ-（Xθ）′KXθ，滞后（X′θ，滞后KXθ，滞后）-1X′θ，最后是lagKXθ，值为Ohm, 轨迹表达式中的矩阵是简单的T.ika，因此轨迹是T k。条件对数概率是：\'\'L（θ）=-T klog（2π）-Tlog（det）(Ohmopt（θ）））-klog（det（λ′λ+Iq））-T kWe注意，这里出现的公式看起来非常像正则回归/协方差公式，但内积由Θ给出-1′TK（θ，T）Θ-1T。让我们讨论连接∑TandΘ的关系-1t和K（θ，T）。这纯粹是一个只涉及θ的代数等式。我们观察到，在纯移动平均的情况下，带有标量θ的似然函数被简化为标量MA（q）模型，并使用Ohm. 比较g（2）在k=1，p=0的情况下，Ohm = σ、 u=0L（θ，X··XT）=-Tlog（2π）-Tlog（det）(Ohm))-对数（det（λ′λ+Iq））-（X′Θ）-1′TK（θ，T）Θ-1TXOhm-1） 对于已知的MA（q）公式，例如（5.5.5）in（Hamilton1994）（注∑Tin），我们的符号是σ-2.Ohm 在那个参考中：L（θ）=-Tlog（2π）-对数（det（∑T））-（X′∑）-1TX）我们得到了所需的方程式。3关于MA（q）的Szeg–o极限对于MA（1）众所周知，det（σ）的大T极限-1T）=det（KT）仅为1- θ. 由于这个决定因素出现在可能性函数中，所以自然要问，一般来说，类似的结果是否成立。结果表明，行列式的大T极限总是非多项式的。强Szeg–o极限定理（Szeg–o 1952；Bingham 2012；Basor和H。

2000）演示了如何计算单位圆盘上某些分析函数产生的截断Toeplitz矩阵行列式的大T极限。对于由有理函数生成的Toeplitz矩阵（基本上是一般的VARMA情况），极限是Knownder状态空间表示，例如（Gohberg、Kaashoek和Schagen 1987）。（Kramer和Rosenblatt 1993）也提到了这个定理。对于MA（q）的情况，表达式非常简单，但我们无法找到一个参考，所以让我们陈述：定理2，如果θ（L）=Qqi=0（1）- λiL）是可逆的→∞det（∑T）-1=（qXi=0θi）（qXi=0(-1） iθi）Y1≤i<j≤问题（1）- λjλj）（35）最后一项是λi中的一个对称多项式，因此它也可以在θi’中表示为一个多项式，适用于函数a（L）=θ（L）θ（L）的Szeg¨o极限定理-1） wehavelimT→∞det（∑T）=exp(∞Xk=1k（log（a））k（log（a））k表示我们取log（a）的洛朗展开式的第k个系数。但是我们有log（a（L））=qXi=1log（1）- λiL）+qXi=1log（1）- λiL-1） 所以我们很容易看到：log（θ（L））k=qXi=1λkikSo指数项是∞Xk=1k（qXi=1λki）=qXi=1Xkλ2kik+2X0≤i<j≤qXk（λiλj）kk与∑TisY（1）的极限-λi）-1Y（1）-λiλj）-2=Y（1）-λi）-1Y（1+λi）-1Y（1）-λiλj）-这就是我们必须证明的。我们可以很容易地计算q（1）的表达式- λiλj）f或小q。下表总结了→∞Σ-1到q=3。q=1- θq=2（1）- θ)(1 - θ） q=3（1）- θ)(1 - θ+ θθ- θ） 对于更高的q，用θ表示的多项式在多项式项上扩展为大量的m，因此用λi表示更简单。我们将在下面看到，它们是对小q强制执行可逆性条件的相同多项式。请注意，“kT”是一个q×q矩阵，因此我们也可以尝试直接计算其极限d。

虽然我们可以在大T极限下证明，\'kt的项是θi中的有理函数，但这些项的显式表达式很复杂，因此行列式是多项式的倒数这一事实很有趣。对于协方差矩阵也可以进行类似的计算。在可逆稳定ARMA（p，q）情况下，结果包括φ（L）/θ（L）的分子和分母的根（我们假设φ（0）=θ（0）=1），用u表示它们-1i和λ-假设|μi |<1，|λj |<1i、 j.极限公式为：limT→∞det∑T=Y（1）- ui）-1Y（1）- uiul）-2Y（1）- λj）-1Y（1）- λjλm）-2Y（1）- uiλj）（36）我们注意到术语q（1）- uiλj）是φ（L）和Lqθ（L）的乘积-1） 达到一定的比例系数。4两个函数L（u，θ，Φ，Ohm) 和L（θ），我们在校准中主要处理第二个。第一种方法的海森误差与回归系数的标准误差有关。我们使用符号国际会计准则简称θi.设T（f，T）为上述Toeplitz映射，将幂级数f映射到其截断的下三角Toeplitz矩阵。

我们有ΘT=T（θ，T）我们有Θ-1T=T（1/θ，T）iΘ-1T=-T（Li/θ，T）Θ*,T-Qθi=Q-i+1，i-1Iq-i+1T-q+i-1.我-1T-q+i-1，q-i+1这里0表示大小为a×b的零矩阵块，因此最后一个等式的右手边有一个大小为q的单位矩阵- 右上角为i+1，其他地方为零。放Θ-尺寸为T×（q）的1英寸块矩阵形式- i+1）和T×（T-q+i- 1):Θ-1T=[（Θ）-1T）T，q-i+1（Θ-1T）T，T-q+i-1] 那么Θ-1TΘ*,T-Qθi=（0T，i）-1, Θ-1T，q-i+1）iλ=我[ΘΘ-1TΘ*,T-q] =(iΘ-1)Θ*,T-q+（0T，i）-1, Θ-1T，q-i+1）i\'K=(iλ）′λ+λ′iλ我知道-1= -“K”-1(好的-1.iXθ=-T（Liθ）-2，T）XiXθ，滞后=-T（Liθ）-1，T）Xθ，对于任意两个矩阵A，B取决于θiAKB=(iA）KB+AK(iB）- A(iλ）`Kλ′B- λ′K(iλ）′B-λ(i′K）λ′B（37）如果A=B，计算会进一步简化，因为结果是对称的，所以我们只需要计算一半的项。应用（37）A，B作为Xθ或Xθ，我们可以计算d（X，θ）=X′θKXθClag（X，θ）=X′θ，lagKXθ，lagBlag（X，θ）=X′θ，lagKXθ的偏导数Ohmopt可以写成D（X，θ）- 布拉格-1拉格布拉格，A，B，C，D由Xθ和Xθ计算得出，因此我们可以借助i（B′C-1B）=(iB′）C-1B+BC-1(iB）- B\'C-1(iC）C-1B。把所有的东西放在一起，我们就得到了偏导数我Ohm选择此外伊洛格（德特）(Ohmopt=Tr(Ohm-1.我Ohm（可选）ilog（det（λ′λ+Iq））=Tr（\'K-1.（好的）i′L（θ）=-T伊洛格（德特）(Ohm（可选）-Kilog（det（λ′λ+Iq）），所以我们有梯度。应用链式规则和各种矩阵导数规则，我们也可以计算‘L’的Hessian。在这里，我们需要我jΘ-1T=2T（Li+j/θ，T）和矩阵乘积规则的导数。这个计算很乏味，但可行。然而，我们不会在这里进行计算。根据最大似然估计的Fisher矩阵的一般理论，L的Hessian与θ和Φ的标准误差估计有关。

我们注意到Klein和Melard（K lein and M’elard 2014）在一般（矩阵θ）VARMAX情况下的工作。从L的表达式来看，Hessian块Hθθ（L）是更复杂的，可以通过FFT卷积来完成，包括卷积1/θ，1/θ，1/θ与X。我们可以用数值计算它。然而，我们注意到，HθΦ（L）和HΦΦ（L）块相当简单，第一个可以从（37）中计算出来，然后用X替换两个Z项中的一个。第二个是VaR情况的直接推广：HΦΦ（L）=Tr（X）Ohm-1′TKOhm-1TXOhm-1） 另一种计算梯度的方法是将L和¨L视为∑T的函数，因此视为函数γi，然后表示θ的γias函数，如（安·德森和塔克穆尔a 1986）所述。他们的计算表明了一个有趣的事实，雅可比矩阵γiθjlooks与行列式K:det的zeg–o极限密切相关(γiθj）=θq+1Y（1）- λ-1i）Y（1+λ）-1）Y（1）- λ-1iλ-1j）最后，我们可以将这种方法推广到θ系数是有限个参数pj的函数的情况。如果我们只处理梯度，我们只需要矩阵(θipj）然后应用链式法则。因此，我们可以应用我们的核心代码来校准更通用的模型。我们将在后面关于扩展模型的章节中讨论这一点。5计算和校准我们已经提到了Ohm-lTS只是1/θL与S的卷积。唯一遇到的长矩阵计算形式是T（1/f，T），其中f是θ，θ，θ之一，a是XorΘ之一*;T-q、 卷积计算可以通过快速傅立叶变换完成。首先是1/f的计算。如果f是一个低阶多项式，我们很可能会遇到，我们可以使用递归算法来计算1/f，或者将f扩展到formc/（1）的部分分数- dL）然后将幂级数展开应用于后者。

另一种方法是使用1/f的快速收敛扩展，例如（Harvey 2011）。为了使用快速傅里叶变换计算卷积，我们假设1/f的系数足够小，可以在c（θ）步后忽略。FFT卷积算法将T分割成可以高效计算FFT的短段，并利用FFT后的偏移卷积转换为互补乘法。例如，使用重叠保存方法将这些部分拼接在一起。实际上，因为我们还需要计算Θ-1TLiX对于i=0，·pit更便于计算-1T+p^X通过FFT卷积。阿斯伯格=十、Xp。。。XT+p写C=Ohm-1T+pand^X，尺寸为p的块体- i、 T和我：Θ-1T+p^X=C[1:（p-i） ，1:（p-i） [0摄氏度[（p-i+1:（T+p）-i） ，1:（p-i） ]Θ-1TC[（T+p-i+1:，1:（p-i） [C]*C**^X[1:（p-i） [LiX^X[（T+p]-i+1）：]我们看到p行的块- i+1至T+p- 我当然-1T+p^X isC[（p-i+1:（T+p）-i） ，1:（p-i） [1:（p-i） ]+Θ-从这里出发-1TLiX是通过将这些行减去firstterm得到的。ΘT+p的子矩阵可以用θ（L）的幂级数展开的分段表示-1.该调整计算的总成本为O（tp（p+1）/2，每个调整计算的成本为O（T×i）。对于大T，调节块的贡献衰减相对较快。我们注意到，如果我们使用FFT，θ的反演是O（T log（T））和卷积isO（kT log（T））。如果Tif T远大于k，则线性回归为O（（kp+1）T）。总体而言，对于常数c，c，似然函数的计算顺序为O（c（kp+1）T+ckT log（T））。这已经比标准卡尔曼滤波计算的O（T）估计有所改进。然而，时间O（T+log（T））的近似值见（Pnevmatikakis等人，2014）。该方法的优点在于校准。我们只需要优化θ参数中的函数L（θ）。

与传统的卡尔曼滤波方法相比，该函数的计算和优化更简单。首先，函数是对称的，我们可以通过应用转置节省一半的计算。这里遇到的矩阵是三角Toeplitz矩阵。其次，K是唯一遇到的大方阵。但我们永远不需要直接计算K。回忆‘K=Iq+λ′λ‘K的大小为q，且对称。因此，我们可以进行Cholesky分解，K=CKC′KHere cke是一个大小为q×q的下三角矩阵。因此，所有后续的计算都基于三角矩阵。因为它- λ′K-1λ′要计算N′KM，N和mh为少量列，T行，我们实际上需要计算N′M，N′λC-1′KandC-1Kλ′M。这里所有的矩阵乘法和求逆都涉及大小为M×T的矩阵<< T下一个问题是如何确定p和q。这里我们回到前面关于N（L）表达式的讨论-1D（L）到标量分母形式，θ（L）=det（N（L））和Φ（L）=NA（L）D（L）。如果N和D来自一个最小的实现，我们看到deg（det（N（L））和deg（NA（L））D（L）都小于McMillan度m。因此我们可以选择在p和q都等于McMillan度的情况下最大化可能性，然后在信息标准的帮助下应用测试来确定哪些更高的阶数可以被消除。这将需要进一步的研究，因为如果q的实际程度远小于m，这一建议可能远远不是最优的。我们希望cannonicalcorrelation分析在这里发挥作用。现在让我们讨论一下似然函数的实际最大化。我们再次回顾，似然公式对复平面上的任意θw根都是有效的。然而，如果我们试图用θ来计算它，在单位d isc内至少有一个r oot，那么两个都是ΘTandOhm将假定一个非常大的值，即使可能性函数仍然是有限的。

第（6）节中讨论的根反转映射在每种情况下都给我们提供了一对（θIR，OhmIR）具有相同的似然函数值，其中θIR的根位于单位圆盘之外。单位卷曲上的根的情况是特殊的。我们将在第6节中详细讨论。我们总是可以选择初始模型进行投资，如果我们想要识别系统，我们必须这样做。θ（L）在单位d isc之外有根的区域是凸连通区域，如第6节所讨论。似然函数不是一般凸函数。事实上，对（Davisand Dunsmuir 1996）中标量MA（1）情况的仔细分析表明，它可能有几个具有不同渐近性的局部极大值。在下一节中，我们将讨论该区域的一些特性，以及如何选择校准的初始点。对于实际校准，读者可以选择他最喜欢的基于梯度的优化器，L-BFGS-B是作者在这种情况下选择的方法。当成本函数在外部时，我们通过给其赋值，迫使优化保持在可逆区域内。考虑到区域的形状以及根反转图下梯度变换的方式，可能会有更好的算法。6根反转图Shansen和Sargent（Hansen和Sargent 1980）（也有（Hamilton 1994）详细阐述）提出了一种方案，将任何Masystem转换为单位圆盘内有根的Masystem。在这种模式下，新模型的自方差母函数与原模型相同。我们在这里证明了这些变换在标量θ为VARMA的情况下有效，并检查它们是否也保持了条件对数似然。

这是标量移动平均情况下相似结果的直接推广。假设方程θ（L）=0有根λ-1, · · · , λ-1q：θ（L）=qYl=1（1）- λlL）回忆Zt=Pqi=0θit-i、 所以Z是一个VMA过程。让γibe在（12）中给出。考虑Z:g（Z；θ；Ohm) = (∞十一=-∞（izi）Ohm =qYl=1（1- λlz）（1- λlz-1)Ohm让我们回忆一下根反转映射是如何构造的。选择与λi···λira对应的指数i···，irc的子集，并考虑多项式θIR（L）=θIRi，··，IR（L）=Yi6∈{i···ir}（1）- λiL）Yi∈{i···ir}（1）- λ-1iL）（38）和协方差矩阵OhmIR=（λi··λIR）Ohm （39）以下定理总结了该映射的一些重要性质：定理3在与i··ir对应的根反转映射下，θIRis不变量的自方差生成函数：g（z；θir；OhmIR）=g（z；θ；Ohm) （40）让ΘIR；Tand′K（θIR），K（θIR，T）是Toeplitz矩阵，K矩阵对应于θIR，我们有：∑IR；T:=ΘIRK-1IR′IR=（λi··λIR）-2∑T（41）det（`K（θIR））=（λi··λIR）-2Tdet（`K（θ））（42）uΦopt（θIR）=uΦopt（θ）（43）Ohmopt（θIR）=（λi···λIR）Ohmopt（θ）（44）条件似然函数在这种变换下也是不变的。L（θIR，u，Φ，OhmIR）=L（θ，u，Φ，Ohm) （45）`L（θIR）=`L（θ）（46）设J[IR]为IR的梯度。L变换的梯度如下：L（θ）= （40）、（41）的证明与（Hamilton 1994）中的标量情况相同。（43）是直接从它的表达。（42）从det（`K）=det（∑T）开始的跟随。最后两个方程来自方向替换。注λ，··，λqare根的Lqθ（L-1） ，我们有时可以称之为根，除非有困惑。可逆性条件为|λi |≤ 1.我们有从根（λ，·λq）到系数（θ，·θq）的Vieta映射v。