ARIMA模型在预测新型冠状病毒传播中的应用

Application of ARIMA Model in Transmission Prediction of the COVID-19

摘要: 目前，新冠病毒仍在全球肆虐。为了更好地模拟病毒的传播情况，本文基于疫情数据，应用ARIMA模型对巴基斯坦的新冠疫情趋势进行预测，对模型的相关性质进行检验。最后分析了该模型在传染病预测中的应用价值，为新冠疫情预测提供实践经验。

Abstract: At present, the COVID-19 is still rampant in the world. In order to better simulate the spread of the virus, based on the epidemic data of Pakistan, the ARIMA model is used to predict the trend of the COVID-19 in this paper, and the relevant properties of the model are tested. Finally, the application value of the model in the prediction of infectious diseases is analyzed, which provides practical ex-perience for the prediction of the COVID-19.

1. 模型背景

1.1. 问题选择

从2019年冬天开始，新型冠状病毒引起的肺炎疫情首先在中国爆发，然后在全球蔓延，对整个世界的人类安全和经济产生了巨大影响。世卫组织报告显示，新冠肺炎的全球死亡率约为6.4%。虽然一些地区疫情得到了有效控制，但全球疫情仍然非常严重。合理有效的数学建模可以为政府的防疫决策、疾病预防和控制提供很大的帮助。

为更好预测疫情的未来发展状况，需要对各类预测模型的优劣进行研究探讨。目前，对新冠疫情的预测主要是应用传染病传播动力学模型，如：SEIR模型、SIR模型等，但传染病动力学模型需要对各种模型参数及其未来的变化趋势有了解，而这些参数很难获取，时间序列模型 [1] [2] [3] [4] 则只需要有病例数的历史序列就可以构建病例数的预测模型，所以我们尝试利用时间序列模型来进行预测。

1.2. 区域选择

巴基斯坦位于南亚次大陆西北部，南濒阿拉伯海，东接印度，东北邻中国，西北与阿富汗交界，西邻伊朗，是我国重要的战略合作伙伴。我国的新型冠状病毒肺炎疫情虽然呈现多点散发态势，但总体已经得到了有效控制，而巴基斯坦每天仍有新增确诊病例，未来疫情发展趋势是大家关注的问题 [5]。

本文基于2021年1月15日至3月15日的巴基斯坦疫情发展相关的时间序列数据，通过数据处理、曲线拟合、参数估计等过程，反复试验，建立ARIMA时间序列预测模型 [2] [3] [4] [5]，对巴基斯坦的新冠疫情趋势进行短期预测，并对模型的合理性进行检验。以此来探讨ARIMA时间序列模型在传染病未来发展走向的预测中的应用方法与研究价值，为新冠疫情预测提供实践经验。

2. 模型理论

2.1. 差分

差分，一般用于以时间为统计维度的分析中，反应了离散量之间的一种变化，它可以减轻数据之间的不规则波动，使其波动曲线更平稳。时间距离为1的两个序列值，做减法运算，可以得到1阶差分，迭代下去，可以得到p阶差分，一般记为： ${\nabla }^p{x}_t={\nabla }^{p-1}{x}_{t-1}$。

2.2. 白噪声的检验

下面来检验序列是否为纯随机序列，也即是进行白噪声检验。理想的结果是观察值序列为非白噪声序列，这样我们就可以使用平稳时间序列对模型进行预测；而残差序列应具有完全随机性，理想的结果是白噪声。

假设条件：H₀： ${\rho }_1={\rho }_2=\cdots ={\rho }_m=0,\forall m\ge 1$ ；

H₁：至少存在某个 ${\rho }_k\ne 0,\forall m\ge 1,k\le m$。

构造统计量：

$Q=n{\displaystyle \underset{k=1}{\overset{m}{\sum }}{\hat{\rho }}_k^2}~{\chi }^2\left(m\right).$

此时可以比较统计量和上分位点 ${\chi }_{1-\alpha }^2\left(m\right)$，来做出拒绝或者接受原假设的判断；或者计算统计量的P值，当P值小于 $\alpha$ 时，能够在 $1-\alpha$ 的置信水平下拒绝原假设，该序列不是白噪声；否则就接受原假设，序列是白噪声。如果样本数量较小，可以使用LB统计量：

$\text{LB}=n\left(n+2\right)\underset{k-1}{\overset{m}{{\displaystyle \sum }}}\left(\frac{{\hat{\rho }}_k^2}{n-k}\right)~{\chi }^2\left(m\right).$

2.3. ARIMA模型

本文选用的主要模型是求和自回归移动平均(ARIMA(p,d,q))模型，该模型能够将非平稳时间序列转化为差分平稳序列，应用差分平稳序列，我们可以进行ARIMA拟合，最后进行预测。

ARIMA(p,d,q)模型的基本形式：

$\{\begin{array}{l}\Phi \left(B\right){\nabla }^d{x}_t=\Theta \left(B\right){\epsilon }_t\\ E\left({\epsilon }_t\right)=0,Var\left({\epsilon }_t\right)={\sigma }_{\epsilon }^2,E\left({\epsilon }_t{\epsilon }_s\right)=0,s\ne t\\ E\left(x_s{\epsilon }_t\right)=0,\forall s<t\end{array}$，

其中， ${\nabla }^d={\left(1-B\right)}^d$ ； $\Phi \left(B\right)=1-{\phi }_{1}B-\cdots -{\phi }_p{B}^p$， $\Theta \left(B\right)=1-{\theta }_{1}B-\cdots -{\theta }_p{B}^p$ 分别是平稳可逆ARMA(p,q)模型的自回归系数多项式和移动平滑多项式； $\left\{{\epsilon }_t\right\}$ 是均值为0的白噪声序列。

ARIMA模型的实质就是差分运算与ARMA模型的组合。这意味着任何非平稳序列如果能通过适当阶数的差分运算变得平稳，就可以对差分后的序列进行ARMA拟合。

ARIMA模型建模过程：

1) 获得观察值序列；

2) 平稳性检验：若序列不平稳，检验差分后的序列是否平稳；

3) 白噪声检验：通过平稳性检验的序列若不是白噪声序列就可以拟合ARMA模型；

4) 利用该模型预测未来趋势。

2.4. 残差自相关检验

下面来检验模型的拟合效果。如果残差序列满足

$E\left({\epsilon }_t,{\epsilon }_{t-j}\right)=0,\forall j\ge 1,$

则说明是白噪声序列，因此模型拟合效果较好。下一步就可以充分提取序列中的相关信息，不需要修正拟合模型；反之，如果残差序列满足

$E\left({\epsilon }_t,{\epsilon }_{t-j}\right)\ne 0,\forall j\ge 1,$

则说明自相关性显著，模型拟合效果不好，需要修正模型再次拟合。具体操作如下：

假设条件：H₀： $E\left({\epsilon }_t,{\epsilon }_{t-1}\right)=0$，即 $\rho =0$ ；H₁： $E\left({\epsilon }_t,{\epsilon }_{t-1}\right)\ne 0$，即 $\rho \ne 0$

DW统计量： $\text{DW}=\frac{{\displaystyle {\sum }_{t=2}^n{\left({\epsilon }_t-{\epsilon }_{t-1}\right)}^2}}{{\displaystyle {\sum }_{t=2}^n{\epsilon }_t^2}}\approx 2\left(1-\frac{{\displaystyle {\sum }_{t=2}^n{\epsilon }_t{\epsilon }_{t-1}}}{{\displaystyle {\sum }_{t=1}^n{\epsilon }_t^2}}\right)$

自相关系数： $\rho =\frac{{\displaystyle {\sum }_{t=2}^n{\epsilon }_t{\epsilon }_{t-1}}}{{\displaystyle {\sum }_{t=1}^n{\epsilon }_t^2}}$

即 $\text{DW}\approx 2\left(1-\rho \right)$

当 $0<\rho \le 1$ 时，序列正相关，且 $\rho \to 1$ 时， $\text{DW}\to 0$ ； $\rho \to 0$ 时， $\text{DW}\to 2$。当 $-1<\rho \le 0$ 时，序列负相关，且 $\rho \to -1$ 时， $\text{DW}\to 4$ ； $\rho \to 0$ 时， $\text{DW}\to 2$。DW值越接近于2， $\rho$ 值越小，拟合效果越好。

2.5. 参数的显著性检验

为了使模型更精简，需要对参数做显著性检验。如果某个参数所对应的自变量对因变量的影响不明显，经过显著性检验之后就可以把对应的自变量剔除，进而实现对模型的精简。

做假设检验：

${\text{H}}_0:{\beta }_j=0\leftrightarrow {\text{H}}_1:{\beta }_j\ne 0,\forall 1\le j\le m.$

构造t检验统计量：

$T=\sqrt{n-m}\frac{{\hat{\beta }}_j}{\sqrt{a_{ij}Q\left(\tilde{\beta }\right)}}~t\left(n-m\right).$

代入样本值，如果满足

$\left|T\right|\ge t_{1-\alpha }\left(n-m\right),$

则拒绝原假设；或者计算检验统计量的P值，当P值小于 $\alpha$ 时，也能够在 $1-\alpha$ 的置信水平下拒绝原假设，则这个参数效果显著，不能剔除。反之，如果这个参数的效果不显著，可以在模型中剔除该参数，重新拟合模型即可。

3. 模型建立与问题求解

3.1. 获得观察值序列

从网易疫情实时数据平台获取巴基斯坦地区2020年3月至2021年6年新冠肺炎确诊人数数据，对该数据进行初步筛选分析后将2021年1月15日至3月15日这60组确诊人数数据作为研究对象。

3.2. 判断序列平稳性

确诊人数序列时序如图1所示。

Figure 1. Time series of confirmed infection

图1. 确诊人数时序图

从图1中可以看出该序列有显著的递增趋势，为典型的非平稳序列。可以考虑对该序列进行差分，获得平稳序列，建立 模型。

3.3. 对原序列进行差分运算

首先尝试对原序列进行一阶差分，一阶差分后序列时序如图2所示：

Figure 2. First order difference time series

图2. 一阶差分时序图

由图2我们发现，差分后序列波动较大，仍不平稳，可以进一步对一阶差分序列再次差分得到二阶差分时序图3。

Figure 3. First order difference time series

图3. 一阶差分时序图

图3显示，二阶差分后序列在均值附近比较平稳地波动，为了进一步确定该时序的平稳性，考察二阶差分后序列的自相关图4。

自相关图4显示序列具有较强的短期相关性，所以可以初步认为二阶差分后序列平稳。

3.4. 对平稳的二阶差分序列进行白噪声检验

白噪声检验结果如表1所示。

Figure 4. Autocorrelogram

图4. 自相关图

Table 1. White noise autocorrelation test

表1. 白噪声自相关检验表格

取显著性水平α = 0.05，从表1可以发现，不管是延迟6阶还是12阶，计算出来的P值均小于0.05，分别为0.0302和0.0168，因此在α = 0.05的显著性水平下，拒绝原假设。因此需要进一步地提取相关信息。

3.5. 对平稳非白噪声差分序列拟合ARIMA模型

观察图5中的二阶差分后序列的自相关图，发现序列的自相关系数具有1阶截尾的特性，下面来计算序列的偏自相关性。

观察偏自相关图5，可以发现序列的截尾性不显著，因此在拟合二阶差分后序列的时候，可以应用MA(q)模型。

3.6. ARIMA模型定阶

为了按照最小信息量准则确定一个相对最优的模型，利用SAS软件分析数据，得到最小信息量表2，结果如下：

Figure 5. Partialautocorrelogram

图5. 偏自相关图

Table 2. Minimum information

表2. 最小信息量表格

根据最小信息量准则，AR(5)模型的结果最优，可以用来拟合差分序列。

Table 3. Conditional least squares estimation and significance test

表3. 条件最小二乘估计以及显著性检验表格

由条件最小二乘估计以及显著性检验表3可以得到拟合模型为：

$\left(1+0.34258B+0.04303B^2+0.31251B^3+0.07688B^4-0.00603B^5\right)\nabla x_t={\epsilon }_t$

3.7. ARIMA模型预测

首先，我们不对参数进行剔除与检验，直接利用该组参数所对应的AR(5)模型进行序列预测，得到的预测值如表4所示：

Table 4. Model AR(5) prediction

表4. AR(5)模型预测表格

为了探究模型是否可得到一定程度的简化，下面对模型不显著的参数进行剔除，只留下通过显著性检验的变量，同时重新估计参数值，得到表5：

Table 5. Parameter estimation and significance test of sparse coefficient model

表5. 疏系数模型参数估计与显著性检验表格

即最终得到的拟合模型为：

$\left(1+0.3116B+0.28172B^3\right)\nabla x_t={\epsilon }_t$

为了判断该模型是否显著，是否还有其他可提取的信息，我们来检验残差序列是否为白噪声序列，结果见表6。从表中可以看出，P值均大于显著性水平0.05，因此该模型是显著的，残差序列具有纯随机性，因此可以应用该模型来进行预测。

Table 6. Residual autocorrelation test

表6. 残差自相关检验表格

利用此拟合模型对确诊人数序列进行7期(一周)预测，预测值以及95%的预测区间如表7所示：

Table 7. Sparse coefficient prediction

表7. 疏系数预测表格

我们将参数剔除前与剔除后的预测值与真实值做对比，并分别计算偏差平方和，以对比剔除前后模型预测效果的优劣，具体结果见表8：

Table 8. Comparison of predicted values

表8. 预测值对比结果

剔除前后两种模型的预测图像如图6所示：

对比可知，剔除前的模型虽然AIC信息量较大，但其拟合的偏差平方和较小，其预测效果相对于剔除后的疏系数模型较为准确。

Figure 6. Predictionresults of AR (5) and sparse coefficient models

图6. AR(5)以及疏系数模型预测图像

4. 模型评价与推广

4.1. 模型评价

本文针对巴基斯坦近两个月新冠肺炎确诊人数非线性非平稳的数据，基于ARIMA模型建立时间序列模型进行预测。首先通过差分得到平稳时间序列模型，基于理论知识与SAS软件的应用对模型的平稳性、参数及残差相关性进行检验。然后用该模型进行短期预测。将预测值与真实值对比发现在只有传染病随时间变化的时间序列历史数据的前提条件下，ARIMA时间序列模型对传染病的未来发展趋势具有较好的预测效果。

但由于时间序列模型对短期预测具有较好的效果这种局限性以及该模型是在传染病的传播发展模式未发生根本性变化的情况下对传染病发展趋势所做的预测，通常难以保证长期预测结果的准确性，该模型只适用于对短期确诊人数的预测。

4.2. 模型推广

在模型后续的改进和优化时可以考虑及时将新的数据加入时间序列，将这种建模及预测过程全程自动化并与病例监测报告系统进行集成，实现实时动态建模和预测，可为今后开展疫情监测提供便捷的手段。

基金项目

中国石油大学(北京)教育教学改革项目(编号XM10720210153、XM10720200035、YJS2020032，YJS2020033)。

参考文献