INAR模型的多变点探测

Multiple Change Point Detection for INAR Model

摘要: 本文考虑了分段平稳的整数值自回归模型的多变点问题。借助最小描述长度函数得到似然比扫描方法，并将似然比扫描方法运用到分段平稳的整数值自回归模型中。此外，还研究了当模型系数之和趋于1时的情况，此时针对似然比扫描方法中的最小描述长度函数进行了调整，以提高变点探测的准确性。然后通过大量的数值模拟，验证了似然比扫描方法在不同的模型参数设置下的有效性，最后并将其运用于一组精神分裂症患者在知觉速度测试中的日常观察得分数据的实证分析之中。

Abstract: This paper considers the multiple change-points problem in piecewise stationary integer-valued autoregressive model. Using the Minimum Description Length function, the likelihood ratio scanning method is obtained and applied to the piecewise stationary integer-value autoregressive model. In addition, when the sum of the model coefficients tends to 1, the Minimum Description Length function in the likelihood ratio scanning method is adjusted to improve the accuracy of the change point detection. Then, through many numerical simulations, the effectiveness of the likelihood ratio scanning method in different model parameter settings was verified. Finally, it was applied to the empirical analysis of the daily observations of the score achieved by a schizophrenic patient on a test of perceptual speed.

1. 引言

在时间序列分析过程中，所关注的某些量可能在某个时刻发生显著变化，这个变化时刻被称作变点。如果在发现变点后仍旧基于之前的模型进行统计分析，就可能得到错误的分析结果，为了降低错误率，提升效益，因此有必要进行变点统计分析。最早的变点研究可以追溯到20世纪50年代，Page [1] 针对质量控制问题发表了一篇关于变点探测的文章，自此变点问题在统计领域得到了广泛的讨论。目前，经过研究人员对变点方法不断地发展与完善，变点的处理方法已经有很多，部分方法也很成熟了。多变点探测常用的方法是通过优化特定的目标函数来寻找变点，如最小二乘法 [2] 、Bayes方法 [3] 等。除此之外还有非参数方法，如U-统计量法 [4] 、Wilcoxon秩和法 [5] 、Kolmogorow检验法 [6] 、累积和(Cumulative Sum, CUSUM)检验法 [7] [8] 。需要注意的是，变点组合的数量会随着样本量的增长而呈指数增长，这使得在优化过程中的计算难度变大。而为了避免优化难题，二分法是一个比较受欢迎的解决办法。二分法最开始是由Vostrikova [9] 和Bai [10] 提出并研究的，二分法更多地应用于均值多变点 [11] 和方差多变点 [12] 的探测中。Fryzlewicz [13] 将二分法扩展为Wild二分法，该拓展方法非常重要，被命名为野生二进制分割(Wild Binary Segmentation, WBS)。除此之外，还有Davis等人 [14] 提出的遗传算法，Killick等人 [15] 提出的截断精确线性时间(Pruned Exact Linear Time, PELT)方法以及Yau和Zhao [16] 提出的似然比扫描方法(Likelihood Ratio Scan Method, LRSM)。与WBS，PELT等方法相比，LRSM不仅降低了计算复杂度，而且提高了变点估计的精度。

近年来兴起了对整数值时间序列的变点问题的研究。随着此类时间序列的几个相应模型的提出，整数值时间序列在各个领域中有着越发广泛的应用，例如在金融(在一个时间周期内事件的发生的概率)、气候学、医学等。关于此类时间序列相关模型的构建，Mckenzie [17] 和Al-osh和Alzaid [18] 提出了一个带有泊松边际分布的INAR(1)过程，即众所周知的PINAR(1)过程，由于其简单的形式非常受欢迎。Alzaid和Alzaid [19] 以及Du和Li [20] 将INAR(1)过程进一步扩展成更一般的INAR(p)过程。尽管INAR模型在实践中得到的应用非常广泛，但在拟合非线性现象时往往会出现不足。例如针对疫情出现拐点的这样的数据，平稳的INAR模型可能无法提供良好的拟合，需要针对变点对其进行分段拟合。因此，对于此类整数值时间序列，研究数据中的变点位置和数量，具有重要的实际意义。目前，据我们所知，关于此类整数值时间序列的变点探测的研究尚不多见，更多的是基于此类序列的变点检验的研究，例如文献 [21] [22] [23] 。

在现有的变点探测方法中，LRSM简单且易于实施，并且当数据量为 $n$ 时，计算复杂度从 $n^2$ 降低到了 $\log \left(n\right)$ ，Yau和Zhao [16] 在模拟研究中将LRSM与其他方法进行了对比，总体而言LRSM的适用性更强并且比其他方法更准确些。因此，本文选择LRSM探测由分段平稳的INAR模型所产生的整数值时间序列数据中的多变点。首先基于分段平稳INAR(1)过程给出变点问题的基本假设和似然比扫描方法的具体实施步骤。其次通过大量的数值模拟并结合实证分析证明LRSM的有效性。

本文的其余部分组织如下。第2节详细介绍了INAR模型以及LRSM的三个步骤。第3节通过大量的数值模拟，说明将似然比扫描方法应用于整数值自回归过程中的有效性。第4节将LRSM方法运用到一个真实案例中。最后，在第5节对全篇内容做出总结，概述本论文所做的工作。

2. INAR模型以及似然比扫描方法

2.1. INAR模型

在许多领域中，例如，经济学、医学、精算统计等，诸多有趣的变量都是非负整数值的，具体的有每次绿灯时路口通过的车辆数量、某城市每日发生抢劫案件的次数、每月某航班的起降次数、某个病人在服药前后的每日发病次数。因此近年来，关于非负整数值时间序列的建模受到越来越多的关注。此类数据相关模型的构建，最常用的是手段就是利用稀疏算子进行建模。Steutal和Van Harn [24] 最早提出了二项稀疏算子“◦”，其定义如下：

$\alpha \circ Y\triangleq {\displaystyle \underset{j=1}{\overset{Y}{\sum }}{z}_j},$

其中 $\alpha \in \left[0,1\right)$ ，Y是非负的整数值随机变量， $\left\{z_j\right\}$ 为独立同分布的以 $\alpha$ 为参数的伯努利随机变量序列，并独立于Y，其分布律为 $P\left\{z_j=1\right\}=1-P\left\{z_j=0\right\}=\alpha$ 。Du和Li [20] 利用算子提出了INAR(p)模型，这也是本文考虑的模型，其具体定义如下：

$Y_t={\alpha }_1\circ Y_{t-1}+\cdots +{\alpha }_p\circ Y_{t-p}+{\epsilon }_t,t\ge p,$

其中 ${\epsilon }_t$ 是独立同分布(i.i.d)的非负整数值随机变量序列，且服从均值为 $\lambda$ 的泊松分布。

关于该模型的的基本性质有：

1) $Y_t$ 的条件均值： $E\left[Y_t|Y_{t-1},\cdots ,Y_{t-p}\right]={\displaystyle \underset{j=1}{\overset{p}{\sum }}{\alpha }_j{Y}_{t-j}}+\lambda$ ；

2) $Y_t$ 的条件方差： $Var\left[Y_t|Y_{t-1},\cdots ,Y_{t-p}\right]={\displaystyle \underset{j=1}{\overset{p}{\sum }}{\alpha }_j\left(1-{\alpha }_j\right)Y_{t-j}}+\lambda$ ；

3) $Y_t$ 的均值： $E\left(Y_t\right)=\lambda /1-{\displaystyle \underset{j=1}{\overset{p}{\sum }}{\alpha }_j}$ ；

4) 当 ${\alpha }_j\in \left[0,1\right)$ ， $j=1,2,\cdots ,p$ ， ${\displaystyle \underset{j=1}{\overset{p}{\sum }}{\alpha }_j}<1$ 时，则模型是平稳遍历的。

2.2. 似然比扫描方法

LRSM相比较其它的变点方法，可以快速有效的探测出时间序列中的变点。此方法主要有三个步骤：1) 扫描由分段平稳的INAR过程所产生的观测序列以获得初始的变点估计；2) 利用MDL函数进行INAR模型的选择过程，得到一致的变点估计；3) 对第二步中得到的变点集合 ${\hat{J}}^{\left(2\right)}$ 中的变点进行最终估计，以获得变点近似阶数达到最优时的变点估计。下面对基于LRSM探测分段平稳INAR过程中变点的三个步骤进行详细介绍。

为更好研究分段平稳INAR过程中的多变点问题，下面先给出一些基本设定。假如有观测序列 ${\left\{X_t\right\}}_{t=1,\cdots ,n}$ 可以被分割成 $m+1$ 段平稳的INAR过程。设 $j=1,\cdots ,m$ ，j表示变点个数，用 ${\tau }_j$ 表示序列从第j段突然变化到第 $j+1$ 段变点的位置(第j个变点的位置)。令 ${\tau }_0\triangleq 0$ ， ${\tau }_{m+1}\triangleq n$ ，则有INAR序列的第j段可以表示为：

$\begin{array}{ccc}{Y}_{t,j}=X_t,& {\tau }_{j-1}<t\le {\tau }_j,& j=1,2,\cdots ,m+1\end{array}$

其中 $\left\{Y_{t,j}\right\}$ 是一个平稳的INAR过程，满足：

$Y_{t,j}={\alpha }_{j,1}\circ Y_{t-1,j}+\cdots +{\alpha }_{j,p_j}\circ Y_{t-p_j,j}+{\epsilon }_{t,j},$

其中算子“◦”为二项稀疏算子， ${\epsilon }_{t,j}$ 为独立同分布的非负整数值随机序列。第j段对应的INAR( $p_j$ )模型的参数向量为 ${\theta }_j=\left({\alpha }_{j,1},\cdots ,{\alpha }_{j,p_j},{\lambda }_j\right)$ 。

假设1 假设 ${\displaystyle \underset{i=1}{\overset{p_j}{\sum }}{\alpha }_{j,i}}<1$ ，以保证每一段INAR过程 ${\left\{Y_{t,j}\right\}}_{j=1,2,\cdots ,m+1}$ 都是平稳的。

假设2 假设观测序列 ${\left\{X_t\right\}}_{t=1,\cdots ,n}$ 被分割每一段INAR过程 ${\left\{Y_{t,j}\right\}}_{j=1,2,\cdots ,m+1}$ 的阶数都是有限的，这里把每一段INAR模型的最大阶数用 $p_{\max }$ 来表示。

假设3 这里用 $J=\left({\tau }_1,\cdots ,{\tau }_m\right)$ ，表示一系列的变点。定义各个变点的相对位置 $\left\{{\varpi }_j\right\}$ ，满足 $0={\varpi }_0<{\varpi }_1<\cdots <{\varpi }_m<{\varpi }_{m+1}=1$ 且 ${\tau }_j=\left[{\varpi }_{j}n\right]$ ，这里 $\left[{\varpi }_{j}n\right]$ 表示 ${\varpi }_{j}n$ 的整数部分。假设存在 ${\epsilon }_{\varpi }={\epsilon }_{\varpi }\left(n\right)>0$ ，使得 ${\min }_{j=0,\cdots ,m-1}\left({\varpi }_{j+1}-{\varpi }_j\right)>{\epsilon }_{\varpi }$ 。同时为了保证模型的可识别性，进一步假设 $n{\epsilon }_{\varpi }>p_{\max }$ 。

接下来详细介绍运用LRSM探测分段平稳INAR过程中多变点的三个步骤。

第一步：扫描来自于分段平稳INAR过程的观测值 ${\left\{X_t\right\}}_{t=1,\cdots ,n}$ 序列，得到最初估计的变点。

定义t时刻的扫描窗口为：

$W_t\left(h\right)=\left\{t-h+1,\cdots ,t+h\right\},$

以及对应的观测值为：

$X_{W_t\left(h\right)}=\left(X_{t-h+1},\cdots ,X_{t+h}\right),$

其中 $t=h,\cdots ,n-h$ ，h为扫描窗口半径。扫描窗口半径h的选择与样本量n以及变点之间的距离有关，关于h的选择，详细可见第3节。

为了获得扫描窗口中最初估计的变点，可以选择似然比统计量。对于一个来自于INAR(p)模型的样本 $D=\left(d_1,\cdots ,d_n\right)$ ，则有INAR模型的条件似然函数如下：

$L\left(\theta \right)=\log {\displaystyle \underset{t=p+1}{\overset{n}{\sum }}\left[P\left(d_t|d_{t-1},\cdots ,d_{t-p}\right)\right]}$ (1)

其中，

$\begin{array}{c}P\left(d_t|d_{t-1},\cdots ,d_{t-p}\right)={\displaystyle \underset{i_1=0}{\overset{\min \left(d_{t-1},d_t\right)}{\sum }}\left(\begin{array}{c}{d}_{t-1}\\ i_1\end{array}\right){\alpha }_1^{i_1}{\left(1-{\alpha }_1\right)}^{d_{t-1}-i_1}}{\displaystyle \underset{i_2=0}{\overset{\min \left(d_{t-2},d_t-i_1\right)}{\sum }}\left(\begin{array}{c}{d}_{t-2}\\ i_2\end{array}\right){\alpha }_2^{i_2}{\left(1-{\alpha }_2\right)}^{d_{t-2}-i_2}}\cdots \\ {\displaystyle \underset{i_p=0}{\overset{\min \left[d_{t-p},d_t-\left(i_1+\cdots +i_{p-1}\right)\right]}{\sum }}\left(\begin{array}{c}{d}_{t-p}\\ i_p\end{array}\right){\alpha }_p^{i_p}}{\left(1-{\alpha }_p\right)}^{d_{t-p}-i_p}\frac{{\text{e}}^{-\lambda }{\lambda }^{d_t-\left(i_1+\cdots +i_p\right)}}{\left[d_t-\left(i_1+\cdots +i_p\right)\right]!},\end{array}$

是在 $d_t$ 之前给定观测值的转移概率。然后扫描窗口 $W_t\left(h\right)$ 的扫描统计量可以定义为：

$S_h\left(t\right)=\frac{1}{h}{L}_{1h}\left(t,{\hat{\theta }}_1\right)+\frac{1}{h}{L}_{2h}\left(t,{\hat{\theta }}_2\right)-\frac{1}{h}{L}_{\cdot h}\left(t,\hat{\theta }\right)$

其中 $L_{1h}\left(t,{\hat{\theta }}_1\right)$ ， $L_{2h}\left(t,{\hat{\theta }}_2\right)$ 和 $L_{\cdot h}\left(t,\hat{\theta }\right)$ 的定义与等式 中的 $L\left(\theta \right)$ 类似，但是使用的观测值是 ${\left\{X_s\right\}}_{t-h+1,\cdots ,t}$ ， ${\left\{X_s\right\}}_{t+1,\cdots ,t+h}$ 以及 ${\left\{X_s\right\}}_{W_t\left(h\right)}$ ， ${\hat{\theta }}_1$ ， ${\hat{\theta }}_2$ 以及 $\hat{\theta }$ 是基于序列 ${\left\{X_s\right\}}_{t-h+1,\cdots ,t}$ ， ${\left\{X_s\right\}}_{t+1,\cdots ,t+h}$ 以及 ${\left\{X_s\right\}}_{W_t\left(h\right)}$ 的参数 $\theta$ 的估计量。接下来，利用 $S_h\left(t\right)$ 扫描观测到的序列，可以得到一系列的似然比扫描统计量值的序列 $\left(S_h\left(h\right),S_h\left(h+1\right),\cdots ,S_h\left(n-h\right)\right)$ ，假如t是变点，则 $S_h\left(t\right)$ 会趋于很大。如果选择的h满足 $2h<n{\epsilon }_{\varpi }$ 以及 $h>p_{\max }$ ，则每个扫面窗口中最多存在一个真实变点。因此，通过对局部变点位置的估计，可以得到潜在变点集合并定义为：

${\hat{J}}^{\left(1\right)}=\left\{m\in \left\{h,h+1,\cdots ,n-h\right\}:S_h\left(m\right)=\underset{t\in \left(m-h,m+h\right]}{\max }{S}_h\left(t\right)\right\},$

其中当 $t<h$ 且 $t>n-h$ 时， $S_h\left(t\right)\triangleq 0$ 。如果以点m为中心的窗口 $\left[m-h+1,m+h\right]$ 中的最大值为 $S_h\left(m\right)$ ，则m为一个局部变点。则有一系列的初始变点 ${\hat{J}}^{\left(1\right)}=\left({\hat{\tau }}_1^{\left(1\right)},{\hat{\tau }}_2^{\left(1\right)},\cdots ,{\hat{\tau }}_{{\hat{m}}^{\left(1\right)}}^{\left(1\right)}\right)$ ，其中 ${\hat{m}}^{\left(1\right)}=\left|{\hat{J}}^{\left(1\right)}\right|$ 。

定理1 记真实变点集合为 $J_0=\left({\tau }_1^0,{\tau }_2^0,\cdots ,{\tau }_{m_0}^0\right)$ ，假设 $2h<n{\epsilon }_{\varpi }$ 且 ${\epsilon }_{\varpi }>c{n}^{-q}$ ，则存在 $q\in \left[0,1\right)$ ， $c>0$ 满足：

$P\left(\underset{{\tau }_i^0\in J_0}{\max }\underset{{\hat{\tau }}_j^{\left(1\right)}\in {\hat{J}}^{\left(1\right)}}{\min }\left|{\tau }_i^0-{\hat{\tau }}_j^{\left(1\right)}\right|<h\right)\to 1.$

定理1表明全部真实变点在初始估计变点集合 ${\hat{J}}^{\left(1\right)}$ 的h-领域内。并且当变点之间的最小距离为 $n{\epsilon }_{\varpi }=O\left(n^{1-q}\right)$ 时，估计的变点个数 ${\hat{m}}^{\left(1\right)}$ 以 $O\left(n^q\right)$ 的速度随着样本量的增加而增加。但此时并不能保证 ${\hat{m}}^{\left(1\right)}$ 等于真实的变点个数 $m_0$ 。也就是说，在第一步中变点的数量可能会被高估。因此接下来通过第二步，即利用合适的信息准则从 ${\hat{J}}^{\left(1\right)}$ 中再次筛选，以获得一致的变点估计。

第二步：优化目标函数，选择INAR模型，得到一致的变点估计。

为了从 ${\hat{J}}^{\left(1\right)}$ 中找出准确的变点，可以使用最小描述长度准则 [25] (简记MDL)进行INAR模型的选择过程。定义MDL如下：

$MDL\left(m,J,P\right)=\log \left(m\right)+\left(m+1\right)\log \left(n\right)+{\displaystyle \underset{j=1}{\overset{m+1}{\sum }}\log \left(p_j\right)}+{\displaystyle \underset{j=1}{\overset{m+1}{\sum }}\frac{p_j+1}{2}\log }\left(n_j\right)-{\displaystyle \underset{j=1}{\overset{m+1}{\sum }}{L}_j}\left({\hat{\theta }}_j\right),$

其中， $J=\left({\tau }_1,\cdots ,{\tau }_m\right)$ 的取值是第一步中得到的初始估计变点集合 ${\hat{J}}^{\left(1\right)}$ ， $\left(n_1,\cdots ,n_{m+1}\right)$ 和 $P=\left(p_1,\cdots ,p_{m+1}\right)$ 分别是所有片段的长度和相应的INAR片段的阶数， $L_j\left({\hat{\theta }}_j\right)$ 是第j个INAR片段的似然统计量。通过MDL准则的优化后，可以得到优化后的变点集合为：

$\left({\hat{m}}^{\left(2\right)},{\hat{J}}^{\left(2\right)},{\hat{P}}^{\left(2\right)}\right)=\underset{\begin{array}{c}m=\left|J\right|,J\subseteq {\hat{J}}^{\left(1\right)}\\ P\in {\left\{1,\cdots ,p_{\max }\right\}}^m\end{array}}{\arg \min }MDL\left(m,J,P\right),$

其中 ${\hat{J}}^{\left(2\right)}=\left({\hat{\tau }}_1^{\left(2\right)},{\hat{\tau }}_2^{\left(2\right)},\cdots ,{\hat{\tau }}_{{\hat{m}}^{\left(2\right)}}^{\left(2\right)}\right)$ ，而 ${\hat{m}}^{\left(2\right)}=\left|{\hat{J}}^{\left(2\right)}\right|$ 。在第二步中得到的变点估计具有如下性质：

定理2 在定理1的成立的条件下，令 $q=0$ ，则有 ${\hat{m}}^{\left(2\right)}\stackrel{p}{\to }{m}_0$ 。并且，当 ${\hat{m}}^{\left(2\right)}=m_0$ ，则有：

$P\left(\underset{j=1,\cdots ,m_0}{\max }\left|{\hat{\tau }}_j^{\left(2\right)}-{\tau }_j^0\right|\right)\to 1,$

$\underset{j=1,\cdots ,m_0+1}{\max }\left|{\hat{p}}_j^{\left(2\right)}-p_j^0\right|\stackrel{p}{\to }0.$

根据定理2可知变点 ${\tau }_j^0\in \left[{\hat{\tau }}_j^{\left(2\right)}-h+1,{\hat{\tau }}_j^{\left(2\right)}+h\right]$ 。由于 $h\ge d\log \left(n\right)\to \infty$ ，因此有 ${\max }_{j=1,\cdots ,{\hat{m}}^{\left(2\right)}}\left|{\hat{\tau }}_j^{\left(2\right)}\right|=O_p\left(h\right)$ 与 $O_p\left(1\right)$ 相比的话得到的变点近似阶数并不是最优的，这里d的取值可参考文献 [16] 。通过运用文献 [26] 中的思想来获得变点位置的最终估计使其近似的阶数达到最优。

第三步：获得使得变点近似阶数达到最优的最终的变点估计。

定义第j个估计变点 ${\hat{\tau }}_j^{\left(2\right)}$ 的拓展的局部窗口和相应的观测值分别为：

$E_j\left(h\right)=\left\{{\hat{\tau }}_j^{\left(2\right)}-2h+1,\cdots ,{\hat{\tau }}_j^{\left(2\right)}+2h\right\},$

$X_{E_j\left(h\right)}=\left(X_{{\hat{\tau }}_j^{\left(2\right)}-2h+1},\cdots ,X_{{\hat{\tau }}_j^{\left(2\right)}+2h}\right),$

令

$L_j\left(\tau ,{\theta }_1,{\theta }_2\right)={\displaystyle \underset{t={\hat{\tau }}_j^{\left(2\right)}-2h+1}{\overset{\tau }{\sum }}{l}_t\left({\theta }_1\right)}+{\displaystyle \underset{t=\tau +1}{\overset{{\hat{\tau }}_j^{\left(2\right)}+2h}{\sum }}{l}_t}\left({\theta }_2\right),$

对于 $j=1,\cdots ,{\hat{m}}^{\left(2\right)}$ ，定义变点最后的估计量为：

${\hat{\tau }}_j^{\left(3\right)}=\underset{\tau \in \left({\hat{\tau }}_j^{\left(2\right)}-h,{\hat{\tau }}_j^{\left(2\right)}+h\right]}{\arg \max }{L}_j\left(\tau ,{\hat{\theta }}_j,{\hat{\theta }}_{j+1}\right),$

其中 ${\hat{\theta }}_j={\hat{\theta }}_j\left(\tau \right)=\arg {\max }_{{\theta }_1}{\displaystyle {\sum }_{t={\hat{\tau }}_j^{\left(2\right)}-2h+1}^{\tau }{l}_t}\left({\theta }_1\right)$ ，同理可以得到 ${\hat{\theta }}_{j+1}$ 的定义。同样地，第j段的INAR模型的阶数取 ${\hat{p}}_j^{\left(2\right)}$ 。

定理3 在定理2成立的条件下，若 $3h<n{\epsilon }_{\varpi }$ ，则有：

${\hat{\tau }}_j^{\left(3\right)}-{\tau }_j^0=O_p\left(1\right),$

且

${\hat{\tau }}_j^{\left(3\right)}-{\tau }_j^0\stackrel{d}{\to }\underset{{\tau }_j}{\arg \max }W\left({\tau }_j,{\theta }_j^0,{\theta }_{j+1}^0\right),$

其中 $W\left({\tau }_j\right):=W\left({\tau }_j,{\theta }_j^0,{\theta }_{j+1}^0\right)$ 是一个双边随机游走如下所示：

$W\left({\tau }_j,{\theta }_j^0,{\theta }_{j+1}^0\right)=\{\begin{array}{ll}{\displaystyle \underset{t={\tau }_j^0+1}{\overset{{\tau }_j^0+{\tau }_j}{\sum }}\left\{l_t\left({\theta }_j^0\right)-l_t\left({\theta }_{j+1}^0\right)\right\},}\hfill & {\tau }_j>0,\hfill \\ 0,\hfill & {\tau }_j=0,\hfill \\ {\displaystyle \underset{t={\tau }_j^0+{\tau }_j}{\overset{{\tau }_j^0-1}{\sum }}\left\{l_t\left({\theta }_{j+1}^0\right)-l_t\left({\theta }_j^0\right)\right\},}\hfill & {\tau }_j<0.\hfill \end{array}$

定理1~3证明过程可参考文献 [16] 。

3. 模拟研究

在本节中，通过大量的数值模拟来评估LRSM在由分段平稳的INAR过程所产生的有限样本下变点识别性能，所有模拟实验均通过R语言程序实现。当样本量小于800时，可以考虑 $h=\max \left\{25,{\left(\log n\right)}^2\right\}$ ，当样本量n大于800时，可以考虑 $h=\max \left\{50,2{\left(\log n\right)}^2\right\}$ 。在本文中，样本量大于800时，取 $h=2{\left(\log n\right)}^2$ ；如果样本量小于800时，取25。所有模拟结果都经过100次循环得到。关于变点数量以及位置 ${\tau }_j=\left[{\varpi }_{j}n\right]$ 的设定，分别设置成没有变点，有一个变点，有两个变点三种情况。在只有一个变点的情况下 ${\varpi }_1$ 取0.5；在有两个变点的情况下， $\left\{{\varpi }_1,{\varpi }_2\right\}=\left\{0.34,0.68\right\}$ 或者 $\left\{{\varpi }_1,{\varpi }_2\right\}=\left\{0.39,0.68\right\}$ 。通过计算频率 $f_0={\displaystyle {\sum }_{i=1}^{100}{I}_{{\hat{m}}_i^{\left(2\right)}=m_0}}$ 来检验LRSM是否可以准确探测出变点数量，其中 $m_0$ 表示真正的变点数量。同时为了评估LRSM估计精度，给出如下定义：如果一个方法若能准确的探测出变点数量，并且当 $n=1024$ 时，每对真实变点和估计变点之间的距离在50以内；当 $n=120$ 时，每对真实变点和估计变点之间的距离在5以内，认为探测出来的变点是有效的。在 ${\hat{m}}_i^{\left(2\right)}=m_0$ 时，还通过计算均值，和均方误差来评估估计精度。

LRSM的三个步骤中都需要求对相应的模型参数的极大似然估计值，但在实际操作中采用的是矩估计值的绝对值，具体原因如下：1) 众所周知，极大似然估计法需要通过数值优化来求解，因此求解过程耗费时间长；矩估计因其有显著的表达式，求解过程简单，耗费时间短。2) 通过模拟研究发现，采用矩估计值几乎不影响LRSM的性能。3) 在模拟研究中，由于生成的随机数，具有一定的随机性，因此有时参数的矩估计值需要调整，而本文采用的调整手段是直接取矩估计值的绝对值。

3.1. INAR(1)的模拟研究

在本小节中，通过数值模拟来评估LRSM在由分段平稳的INAR(1)过程所产生的有限样本下变点识别性能，所考虑的模型有：模型A-D。模型A-D的样本路径图见图1，模拟结果见表1。

Table 1. Simulation results from Model A to Model D

表1. 从模型A到模型D的模拟结果

¹正确识别变点数量的次数；²有效识别变点的次数；³真实变点的位置；⁴有效识别到的变点的平均值；⁵有效识别到的变点的均方误差。

模型A：没有变点的平稳INAR(1)过程

模型A用来评估LRSM在没有变点时的性能，即当观测值中无变点时，该方法能否准确识别。设样本量 $n=1024$ ，并采用多种参数组合的INAR(1)过程： $X_t=\alpha \circ X_{t-1}+{\epsilon }_t$ 生成观测值，其中 ${\epsilon }_t$ 服从均值为 $\lambda$ 的泊松分布， $\alpha$ 和 $\lambda$ 的取值见表1。从表1中可以看出，在 $\alpha =0.3$ ， $\lambda =3$ 以及 $\alpha =0.9$ ， $\lambda =1$ 时正确识别变点数量的次数只有98次，而其它参数组合的情况下，正确识别变点的次数都是100次。由于随机因数导致LRSM偶尔会有一两次错误的识别到变点。总的来看在没有变点的情况下，LRSM是非常准确有效的。

模型B：有一个变点的分段平稳INAR(1)过程

$X_t=\{\begin{array}{ll}0.5\circ X_{t-1}+{\epsilon }_{t,1}\left(1\right)\hfill & 1\le t\le 512,\hfill \\ 0.5\circ X_{t-1}+{\epsilon }_{t,2}\left(4\right)\hfill & 513\le t\le 1024.\hfill \end{array}$

其中 ${\epsilon }_{t,j}\left({\lambda }_j\right)$ 表示分段平稳的INAR模型的新息过程 ${\epsilon }_{t,j}$ 服从均值为 ${\lambda }_j$ 的泊松分布。下面模型C-D的 ${\epsilon }_{t,j}\left({\lambda }_j\right)$ 的含义与这里的 ${\epsilon }_{t,j}$ 的一致。

模型B在512处有一个变点。从表1可以看出，LRSM正确识别变点数量的频率高达99次，且全都满足了估计变点与真实变点之间的距离均小于50；而估计变点的均值约为512，且均方误差也是合理的。所以针对模型B这种情况，LRSM的估计效果也很理想。

模型A ( $\alpha =0.4$ ) 模型B模型C 模型D

Figure 1. Sample path from Model A to Model D

图1. 模型A到模型D的样本路径图

模型C：有两个变点的分段平稳INAR(1)过程

$X_t=\{\begin{array}{ll}0.3\circ X_{t-1}+{\epsilon }_{t,1}\left(1\right)\hfill & 1\le t\le 400,\hfill \\ 0.6\circ X_{t-1}+{\epsilon }_{t,2}\left(2\right)\hfill & 401\le t\le 612,\hfill \\ 0.8\circ X_{t-1}+{\epsilon }_{t,3}\left(3\right)\hfill & 613\le t\le 1024.\hfill \end{array}$

模型C中有两个变点，分别在400和612。正确识别变点个数的次数有98次，而98次的结果都是有效的。探测出来的两个变点的均值皆接近真实值，均方误差也不是太大。但是随着变点个数的增加，变点识别难度也会提升，LRSM变点识别准确率有所降低，总的来看估计效果还是比较理想。

模型D：在小样本情况下，有一个变点的分段平稳INAR(1)过程

$X_t=\{\begin{array}{ll}0.455\circ X_{t-1}+{\epsilon }_{t,1}\left(38.14\right)\hfill & 1\le t\le 69,\hfill \\ 0.64\circ X_{t-1}+{\epsilon }_{t,2}\left(13.62\right)\hfill & 70\le t\le 120.\hfill \end{array}$

模型D用于评估当样本量较小时，LRSM的性能。所有参数以及变点位置均根据第4节中的真实案例分析设定。从表1中可以看出，无论是在准确探测变点数量上还是有效探测变点位置上，LRSM都表现良好。

综上所述，LRSM能较为高效准确的探测出INAR(1)模型中的变点个数和位置。

3.2. INAR(2)的模拟研究

在本小节中，通过数值模拟来评估LRSM在由高阶的分段平稳INAR过程所产生的有限样本下变点识别性能，所考虑的模型有：模型E~G。模型E~G的样本路径图见图2，模拟结果见表2。

Table 2. Simulation results from Model E to Model G

表2. 从模型E到模型G的模拟结果

¹正确识别变点数量的次数；²有效识别变点的次数；³真实变点的位置；⁴有效识别到的变点的平均值；⁵有效识别到的变点的均方误差。

模型E：无变点的平稳INAR(2)过程

$\begin{array}{cc}{X}_t=0.5\circ X_{t-1}+0.2\circ X_{t-2}+{\epsilon }_{t,1}\left(1\right)& 1\le t\le 1024.\end{array}$

模型F：具有一个变点的分段平稳INAR(2)过程

$X_t=\{\begin{array}{ll}0.5\circ X_{t-1}+0.1\circ X_{t-2}+{\epsilon }_{t,1}\left(1\right)\hfill & 1\le t\le 512,\hfill \\ 0.3\circ X_{t-1}+0.1\circ X_{t-2}+{\epsilon }_{t,2}\left(4\right)\hfill & 513\le t\le 1024.\hfill \end{array}$

模型G：具有两个变点的分段平稳INAR(2)过程

$X_t=\{\begin{array}{ll}0.25\circ X_{t-1}+0.2\circ X_{t-2}+{\epsilon }_{t,1}\left(1\right)\hfill & 1\le t\le 350,\hfill \\ 0.25\circ X_{t-1}+0.2\circ X_{t-2}+{\epsilon }_{t,2}\left(4\right)\hfill & 351\le t\le 700,\hfill \\ 0.25\circ X_{t-1}+0.2\circ X_{t-2}+{\epsilon }_{t,3}\left(1\right)\hfill & 701\le t\le 1024.\hfill \end{array}$

其中 ${\epsilon }_{t,j}\left({\lambda }_j\right)$ 表示分段平稳的INAR模型的新息过程 ${\epsilon }_{t,j}$ 服从均值为 ${\lambda }_j$ 的泊松分布。

由于在模型INAR(2)中LRSM变点识别性能与在模型INAR(1)中表现相似，因此这里只简单阐述上述三种模型的结果。从表2中可以看出，在模型E中，也就是在没有变点的情况下，LRSM是完全准确和有效的。在模型F和模型G中，LRSM有97次准确的识别出变点的数量，并且均满足估计出的变点与真实变点的距离小于50。其中有效估计出的变点中，其变点位置的均值与真实变点之间均相差1.5以内。综上所述，LRSM可以有效准确地探测出INAR(2)模型变点的数量以及位置。

模型E模型F

模型G模型H模型I模型J

Figure 2. Sample path from Model E to Model J

图2. 模型E到模型J的样本路径图

3.3. INAR模型的系数之和接近于1时的模拟研究

当使用LRSM处理INAR模型中 ${\displaystyle \underset{i=1}{\overset{p_j}{\sum }}{\alpha }_{j,i}}\to 1$ 的变点问题时，发现LRSM在第二步偶尔会保留非变点的情况，使得最终估计的变点数量很容易大于真实的变点数量。第二步中使用的MDL准则来源于计算理论，是Rissanen [27] 研究通用编码时提出的，并且证明了残差的编码长度由拟合模型的对数似然的负数给出。因此，本文对MDL准则进行了调整。也就是说，残差的编码长度由拟合模型的对数似然的负数给出，而不是条件对数似然的负数。调整后的MDL函数如下：

$\text{MDL}\left(m,J,P\right)=\log \left(m\right)+\left(m+1\right)\log \left(n\right)+{\displaystyle \underset{j=1}{\overset{m+1}{\sum }}\log \left(p_j\right)}+{\displaystyle \underset{j=1}{\overset{m+1}{\sum }}\frac{p_j+1}{2}\log }\left(n_j\right)-{\displaystyle \underset{j=1}{\overset{m+1}{\sum }}{L}_j\left({\hat{\theta }}_j\right)}-{\displaystyle \underset{j=1}{\overset{m+1}{\sum }}{\stackrel{⌢}{L}}_j}\left({\hat{\theta }}_j\right),$

其中 $\stackrel{⌢}{L}\left(\theta \right)=\log P\left(d_1,d_2,\cdots ,d_p\right)+L\left(\theta \right)$ ， $L\left(\theta \right)$ 的定义见式子。为了便于表述，调整后的LRSM被记录为全似然比扫描方法(简记FLRSM)。在本小节中，通过把FLRSM用于由模型H~J生成的分段平稳整数值自回归多变点过程来评估FLRSM的性能，同时将其与以前的LRSM进行比较。模型H~J的样本路径图见图2，模拟结果见表3。

Table 3. Simulation results from Model H to Model J

表3. 从模型H到模型J的模拟结果

¹正确识别变点数量的次数；²有效识别变点的次数；³真实变点的位置；⁴有效识别到的变点的平均值；⁵有效识别到的变点的均方误差。

模型H：具有一个变点的分段平稳INAR(1)过程

$X_t=\{\begin{array}{ll}0.9\circ X_{t-1}+{\epsilon }_{t,1}\left(1\right)\hfill & 1\le t\le 512,\hfill \\ 0.9\circ X_{t-1}+{\epsilon }_{t,2}\left(2\right)\hfill & 513\le t\le 1024.\hfill \end{array}$

模型I：具有两个变点的分段平稳INAR(1)过程

$X_t=\{\begin{array}{ll}0.85\circ X_{t-1}+{\epsilon }_{t,1}\left(1\right)\hfill & 1\le t\le 400,\hfill \\ 0.5\circ X_{t-1}+{\epsilon }_{t,2}\left(2\right)\hfill & 401\le t\le 612,\hfill \\ 0.85\circ X_{t-1}+{\epsilon }_{t,3}\left(1\right)\hfill & 613\le t\le 1024.\hfill \end{array}$

模型J：具有两个变点的分段平稳INAR(2)过程

$X_t=\{\begin{array}{ll}0.83\circ X_{t-1}+0.1\circ X_{t-2}+{\epsilon }_{t,1}\left(1\right)\hfill & 1\le t\le 350,\hfill \\ 0.49\circ X_{t-1}+0.1\circ X_{t-2}+{\epsilon }_{t,2}\left(2\right)\hfill & 351\le t\le 700,\hfill \\ 0.25\circ X_{t-1}+0.1\circ X_{t-2}+{\epsilon }_{t,3}\left(1\right)\hfill & 701\le t\le 1024.\hfill \end{array}$

从表3可以看出，数据波动越大， ${\displaystyle \underset{i=1}{\overset{p_j}{\sum }}{\alpha }_{j,i}}$ 越接近于1，探测变点的难度就越大。在探测变点的数量和位置上，FLRSM的性能优于LRSM。但两种方法得到的均方误差都比较大，估计精度还有待提高。

4. 实证分析

在本节中将LRSM运用到真实的整数值时间序列数据中。在拟合具体模型时需要估计相应的参数值，主要考虑矩估计值和最大似然估计值。通过计算拟合模型预测序列的均方误差，本文最终选择了均方误差偏小的矩估计值。

这个例子以日常观察精神分裂症患者在感知速度测试中获得的分数为观测值，见McCleary和Hay [28] 数据如图3所示。数据由120个连续的日常得分组成。然而从第61天起，患者开始注射一种强效镇静剂(氯丙嗪)，此药剂有望降低感知速度。Neal和Subba [29] 猜测第60天是变点，因此通过两个不同的INAR模型拟合数据的两个部分。然而，他们没有提供任何正式估计变点的方法。Kashikar等人 [30] 通过马尔可夫链蒙特卡罗(MCMC)程序，估计的变点位置在第80天。

Figure 3. Schizophrenic actual data and fitted data

图3. 精神分裂症的实际数据和拟合数据

针对该数据，首先运用LRSM探测出变点位置在第69天，并进行分段INAR(1)模型拟合，最后估计了拟合模型的参数。同时注意到，在服用镇静剂后，患者评分开始下降，这是Neal和Subba猜测变点在第60天的基础。从评分图(图3)可以看出，患者在69天后的评分明显低于69天前的患者评分。80天后，患者评分下降更为明显。其中图中的黑色曲线表示实际数据，红色曲线表示模型拟合值。

这意味着镇静剂从使用的第一天开始起作用，到服药后第9天疗效更明显，20天后，即第80天，药效才得到充分发挥。因此，把第60天、第69天、第80天看作变点是合理的。为此，本文对三个疑似变点分别进行分段的INAR(1)模型拟合，并利用预测值序列 $\left\{{\hat{X}}_t\right\}$ 的均方误差来选择合适的变点，均方误差计算公式如下：

$MSE=\frac{1}{120}{\displaystyle \underset{t=1}{\overset{120}{\sum }}{\left(X_t-{\hat{X}}_t\right)}^2},$

其中预测值序列 ${\hat{X}}_t=E\left(X_t|X_{t-1}\right)$ 。通过计算得到，当变点为60天时，此时预测值的均方误差为137.31；变点为第69天时相应的均方误差为110.31；而变点为第80天时，模型预测值的均方误差为111.81。这清楚地表明，当变点为第69天时，拟合的INAR(1)模型的性能更好。同时也可以观察到图3所示的模型的拟合曲线支持变点为第69天时拟合的INAR(1)模型。此时相应拟合模型为：

$X_t=\{\begin{array}{ll}0.455\circ X_{t-1}+{\epsilon }_{t,1}\left(38.14\right)\hfill & 1\le t\le 69,\hfill \\ 0.64\circ X_{t-1}+{\epsilon }_{t,2}\left(13.62\right)\hfill & 70\le t\le 120.\hfill \end{array}$

5. 小结

本文主要运用LRSM探测分段平稳的整数值时间序列中变点的数量和位置。此外，针对模型系数之和趋近于1的情况，提出了FLRSM，有效提高了变点探测精度。最后，将LRSM应用于一个真实的生物特征数据集，即精神分裂症患者数据。该数据集里面存在一个变点，而LRSM能有效地探测出该变点位置。通过模拟研究和实际数据分析的结果表明：LRSM能够有效地探测出由INAR模型生成的分段平稳整数值自回归过程中的多变点。

基金项目

辽宁科技大学博士启动资金(601010391)。

NOTES

^*通讯作者。