基于经验似然INAR(1)模型离群点的检测与估计

Detection and Estimation of Outliers in the Empirical Likelihood INAR(1) Model

摘要: 整值时间序列可出现在教育，金融，医疗，交通等诸多领域，本文旨在研究利用经验似然方法对整值时间序列中的加性离群点与新息离群点进行检测与估计，并针对凸包问题进行了详细讨论，最后通过数值模拟充分验证了经验似然方法检测离群点的有效性。仿真实验结果表明，经验似然方法可以有效检测与估计出不同新息分布下一阶整值时间序列模型中的离群点。

Abstract: Integer-valued time series can appear in various fields such as education, finance, healthcare, and transportation. This paper aims to investigate the detection and estimation of additive outliers and innovation outliers in integer-valued time series based on the empirical likelihood method. Additionally, the convex hull problem is discussed in detail. Finally, numerical simulations are conducted to fully verify the effectiveness of the empirical likelihood method in detecting outliers. The simulation results show that the empirical likelihood method can effectively detect and estimate outliers in first-order integer-valued time series models with different innovation distributions.

1. 引言

经验似然是一种众所周知的方法，通过经验似然方法获得的置信区间没有对称性约束，因此可以更好地描述分布的真实形状。它有类似于Bootstrap的抽样特性，它允许在不指定数据概率分布的情况下对观察样本进行统计推断。它基于样本数据，通过估计其概率分布的最大似然参数来推断总体参数的值。具体来说，经验似然方法将样本数据看作样本的概率密度函数，并通过最大化该密度函数来估计总体参数。经验似然方法常用于建立概率模型、参数优化、假设检验等问题中。它是一种基本的统计方法，具有广泛的应用和重要的理论意义。

自从Owen [1] [2]的开创性成果以来，经验似然方法受到学术界广泛的关注。在理论研究中，Owen [3]首先提出经验似然方法并应用到线性模型中，Qin和Lawless [4]基于无偏估计函数推进发展了经验似然方法，并将经验似然方法应用于广义估计方程。Tang和Leng [5]提出了一种利用经验似然的方法来同时进行参数估计和变量选择。这种方法通过优化经验似然函数来达到目的，并引入了惩罚项来约束参数的数量和复杂度。Shi和Lau [6]、Wei [7]等人以及Liu和Yuan [8]等学者在各自的研究中将经验似然方法应用于不同的领域，展示了其在实际应用中的有效性。与此同时，Kitamura [9]和Monti [10]等人也利用经验似然方法对相依数据进行了分析和预测工作。他们的研究表明，经验似然方法在数据科学领域的广泛应用，为研究人员提供了一种有效的工具来处理复杂的数据模型，帮助他们更好地理解并预测数据的行为。

在时间序列领域中，Chuang和Chan [11]以及Chan和Ling [12]在其研究中采用经验似然方法，针对不同的假设条件下自回归模型进行了建模工作。其中，他们建立了在平稳和非平稳情况下的对数经验似然比统计量的极限分布。在研究中，Zhao和Wang [13]应用了经验似然方法，提供了自回归模型中回归参数的置信区间估计。他们的研究结果表明，对数经验似然比统计量在一定条件下收敛于卡方分布，为在自回归模型中对参数估计提供了理论依据。Nordman和Lahiri [14]给出了全面系统性的回顾和总结。

检测数据集中离群值的位置并估计其大小是统计学中比较常见的主题之一。Fox [15]首先在时间序列模型中提出加性离群点(Additive Outlier, AO)以及新息离群点(innovation Outlier, IO)，在此之后的时间里，对于离群点检测的方法研究受到了学术界广泛关注。如Tsay [16]提出了一个迭代方法来识别离群点，以及消除它们对时间序列的影响。近几年研究人员的注意力转移到整值计数时间序列中离群值的检测上。Fokianos和Fried [17]提出了整数值广义条件异方差(INGARCH)模型框架下两类离群点的检测和估计方法。Fried [18]等人提出了在存在离群点和干预效应的情况下，采用广义线性模型对计数时间序列进行稳健拟合。INAR(1)是现有文献中提出的最成功的整值时间序列模型之一，Barczy [19] [20]分别考虑了在已知时间段内受创新离群点(IO)和加性离群点(AO)污染的INAR(1)模型参数的条件经验似然估计。Silva [21]研究了INAR(1)模型中含有加性离群点时的检测问题，并对加性离群点(AO)进行检测和定位，并同时对离群值大小进行估计。Baragona等人[22]提出了使用经验似然方法检验AR(p)模型中的离群点Hua Shang [23]提出了吉布斯抽样算法来估计在泊松INAR(1)模型下的参数和离群点的大小。本研究首次将经验似然方法系统性地引入一阶整值自回归(INAR(1))模型的离群点分析领域，以及首次提出了将平衡经验似然用于INAR(1)模型中进行优化，经验似然方法创新性地构建了基于约束优化与非参数似然比统计量的离群点联合检测框架，为后续离群点检测研究提供了理论支持。

本文的结构是，在第二节中，详细介绍了经验似然方法，第三节给出了含有加性离群点与新息离群点的INAR(1)模型，并给出了估计统计量的渐近结果以及经验似然的改进方法——平衡经验似然。在第四节中通过模拟实验及其分析验证了本文所提经验似然方法的有效性，第五节介绍了经验似然方法检测离群点在实际生活中的应用，第六节为结论与展望，对全文进行总结，并提出未来可研究的方向。

2. 经验似然

经验似然(EL)方法是由Owen [2]最早提出的一种非参数统计推断方法，与传统的参数统计方法不同，它无需事先假设数据服从某个特定的分布族。这种特性使得EL方法在处理实际问题时更具灵活性和广泛的适用性。通过直接计算总体分布的各阶矩，EL方法避免了对参数渐近方差的估计，从而能够提供更为准确的统计推断结果。

设随机变量序列$X_1,X_2,\cdots ,X_n$ 具有独立同分布(i.i.d.)的特性，其具有共同且未知的概率分布$F_{\theta }$ ，感兴趣的未知参数记为$\theta ={\left({\theta }_1,\cdots ,{\theta }_p\right)}^{\text{T}}$ ，那么此时经验似然定义为：

$L\left(F_{\theta }\right)={\displaystyle \prod _{i=1}^{n}d}{F}_{\theta }\left(X_i\right)={\displaystyle \prod _{i=1}^n{p}_i},$ (1)

其中，$p_i=d{F}_{\theta }\left(X_i\right)=\mathrm{Pr}\left(X=X_i\right)$ ，$0\le p_i\le 1$ ，${\displaystyle \sum _{i=1}^n{p}_i}=1$ 。

而经验累积分布函数$F_n\left(x\right)=n^{-1}{\displaystyle \sum _{i=1}^{n}I}\left(X_i<x\right)$ 是$F_{\theta }$ 的非参数极大似然估计，经验似然比(ELR)函数定义为：

$R\left(F_{\theta }\right)=\frac{L\left(F_{\theta }\right)}{L\left(F_n\right)}={\displaystyle \prod _{i=1}^{n}n}{p}_i.$ (2)

若对参数$\theta$ 进行估计推断，通常需构造r个估计方程满足：$E_F\left\{g_j\left(X,\theta \right)\right\}=0$ ，$g_j\left(X,\theta \right)$ ，$j=1,2,\cdots ,r$ ，$r\ge p$ ，将向量表达式设为：

$g\left(X,\theta \right)={\left(g_1\left(X,\theta \right),\cdots ,g_r\left(X,\theta \right)\right)}^{\text{T}},$

此时经验似然比函数为：

$R\left(\theta \right)=\sup \left\{{\displaystyle \prod _{i=1}^{n}n}{p}_i|p_i\ge 0,{\displaystyle \sum _i^n{p}_i}=1,{\displaystyle \sum _{i=1}^n{p}_i}g\left(X,\theta \right)=0\right\}.$ (3)

由拉格朗日乘子法，可作辅助函数：

$H={\displaystyle \sum _i^n\log p_i}+\lambda \left(1-{\displaystyle \sum _i^n{p}_i}\right)-n{t}^{\text{T}}{\displaystyle \sum _i^n{p}_i}g\left(X_i,\theta \right),$ (4)

其中，$\lambda$ 和$t={\left(t_1,t_2,\cdots ,t_r\right)}^{\text{T}}$ 均为拉格朗日乘子。将辅助函数(4)对其$p_i$ 求偏导，得到$p_i=\frac{1}{n}\cdot \frac{1}{1+t^{\text{T}}g\left(X_i,\theta \right)}$ ，并且$t$ 满足${\displaystyle \sum _{i=1}^n\frac{1}{n}}\cdot \frac{g\left(X_i,\theta \right)}{1+t^{\text{T}}g\left(X_i,\theta \right)}=0$ 。

由此，可得到对数经验似然比为：

$l_E\left(\theta \right)=-2\log R_E\left(\theta \right)=2{\displaystyle \sum _{i=1}^n\log }\left(1+t^{\text{T}}g\left(X_i,\theta \right)\right).$ (5)

Owen (1991)证明了，当$n$ 趋于无穷大时，若$\theta ={\theta }_0$ ，上述对数经验似然比统计量$l_E\left(\theta \right)$ 渐近服从自由度为$p$ 的${\chi }^2$ 分布。

3. 基于经验似然INAR(1)模型的离群点检测

本章针对INAR(1)模型的离群点检测问题。INAR(1)模型作为离散时间序列的一种特殊形式，广泛应用于描述计数数据(如事件发生次数、事故发生频率等)。在该模型中，离群点的存在可能会影响自回归参数的估计，进而导致不准确的预测和决策。因此，本章提出了一种基于经验似然方法对加性离群点(AO)与新息离群点(IO)的检测方案。

3.1. 含有离群点的INAR(1)模型及经验似然

首先介绍传统INAR模型，INAR模型是一种常见的离散时间模型，广泛应用于金融、生物医学、通信等领域。其研究方法是通过稀疏算子构建模型，其中基于二项稀疏算子的整值一阶自回归INAR(1)模型是最为经典的模型之一。其定义如下：

$Z_t=\alpha \circ Z_{t-1}+{\epsilon }_t,t=0,\pm 1,\pm 2,\cdots$ (6)

其中，$Z$ 是非负的整数值随机变量，$\alpha \in \left[0,1\right)$ ，$\left\{{\epsilon }_t\right\}$ 为独立同分布且均值为$\mu$ ，方差为${\sigma }^2$ 的取值为整数值的随机变量序列，“$\circ$ ”为二项稀疏算子，其被定义为$\alpha \circ Z_{t-1}={\displaystyle \sum _{i=1}^{Z_{t-1}}{Y}_i}$ ，这里$\left\{Y_i\right\}$ 是一列独立同分布的以$\alpha$ 为参数的伯努利随机变量，并独立于$Z$ ，其分布律为$P\left\{Y_i=1\right\}=1-P\left\{Y_i=0\right\}=\alpha$ 。

在时间序列中，离群点的种类繁多，一般分为四类：加性离群点(AO 类)、新息离群点(IO类)、水平移位离群点(LS类)和暂时变更离群点(TC类)。

其中，加性离群点(AO)是指在某个特定时间点发生的外部误差或外生变化，这种误差或变化只影响到该时间点的观测值，而不对其他时间点的观测值产生影响，这种离群点通常是由于观测或记录过程中的错误或误差所导致。新息离群点(IO)是由内部变化或噪声过程中的内源性因素引起的，其影响会传播到所有后续的观测值，并且其影响方式与系统的动态模型密切相关。水平移位离群点(LS)的出现会导致序列均值的偏移，从而影响后续所有观测值。这种离群点的产生通常是由于干扰造成的，进而改变了系统的整体结构，进而对时间序列产生持久的影响。暂时变更离群点(TC)的出现不仅会影响当前的观测值，还会对后续所有观测值产生影响，不过这种影响的程度会随着时间的推移而呈指数衰减。本文主要研究AO与IO两种离群点形式的检测。

定义1 随机过程${\left(Y_n\right)}_{n\in {\mathbb{Z}}_{+}}$ 称为含有有限多个加性离群点(AO)的INAR(1)时间序列模型，若满足：

$Y_k=Z_k+{\displaystyle \sum _{i=1}^I{c}_{k,s_i}}{\omega }_i,k\in {\mathbb{Z}}_{+}$ (7)

其中，${\left(Z_n\right)}_{n\in {\mathbb{Z}}_{+}}$ 为满足INAR(1)模型的序列，${\omega }_i$ 为离群点的大小；$c_{k,s_i}$ 为一个脉冲，当$k=s_i$ 时取值为1，否则为0。此模型假设只有在特定的时间点$s_i$ 发生离群扰动，并且扰动的大小由${\omega }_i$ 决定

若单个AO位置任意，假设($I=1$ ，$s_1:=s$ )，由公式(7)可以得出：

$\begin{array}{c}\text{E}\left(Y_k|{ℱ}_{k-1}^Y\right)=\alpha X_{k-1}+{\mu }_{\epsilon }+c_{k,s}\omega \\ =\alpha \left(Y_{k-1}-c_{k-1,s}\omega \right)+c_{k,s}\omega \\ =\alpha Y_{k-1}+{\mu }_{\epsilon }+\left(-\alpha c_{k-1,s}+c_{k,s}\right)\omega \\ =\{\begin{array}{ll}\alpha Y_{k-1}+{\mu }_{\epsilon }\hfill & \text{if}k=1,\cdots ,s-1\hfill \\ \alpha Y_{k-1}+{\mu }_{\epsilon }+\omega \hfill & \text{if}k=s\hfill \\ \alpha Y_{k-1}+{\mu }_{\epsilon }-\alpha \omega \hfill & \text{if}k=s+1\hfill \\ \alpha Y_{k-1}+{\mu }_{\epsilon }\hfill & \text{if}k\ge s+2\hfill \end{array}.\end{array}$ (8)

(8)公式表明，在时间点$s$ 发生离群扰动时，模型的期望值会在其前后位置发生变化。

因此，计算其残差平方和得到：

$\begin{array}{c}{\displaystyle \sum _{k=1}^n{\left(Y_k-\text{E}\left(Y_k|{ℱ}_{k-1}^Y\right)\right)}^2}={\displaystyle \sum _{k=1}^n{}^{\left(s,s+1\right)}}{\left(Y_k-\alpha Y_{k-1}-{\mu }_{\epsilon }\right)}^2+{\left(Y_s-\alpha Y_{s-1}-{\mu }_{\epsilon }-\omega \right)}^2\\ +{\left(Y_{s+1}-\alpha Y_s+{\mu }_{\epsilon }+\alpha \omega \right)}^2.\end{array}$ (9)

为了对模型中的参数$\alpha$ ，$\mu$ ，$\omega$ 进行最优化估计，可以使用最小二乘法(LS)来获得。针对模型(7)，可以通过一组估计方程来获得参数的最小二乘估计。

最小二乘估计可以写成如下形式：

$\begin{array}{l}{g}_1\left(Y_k,\alpha ,\mu ,\omega \right)=-2{\displaystyle \sum _{k=1}^n{Y}_{k-1}}\left(Y_k-\alpha Y_{k-1}-{\mu }_{ϵ}\right)-2{\displaystyle \sum _{k=s+1}^n{Y}_{k-1}}\left(Y_k-\alpha Y_{k-1}-{\mu }_{ϵ}\right)\\ +2\left(Y_{s+1}-\alpha Y_s-{\mu }_{ϵ}+\alpha \omega \right)\left(-Y_s+\omega \right)\\ g_2\left(Y_k,\alpha ,\mu ,\omega \right)=-2{\displaystyle \sum _{k=1}^n\left(Y_k-\alpha Y_{k-1}-{\mu }_{ϵ}\right)}-2{\displaystyle \sum _{k=s+1}^n\left(Y_k-\alpha Y_{k-1}-{\mu }_{ϵ}\right)}\\ -2\left(Y_{s+1}-\alpha Y_s-{\mu }_{ϵ}+\alpha \omega \right)\\ g_3\left(Y_k,\alpha ,\mu ,\omega \right)=-2\alpha \left(Y_{s+1}-\alpha Y_s-{\mu }_{ϵ}+\alpha \omega \right).\end{array}$ (10)

定义2 令${\left({\epsilon }_k\right)}_{k\in \mathbb{N}}$ 为独立同分布(i.i.d)非负整值随机变量序列，随机过程${\left(Y_n\right)}_{n\in {\mathbb{Z}}_{+}}$ 称为含有有限多个新息离群点(IO)的INAR(1)时间序列模型，若满足

$Y_k={\displaystyle \sum _{j=1}^{Y_{k-1}}{\xi }_{k,j}}+{\eta }_k,k\in \mathbb{N}.$ (11)

其中，${\left({\xi }_{k,j}\right)}_{j\in \mathbb{N}}$ 为独立同分布的均值$\alpha \in \left(0,1\right)$ 的伯努利随机变量，且与${\left({\epsilon }_l\right)}_{l\in \mathbb{N}}$ 相互独立，$Y_0$ 为独立于${\left({\xi }_{k,j}\right)}_{j\in \mathbb{N}},k\in \mathbb{N}$ 以及${\left({\epsilon }_l\right)}_{l\in \mathbb{N}}$ 的非负整值随机变量，且

${\eta }_k={\epsilon }_k+{\displaystyle \sum _{i=1}^I{c}_{k,s_i}}{\omega }_i,k\in {\mathbb{Z}}_{+}.$

其中，$I\in \mathbb{N},s_i,{\omega }_i\in \mathbb{N},i=1,\cdots ,I$ ；${\omega }_i$ 为离群点大小；$c_{k,s_i}$ 为一个脉冲，当$k=s_i$ 时取值为1，否则为0。

由此可见，若单个IO位置任意，假设($I=1$ ，$s_1:=s$ )，由公式(11)，可以得出

$\text{E}\left(Y_k|{ℱ}_{k-1}^Y\right)=\alpha Y_{k-1}+\text{E}{\eta }_k=\alpha Y_{k-1}+{\mu }_{\epsilon }+c_{k,s}\omega ,k\in \mathbb{N}$ .

因此，对于所有，$n\ge \max \left(3,s+1\right)$ ，都有

${\displaystyle \sum _{k=1}^n{\left(Y_k-\text{E}\left(Y_k|{ℱ}_{k-1}^Y\right)\right)}^2}={\displaystyle \sum _{k=1}^n{\left(Y_k-\alpha Y_{k-1}-{\mu }_{\epsilon }\right)}^2}+{\left(Y_s-\alpha Y_{s-1}-{\mu }_{\epsilon }-\omega \right)}^2.$ (12)

通过分别对其$\alpha$ ，$\mu$ ，$\omega$ 求偏导得出最小二乘估计方程为：

$\begin{array}{l}{g}_1\left(Y_k,\alpha ,\mu ,\omega \right)=-2{\displaystyle \sum _{k=1}^n{Y}_{k-1}}\left(Y_k-\alpha Y_{k-1}-{\mu }_{ϵ}\right)-2{\displaystyle \sum _{s=1}^n{Y}_{s-1}}\left(Y_s-\alpha Y_{s-1}-{\mu }_{ϵ}-\omega \right)\\ g_2\left(Y_k,\alpha ,\mu ,\omega \right)=-2{\displaystyle \sum _{k=1}^n\left(Y_k-\alpha Y_{k-1}-{\mu }_{ϵ}\right)}-2{\displaystyle \sum _{s=1}^n\left(Y_s-\alpha Y_{s-1}-{\mu }_{ϵ}-\omega \right)}\\ g_3\left(Y_k,\alpha ,\mu ,\omega \right)=-2{\displaystyle \sum _{s=1}^n\left(Y_s-\alpha Y_{s-1}-{\mu }_{ϵ}-\omega \right)}.\end{array}$ (13)

此时，针对两种不同类型的离群点构造，经验似然比(ELR)为：

$\ell \left(\alpha ,\mu ,\omega \right)=-2\underset{k}{\sup }\{{\displaystyle \sum _{k=2}^N\log }\left(N{w}_k\right):{\displaystyle \sum _{k=2}^N{w}_k}=1,{\displaystyle \sum _{k=2}^N{w}_k}{g}_j\left(Y_k,\alpha ,\mu ,\omega \right)=0,j=1,2,3\}.$ (14)

定理1 若$\left({\alpha }_0,{\mu }_0,{\omega }_0\right)$ 为模型参数的真值，那么当样本量趋于无穷时，统计量$\ell \left({\alpha }_0,{\mu }_0,{\omega }_0\right)$ 渐近服从以自由度为3的中心${\chi }^2$ 分布。

通过定理1对$\ell \left({\alpha }_0,{\mu }_0,{\omega }_0\right)$ 的渐近分布，很容易导出其置信区域，此外，还可以根据这个分布推导出参数检验，其定理的证明可以仿造Baragona [22]中的证明步骤，在此不给出具体的证明。若有多个离群值，可以仿造上述思路构建经验似然比统计量。

3.2. 平衡经验似然(BEL)

在经验似然方法中，解决凸包问题一直是个重要的研究方向，本文参考了文献中提出的几种解决方案(Chen等人，Emerson和Owen以及Tsao) [24]-[26]，研究表明，调节常数的有限样本校准是影响方法性能的关键因素。基于此，可以采用平衡经验似然(BEL)的方法，其优势在于可通过解析准则选择调节常数，在保持经验似然渐近性质的同时，减少了计算复杂性，并提高了对小样本数据的适应性。

在检测离群点问题中，估计方程约束系统的不可行性通常仅发生于最后一个矩条件。基于此特性，可以提出假设：${\overline{g}}_j\left(\alpha ,\mu ,\omega \right)\simeq 0$ ，$j=1,2,3$ 。此时，BEL方法通过加入两个人工观测值来简化约束系统，从而避免无解的情况：

$g\left(y_{N+1},\alpha ,\mu ,\omega \right)={\left[0,0,\cdots ,0,-\delta {\overline{g}}_2\left(\alpha ,\mu ,\omega \right)\right]}^{\text{T}}$

$g\left(y_{N+2},\alpha ,\mu ,\omega \right)={\left[0,0,\cdots ,0,\left(2+\delta \right){\overline{g}}_2\left(\alpha ,\mu ,\omega \right)\right]}^{\text{T}}.$

在该方法中，引入调节参数$\delta >0$ 的目的是确保所添加的人工数据不会对最终的求解结果造成显著的干扰。借助这种处理方式，原本因约束条件无法满足而无解的问题得以成功解决。

此外，考虑到${\overline{g}}_2\left(\alpha ,\mu ,\omega \right)$ 的量级相对于调节参数$\delta$ 来说可以忽略，因此，可以假设：

$\delta {\overline{g}}_2\left(\alpha ,\mu ,\omega \right)=\left(2+\delta \right){\overline{g}}_2\left(\alpha ,\mu ,\omega \right)=\Delta ,$

基于上述，BEL方法中加入的两个人工观测值可以写成：

$\begin{array}{l}g\left(y_{N+1},\alpha ,\mu ,\omega \right)={\left[0,0,\cdots ,0,-\Delta \right]}^{\text{T}}\\ g\left(y_{N+2},\alpha ,\mu ,\omega \right)={\left[0,0,\cdots ,0,\Delta \right]}^{\text{T}},\end{array}$

那么，平衡经验似然比(BELR)为：

$\begin{array}{l}{\ell }_{\Delta }^{*}\left(\alpha ,\mu ,\omega \right)=-2\underset{w}{\sup }\{{\displaystyle \sum _{k=2}^{N+2}\log }\left[\left(N+1\right)w_k\right]:{\displaystyle \sum _{k=2}^{N+2}{w}_k}=1;\\ {\displaystyle \sum _{k=2}^N{w}_k}{g}_j\left(y_k,\alpha ,\mu ,\omega \right)=0,j=1,2,3;\\ {\displaystyle \sum _{k=2}^N{w}_k}{g}_2\left(y_k,\alpha ,\mu ,\omega \right)-w_{N+1}\Delta +w_{N+2}\Delta =0\}.\end{array}$ (15)

常数$\Delta$ 的选取对BELR计算结果具有重要影响，如果选择$\Delta =o\left(N^{1/2}\right)$ 那么根据定理2的结果，${\ell }_{\Delta }^{*}\left({\alpha }_0,{\mu }_0,{\omega }_0\right)$ 在渐近意义下将与经验似然比(ELR)(14)相等。

定理2 当$\Delta =o\left(N^{1/2}\right)$ 时，ELR(14)与BELR(15)的差别随着样本量$N$的增大依概率收敛为0。

最后注意到，

${\displaystyle \sum _{k=2}^{N+2}{g}_2}\left(y_k,\alpha ,\mu ,\omega \right)={\displaystyle \sum _{k=2}^N{g}_2}\left(y_k,\alpha ,\mu ,\omega \right).$

这表明在计算平衡经验似然比(BELR)时，加入的人工观测值对整体结果没有影响。由此可知，参数的最大平衡经验似然估计量(MBELE)与最大经验似然估计量(MELE)以及联合最小二乘估计是等价的。

当BELR在与真实参数不同的点上计算时，它往往会变得任意大，如定理3所示。

定理3 当$\omega \ne {\omega }_0$ ，随着N逐渐增大，BELR具有发散特性：

$N^{-1/3}{\ell }_{\Delta }^{*}\left({\alpha }_0,{\mu }_0,\omega \right)\to \infty$

定理2与3的证明可以仿造Baragona [22]中的证明步骤，在此不给出具体的证明。

这一特性揭示了在计算平衡经验似然比(BELR)时，若所采用的参数与真实值存在偏差，哪怕这种偏差微乎其微，BELR的值也会显著上升。由此可知，该方法对于参数估计的准确性极为敏感，一旦参数估计偏离真实值，BELR的显著变化即可作为模型误设检测的有力理论依据。

4. 模拟实验及其分析

在本节中，进行了一项模拟研究，生成1000组长度为$N=100$ 的模拟时间序列。所使用的一阶整值自回归模型如下：

模型1：一阶整值自回归模型$y_k=0.2\circ y_{k-1}+{\epsilon }_t$ 。

模型2：一阶整值自回归模型$y_k=0.5\circ y_{k-1}+{\epsilon }_t$ 。

模型3：一阶整值自回归模型$y_k=0.8\circ y_{k-1}+{\epsilon }_t$ 。

其中，$\left\{{\epsilon }_t\right\}$ 是独立同分布(i.i.d.)的取值为整数值的随机变量序列，且根据以下三种不同的概率分布生成：泊松分布(Poisson)，几何分布(Geometric)，二项分布(Binomial)。

本模拟采用了三种不同的方法来进行离群点检测，每种方法如下：

方法1：使用公式(14)中的ELR，其检验统计量为$\underset{\alpha ,\mu }{\inf }\ell \left(\alpha ,\mu ,0\right)$ 。

方法2：使用公式(15)中的BELR，其检验统计量为$\underset{\alpha ,\mu }{\inf }{\ell }_{\Delta }^{*}\left(\alpha ,\mu ,0\right)$ 。

方法3：使用贝叶斯检测法，参考Silva等人[21]的计算步骤。

在本研究中，采用一种顺序检测的程序，该程序在每个步骤中识别最可疑的观测值，并去除其对剩余数据的影响。具体程序如下：

对每个时间点$s$ ($2\le s\le N$ )计算了相应的统计量，并计算所有时间点统计量的最大值。根据卡方分布的$1-\alpha /\left(N-p\right)$ 分位数进行假设检验，如果统计量的最大值超过了该分位数，则拒绝原假设，即认为数据中存在离群点，离群点发生的时间被假定为统计量最大值对应的时间点，并在此时间点处估计离群值的大小再进行调整去除，在调整后的数据序列$\left\{{y^{\prime }}_k\right\}$ 上重复上述检测程序，直到没有显著的统计量为止。

首先，考虑泊松分布生成的模拟时间序列。在1000次重复实验中，对于检验无离群点的假设与存在加性离群点的假设之间的对比，理论显著性水平$\alpha =0.05$ ，每种模型和方法的观察频率如表1所示。可以看出，观察到的大小总是非常接近名义值，贝叶斯方法相较于其他两种方法略差一些，对于方法2 (BELR)，其结果稍微更准确。

Table 1. Observed relative frequency of rejection of the hypothesis of the absence of outliers on 1000 series generated according to Poisson innovations

表1. 基于泊松创新生成的1000组序列中，拒绝无离群点假设的观测相对频率

模型	方法
	1	2	3
1	0.078	0.064	0.073
2	0.084	0.067	0.092
3	0.072	0.078	0.081

随后，在每个模拟系列的时间点$q=50$ 插入一个加性离群点，大小为${\omega }_0=3.5$ ，并重复检验，得到了表2所示的观察功效。并且区分了检测到的离群点是否正确定位在时间50。

Table 2. Observed relative frequency of rejection of the hypothesis of the absence of outliers on 1000 series generated according to Poisson innovations ($q=50$ , ${\omega }_0=3.5$ )

表2. 基于泊松序列创新生成1000组包含加性离群点的序列中，无离群点的假设被拒绝的观测相对频率，其中$q=50$ ，${\omega }_0=3.5$

模型	时间位置	方法
		1	2	3
1	准确	0.643	0.818	0.532
	错误	0.034	0.020	0.044
2	准确	0.382	0.671	0.431
	错误	0.043	0.029	0.072
3	准确	0.552	0.842	0.623
	错误	0.048	0.035	0.423

首先，我们注意到，使用方法1与方法2 (即ELR与BELR)的结果差距较为明显，并且方法2 (BELR)错误定位的频率相对较低，效果更为准确。而方法3的检测功效相对其他两种方法相对较低。

接下来，考虑非泊松创新的情况。在模拟实验中，使用了模型1，并通过几何分布(Geometric)和二项分布(Binomial)生成创新项。我们应用上述三种方法进行检验，并将结果报告在表3中。

可以观察到，对于Geometric分布生成的创新项，方法2 (BELR)观测到的检验大小方与泊松创新的情况没有显著变化。然而，方法1 (ELR)与方法3观察到的检验大小明显增加。对于Binomial分布生成的创新项，结果则显示出更加严重的问题：除了方法2的观测大小依然膨胀外，方法1观察到的检验大小是原来的六倍，方法3观测到的检验大小也出现了显著膨胀，这表明这些方法在面对Binomial分布的创新项时，均存在显著的偏差。

若考虑在非泊松创新的情况下，依旧向每个模拟序列的时间点$q=50$ 插入大小为${\omega }_0=3.5$ 的加性离群点，并重复进行检测，所得到的观察功效如表4所示，可以得知，对于Geometric分布生成的创新项，方法1，方法2与方法3的观察功效与泊松创新的情形相似，增加离群值的大小并未导致功效有显著提升，对其产生的影响是有限的，但方法3的功效略差于其他两种方法，方法2错误定位的频率要明显低于方法1与方法3，在检测到离群点时能更准确地定位到正确的时间点。对于Binomial分布生成的创新项，三种方法均存在显著偏差，但方法1和方法2的效果更优于方法3。

Table 3. Observed relative frequency of rejection of the hypothesis of the absence of outliers on 1000 series generated according to model 1 and non-Poisson innovations

表3. 基于模型1和非泊松创新生成的1000组序列中，拒绝无离群点假设的观测相对频率

创新	方法
	1	2	3
Geometric	0.254	0.078	0.332
Binomial	0.478	0.247	0.572

Table 4. Observed relative frequency of rejection of the hypothesis of the absence of outliers on 1000 series with an additive outlier at $q=50$ and size ${\omega }_0=3.5$ generated according to model 1 and non-Poisson innovations

表4. 基于模型1和非泊松创新生成的1000组序列中，假设无离群点的假设被拒绝的观测相对频率，其中$q=50$ ，${\omega }_0=3.5$

创新	时间位置	方法
		1	2	3
Geometric	准确	0.684	0.859	0.578
	错误	0.041	0.032	0.052
Binomial	准确	0.287	0.543	0.256
	错误	0.032	0.029	0.042

在接下来的模拟实验中，继续考虑非泊松分布情况。具体而言，我们采用了模型2，并利用几何分布(Geometric)和二项分布(Binomial)来生成创新项。随后，我们运用上述三种方法对这些数据进行了检验，并将检验结果汇总在表5中。

Table 5. Observed relative frequency of rejection of the hypothesis of the absence of outliers on 1000 series generated according to model 2 and non-Poisson innovations

表5. 基于模型2和非泊松创新生成的1000组序列中，拒绝无离群点假设的观测相对频率

创新	方法
	1	2	3
Geometric	0.342	0.087	0.432
Binomial	0.562	0.278	0.572

我们观察到，当创新项由几何分布(Geometric)生成时，方法2 (BELR)的检验结果与泊松分布创新项的情况基本一致，检验大小仅稍微增大。而方法1 (ELR)与方法3的检验大小却明显增大。

当创新项由二项分布(Binomial)生成时，问题更为严重。方法2的检验大小仍然存在膨胀现象，而方法1与方法3的检验大小更是可以观察到显著增加。这说明在处理二项分布(Binomial)创新项时，这三种方法都存在明显的偏差。

在非泊松创新的场景下，我们对每个模拟序列的时间点$q=50$ 插入大小为${\omega }_0=3.5$ 的加性离群值，并重复进行检测。检测结果如表6所示。分析发现：对于由几何分布(Geometric)生成的创新项，方法1，方法2与方法3的检测功效与泊松分布创新项的情况相近，方法3的功效略低于其他两种方法，可以表明尽管增加离群值的大小，但检测功效并未显著提升。此外，方法2在检测到离群点时，错误定位的频率明显低于方法1与方法3，能够更准确地定位到离群值的正确时间点。对于由二项分布(Binomial)生成的创新项，三种方法的检测结果都存在显著偏差。

Table 6. Observed relative frequency of rejection of the hypothesis of the absence of outliers on 1000 series with an additive outlier at $q=50$ and size ${\omega }_0=3.5$ generated according to model 2 and non-Poisson innovations

表6. 基于模型2和非泊松创新生成的1000组序列中，假设无离群点的假设被拒绝的观测相对频率，其中$q=50$ ，${\omega }_0=3.5$

创新	时间位置	方法
		1	2	3
Geometric	准确	0.375	0.685	0.411
	错误	0.052	0.025	0.068
Binomial	准确	0.266	0.382	0.213
	错误	0.043	0.033	0.051

对于模型3，我们依然考虑非泊松分布情况，并利用几何分布(Geometric)和二项分布(Binomial)来生成创新项，继续运用上述三种方法对这些数据进行了检验，并将检验结果汇总在表7中。我们观察到，当创新项由几何分布(Geometric)生成时，方法2 (BELR)的检验结果与泊松分布创新项的情况基本一致，而方法1 (ELR)与方法3的检验大小显著增大。当创新项由二项分布(Binomial)生成时，问题依旧更为严重。方法1，方法2与方法3均显著膨胀。

Table 7. Observed relative frequency of rejection of the hypothesis of the absence of outliers on 1000 series generated according to model 3 and non-Poisson innovations

表7. 基于模型3和非泊松创新生成的1000组序列中，拒绝无离群点假设的观测相对频率

创新	方法
	1	2	3
Geometric	0.314	0.091	0.472
Binomial	0.563	0.315	0.682

在非泊松创新的场景下，我们对每个模拟序列的时间点$q=50$ 插入大小为${\omega }_0=3.5$ 的加性离群点，并重复进行检测。检测结果如表8所示。分析可观察到，对于由几何分布(Geometric)生成的创新项，方法1和方法2的检测功效与泊松分布创新项的情况相近，方法3略差一些。此外，方法2在检测到离群点时，错误定位的频率显著低于方法1与方法3。对于由二项分布(Binomial)生成的创新项，三种方法的检测结果都存在显著偏差。

综合来看，对于泊松数据，BELR的检验方法更加有效，而在非泊松创新的情况下，BELR的方法效果依旧比ELR与贝叶斯方法显著，此外，也充分证明了在INAR(1)模型中对离群点的检测与估计使用经验似然方法是有效的。

Table 8. Observed relative frequency of rejection of the hypothesis of the absence of outliers on 1000 series with an additive outlier at $q=50$ and size ${\omega }_0=3.5$ generated according to model 3 and non-Poisson innovations

表8. 基于模型3和非泊松创新生成的1000组序列中，假设无离群点的假设被拒绝的观测相对频率，其中$q=50$ ，${\omega }_0=3.5$

创新	时间位置	方法
		1	2	3
Geometric	准确	0.674	0.836	0.534
	错误	0.051	0.031	0.055
Binomial	准确	0.316	0.427	0.312
	错误	0.047	0.031	0.062

5. 实例分析

曲杆菌病是由弯曲杆菌属病原体引发的急性细菌感染。该病原体主要攻击人体消化系统，是全球细菌性胃肠炎最主要的致病菌之一。临床研究显示，感染后潜伏期存在个体差异，多数患者在暴露后2~5日出现典型症状，但存在1~10日的波动区间，这与宿主免疫状态及感染剂量密切相关。

值得注意的是，该菌株具有血行播散潜能，可能导致菌血症等严重并发症。流行病学数据显示，婴幼儿、老年群体及免疫功能受损患者的致死风险显著升高，病死率可达普通人群的3~5倍。其病例数据常表现为整数值时间序列，适合使用一阶整值自回归模型(INAR(1))进行建模。本节针对加拿大魁北克省北部1990~2000年弯曲杆菌病例数据[27] (每28天记录一次，取120个观测值)，其样本均值为$\overline{X}=10.1833$ ，样本方差为$s^2=39.6933$ ，绘制时间序列图，如图1所示，红色三角形标记检测到的离群点，蓝色虚线标识真实异常。

Figure 1. Campylobacter data

图1. 弯曲杆菌数据

由图1可以看出在时间$t=70$ 前病例数围绕某一基线水平(约8~12例/周期)波动，表现出弱平稳性，标准差约为2.5例，表明感染传播处于相对可控状态。在时间$t=70$ 后病例数出现显著爆发式增长，峰值达25例($t=73$ )，较基线水平上升约150%，且后续持续高位震荡(15~20例/周期)，表明可能存在离群点，此现象可能由于外部系统性因素导致，如魁北克北部冬季寒冷，可能导致饮用水管道冻裂，引发水源性污染事件，短期内病例数激增，或宿主免疫波动，导致病例异常增加。离群点的存在会导致模型中参数估计存在偏差，以及预测性能下降。

接下来运用INAR(1)与该数据进行了拟合，并使用经验似然方法检测离群点，展示了检测结果与病例数对比，可以观察到在$t=70$ 处(对应真实病例激增位置)，经验似然方法成功检测到离群点，除真实离群点外，由于强离群点可能存在掩盖效应，仅在$t=85$ 被误判为离群点，以上充分说明了经验似然方法对检测离群点的有效性，通过非参数框架、动态约束整合和稳健权重分配，实现了离群点检测的高精度与强解释性。

6. 结论与展望

本研究针对一阶整值自回归模型(INAR(1))的离群点检测与参数估计问题，提出了一套基于改进经验似然法的创新解决方案。针对传统经验似然方法在离散时间序列分析中存在的凸包问题，创新性地引入平衡经验似然(BEL)方法，通过构造对称约束条件有效解决了有限样本下矩方程无解的数值稳定性难题。

通过大量的模拟实验结果表明，经验似然方法可以有效地检测出INAR(1)模型中的离群点，并且相对于传统的检测方法，经验似然方法对离群点的检测效果更加显著。

然而，当前研究仍存在一些局限性，未来的研究可以从以下几个方向进行拓展和深化：

• 本文主要关注一阶整值时间序列模型，但在实际应用中，许多时间序列数据可能具有更高阶的自相关性或更复杂的动态结构。因此，未来的研究可以考虑将经验似然方法应用于更高阶的整值时间序列模型，如二阶或更高阶的自回归模型(如INAR(p)模型)，以提高离群值检测的适用性和准确性。

• 随着数据分析的复杂性增加，多变量时间序列数据在金融、气象、生物医学等领域中越来越常见。将经验似然方法扩展到多变量时间序列的离群值检测中，将是一个具有挑战性和应用价值的研究方向。这需要考虑变量之间的相关性以及离群值在多维空间中的表现形式。

• 尽管经验似然方法在离群值检测中表现良好，但其计算复杂度可能在处理大规模数据时成为瓶颈。未来的研究可以探索更高效的算法实现，例如通过并行计算、近似算法或利用机器学习技术来加速经验似然的计算过程，从而提高其在实际应用中的可行性。

• 离群值检测是一个多学科交叉的研究领域，经验似然方法可以与其他先进的统计或机器学习方法相结合，以进一步提高检测性能。例如，将经验似然与深度学习模型(如自编码器、神经网络)相结合，可能会在复杂时间序列数据中发现更隐蔽的离群值模式。

基金项目

辽宁科技大学博士启动基金(6003000310)。

NOTES

^*通讯作者。