缺失数据重尾AR模型的参数估计

（18） 然后，对完整数据日志可能性的期望可以表示为qθ|θ（k）=ZZlog（L（θ；y，τ））pym，τ| yo；θ（k）dymdτ=ZZlogh（y，τ）exp-ψ（θ）+Ds（yo，ym，τ），φ（θ）E×pym，τ| yo；θ（k）dymdτ=ZZlog（h（y，τ））pym，τ| yo；θ（k）dymdτ- ψ（θ）+DZZs（yo，ym，τ）pym，τ| yo；θ（k）dymdτ，φ（θ）E=-ψ（θ）+D′sθ（k）, φ（θ）E+常数。，（19） 在哪里θ（k）=ZZs（yo，ym，τ）pym，τ| yo；θ（k）dymdτ。（20） 通过利用指数族的特性，可以方便地简化EM算法。EMalgorithm的步骤简化为计算经验最小有效统计数据θ（k）, 并且M阶跃减小到函数（19）的最大值。A、 步骤给定yo和θ的ymandτ的条件分布为：p（ym，τ| yo；θ）=p（y，τ；θ）p（yo；θ）=p（y，τ；θ）RRp（y，τ；θ）dymdτ∝ p（y，τ；θ）=TYt=2νντν-1tΓν√2πσexp-τt2σ（yt- φ- ^1yt-1)-ντt∝TYt=2τν-1试验-τt2σ（yt- φ- ^1yt-1)-ντt. （21）由于积分rrp（y，τ；θ）dymdτ没有表达式，我们只知道p（ym，τ| yo；θ）u p to a标量。此外，比例表达式是复杂的，我们无法得到条件期望的闭式表达式θ（k）或Qθ|θ（k）. 因此，我们研究了SAEM-MCMC算法，该算法根据马尔可夫链过程中的条件分布生成样本，并逼近期望值θ（k）和Qθ|θ（k）通过随机逼近。我们建议使用吉布斯抽样方法来生成马尔可夫链。吉布斯采样器将最新的t变量（ym，τ）分为两个块τ和ym，然后通过从分布p（ym，τ| yo；θ）的条件分布p（τ| ym，yo；θ）和p（ym |τ，yo；θ）交替绘制实现，从分布p（ym，τ| yo；θ）生成样本s的阿马尔科夫链。

更具体地说，在迭代k时，给定当前估计θ（k），吉布斯抽样器开始于τ（k-1，l），y（k-1，l）m（l=1，2…，l）并生成下一个样本τ（k，l），y（k，l）m通过以下方案：o来自p的样本τ（k，l）τ| y（k-1，l）m，yo；θ（k）,o 样品y（k，l）mfrom pym |τ（k，l），yo；θ（k）.然后是预期的最低充分统计数据θ（k）和预期的完全数据可能性Qθ|θ（k）近似为^s（k）=^s（k-1） +γ（k）LLXl=1syo，y（k，l）m，τ（k，l）-^s（k-1)!,（22）^Qθ、 ^s（k）= -ψ（θ）+D^s（k），φ（θ）E+常数。（23）引理1和2给出了两个条件分布p（τ| ym，yo；θ）和p（ym |τ，yo；θ）。基本上，要从mthem中采样，我们只需要从某些高斯分布和伽马分布中提取实现，这很简单。基于上述抽样方案，我们可以得到马尔可夫链的转移概率密度函数如下：∏θ（ym，τ，y′m，τ′）=p（τ′ym，yo；θ）p（y′m |τ′，yo；θ）。（24）引理1。给定ym、yo和θ，混合物重量{τt}彼此独立，即p（τ| ym，yo；θ）=TYt=2p（τt | ym，yo；θ）。（25）此外，τt遵循伽马分布：τt | ym，yo；θ~ γν+1，（yt- φ- ^1yt-1)/σ+ ν!.（26）证明：见Ap pendix A-A.Lemma 2。给定τ、yo和θ，缺失块yd=【ytd+1，ytd+2，…，ytd+nd】T，其中d=1，2，D、 相互独立，即p（ym |τ，yo；θ）=DYd=1p（yd |τ，yo；θ）。

（27）此外，ydonly的条件分布依赖于两个最近的观测样品YtAnd和ytd+nd+1以及Yd |τ，yo；θ ~ N（ud，∑d），（28），其中udud（i）的第i分量=i-1Xq=0ДqД+Дiytd+Piq=1Дi-2qτtd+qPnd+1q=1хnd+1-2qτtd+q×ytd+nd+1-ndXq=0хqх- ^1nd+1年初至今！，（29）和∑d∑d（i，j）的第i列和第j行中的分量=最小（i，j）Xq=1хi+j-2qτtd+q-Piq=1хi-2qτtd+qPjq=1Дj-2qτtd+qPnd+1q=1х-2qτtd+qσ、 （30）其中ud（i）中的几何级数之和可以简化为asi-1Xq=0ДqД=（iД，Д=1，Д（Дi-1)φ-1，Д6=1，（31）和ndxq=0Дq=（nd+1）Д，Д=1，Д^1nd+1-1.φ-1, φ6= 1.（32）证明：获得近似值^Q后，参见Ap pendix A-B.B.M步骤θ、 ^s（k）在（23）中，我们需要将其最大化以更新估计值。函数^Qθ、 ^s（k）可以重写为^Qθ、 ^s（k）= -ψ（θ）+D^s（k），φ（θ）E+常数=（T- 1)νlogν- 日志Γν-日志σ+ν^s（k）-^s（k）2σ-Д^s（k）2σ-^s（k）2σ+^s（k）σ+^s（k）σ-И^s（k）σ+常数，（33），其中^s（k）i（i=1，2，…，7）是^s（k）的第i个分量。ν、Д和σ的优化与ν的优化分离。设置^Q的导数θ、 ^s（k）关于Д、Д和σ到0，给出了Д（k+1）=^s（k）- Д（k+1）^s（k）^s（k），（34）Д（k+1）=^s（k）^s（k）- ^s（k）^s（k）^s（k）^s（k）-^s（k）, （35）和σ（k+1）=T- 1.^s（k）+^1（k+1）^s（k）+^1（k+1）^s（k）- 2х（k+1）^s（k）- 2Д（k+1）^s（k）+2Д（k+1）Д（k+1）^s（k）.（36）ν（k+1）可通过以下公式找到：ν（k+1）=arg maxν>0fν、 ^s（k）（37）带fν、 ^s（k）=νlogν- 日志Γν+ν^s（k）2（T-1). 根据文献[30]中的命题1，ν（k+1）总是存在且不唯一。如文献[30]所述，最大化子ν（k+1）可以通过一维搜索获得，如半间隔法[31]。由此产生的SAEM-MCMC算法是经过归纳的inAlgorithm 1。C、 特殊情况下，在θ中的某些参数已知的情况下，我们只需要相应地改变M步中的upd，并且模拟和近似步骤保持不变。

例如，如果我们知道时间序列是零均值[1]，[12]，即φ=0，则应将φ（k+1）和φ（k+1）的更新替换为φ（k+1）=0，（38）算法1 SAEM-MCMC算法，用于学生的t A R（1）1：初始化θ（0）∈ Θ，^s（0）=0，k=0，y（0，l）mfor l=1，2，L2： 对于k=1，2。do3：模拟：4：对于l=1，2，L do5：来自p的样本τ（k，L）τ| y（k-1，l）m，yo；θ（k）usingLemma 1,6：样品y（k，l）mfor pym |τ（k，l），yo；θ（k）usingLemma 2.7：结束8：随机预测：评估^s（k）和^Qθ、 ^s（k）分别如（22）和（23）所示。9： 最大化：更新θ（k+1），如（3-4）、（35）、（36）和（37）所示。10： 如果满足停止c标准，则11：终止循环12：结束if13：结束forandД（k+1）=^s（k）^s（k），（39）如果已知时间序列遵循随机游走模型[14]，这是AR（1）模型的特殊情况，且Д=1，则应将Д（k+1）和Д（k+1）的更新替换为Д（k+1）=^s（k）- ^s（k）^s（k），（40）和Д（k+1）=1。（41）D.AR（p）的泛化上述ML估计方法可立即泛化为学生的t AR（p）模型：yt=Д+pXi=1Дiyt-i+εt，（42），其中εti。i、 d。~ t型0, σ, ν. 同样，我们可以应用SAEMMCM算法，通过考虑τandymas潜在数据和yoas观测数据来获得估计。在每次迭代中，我们根据条件分布p得出τ和ym的一些实现ym，τ| yo；θ（k）近似经验函数Qθ|θ（k）, 最大化近似值^Qθ|θ（k）更新估算。主要区别在于，AR（p）的条件分布将比AR（1）的条件分布更加复杂，因为AR（p）的每个样本都更依赖于之前的样本。

为了应对这一挑战，在应用吉布斯抽样时，我们可以将潜在数据（ym，τ）分成更多的块，τ作为一个块，每个块∈CMA是一个块，因此，以其他潜在变量为条件的每个潜在变量块的分布将很容易获得并从中采样。由于空间的限制，我们在这里不详细介绍，我们将在未来的工作中考虑这一点。五、 收敛在本节中，我们为所提出的算法m的收敛性提供了理论保证。简单确定性EM算法的收敛性已由许多不同的作者解决，从[23]中的基本工作开始，到[24]中的一个更一般的考虑。然而，由于采样的随机性，EMalgorithm的随机变量（如MCEM、SAEM和SAEM-MCMC算法）的收敛性分析具有挑战性。有关这些随机EM算法及其收敛性分析的更一般概述，请参见[18]、[19]、[32]–[35]。出于特殊的兴趣，作者在[18]中介绍了SAEMalgorithm，并建立了在mildaditional条件下几乎肯定会收敛到观测数据（如od）的静态点。文[19]中的作者将SAEM框架与MCMC过程耦合，并给出了SAEM-MCMC算法在完全数据似然为指数函数时的收敛条件。

在我们的例子中，给定的一组条件如下。（M1）对于任何θ∈ Θ，Z Zks（yo，ym，τ）kp（ym，τ| yo；θ）dymdτ<∞.（43）（M2）ψ（θ）和φ（θ）在Θ上是两次连续可微的。（M3）函数\'s（θ）=Z Zs（yo，ym，τ）p（ym，τ| yo；θ）dymdτ（44）在920;上连续可微。（M4）目标函数l（θ；yo）=logZ Zp（y，τ；θ）dymdτ（45）在Θ上连续可微，且θZ Zp（y，τ；θ）dymdτ=Z Zθp（y，τ；θ）dymdτ。（46）（M5）对于Q（θ，’s）=-ψ（θ）+h's，φ（θ）i+常数。，存在一个函数|θ（(R)s），使得\'s和θ ∈ Θ，Q°θ（\'s），\'s≥ Q（θ，’s）。此外，函数θ（(R)s）是连续可微的。（SAEM1）对于所有k，γ（k）∈ [0，1]，P∞k=1γ（k）=∞ 存在<λ≤ 1此类P∞k=1γ（k）1+λ<∞.（SAEM2）l（θ；yo）在Θ上是d次可微的，其中d=7是s（yo，ym，τ）的维数，而∧θ（s）是d次可微的。（SAEM3）1）链在紧集中获取其值Ohm.2） s（yo，ym，τ）有界于Ohm, 和序列^s（k）在com pactsubset中获取其值。3） 对于Θ的任何紧子集V，都存在面积常数L，因此对于任何θ, θ′在Vsup中（ym，τ，y′m，τ′）∈Ohmπθ（ym，τ，y′m，τ′）- πθ′（ym，τ，y′m，τ′）≤ L |θ- θ′|.（47）4）转移概率∏θ生成一个多遍历链，其不变量概率为条件分布p（ym，τ| yo；θ）。总之，条件（M1）-（M5）都是关于模型的，并且是de TerministicEM算法收敛的条件。条件（M1）和（M3）要求充分统计量的期望值有界且连续可微。条件（m2）和（M4）保证了完全数据似然l（θ；y，τ）的连续可微性，以及完全数据似然Q的期望θ|θ（k）, 观测数据的对数似然度（θ；yo）。

条件（M5）表明Q（θ，’s）存在全局最大化子。条件（SAEM1）-（SAEM3）是SAEM-MCMC收敛的附加要求。条件（SAEM1）与步长有关γ（k）. 通过选择合适的步长，可以很容易地满足这一条件。建议将γ（k）=1设为1≤ k≤ K和γ（K）=K-k或k≥ K+1，其中K是一个正整数，因为初始猜测θ（0）可能远离我们正在寻找的ML估计，选择第一个K步长等于1，则序列nθ（K）oto会有很大的变化，然后收敛到最大似然的邻域[27]。条件（SAEM2）要求d=l（θ；yo）和^θ的7倍可微性^s（k）. 条件（SAEM3）对生成的马尔可夫链施加了一些约束。在[19]中，作者建立了AEM MCMC算法对平稳点的c收敛性。然而，他们的分析假设，完整的数据可能性属于曲线指数族，所有这些条件（M1）（M5）和d（SAEM1）-（SAEM3）都是满足的。这些假设是非常特殊的问题，对于我们的情况并不简单，因为我们的潜变量的条件分布非常复杂。为了评价我们提出的算法的收敛性，我们需要逐个建立条件（M1）（M5）和d（SAEM1）-（SAEM3）。最后，我们得到了我们所提出算法的收敛结果，总结如下定理。定理1。

假设参数空间Θ被设置为一个很大的边界集，参数ν>2，由（25）和（27）生成的马尔可夫链取紧集中的值，由算法生成的序列nθ（k）这意味着（33）的无约束最大化子（由（34）、（35）、（36）和（37）给出）位于该边界集中。理论上，由（25）和（27）生成的马尔可夫链取无界集中的值。然而，在实践中，链不会取非常大的值，我们可以考虑链取非常大的紧集中的值【19】，【27】。1具有以下渐近性质：概率为1，limk→+∞dθ（k），L= 0，其中dθ（k），L表示从θ（k）到观测数据的平稳点集的距离，类似于对数d L=nθ∈ Θ ,l（θ；yo）θ=0°。证明：有关条件（M1）-（M5）和（SAEM2）-（SAEM3）的证明，请参阅附录B。如前所述，通过正确选择步长，可以轻松满足条件（SAEM1）。在建立这些条件后，这个定理的证明遵循了[19]中对工作的分析。六、 仿真在本节中，我们对所提出的ML估计器的性能和所提出算法的收敛性进行了仿真研究。首先，我们表明，所提出的估计器能够从不完整的时间序列中很好地估计参数，这些时间序列已被系统化以拟合模型。其次，我们展示了它对外部创新的鲁棒性。最后，我们用一个真实的金融时间序列，即恒生指数对其进行了检验。A、 参数估计在这一小节中，我们展示了ProposedSAEM MCMC算法的收敛性以及propo-SED估计器在不同样本数和缺失百分比的不完全Student t AR（1）时间序列上的性能。

估计误差由平方误差（MSE）测量：MSE（θ）：=E^θ - θ真,其中，^θ是参数θ的估计值，θtrueis是其真值。参数θ可以是Д、Д、σ和ν。使用100个独立的不完整时间序列，通过蒙特卡罗模拟来近似期望值。我们将Дtrue=1，Дtrue=0.5，（σtrue）=0.01，νtrue=2.5。对于每个不完整的数据集yo，我们首先根据Student的T AR（1）模型生成一个具有T个样本的完整时间序列{yt}。然后随机删除nmisnumber样本以获得不完整的时间序列。不完整时间序列的缺失百分比ρ：=nmisT×100%。在第五节中，我们建立了所提出的SAEM-MCMC算法对观测数据似然平稳点的收敛性。然而，据观察，由于多个平稳点的存在，由算法m获得的估计结果可能对初始化敏感。这是一个不可避免的问题，因为它是一个非对流优化问题。有趣的是，我们还观察到，与随机初始化相比，当我们使用ML估计对G a ussian AR（1）模型进行初始化时，最终估计得到了显著改善。文献[13]介绍了不完全数据中高斯AR模型的最大似然估计，通过确定性EM算法可以很容易地得到估计值。我们初始化Д（0）、Д（0）和σ(0)使用1.00 1.04 1.08 1.12φ00.44 0.46 0.48 0.50φ10.010 0.020 0.030σ20 40 60 80 1002迭代kν图1中的估计值。

估算与迭代。表I未完成学生t AR的估计结果（1）。（σ）真值1.000 0.500 0.010 2.5高斯AR（1）1.119 0.442 0.033+∞学生t AR（1）0.989 0.501 0.009 2.234高斯AR（1）模型（Д）g，（Д）g，和σg、 利用条件分布的me a n初始化y（0，l）ym；yo，（Д）g，（Д）g，σg级, 这是一个高斯分布。参数ν（0）初始化为随机正枚举。在每次迭代中，我们绘制L=10个样本。对于步骤大小，我们将γ（k）=1设为1≤ k≤ 30和dγ（k）=k-30对于k≥ 图1给出了应用所提出的SAEM-MCMC算法估计合成AR（1）数据集（T=300，缺失百分比ρ=10%）参数的示例。我们可以看到，该算法的收敛时间小于10次迭代，其中每次迭代只需要sL=10次吉布斯采样，而且最终估计误差很小。表一比较了Student的t A R模型和Gaussian AR模型的估计结果。这证明了我们的观点，即对于不完全重尾数据，传统的不完全高斯AR时间序列方法太不有效，并且可以通过在重尾模型下设计算法来实现显著的性能增益。图2显示了φ0ρ=40%ρ=30%ρ=20%ρ=10%ρ=00.000 0.004 0.008MSEφ10.0e+00 4.0e的数字0.00 0.01 0.02 0.03 0.04MSE的估算结果-06 8.0e-06 1.2e-σ2100 200 300 400 5000.0 0.2 0.4 0.6 0.8 1.0的05MSEν图2的样本数TMSE。具有不同样本数和缺失百分比的不完整时间序列的均方误差。样品T=100、200、300、400、500，缺失百分比ρ=10%、20%、30%、40%。