Linear Time Series Analysis: Stationarity, ARMA, and ARIMA

引论与一些基本概念

基本定义

时间序列是一系列时间点上的数据 这些时间点可以有序排列

观测一系列时间点上的数据在研究中非常常见 所以时间序列分析在不同的方面应用非常的广泛

时间序列分析的研究目的主要有两个方面

• 研究时间序列的产生机制 了解过去
• 基于历史数据和其他相关因素 对未来的可能性做出预测 预测未来

为了实现时间序列分析 我们往往需要要求时间序列存在一些内在的结构 如果他就是完全随机的 我们往往缺少研究时间序列的必要

时间序列的一个举例

如下图所示 这就是一个时间序列 它包含的一个时间（可以是年份 月份 小时等）作为横坐标 一个观测值作为纵坐标 这种图就是表示时间序列最基本的图形

另一种常见的时间序列图 它使用上一年的数据作为横坐标 当年的数据作为纵坐标 ，这种图可以让我们研究两年的观测值之间是否存在一定的影响

时间序列分析对我们研究图的能力要求比其他领域更高，我们需要根据研究目标作图 分析图中的信息 培养对图表的理解能力非常重要 我们往往会根据图反应的信息来选取合适的时间序列分析模型

线性回归与时间序列分析

时间序列不能等同于回归分析

时间序列数据的特征是：

• 自相关性（线形or非线性）
• 不可交换性（样本顺序不能交换）

回归分析的应用场景中最简单的数据结构是可交换独立同分布数据，当时间序列数据满足一些条件的时候 可以使用回归分析的手法来处理 （比如$AR$模型和普通的多元线性回归模型有着相同的形式）但是我们不能认为他们相同

时间序列分析中有很多独有的方法 他们和回归分析并没有什么联系 所以时间序列不能等同于回归分析

随机过程与时间序列分析

时间序列是来自随机过程的一组观测：我们把它对应的随机过程称为时间序列过程

时间序列过程是一种特殊的随机过程 他的下标是时间变量

时间序列的数字特征

比较常用的如下所示 均值函数定义为

$$\mu_{t}=E(Yt)$$

自协方差函数定义为

$$Cov(X_s,X_t)=\operatorname{E}\left[(X_s-m_X(s))(X_t-m_X(t))\right]$$

自相关函数定义为

$$p_{s,t}=\mathrm{Corr}(Y_{t},Y_{s})=\frac{\mathrm{Cov}(Y_{t},Y_{s})}{\sqrt{\mathrm{Var}(Y_{t})\mathrm{Var}({Y}_{s})}}=\frac{\gamma_{t,s}}{\sqrt{\gamma_{t,t}\gamma_{s,s}}}$$

他们有一些基础的性质为

$$\begin{aligned}\gamma_{\iota.\iota}&=\operatorname{Var}(Y_{\iota})&\rho_{\iota.\iota}&=1\\\gamma_{t.s}&=\gamma_{s.t}&\rho_{t,s}&=\rho_{s,t}\\\mid\gamma_{\iota,s}|&\leqslant\sqrt{\gamma_{\iota.t}\gamma_{s,s}}&\mid\rho_{t.s}\mid\leqslant1\end{aligned}$$

我们这里给出一个有用的定理

$$\mathrm{Cov}\bigg[\sum_{i=1}^{m}c_{i}Y_{t_{i}},\sum_{j=1}^{n}d_{j}Y_{s_{j}}\bigg]=\sum_{i=1}^{m}\sum_{j=1}^{n}c_{i}d_{j}\mathrm{Cov}(Y_{t_{i}},Y_{s_{j}})$$

我们在这里再简单复习一下随机过程平稳性的概念 他在时间序列分析中确实有点太常用了 随机过程基础 的“平稳过程的定义”一节

时间序列的分解

时间序列需要有内在的结构才可以被我们分析 下面是一种比较完全的结构形式

$$X_t=T_t+S_t+R_t,t=1,2,\ldots$$

• 趋势项
• 季节项
• 随机项 至于如何估计这些项 我们有很多自然的想法 比如用回归拟合趋势 剩余的项中用季度平均捕捉季节等 我们后面会展开详细的研究

趋势项和季节项都可以被当做非随机的时间序列处理, 它们的预测问题往往是简单的。 随机项一般是平稳序列。

常见时间序列举例

随机游走

令$e_1,e_2,...$为均值为0，方差是$\sigma^2$。的独立同分布的随机变量序列，观测时间序列 $\{Y_t:t=$ $1,\:2,\:...\}$ 构造如下：

$$\left.\left.\begin{array}{rl}Y_1&=e_1\\Y_2&=e_1+e_2\\&\vdots\\Y_t&=e_1+e_2+\cdots+e_t\end{array}\right.\right\}$$

我们可以轻松的计算到 均值函数$\mu_{t}=0$ 方差函数 $\mathrm{Var}(Y_{\imath})=t\sigma_{e}^{2}$ 自相关函数 $\rho_{t,s}=\frac{\gamma_{t,s}}{\sqrt{\gamma_{t,t}\gamma_{s,s}}}={\sqrt{\frac{t}{s}}}$

基于此我们可以给出一些关于随机游动过程的理解

随着时间的推移，相邻时点上$Y$值的正相关程度越来越强，而另一方面，对时点相距遥远的$Y$值，其相关程度越来越弱，对单位根过程的样本作ACF图， 其衰减速度很慢很慢

虽然理论均值为0 但是方差随着时间增长 所以预期过程会在远离零的位置摆动，也说明了此模型不可预测

滑动平均

还是前面的假设 构造

$$Y_t=\frac{e_t+e_{t-1}}2$$

均值函数$\mu_{t}=0$ 方差函数 $\mathrm{Var}(Y_{\imath})=0.5\sigma_{e}^{2}$ 自相关函数

$$\rho_{t,s}=\begin{cases}1&\mid t-s\mid=0\\0.5&\mid t-s\mid=1\\0&\mid t-s\mid>1\end{cases}$$

滑动平均过程的是一个平稳过程 一般用来作为引入的例子

白噪声

平稳过程中一个很重要的例子是所谓的白噪声过程，定义为同分布的随机变量序列$\langle e_i\rangle$. 其重要性并非源自它是有趣的模型，而是因为许多有用的过程可由白噪声过程构造出来

对于白噪声过程有 $E[e_t]=\mu$ $cov(e_t,e_s)=\begin{cases}\sigma^2&t=s\\0&else\end{cases}$ 如果 随机变量序列相互独立 则称为独立白噪声 如果此时均值为0 则称为零均值独立白噪声 此时如果方差为1 则称为标准独立白噪声

我们一般只研究最一般的标准独立白噪声情况 他明显是一个严平稳过程 我们能看出 随机游走和滑动平均都是根据白噪声过程构造的

随机余弦波

$$Y_{t}=\cos\left[2\pi(\frac{t}{12}+\Phi)\right]$$

其中 $\Phi\sim U(0,1)$ 我们能看出 这个随机过程有着很强的确定性 它具有周期性 他唯一的随机性在于我们的初相如何进行选取 我们可以计算出来 均值函数 $\mu_{t}=0$ 自相关函数 $\rho_k=\cos(2\pi\frac k{12})$ 也就是 我们能判断 这是一个平稳的时间序列 观测滑动平均和随机余弦波的模拟序列图能看出 仅仅依靠我们的观测时间序列图 想要判断时间序列是否平稳是不现实的 我们需要再后面需要其他的处理方法

差分平稳化

我们知道随机游走序列不是平稳的 但是我们我们不单独考虑 $Y_t$而是考虑他的差分 也就是

$$Z_t=Y_t-Y_{t-1}$$

容易看出 差分后的序列是平稳的 也就是我们可以依靠差分这样简单的技巧 让原本并不平稳的序列拥有统计学上的平稳状态

趋势

这是对本文“时间序列的分解”一节的一些解释 解释趋势项的含义

确定性趋势和随机趋势

趋势是对现实情况的观测得到的产物，但是我们知道，时间序列具有随机性，我们所看到的趋势不一定是序列真实具有的特征，因此我们需要多次的模拟，从多个表面趋势中寻找真正的内在的信息

• 随机趋势：多次模拟的结果表现趋势完全不同的时间序列
• 确定性趋势：多次模拟的结果表现的趋势接近的时间序列

非常明显的是：如果我们只有一次观测的机会，那么完全没有办法判断趋势是随机还是确定性的，我们在后面的研究中将会以确定性趋势为主导

常数均值

作为最简单的情况，也是所谓平稳的时间序列 我们假设均值函数是常数

$$Y_t=\mu+X_t$$

其中随机扰动$X_t$有$EX=0$ 当不增加任何其他假设的情况下 我们可以对均值作下面的估计

$$\overline{Y}=\mu$$

为了研究这种估计的精度 我们需要对$X_t$增加更多的假设 其中不同的假设会有着不同的精度估计，这里对此不多做研究了

非常数均值

研究那些并不平稳的序列 我们对其中均值的变化要作为一系列假设才可以具体的分析问题

线性趋势

考虑

$$\mu_t=\beta_0+\beta_1t$$

想要对这种情况进行估计 我们需要在回归分析中最为常用的最小二乘法

季节性趋势

还是我们最基础的模型

$$Y_t=\mu_t+X_t$$

一种比较常用的季节模型的思路是按照每个月份来

$$\mu_t=\begin{cases}\beta_1&t=1,13,25,\cdots\\\beta_2&t=2,14,26,\cdots\\\vdots\\\beta_{12}&t=12,24,36,\cdots\end{cases}$$

对于这种趋势的估计依赖经验和一些我们没有进行过介绍统计学方法来处理，也需要统计软件的使用，我们知道这样可以分析就够了

余弦趋势

季节均值模型包含了很多独立的参数，但是它对季节趋势的形状是没有考虑的，只有很少的时间序列模型在接近的时候反而是没有什么相关性的，所以我们引入了余弦趋势模型，它在存在周期性的时候非常常用 考虑下面的均值模型

$$\mu_{t}=\beta\mathrm{cos}(2\pi ft+\Phi)$$

我们可以把它进行下面的变形，方便使用回归来处理这个问题

$$\mu_{t}=\beta_{0}+\beta_{1}\cos(2\pi ft)+\beta_{2}\sin(2\pi ft)$$

平稳随机时间序列ARMA

我们来讨论一大类参数时间序列模型 自回归滑动平均模型（ARMA）的基本概念，他们在真实过程的建模中有着重要的作用，他们的核心特征是平稳

一般线性过程

我么认为$Y_t$是我们观测到的时间序列，认为$e_t$是没有观测到的白噪声序列； 不妨设这个白噪声序列的均值为0，也就是抹除了均值信息，在真实的建模中，可以认为是我们减掉了均值 那么一般线性过程就可以看成过去和现在的白噪声的加权线性组合

$$Y_{t}=e_{t}+\psi_{1}e_{t-1}+\psi_{2}e_{t-2}+\cdots$$

右边是一个无穷级数，我们后面会对它进行更多的限制来得到我们想要研究的目标

滑动平均过程（MA）

一个非常自然的思想是，离的非常远的噪声对现在的影响可以忽略不计，所以我们可以把一般线性过程改造为

$$Y_t=e_t-\theta_1e_{t+1}-\theta_2e_{t-2}-\cdots-\theta_qe_{t-q}$$

系数就是随便取的，取成正的也可以 我们称该方程为$q$阶滑动平均过程 记为$MA(q)$ 所谓滑动平均 是指根据权数进行平均，然后滑动一次重新进行平均，以此类推

MA（1）

模型表达为

$$Y_{t}=e_{t}-\theta e_{t-1}.$$

容易计算得到

$$E(Y_{t})=0$$ $$\mathrm{Cov}(Y_t,Y_{t-1})=\mathrm{Cov}(e_t-\theta e_{t-1},e_{t-1}-\theta e_{t-2})=\mathrm{Cov}(-\theta e_{t-1},e_{t-1})=-\theta \sigma_t^2$$ $$\mathrm{Cov}(Y_t,Y_{t-2})=\mathrm{Cov}(e_t-\theta e_{t-1},e_{t-2}-\theta e_{t-3})=0$$

能看出，过程大于1阶滞后以后，不存在自相关 我们可以根据具体的$\theta$的取值来分析相关系数的大小，也可以根据散点图的表现分析相关性，正如我们在开篇中展现的一样 本文“时间序列举例”一节

MA（2）

模型为

$$Y_{t}=e_{t}-\theta_{1}e_{t-1}-\theta_{2}e_{t-2}$$

计算协方差函数有

$$\gamma_0=\mathrm{Var}(Y_t)=\mathrm{Var}(e_t-\theta_1e_{t-1}-\theta_2e_{t-2})=(1+\theta_1^2+\theta_2^2)\sigma_{t}^{2}$$ $$\begin{aligned} \gamma_{1}& =\mathrm{Cov}(Y_{i},Y_{i-1})=\mathrm{Cov}(e_{i}-\theta_{1}e_{i-1}-\theta_{2}e_{i-2},e_{i-1}-\theta_{1}e_{i-2}-\theta_{2}e_{i-3}) \\ &=\mathrm{Cov}(-\theta_1e_{t-1},e_{t-1})+\mathrm{Cov}(-\theta_1e_{t-2},-\theta_2e_{t-2}) \\ &=[-\theta_1+(-\theta_1)(-\theta_2)]\sigma_e^2=(-\theta_1+\theta_1\theta_2)\sigma_e^2 \end{aligned}$$ $$\begin{gathered} \gamma_{2} =-\theta_{2}\sigma_{e}^{2} \end{gathered}$$

也就是$MA(2)$ 模型大于二阶滞后不存在自相关

MA（q）

我们直接给出关于相关函数的结论

$$\rho_k=\begin{cases}\frac{-\theta_k+\theta_1\theta_{k+1}+\theta_2\theta_{k+2}+\cdots+\theta_{q-k}\theta_q}{1+\theta_1^2+\theta_2^2+\cdots+\theta_q^2}&\quad k=1,2,\cdots,q\\0&\quad k>q\end{cases}$$

结论已经很明显了

$MA(q)$过程在大于$q$阶滞后的时候不存在自相关

关于滑动平均过程我们研究这些已经完全足够了，下面我们来介绍另一种重要的模型

自回归过程（AR）

顾名思义，自回归指的是用自身作为回归变量。如下形式，$p$阶自回归过程满足方程

$$Y_t=\phi_1Y_{t-1}+\phi_2Y_{t-2}+\cdotp\cdotp\cdotp+\phi_pY_{t-p}+e_t$$

也就是$p$个滞后项和新的信息项$e_t$

AR（1）

还是从简单的形式开始研究 模型形式为

$$Y_{t}=\phi Y_{t-1}+e_{t}$$

我们还是假设均值为0，以后它对应的是平稳的条件 我们计算

$$\gamma_0=\frac{\sigma_e^2}{1-\phi^2}$$ $$\gamma_k=\phi^k\frac{\sigma_e^2}{1-\phi^2}$$ $$\rho_k=\frac{\gamma_k}{\gamma_0}=\phi^k$$

由于方差必须为正 所以$-1<\phi<1$ 因此 我们知道 自相关函数呈现指数递减的特征 根据正负的不同会有锯齿状了平滑状两种 我们不妨取$\phi=0.9$ 能看出 哪怕是三阶滞后都存在比较强的自相关性 我们很容易证明 当且仅当$-1<\phi<1$ 的时候 $AR(1)$过程是平稳的

AR（2）

模型的形式为

$$Y_t=\phi_1Y_{t-1}+\phi_2Y_{t-2}+e_t$$

为了满足模型是平稳的 我们需要满足条件

$$\phi_1+\phi_2<1,\quad\phi_2-\phi_1<1\text{,}\quad|\phi_2|<1$$

为了研究自相关函数 我们沿用上面的思路得到下面的递推方程 Yule-Walker 方程

$$\rho_{k}=\phi_{1}\rho_{k-1}+\phi_{2}\rho_{k-2},\quad k=1,2,3,\cdotp\cdotp\cdotp$$

其中为了进行递推有

$$\rho_1=\frac{\phi_1}{1-\phi_2}$$ $$\rho_2=\phi_1\rho_1+\phi_2\rho_0=\frac{\phi_2(1-\phi_2)+\phi_1^2}{1-\phi_2}$$

我们研究自相关函数的递推公式可以分析自相关函数的性质，随着滞后阶数的增加自相关系数指数递减，正负变动都有可能 这个递减可能是指数递减，也可能是阻尼正弦波动

AR（p）

考虑模型为

$$Y_t=\phi_1Y_{t-1}+\phi_2Y_{t-2}+\cdotp\cdotp\cdotp+\phi_pY_{t-p}+e_t$$

给出平稳性的必要条件为

$$\begin{aligned}\phi_1+\phi_2+\cdots+\phi_p<1\\|\phi_p|<1\end{aligned}$$

给出Yule-Walker方程为

$$\begin{aligned}\rho_1&=\phi_1+\phi_2\rho_1+\phi_3\rho_2+\cdots+\phi_p\rho_{p-1}\\\rho_2&=\phi_1\rho_1+\phi_2+\phi_3\rho_1+\cdots+\phi_{p,p-2}\\&\vdots\\\rho_p&=\phi_1\rho_{p-1}+\phi_2\rho_{p-2}+\phi_3\rho_{p-3}+\cdots+\phi_p\end{aligned}$$

只要给定具体的系数值 我们就可以通过解方程的思路求解相关系数 给出相关系数的性质有：相关系数是一些阻尼递减项和一些阻尼正弦波动项的线性组合

自回归滑动平均混合过程（ARMA）

如果模型中一部分是自回归的，另一部分是滑动平均的 则有一个比较普遍的模型形式（还是平稳的）

$$Y_t=\phi_1Y_{t-1}+\phi_2Y_{t-2}+\cdots+\phi_pY_{t-p}+e_t-\theta_1e_{t-1}-\theta_2e_{t-2}-\cdots-\theta_qe_{t-q}$$

我们称为$ARMA(p,q)$ 我们下面只介绍一种比较简单的形式

ARMA（1）

模型定义为

$$Y_t=\phi Y_{t-1}+e_t-\theta e_{t-1}$$

由于含有$AR$成分 我们需要考虑平稳性条件有

$$\mid\phi\mid<1$$

通过一系列运算有 自相关函数为

$$\rho_k=\frac{(1-\theta\phi)(\phi-\theta)}{1-2\theta\phi+\theta^2}\phi^{k-1},\quad k\geqslant1$$

他也是一种指数递减的形式 递减的初始值依赖于$p_0$ (依赖于$\theta$) 这就和$AR(1),MA(1)$ 都不一样了 他们一个是一阶滞后依赖于$\theta$ 多阶滞后相关为0 另一个是指数递减但是从1开始

对于一般化的ARMA模型 我们在满足平稳性条件

$$\text{当月仅当 AR 特征方程 }\phi(x)=0\text{ 的根的模大于 }1$$

此时有自相关函数满足

$$\rho_k=\phi_1\rho_{k-1}+\phi_2\rho_{k-2}+\cdots+\phi_p\rho_{k-p},\quad k>q$$

一个类似Yule-Walker方程的形式 当$k<q$的时候 自相关函数会含有$\theta$成分

可逆性

我们前面其实能看到 对于$MA$模型 我们用$MA(1)$ 举例有

$$\mathrm{Cov}(Y_t,Y_{t-1})=\mathrm{Cov}(e_t-\theta e_{t-1},e_{t-1}-\theta e_{t-2})=\mathrm{Cov}(-\theta e_{t-1},e_{t-1})=-\theta \sigma_t^2$$

然后有

$$\rho_{1}=(-\theta)/(1+\theta^{2})$$

容易发现 代入$\theta,\frac{1}{\theta}$ 得到的相关系数是相同的

也就是哪怕我们代入一个已经知道的相关系数的值，得到的系数并不唯一 对应的模型不唯一，这个问题必须解决

我们知道（平稳的）AR过程可以被表示为一般线性过程。那么MA过程呢 我们考虑

$$Y_t=e_t-\theta e_{t-1}$$

转化为

$$e_{t}=Y_{t}+\theta e_{t-1}$$

然后不断的代换其中的$e_{t-1}$ 有

$$e_t=Y_t+\theta Y_{t-1}+\theta^2Y_{t-2}+\cdotp\cdotp\cdotp$$

所以有

$$Y_t=(-\theta Y_{t-1}-\theta^2Y_{t-2}-\theta^3Y_{t-3}-\cdots)+e_t$$

也就是 如果我们有 $|\theta|<1$ 则$MA(1)$ 可以转化为一个自回归模型 此时我们称其为可逆的$MA$ 模型 对于一般的$MA,ARMA$ 模型 他们可逆的条件是：特征方程的根的模大于1 我们容易证明：对于可逆的MA过程，给定自相关函数的情况下可以得到一组唯一的参数 我们一般研究的$ARMA$模型都是同时满足平稳和可逆的

非平稳随机时间序列ARIMA

并不是所有事件序列模型都是平稳的；事实上 现实世界中大部分的时间序列模型都是非平稳的 强行使用本文“平稳随机时间序列ARMA”一节中的方法建模只会产生一些荒谬的结论；

幸运的是 我们只需要一些很简单的手段就可以研究非平稳的问题；

AR模型涉及平稳的问题（系数） MA模型涉及可逆的问题（还是和系数有关系） 最后导致ARMA模型 ARIMA模型将会同时研究他们

差分平稳化

我们考虑一个$AR(1)$ 模型

$$Y_t=3Y_{t-1}+e_t$$

这不符合前面介绍的$AR(1)$ 模型所研究的平稳化条件 事实上

$$\mathrm{Var}(Y_{t})=\frac{1}{8}(9^{t}-1)\sigma_{e}^{2}$$

方差有个指数爆炸的趋势

$$\mathrm{Corr}(Y_{t},Y_{t-k})=3^{k}\sqrt{\frac{9^{t-k}-1}{9^{t}-1}}\approx1.$$

对于较大的$t$ 和中等强度的$k$ 指数爆炸的原因就是这种正相关 也就是无论如何 这个时间序列会呈现指数爆炸（可能正也可能负）

我们考虑

$$Y_t=Y_{t-1}+e_t$$

计算它的一阶差分

$$\nabla Y_t=e_t$$

容易看出 差分可以让原本不平稳的模型实现平稳化 前面所述的现象是广泛存在的，差分有利于将不平稳的模型平稳化，如果一阶差分不够那就再进行差分，在经验上 3阶以内的差分就已经足够了

ARIMA模型

如果一个时间序列模型$\{Y_t\}$ 的$d$次差分是一个平稳的ARMA模型 则称其为一个ARIMA模型 如果差分服从$ARMA(p,d)$ 则称 $Y_t$ 服从$ARIMA(p,d,q)$ 下面考虑一个$ARIMA(p,1,q)$ 令$W_t=Y_t-Y_{t-1}$

$$W_t=\phi_1W_{t-1}+\phi_2W_{t-2}+\cdots+\phi_pW_{t-p}+e_t-\theta_1e_{t-1}-\theta_2e_{t-2}-\cdots-\theta_qe_{t-q}$$

或者表示为

$$\begin{aligned}Y_t-Y_{t-1}=&\phi_1(Y_{t-1}-Y_{t-2})+\phi_2(Y_{t-2}-Y_{t-3})+\cdots+\phi_p(Y_{t-p}-Y_{t-p-1})\\&+e_t-\theta_1e_{t-1}-\theta_2e_{t-2}-\cdots-\theta_qe_{t-q}\end{aligned}$$

可以改写为

$$\begin{aligned}Y_t&=(1+\phi_1)Y_{t-1}+(\phi_2-\phi_1)Y_{t-2}+(\phi_3-\phi_2)Y_{t-3}+\cdots\\&+(\phi_p-\phi_{p-1})Y_{t-p}-\phi_pY_{t-p-1}+e_t-\theta_1e_{t-1}-\theta_2e_{t-2}-\cdots-\theta_qe_{t-q}\end{aligned}$$

要理解这三个形式 第一个体现为了一个$ARMA(p,q)$ 第二个知只是进行了代换 第三个则体现为了一个 $ARMA(p+1,q)$

IMA(1,1)

我们还是从最简单的形式考虑

$$Y_{t}=Y_{t-1}+e_{t}-\theta e_{t-1}$$

计算模型的方差和相关系数有

$$\left.\mathrm{Var}(Y_{i})=\left[\begin{matrix}{1+\theta^{2}+(1-\theta)^{2}(t+m)}\\\end{matrix}\right.\right]\sigma_{\epsilon}^{2}$$

呈现爆炸

$$\begin{aligned} \mathrm{Corr}(Y_{\iota},Y_{t-k})& =\frac{1-\theta+\theta^2+(1-\theta)^2(t+m-k)}{[\operatorname{Var}(Y_t)\operatorname{Var}(Y_{t-k})]^{1/2}}\approx\sqrt{\frac{t+m-k}{t+m}} \\ &\approx1 \end{aligned}$$

爆炸的原因是强相关性

IMA(2,2)

模型形式为

$$\nabla^2Y_t=e_t-\theta_1e_{t-1}-\theta_2e_{t-2}$$

我们不加计算 直接回答 方差和相关系数的呈现情况同本文“IMA(1,1)”一节

ARIMA中的常数项

在本文“平稳随机时间序列ARMA”一节中 常数项是不影响我们的研究的 减去均值以后针对0均值展开研究 最后所有的结果再加上均值就可以了； 但是在ARIMA中 就没有这么简单了 把模型引入差分后的模型是很容易的 如下

$$\begin{aligned} W_{t}-\mu=& \phi_{1}(W_{t-1}-\mu)+\phi_{2}(W_{t-2}-\mu)+\cdots+\phi_{p}(W_{t-p}-\mu) \\ &+e_{t}-\theta_{1}e_{t-1}-\theta_{2}e_{t-2}-\cdots-\theta_{q}e_{t-q} \end{aligned}$$

或者

$$W_t=\theta_0+\phi_1W_{t-1}+\phi_2W_{t-2}+\cdots+\phi_pW_{t-p}+e_t-\theta_1e_{t-1}-\theta_2e_{t-2}-\cdots-\theta_qe_{t-q}$$

其中

$$\theta_{0}=\mu(1-\phi_{1}-\phi_{2}-\cdots-\phi_{p})$$

两者其实是等价的

下面我们需要考虑这样的非0均值对ARIMA模型的原始模型的影响 考虑 IMA(1,1). 迭代$Y_t$有

$$\begin{aligned}Y_t=&e_t+(1-\theta)e_{t-1}+(1-\theta)e_{t-2}+\cdots+(1-\theta)e_{-m}-\theta e_{-m-1}\\&+(t+m+1)\theta_0\end{aligned}$$

也就是一个非0的均值 在IMA(1,1)中 导致了一个线性于时间的趋势 在$d=2$ 的时候 非0均值导致二次时间趋势项

时间序列分析中的数据变换

在很多的现实情况中 我们会呈现出一个百分比增长的趋势 尤其是在经济数据和生物学数据中，这种百分比增长体现在时间序列中实际上是一个指数的增长形式

我们进行对数变换在此时非常有效的 而对数变化就属于一个非常重要的变换族

BoxCox当然不仅仅用于回归分析与方差分析，有许多辅助我们选择具体的$\lambda$的方法，他们一般都是在正态性和方差齐性上研究

数据变换也数据分析中非常常见的一种操作，不和TSA强制绑定，在什么时候怎么选择变换我们也不进行研究，仅仅作为一种提示：时间序列分析中的数据变换也是非常有用的一个手法

模型识别

我们迄今为止已经开发出了一大类时间序列模型 ARIMA 模型；现在我们要从前面的研究开始学习进行统计推断了，下面我们的工作主要分为四个部分

• 对于给定的时间序列数据 选定合适的$p,d,q$值
• 估计选定的模型的参数
• 检验拟合的模型并且进行改进
• 使用前面确定的模型对未来的数据进行预测 也就是模型识别 参数估计 模型检验 模型预测四个部分 我们会在包括本章在内的四章介绍这些问题；

MA模型的识别

观测样本序列 我们可以估计样本的自相关函数有

$$r_{k}=\frac{\sum_{t=k+1}^{n}\left(Y_{t}-\overline{Y}\right)\left(Y_{t-k}-\overline{Y}\right)}{\sum_{t=1}^{n}\left(Y_{t}-\overline{Y}\right)^{2}},\quad k=1,2\cdots$$

我们目的就是从样本的自相关函数来识别出$r_k$的模式 来比较ARIMA模型中已经知道的模式 实现选定合适的$p,d,q$

对于$MA(q)$ 模型 当滞后阶数超过$q$以后 自相关函数为0 这意味着样本自相关函数已经是$MA$过程的一个良好指示器了

使用样本ACF识别

AR模型的识别

但是对于$AR(p)$模型 自相关函数呈现逐渐衰减的特点 我们还不能用样本自相关函数指示$AR(p)$ 模型 我们引入 $k$阶滞后的偏自相关函数形式为

$$\begin{aligned}\phi_{kk}=\mathrm{Corr}(Y_t-\beta_1Y_{t-1}-\beta_2Y_{t-2}-\cdots-\beta_{k-1}Y_{t-k+1},\\Y_{t-k}-\beta_1Y_{t-k+1}-\beta_2Y_{t-k+2}-\cdots-\beta_{k-1}Y_{t-1})\end{aligned}$$

样本的偏自相关函数可以用下面的形式计算

$$\phi_{kk}=\frac{\rho_k-\sum_{j=1}^{k-1}\phi_{k-1,j}\rho_{k-j}}{1-\sum_{j=1}^{k-1}\phi_{k-1,j}\rho_j}.$$

而偏自相关函数有什么特点呢？ 对于$AR(p)$模型 则

$$\phi_{kk}=0,\quad k>p$$

也就是 $AR(p)$ 模型会被偏自相关函数良好指示 $MA(q)$模型没有关于偏自相关函数的指示性，他会指数衰减而不是归零

使用样本PACF识别

ARMA模型的识别

ARMA(p,q)对于两种指示函数都不存在截尾的性质 我们需要提出新的方法来实现这个目标 有不止一种方法被我们提出用来解决ARMA模型的识别问题，我们介绍其中的EACF法（扩展自相关函数）这种方法目前被模拟认为是具有比较优良的性质

EACF法核心思想是：如果一个模型的AR事已知的，则观测序列滤除自回归部分将会得到一个纯MA模型，这个MA模型可以用ACF确定阶数

考虑$ARMA(1,1)$来介绍EACF的使用

$$Y_t=\phi Y_{t-1}+e_t-\theta e_{t-1}$$

此时 应用$Y_t$对$Y_{t-1}$ 做简单的线性回归可以得到$\phi$的不一致估计量（存在系统性偏差，因为回归系数估计的量本质含有$\theta$）但是这个回归的残差确实可以帮助我们分析 第二次回归用$Y_t$ 对第一次回归残差的一阶滞后进行回归 得到的系数$\tilde{\phi}$ 就是$\phi$的一致估计量 也就是

$$W_{t}=Y_{t}-\tilde{\phi}Y_{t-1}$$是一个$MA(1)$过程

对于ARMA模型，考虑EACF和对应的0锐角定点是一个不错的手段 考虑其中的零三角能看出 考虑$MA(1),AR(1or2)$ 都是可以接受的

ARIMA模型识别

模型的识别思路

我们在本文“非平稳随机时间序列ARIMA”一节中介绍了可以被ARIMA模型解释的非平稳性； 所有使用非平稳的序列计算出的ACF不表示任何东西（ACF计算本身假定平稳）

但是我们在大量的研究中发现了性质：非平稳序列的ACF呈现非指数下降的缓慢下降趋势

这一点是MA AR ARMA 模型都不具有的 可以认为是判断ARIMA模型的良好方法

在判断出ARIMA模型后 差分就是我们唯一的选择

关于过度差分

我们在本文“非平稳随机时间序列ARIMA”一节介绍了 平稳序列的差分依旧是平稳的；但是过度差分是不可取的

当我们多进行一次差分 就会需要多估计一次根本不存在的$\theta$值 他还会严重影响参数估计工作

模型的建立应当尽可能遵循简洁的原则，过度差分是不可取的 不过明显没有差分到位的情况 我们也要果断的差分

其他模型识别方法

我们有很多基于纯数值的方法 Dickey - Fuller 单位根检验 使用AIC BIC 等信息量准则选择最好的模型

参数估计

前面我们已经解决了模型的选择的问题 确定了模型的$p,d,q$ 现在我们确定模型中的那些参数了，我们只用考虑ARMA模型 ARIMA模型差分后转化为ARMA模型考虑 对于那些非零均值的ARMA模型 减去他们的均值后再进行分析

矩估计

矩估计并不为认为是最有效的方法 但是他是最简单的方法 还是样本矩等于理论矩的方法 它组成的方程组可以求任意未知参数的估计，这个方法最经典的例子是通过样本均值来估计总体均值

自回归模型AR

对于$AR(1)$模型 我们有

$$\rho_{1}=\phi.$$

所以我们计算样本自相关函数$r_k$有 $r_1=\phi$ 对于$AR(2)$模型

$$r_1=\phi_1+r_1\phi_2,\quad r_2=r_1\phi_1+\phi_2$$

更高阶的$AR$模型也是一样的原理 用Yule-Walker 模型计算就可以了

滑动平均模型MA

对于滑动平均模型 矩估计效果并不是很好 考虑$MA(1)$有

$$\rho_1=-\frac\theta{1+\theta^2}$$

代入我们的$r_1$ 我们是在处理一个二次方程

$$\hat{\theta}=\frac{-1+\sqrt{1-4r_1^2}}{2r_1}$$

这是在$|r_1|<0.5$的情况下才可以解出来的 否则这个方程没有实数解； 对于阶数更高的MA模型 会迅速的更加复杂化 我们只能考虑一些数值求解方法来

ARMA

我们还是考虑$ARMA(1,1)$的情况

$$\hat{\phi}=\frac{r_{2}}{r_{1}}$$ $$r_1=\frac{(1-\theta\hat{\phi})(\hat{\phi}-\theta)}{1-2\theta\hat{\phi}+\theta^2}$$

我们还是需要处理一个二次方程 在它存在解的时候给出 比较麻烦

噪声方差

我们最后一个需要估计的量是噪声方差$\sigma_e^2$ 首先我们知道可以用样本方差估计序列方差

$$s^{2}=\frac{1}{n-1}\sum_{\iota=1}^{n}{(Y_{t}-\overline{Y})^{2}}$$

然后我们使用前面在平稳模型中介绍过的关于方差的内容来估计噪声方差

$AR(p)$ $$\hat{\sigma}_{e}^{2}=(1-\hat{\phi}_{1}r_{1}-\hat{\phi}_{2}r_{2}-\cdots-\hat{\phi}_{p}r_{p})s^{2}$$ $MA(q)$ $$\hat{\sigma}_{\epsilon}^{2}=\frac{s^{2}}{1+\hat{\theta}_{1}^{2}+\hat{\theta}_{2}^{2}+\cdots+\hat{\theta}_{q}^{2}}$$ $ARMA(1,1)$ $$\hat{\sigma}_e^2=\frac{1-\hat{\phi}^2}{1-2\hat{\phi}\hat{\theta}+\hat{\theta}^2}s^2$$

总结

根据前面计算的结果 我们可以很容易的给出下面的结论

• 自回归模型的矩估计结果都是可以接受的
• 滑动平均模型的矩估计结果难以接受
• 混合模型的矩估计结果难以接受
• 也就是 MA 成分导致矩估计结果效果很差

最小二乘估计

此时 我们在平稳模型中引入一个均值$\mu$ 把它作为一个需要估计参数 从这点上看 最小二乘估计是一个还算不错的方法

自回归模型AR

考虑$AR(1)$的情况

$$Y_t-\mu=\phi(Y_{t-1}-\mu)+e_t$$

我们可以看成用$Y_t$作因变量 $Y_{t-1}$ 作为自变量的回归模型 最小二乘法研究的是偏差平方和的最小化 也就是

$$S_{\epsilon}(\phi,\mu)=\sum_{\iota=2}^{\pi}\bigl[(Y_{\iota}-\mu)-\phi(Y_{t-1}-\mu)\bigr]^{2}$$

的最小化 根据最小二乘法的相关计算方法 我们可以给出估计

$$\mu=\overline{Y}$$ $$\hat{\phi}=\frac{\sum_{t=2}^n{(Y_t-\overline{Y})(Y_{t-1}-\overline{Y})}}{\sum_{t=2}^n{(Y_{t-1}-\overline{Y})^2}}\quad.$$

他们并不是精确估计 但是就平稳过程而言 缺项导致的误差可以忽略 我们可以轻松的推广这些结果到更高阶的AR模型中 最小二乘法给出的结果在AR模型和矩估计区别不是很大

滑动平均模型MA

考虑最简单的$MA(1)$情况

$$Y_t=e_t-\theta e_{t-1}$$

这看起来并不是可以应用最小二乘法的样子 但是我们可以把MA模型表示成一个近似的自回归模型的性质（如果可逆） 如下

$$Y_t=-\theta Y_{t-1}-\theta^2Y_{t-2}-\theta^3Y_{t-3}-\cdots+e_t$$

由于其中需要求解的参数$\theta$存在非线性的存在 所以我们需要使用一些数值求解方法 针对更高阶的情况 迭代的求解是可行的

混合模型

考虑$ARMA(1,1)$

$$Y_t=\phi Y_{t-1}+e_t-\theta e_{t-1}$$

我们还是考虑残差平方和的最小化

$$e_t=Y_t-\phi Y_{t-1}+\theta e_{t-1}$$

此时 我们涉及到选取$e_i$的初始值 不过我们可以自由选取 在大样本情况下它对最终结果基本没有影响

极大似然估计

对于长度适中的序列和随机季节模型来说 初始值的选取对参数最后的估计结果影响很大 因此我们引入了最好用的估计方法———MLE

似然函数$L$定义为 为获取实际观测数据的概率密度函数 也被看作是当观测数据固定的时候 模型的未知参数的函数；对于$ARIMA$模型 给定观测值后的$L$是模型参数的函数 我们把它最大化的结果就是极大似然估计量

具体的估计方法这里就略去了 调用函数就可以解决问题了

模型诊断

现在 我们来考虑检验模型的拟合优度 分析拟合的残差 分析过度参数化的模型（比目标更一般的模型） 他们也是回归方法中诊断的常用思想（残差需要预测才能得到，不过我们把诊断拿到这里研究了）

残差分析

残差分析是涉及拟合分析问题中最广泛的一个分析方法 也是整个模型诊断中最基础的部分 我们在线性回归中就介绍过它 线性回归基础 的“回归诊断”一节 很多模型诊断方法都是基于残差进行的 这里我们介绍一些时间序列分析中的残差分析方法 残差的定义式还是以前的样子

$$\text{残差}=\text{实际值}-\text{预测值}$$

残差的时间序列图

一个理想的 没有任何模式的残差图 应该呈现出一个围绕零水平线 没有趋势的长方形散点图 如下 这就是一个基本理想的残差时间序列图 残差的时间序列图一般还会用于异常值的检验

残差的正态性检验

在模型拟合的过程中 我们假设了残差具有正态分析 此时进行残差的正态性检验非常的合适 比较常用的正态性检验方法有

• QQ图
• Shapiro-Wilk 正态性检验等非参数检验方法

残差的自相关性

我们在研究时间序列模型的时候 要求噪声项是独立的 这个噪声项在样本中被表示为残差项； 对于真正的白噪声（白噪声序列） 和 比较大的$n$ 样本自相关函数近似无关且有着0均值的正态分布； 不过 哪怕对于参数估计有效的正确识别的模型 残差也总是有着不同的特性； 一般的 残差近似服从0均值的正态怒分布 对于较小的滞后 方差远小于$\frac{1}{n}$ 对于较大的滞后 才可以用近似方差$\frac{1}{n}$

直观方法

现在我们需要结合具体的问题来研究了；我们可以绘制样本残差的ACF曲线 如果它满足小于临界虚线（用方差研究的） 那么就可以基本认为残差没有自相关性

需要特别注意的是，对于季节模型，滞后4 （季度数据）12（月度数据）的残差ACF很可能出现超过临界的情况，这点我们在季节模型的位置再详细的介绍

Ljung-Box检验

我们前面研究的是单独滞后的残差的自相关系数，事实上我们把这些系数放在一起研究也是很有意义的 类似于数理统计 的“复相关系数”一节中的研究 我们提出了统计量

$$Q=n(\hat{r}_{1}^{2}+\hat{r}_{2}^{2}+\cdots+\hat{r}_{k}^{2})$$

如果估计是正确的$ARMA(p,q)$模型 在大样本情况下，它近似服从$\chi^2(K-p-q)$

可惜 上面的这个结论对于较小的样本效果并不好 于是Ljung 和 Box提出了修正 对于典型的样本容量而言

$$Q_{\star}=n(n+2)\Big(\frac{\hat{r}_1^2}{n-1}+\frac{\hat{r}_2^2}{n-2}+\cdots+\frac{\hat{r}_k^2}{n-K}\Big)$$

比上面的统计量近似卡方分布的效果更好

在实际使用中，我们会根据自由参数的数量来降低自由度，调整Ljung-Box检验的参数

过度拟合和参数冗余

这个另一个在时间序列分析中比较重要的诊断方式；我们希望结局的问题是 在部分情况下 如$AR(2)$模型这种小参数量的模型效果已经很好了 但是我们在前面不小心拟合了如$AR(3)$这样的更多的参数的模型 并且进行了参数估计 现在就是我们修正这个错误的时候 在条件允许的情况下，我们需要尽可能的减少模型的复杂程度

在如下的情况下 我们认为存在了过度拟合

• 额外的参数不显著不为0
• 共有参数和原始估计相比没有什么变换 在部分的情况下 我们会参考模型的拟合效果量 如AIC值来辅助一些这里的判断

为了避免过度拟合的问题产生 我们也有一些在模型设计的阶段的建议

• 在条件允许的情况下尽可能使用简单的模型，如果拿不准，就同时考虑进行如残差分析之类的手段再处理
• 不要同时增加 AR MA I 三个部分的阶数

模型预测

事实上 预测才是时间序列建模的目的 时间序列分析根本就不像回归分析一样有系数可以分析 因此 进行预测 并且评估预测的精度（预测不确定性的度量）是整个时间序列分析的最后部分

最小均方误差预测

序列可以获得直到时间$t$的数据 也就是$Y_1,...Y_t$ 我们希望预测未来$l$期的数据 也就是$Y_{t+l}$ 我们称时间$t$为预测起点 $l$为预测步数或者预测前置时间 在我们省略的一些证明中 我们给出结论有 第$l$步的最小均方误差预测为

$$\hat{Y}_t(\ell)=E(Y_{t+\ell}|Y_1,Y_2,\cdotp\cdotp\cdotp,Y_t)$$

这是一个非常重要的结论 我们后面采取的所有预测都是最小MLE预测 这个条件期望将为我们简化很多问题

确定性趋势

这里我们的确定性趋势指的并不是观测代表着真实情况的含义 而是说我们已知模型内在机理 这里的例子希望我们能够理解我们预测的基本方法 考虑

$$Y_t=\mu_t+X_t$$

其中$X_t$是零均值的已知方差白噪声 则

$$\begin{aligned} \hat{Y}_{t}(\ell)& =E(\mu_{t+\ell}+X_{t+\ell}|Y_{1},Y_{2},\cdotp\cdotp\cdotp,Y_{\iota}) \\ &=E(\mu_{t+\ell}\mid Y_1,Y_2,\cdots,Y_t)+E(X_{t+\ell}\mid Y_1,Y_2,\cdots,Y_t) \\ &=\mu_{t+\ell}+E(X_{t+\ell}) \end{aligned}$$

也就是

$$\hat{Y}_{\iota}(\ell)=\mu_{\iota+\ell}$$

预测误差为

$$e_{t}(\ell)=Y_{t+\ell}-\hat{Y}_{t}(\ell)=\mu_{t+\ell}+X_{t+\ell}-\mu_{t+\ell}=X_{t+\ell}$$

研究预测误差有

$$E(e_{t}(\ell))=E(X_{t+\ell})=0$$ $$\mathrm{Var}(e_{t}(\ell))=\mathrm{Var}(X_{t+\ell})=\gamma_{0}$$

也就是预测是无偏的且方差为确定的数

ARIMA预测

这是承接我们的确定性趋势的 现在我们是从数据来推测模型 然后估计参数 最后进行预测 各种不同的模型有着不太一样的方法 我们还需要添加一些承接性的小节 后面的叙述将是一个整体 要注意整体的理解

AR模型

我们从非零均值的$AR(1)$模型开始 然后自然的推广理解整个AR模型 模型形式

$$Y_{t}-\mu=\phi(Y_{t-1}-\mu)+e_{t}$$

我们进行一步预测有

$$Y_{t+1}-\mu=\phi(Y_t-\mu)+e_{t+1}$$

根据最小MSE预测的要求 我们取条件期望有

$$\hat{Y}_{t}(1)-\mu=\phi[E(Y_{t}\mid Y_{1},Y_{2},\cdots,Y_{t})-\mu]+E(e_{t+1}\mid Y_{1},Y_{2},\cdots,Y_{t})$$

根据条件期望的性质我们可以实现化简

$$\hat{Y}_{t}(1)=\mu+\phi(Y_{t}-\mu)$$

能看出 一阶的AR模型本质上有压缩的现象（系数绝对值小于1） 最后模型会收敛到均值 想要预测更多的阶数指数迭代我们的预测结果就可以了

模型的预测误差也是非常好研究的 一步向前预测误差为

$$e_{t}(1)=Y_{t+1}-\hat{Y}_{t}(1)=\left[\phi(Y_{t}-\mu)+\mu+e_{t+1}\right]-\left[\phi(Y_{t}-\mu)+\mu\right]$$

也就是

$$e_t(1)=e_{t+1}$$

一步向前预测误差是AR噪声 我们也可以轻松的给出预测误差的方差（期望为0必然）

$$\mathrm{Var}(e_{t}(1))=\sigma_{e}^{2}$$

研究多步预测的误差有

$$e_t(\ell)=e_{t+\ell}+\psi_1e_{t+\ell-1}+\psi_2e_{t+\ell-2}+\cdots+\psi_{t-1}e_{t+1}$$

其中的$\psi$是原本系数变形的形式 这个形式对ARIMA模型都成立 容易判断 这个误差的期望为0 计算方差有

$$\mathrm{Var}(e_t(\ell))=\sigma_e^2(1+\psi_1^2+\psi_2^2+\cdotp\cdotp\cdotp+\psi_{\ell-1}^2)$$

这个形式也对所有的ARIMA模型成立

至于这个值的其他性质 我们直接给出结论有 对于所有的平稳的ARMA模型 有

$$\mathrm{Var}(e_{\iota}(\ell))\approx\mathrm{Var}(Y_{\iota})=\gamma_{0},\text{对较大的 }\ell$$

本小节 我们给出了三个不止局限于AR模型的结论

MA模型

我们这里考虑对于滑动平均成分应该怎么处理 考虑$MA(1)$ 有

$$Y_t=\mu+e_t-\theta e_{t-1}$$

进行一步预测 同时取条件期望有

$$\hat{Y}_{\iota}(1)=\mu-\theta E(e_{t}|Y_{t},Y_{2},\cdots,Y_{t})$$

然而我们知道（这是一个在$t$较大的时候成立的逼近结论）

$$E(e_t\mid Y_1,Y_2,\cdotp\cdotp\cdotp,Y_t)=e_t$$

因此 一步预测的表达式为

$$\hat{Y}_t(1)=\mu-\theta e_t$$

其中的$e_t$作为第$t$步的残差 在模型拟合的时候就被我们确定了

MA模型的多步预测稍微出现了一些变化 如下

$$\hat{Y}_t(\ell)=\mu+E(e_{t+\ell}\mid Y_1,Y_2,\cdotp\cdotp\cdotp,Y_t)-\theta e(e_{t+\ell-1}\mid Y_1,Y_2,\cdotp\cdotp\cdotp,Y_t)$$

我们发现 我们根本对超过$t$步的残差没有任何了解 根本没有真实值当然不用提残差 不过我们知道 残差$e_{t+l}$ 和$Y_i$无关 因此条件期望就是0 所以MA模型的多步预测为

$$\hat{Y}_{\iota}(\ell)=\mu,\ell>1$$

MA的预测误差我们不进行介绍，沿用AR给出的三个结论就好了

ARMA模型

我们直接给出预测的形式有

$$\begin{aligned}\hat{Y}_i(\ell)&=\phi_1\hat{Y}_i(\ell-1)+\phi_2\hat{Y}_i(\ell-2)+\cdots+\phi_p\hat{Y}_t(\ell-p)+\theta_0-\theta_1E(e_{t+t-1}\mid Y_1,Y_2,\cdot\cdot\cdot,Y_t)\\&-\theta_2E(e_{t+t-2}\mid Y_1,Y_2,\cdot\cdot\cdot,Y_t)-\cdot\cdot\cdot-\theta_eE(e_{t+t-q}\mid Y_1,Y_2,\cdot\cdot\cdot,Y_t)\end{aligned}$$

其中有

$$E(e_{t+j}\mid Y_1,Y_2,\cdots,Y_t)=\begin{cases}0&j>0\\e_{t+j}&j\leqslant0\end{cases}$$

这个关于误差的条件期望的取值情况是很合理的 在$j>0$的时候 我们根本没有残差可以使用 误差的期望当然是0 这意味着 在步数足够大的时候，噪声项的影响将会变得微弱，主要被自回归参数影响 综上所述 当预测步数小于$q$的时候 噪声项还可以直接参与我们的预测模型 否则 噪声项通过影响前面的自回归项取值情况来影响后面的模型

根据下面的公式

$$\begin{aligned}\hat{Y}_t(\ell)-\mu&=\phi_1[\hat{Y}_t(\ell-1)-\mu]+\phi_2[\hat{Y}_t(\ell-2)-\mu]+\cdots\\&+\phi_p[\hat{Y}_t(\ell-p)-\mu],\quad\ell>q\end{aligned}$$

类似于$ARMA$模型的$p_k$ 我们根据Yule-Walker递推可以告知结论 前面的式子将会指数衰减配合阻尼正弦快速衰减到0 也就是 平稳的ARMA模型长期预测会收敛于均值 这也回应了AR MA两个模型的相关性质

我们在这里继续复述前面已经给出的关于残差的结论有

$$\mathrm{Var}(e_{\iota}(\ell))\approx\mathrm{Var}(Y_{\iota})=\gamma_{0},\text{对较大的 }\ell$$

含有漂移的随机游动

为了处理ARIMA模型 我们先在这里进行一些必要的介绍 考虑模型

$$Y_t=Y_{t-1}+\theta_0+e_t$$

此时进行一步预测 进行条件期望有

$$\hat{Y}_\iota(1)=Y_\iota+\theta_0$$

不断差分就可以进行多步预测 能看出 如果$\theta_{0}\ne 0$ 那么无论多少步预测 都不会收敛 而是沿着直线前进

由于常数项会明显的改变预测的性质 因此非平稳的ARIMA模型在差分后研究的时候应该尽量避免常数项的存在 除非我们明显的发现 差分后序列的均值部位0

它的预测误差我们不单独研究了 后面在ARIMA模型中再进行介绍

ARIMA模型

ARIMA模型的预测并不是什么难以处理的问题 它的核心问题有两个

非稳定趋势

如果差分阶数不为0 并且ARMA模型均值不为0 则ARIMA模型则有差分阶数的变化趋势 而不是平稳的趋势

关于误差

直接给出关于误差的结论

$$E(e_{t}(\ell))=0,\ell\geqslant1$$ $$\mathrm{Var}(e_t(\ell))=\sigma_e^2\sum_{j=0}^{\ell-1}\psi_i^2,\ell\geqslant1$$

其中后者是一个不收敛的级数 也就是 非平稳的预测误差的方差会不断增大，没有上线 这点非常的合理 毕竟非平稳序列的未来 相当难以确定

变换序列后的预测

差分

我们有着非常自然的想法 预测差分后的平稳序列，然后加总得到原序列的值 事实上这种方法的效果非常的不错

如Box Cox变换之类的变换

一种非常自然的思想是 我们对变换后的序列建模 然后把预测的结果逆变换回去 非常遗憾的是

$$E(Y_{t+t}|Y_{t},Y_{t-1},\cdotp\cdotp\cdotp,Y_{1})\geqslant\exp[E(Z_{t+t}|Z_{t},Z_{t-1},\cdotp\cdotp\cdotp,Z_{l})]$$

也就是上面的思路 无法保证最小的MSE

好在我们需要额外进行的工作并不多 根据矩母函数相关结论有 如果$X$服从正态分布 则

$$E[\exp(X)]=\exp[\mu+\frac{\sigma^2}{2}]$$

因此原始序列的最小MSE预测为

$$\exp\left\{\hat{Z}_{t}(\ell)+\frac{1}{2}\mathrm{Var}[e_{t}(\ell)]\right\}$$

其中后者是预测的误差的方差 不太要求最小MSE的情况直接用前面的直接变换思路就可以了

单位根过程

典型的非平稳时间序列模型是单位根(unit root)非平稳时间序列

随机游走

参见本文“随机游走”一节。

含有漂移项的随机游走

我们考虑随机游走模型的变形 有

$$p_t=\mu+p_{t-1}+\varepsilon_t,t=1,2,\ldots$$

考虑我们前面考虑的特征有

• 方差不发生变化
• 均值增加了一个趋势项 模型依旧不可预测，理论上在直线附近摆动的序列因为过大的方差而失去规律性（预测的MSE接近无穷）

带漂移的随机游动$p_t$,可以分解为两部分：

$$p_t=(p_0+\mu t)+p_t^*.$$

其中$p_t^*=\sum_{j=1}^t\varepsilon_t$是从0出发的不带漂移的随机游动，$p_0+\mu t$是一个非随机的线性趋势。

固定趋势模型

固定趋势模型的介绍可以参考 本文“趋势”一节

固定趋势模型与随机游走的联系

随机游动$p_t=p_{t+1}+\varepsilon_t$与固定趋势加扰动 $Y_t=a+bt+X_t$(其中$\{X_t\}$平稳)都能呈现出缓慢的趋势变化。

区别在于：

• 随机游动的方差是线性增长的，固定趋势的观测值方差不变；
• 随机游动的扰动的影响是永久的，固定趋势的扰动的影响仅在一个时刻(如果扰动$X_t$是白噪声)或者很短时间 (如果是扰动$X_t$是线性时间序列) ;
• 随机游动的趋势没有固定方向，固定趋势的变化形状是固定的；
• 固定趋势模型$Y_t$减去一个固定的回归函数$Y=a+bt$就可以变成平稳列，随机游动减去任意的非随机函数都不能变平稳，可以用差分运算变成平稳。

ARIMA模型

ARMA模型的叙述可以参考 本文“平稳随机时间序列ARMA”一节 ARIMA模型只是在ARMA模型的基础上增加差分阶数 参考 本文“非平稳随机时间序列ARIMA”一节

如果$Y_t$本身已经是弱平稳列，则不应对$Y_{t}$进行差分。（非常自然） 如果$Y_t$是非随机的线性趋势加平稳列，虽然差分能将其变成平稳列，但是也不应该使用差分来做而是应该用回归来做，用差分来做会在ARMA模型的MA部分引入不必要的单位根。

指数平滑模型

指数平滑最早是来自一种简单的预测方法： 用历史数据的线性组合预测下一时间点的值， 线性组合系数随距离变远而按负指数（几何级数）衰减

$$\hat{x}_h(1)\approx wx_h+w^2x_{h-1}+\cdots=\sum_{j=1}^\infty w^jx_{h+1-j}$$

加权平均需要满足权重和为1 则

$$\hat{x}_h(1)=(1-w)(x_h+wx_{h-1}+w^2x_{h-2}+\ldots)=(1-w)\sum_{j=0}^\infty w^jx_{h-j}$$

能发现 这和ARIMA(0,1,1) 有着一样的形式 因此对于指数平滑模型 我们可以用ARIMA模型的方法来研究指数平滑模型 不需要浪费时间单独进行建模

单位根检验

我们可以使用单位根检验来判断一个过程是不是单位根过程；它假设一个过程是单位根过程（有单位根） 不拒绝原假设则可以认为是单位根过程

顺利判断单位根过程可以认定模型非平稳

我们可以选定单位根过程的趋势形式 它将在拟合趋势后再判断剩余残差是否支持非平稳

季节模型

季节数据 或者说存在周期性的数据在时间序列分析中应该也是很常见的；我们前面介绍的模型无法解释这些数据

我们会发现 简单拟合后的残差在许多滞后还是会出现高度的自相关（ARIMA当然不能接受这些）,我们在下面介绍的随机季节模型却能够很好的拟合这样的序列

季节 ARMA 模型

季节MA模型

我们还是从平稳模型的情况开始 用$s$表示已知的季节周期 月度数据一般$s=12$ 季度数据一般$s=4$ 考虑下面形式的模型

$$Y_t=e_t-\Theta e_{t-12}$$

我们容易验证：

$$\mathrm{Cov}(Y_{t},Y_{t-1})=\mathrm{Cov}(e_{t}-\Theta e_{t-12},e_{t-1}-\Theta e_{t-13})=0$$ $$\mathrm{Cov}(Y_{t},Y_{t-12})=\mathrm{Cov}(e_{t}-\Theta e_{t-12},e_{t-12}-\Theta e_{t-24})=-\Theta\sigma_{e}^{2}$$

容易看出 这个序列是平稳的 并且只有12阶滞后才有自相关性

根据上面的想法 我们定义季节周期为$s$的$Q$阶MA模型 $MA(Q)$

$$Y_t=e_t-\Theta_1e_{t-s}-\Theta_2e_{t-2s}-\cdots-\Theta_Qe_{t-Qs}$$

关于可逆性的条件和前面还是一样的 它的自相关函数还是和前面的一样 在若干阶后变为0

季节AR模型

非常自然的可以定义季节AR模型 我们取$AR(P)$

$$Y_{\iota}=\Phi_{1}Y_{t-s}+\Phi_{2}Y_{t-2s}+\cdots+\Phi_{P}Y_{t-Ps}+e_{t}$$

我们直接给出结论

• 关于平稳性的条件不变
• 自相关函数是指数衰减和阻尼正弦的组合

乘法季节ARMA模型

考虑前面那些只有季节滞后上包含自相关性的模型是没有价值的，分析起来也是和ARIMA完全一样 其实没有什么意思；

现在我们希望结合季节模型ARIMA和前面研究ARIMA上的思考 研究那些不仅在季节滞后上包含相关性 还在临近上包含相关性的模型

我们给出两个例子

$$Y_t=e_t-\theta e_{t-1}-\Theta e_{t-12}+\theta\Theta e_{t-13}$$ $$Y_{t}=\Phi Y_{t-12}+e_{t}-\theta e_{t-1}$$

同时给出两个典型的ACF曲线 他们都是典型的 符合我们前面介绍的模型

我们把此时我们介绍的模型记为

$$ARMA(p,q)\times(P,Q)$$

它的含义是

• 非季节部分 含有$ARMA(p,q)$
• 季节部分本身 符合$ARMA(P,Q)$

非平稳季节ARIMA模型

差分 还是研究非平稳的核心步骤 我们有周期为$s$的季节差分为

$$\nabla_{s}Y_{t}=Y_{t}-Y_{t-s}$$

再结合本身模型的差分 就可以得到乘法季节ARIMA模型为 它有着

• 季节周期$s$
• 非季节的阶数$ARIMA(p,d,q)$
• 季节的阶数$ARIMA(P,D,Q)$ 它是一大类模型

季节ARIMA模型识别 拟合 检验 预测

核心的方法我们已经在本文“模型识别”一节本文“参数估计”一节本文“模型诊断”一节中介绍过了 这里针对季节模型进行一些单独的介绍

季节ARIMA模型识别

• 研究时间序列图 ACF PACF 判断平稳性和周期性
• 考虑普通差分 尝试捕捉平稳性
• 考虑季节差分 尝试消除周期性
• 研究样本ACF 是否消除自相关（前面的差分是在消除所有自相关）

季节ARIMA模型经过一阶差分和季节差分后的序列常常仍在滞后1、4、5这些位置呈现出自相关 （如果是月度数据，则为滞后1、12、13）要综合ACF PACF 时间序列图判断是否已经消除了平稳与周期性

季节ARIMA模型拟合

我们首选MLE的拟合方式 代码可以很好的帮助我们计算这些结果了

季节ARIMA模型诊断

还是研究残差 具体方法没有任何变换

季节ARIMA模型预测

同预期的一样 季节模型的最好的预测方法就是运用差分进行预测 考虑 $ARIMA(0,1,1)\times(1,0,1)_{12}$ 有（想理解差分，最好从差分后的ARMA入手）

$$Y_t-Y_{t-1}=\Phi(Y_{t-12}-Y_{t-13})+e_t-\theta e_{t-1}-\Theta e_{t-12}+\theta\Theta e_{t-13}$$

一步向前预测就有

$$\hat{Y}_t(1)=Y_t+\Phi Y_{t-11}-\Phi Y_{t-12}-\theta e_t-\Theta e_{t-11}+\theta\Theta e_{t-12}$$

更多的步数同理 我们还是需要考虑噪声项有时候以残差形式纳入预测 有时候以自回归的形式纳入的问题 季节的ARIMA模型含有两个部分，一个是前置时间的趋势，一个是周期性部分的趋势，他们的和，就是我们的季节ARIMA

季节虚拟变量

另一种表示季节性的方法是用非随机的回归项表示固定的季节模式。 这样的模式虽然也可以通过季节差分消除， 但是与动态模型和非随机线性趋势模型的关系类似， 固定的季节模型不应该用季节差分处理。 如果能找到确定的趋势 减除趋势就可以平稳化，那么就不需要采用差分处理

非随机的季节因素用回归哑变量表示。$s=4$时，用3个哑变量就可以表示4个不同季节的固定水平。为了判断非随机季节模型是否使用，可以先拟合动态的季节ARIMA$(1,0,1)(1,0,1)_s$模型，当发现其中的季节因素可以忽略时，就可考虑采用非随机的季节模型。下面举例说明。

给出序列图 不能观测出明显的周期性 作ACF图 12阶滞后明显不为0 体现了周期性 拟合动态的季节ARIMA$(1,0,1)(1,0,1)_{12}$ 发现 sar1 = 0.9882,sma1 = -0.9142 写成实际模型

$$(1+0.0639B)(1-0.9882B^{12})(X_t-0.0117)=(1+0.2508B)(1-0.9142B^{12})\varepsilon_t$$

他们可以近似约去 这意味着可以考虑采用季节 dummy 变量的回归模型 这本质上就是在考虑一个回归问题了 我们只是在回归从而计算趋势项 时间序列的特征此时不考虑了

这样的操作看似非常的可信，但是实际上，时间序列数据存在很强的自相关性，我们对回归的残差分析就会发现这一点，在时间序列中直接考察趋势进行回归就结束分析是草率且不负责的

含有时间序列误差的回归模型

在统计学的数据分析中， 线性回归分析是最常用的分析工具之一 一元线性回归如下

$$Y_t=\beta_0+\beta_1X_t+e_t,t=1,2,\ldots,T$$

我们在回归模型的时候要求残差项$e_t$ 独立同正态分布

但是在金融涉及的回归分析中 时间序列是频繁存在的 包括我们在去趋势的时候采用的回归分析 他们都在存在并不独立的残差项 此时基于残差项$e_t$ 独立同正态分布作出的判断 比如估计量的标准误差估计、假设检验都不再成立 回归系数的估计依旧是可信的

有研究指出 当残差项存在正相关的时候 回归系数的标准误差估计偏低， 使得相应的$t$和$F$检验统计量的绝对值偏大

当$\{e_t\}$是平稳可逆ARMA序列时，可以将线性回归模型与平稳可逆ARMA序列同时估计，可以得到需要标准误差估计、假设检验和预测。 arima()函数提供了一个 xreg= 用来引入回归自变量。

如果不关心$\{e_t\}$的具体模型，而只关心对回归系数的SE的正确估计以及假设检验的正确性，可以仅假设$\{e_t\}$的协方差结构而不考虑$\{e_t\}$的建模。如线性回归基础 的“广义最小二乘估计（加权OLS）”一节

建模的基本步骤为

1. 拟合一个线性回归模型，并检验残差的序列相关性
2. 如果残差序列是单位根非平稳的，则对因变量和自变量都做一阶差分。 然后对差分后的序列再进行第一步。 如果残差序列是平稳的， 则对残差序列识别一个ARMA模型， 并相应地修改线性回归模型。
3. 用最大似然估计法对回归模型与ARMA模型进行联合估计， 并对模型进行检验， 看是否需要改进。 主要可以使用Ljung-Box对残差进行白噪声检验。

长记忆模型

ACF是时间序列建模的重要参考。

• 对于ARMA序列，当滞后$k\to\infty$时其样本ACF是负指数速度趋于零的。
• 对于单位根非平稳列，其理论ACF无定义(因为自协方差是针对弱平稳列定义的),其样本ACF在样本量$T\to\infty$时每个$\hat{\rho}_k$都趋于1 $( k> 0)$ 。

有一些平稳时间序列的ACF（PACF）虽然也随滞后$k\to\infty$趋于零，但是收敛到零的速度比较慢，只有负幂次$k^{-\alpha}$ 这样的速度。这代表着序列的自相关性随着距离变远而减小得比较慢，称这样的序列是长记忆时间序列。

注意，长记忆时间序列仍是弱平稳的，单位根非平稳列虽然远距离的自相关性很强但不称为长记忆。

在金融时间序列建模中， 如果样本ACF数值不大但是衰减特别缓慢， 可以考虑长记忆模型。 如果数值很大同时衰减慢则可能是单位根非平稳， 或者具有很接近1的特征根的ARMA序列。

长记忆时间序列的典型模型是分数差分弱平稳列，模型为

$$(1-B)^dX_t=\xi_t,\:-0.5<d<0.5$$

其中$\{\xi_t\}$是零均值独立同分布白噪声

如果我们讲这个白噪声推广到可以取$ARMA(p,q)$序列 这样的模型称为$ARFIMA(p,d,q)$ 模型 称为分数阶差分ARMA模型

突变导致的非平稳性

非平稳性产生的第二种原因是总体回归函数在样本期内发生了变化。在经济学中， 很多原因会导致上述现象的出现，如经济政策的变更、经济结构的变化、创新导致某个行业突然的变化。如果这些变化或“突变”发生，但在模型中没有考虑到这些因素，就会动摇我们进行预测和推断的基础。

本节介绍了两个检测时间序列回归模型中突变的思路。

第一个思路是从假设检验的角度来寻找可能的突变点并通过$F$ 统计量来检验回归系数变化的显著性。

第二个思路是从预测的角度来寻找可能的突变点：假设样本在样本期间实际结束前就已经结束了，在此基础上进行预测，并对该预测结果进行评价。如果预测结果的精度大幅下降，就推断出现了突变。

我们前面仅仅介绍了趋势产生的非平稳性，而没有考虑突变的问题。

什么是突变

突变可能产生于总体回归函数在某一个特定时间的突然变动，也可能产生于长期的演变过程之中。

宏观经济数据的突然变动可能产生于宏观政策的大规模变动。

突变也可能产生于总体回归函数随时间的演变过程中，例如经济政策的缓慢改革和经济结构的逐步变革。

本节所介绍的检测突变的方法既可以用来检验突然的突变，也可以用来检验长时间演变中所造成的突变。

突变的检验

探测突变的方法之一是检验回归系数的离散变化或突变。具体应该如何来检验则取决于突变点发生的时间是否可知。

已知时间

对于已知时间突变点的检测。有些情况下，你可能怀疑在某一个特定的时间点发生了突变。

如果可能发生突变的日期是已知的，则可以运用二元变量交叉回归模型来对零假设进行检验，此时零假设应为没有突变。为了简单起见，我们下面考虑 ADL(1,1)模型，模型中的解释变量包括截距项、$Y_t$的一阶滞后、 $X_{\iota}$的一阶滞后。用$\tau$代表假设的突变点发生日期，$D_\iota(\tau)$ 为二元变量，在突变期之前其取值为 0,在突变期之后，其取值为 1,即当$t\leqslant\tau$时，$D_t(\tau)=0$,当$t>\tau$时， $D_{\iota}(\tau)=1$。包含二元突变指示变量和所有交叉项的回归方程为：

$$Y_t=\beta_0+\beta_1Y_{t-1}+\delta_1X_{t-1}+\gamma_0D_t(\tau)+\gamma_1\begin{bmatrix}D_t(\tau)\times Y_{t-1}\end{bmatrix}+\gamma_2\begin{bmatrix}D_t(\tau)\times X_{t-1}\end{bmatrix}+u_t$$

如果样本期内没有突变发生，则总体回归函数在两个阶段应该是相同的，即所有包含着$D_t(\tau)$的项的系数都应为 0。也就是说，零假设应为样本期内无突变，即$\gamma_0=\gamma_1=$ $\gamma_2=0$。备择假设为样本期内存在突变，这意味着总体回归函数在突变点$\tau$前后有所不同，即$\gamma_0$、$\gamma_1$、$\gamma_2$至少有一个不为零。所以样本期内是否含有突变，可以通过$F$统计量进行检验

未知时间

如果我们未知具体突变发生的具体时间点，那么可以考虑进行多次已知时间的突变检验，幸运的是，已经有研究研究了其整合的形式

这种被改进了的邹氏检验通常被称作匡特似然比检验Quandt Likelihood Ratio(QLR) Statistic(在下文我们将运用该术语表示此检验),或者称为 sup-Wald 统计量。

伪样本外预测检验

对于模型预测精度的检验最终还是要看其样本外预测的能力，即要看在模型估计完成后，其在“实际预测区间”上的预测能力。

伪样本外预测(Pseudo Out-of-sample Forecasting) 是一个用来模拟预测模型在实际预测区间预测表现的方法。它的思路很简单：选择样本区间末端的一个时间点，运用该时间点以前的数据对模型进行估计，然后运用估计出来的模型对样本末端的观测值进行预测。在样本期末端的多个时间点上重复上述步骤可以得到多个伪预测值和多个伪预测误差。这些误差可以用来检验在预测关系在稳定假设情况下是否让我们满意

这种方式也可以帮助我们判断是否有突变这种非平稳性产生

突变的处理

突变的处理是一个非常复杂的问题，我们这里不介绍如何处理突变