时间序列与R语言应用(part2)--ADF单位根检验

学习笔记，仅供参考，有错必纠
博客阅读索引：博客阅读及知识获取指南

参考文献：《计量经济学》、《时间序列分析及应用R语言》
上一篇：时间序列基本概念
下一篇：趋势平稳与差分平稳

---

---

ADF单位根检验

学习单位根检验之前，先来几道开胃菜，看几个知识点。

开胃菜

自回归过程AR

顾名思义，自回归过程是用自身做回归变量，具体来说，$p$阶自回归过程{$Y_t$}满足方程：
$$Y_t={\varphi }_1{Y}_{t-1}+{\varphi }_2{Y}_{t-2}+...+{\varphi }_p{Y}_{t-p}+e_t\text{\hspace{1em}(1)}$$
序列$Y_t$的当期值是自身$p$阶滞后项和新息项$e_t$的线性组合，其中$e_t$包括了序列在$t$期无法用过去值来解释的所有新信息，因此，对于每一个$t$,假设$e_t$独立于$Y_{t-1},Y_{t-2}{Y}_{t-3}...$

特征根

为了学习特征根，我们以二阶自回归过程为例。

我们考虑以下二阶自回归过程：
$$Y_t={\varphi }_1{Y}_{t-1}+{\varphi }_2{Y}_{t-2}+e_t\text{\hspace{1em}(2)}$$
为了讨论平稳性，引入AR特征多项式：
$$\varphi (x)=1-{\varphi }_{1}x-{\varphi }_2{x}^2\text{\hspace{1em}(3)}$$
和相应的AR特征方程:
$$1-{\varphi }_{1}x-{\varphi }_2{x}^2=0\text{\hspace{1em}(4)}$$
二次方程总是有两个根(可能有复数根)，这两个根就叫特征根。

可以证明，当$e_t$独立于$Y_{t-1},Y_{t-2}{Y}_{t-3}...$的条件下，当且仅当AR特征方程的根的绝对值(模)大于1时，方程(2)才存在平稳解。

当方程(4)存在等于1的特征根时，我们称这个特征根为单位根，此时该二阶自回归过程不满足平稳条件。

趋势

• 随机趋势

首先我们观察下图，它是否有趋势？

可以看出，该序列有一个普通的上升趋势，随着时间的增加，观测值也增大。但是我们并不能说，这个序列存在一个确定性的因素使其观测值不断增大，因为这是一个随机游动过程。我们知道，随机游动过程在任何时间上都有零均值，在相邻时点上，序列值存在很强的正相关，并且随机过程的方差随时间的增加而增加，它表面上趋势是由此得到的产物，有些作者把这样的趋势称为随机趋势。

• 确定性趋势

若有如下公式：
$$Y_t={\mu }_t+X_t\text{\hspace{1em}(5)}$$
我们可以假设$X_t$为对所有$t$都零均值的扰动，${\mu }_t$是观测序列$Y_t$的均值函数，与随机趋势相比，我们可以把这种模型描述成具有确定趋势。我们可以假设确定性趋势关于时间是线性的，即${\mu }_t={\beta }_0+{\beta }_{1}t$,或者具有二次时间趋势，即${\mu }_t={\beta }_0+{\beta }_{1}t+{\beta }_2{t}^2$

非平稳时间序列

具有时变均值的任何时间序列都是非平稳的，形如：
$$Y_t={\mu }_t+X_t\text{\hspace{1em}(6)}$$
的模型，其中${\mu }_t$是时变均值函数，$X_t$为零均值平稳序列，可以说这个序列是趋势非平稳的。

当我们观察一个时序图时，仅因序列的一部分看似线性递增(或递减)，不足以确信过程具有内在的线性特征，并长期保持下去。比如，下面这个原油月度价格图：

我们看到这个序列图，变动较大，由此可以判断该过程不是平稳的。貌似也没有适合此序列的确定性趋势模型，故一个含有随机趋势的非平稳模型才是适当的。

利用ADF检验区分趋势非平稳和差分非平稳

需要注意的是，趋势非平稳和差分非平稳是不同的，趋势非平稳在某种意义上指，序列具有一个确定性趋势(比如，线性趋势)，否则则是平稳的。我们可以利用单位根检验来区分趋势非平稳和差分非平稳。比如，可以通过对去趋势后的序列进行ADF检验。

单位根检验

DF检验

我们己知道，随机游走序列：
$$X_t=X_{t-1}+{\mu }_t\text{\hspace{1em}(7)}$$

是非平稳的,其中${\mu }_t$是白噪声.而该序列可看成是随机模型：
$$X_t=\rho X_{t-1}+{\mu }_t\text{\hspace{1em}(8)}$$
当参数$\rho =1$时的情形。

也就是说，对(8)式做回归，如果发现$\rho =1$,则称随机变量$X_t$有一个单位根。显然，一个有单位根的时间序列就是随机游走序列，而随机游走序列是非平稳的。因此，要判断某时间序列是否是平稳的，可判断(8)式是否有单位根。

(8)式可变形成差分形式：
$$\begin{gathered} \Delta X_t=(\rho -1)X_{t-1}+{\mu }_t \\ =\delta X_{t-1}+{\mu }_t\text{\hspace{1em}(9)} \end{gathered}$$
检验(8)式是否存在单位根，也可通过(9)式判断是否存在$\delta =0$

一般地，检验一个时间序列$X_t$的平稳性，可通过检验带有截距项的一阶自回归模型：
$$X_t=\alpha +\rho X_{t-1}+{\mu }_t\text{\hspace{1em}(10)}$$
中的参数$\rho$是否小于1，或者说检验其等价变形式:
$$\Delta X_t=\alpha +\delta X_{t-1}+{\mu }_t\text{\hspace{1em}(11)}$$
中的参数$\delta$是否小于0.

若(8)式中$\rho$大于或等于1，则序列非平稳,证明如下:

有(8)式的特征多项式为：
$$\rho (x)=1-\rho x$$
则相应的特征方程为：
$$1-\rho x=0$$
其唯一特征根为：
$$x=\frac{1}{\rho }$$
由于$\rho$大于等于1，则特征根小于等于1，则方程(8)不存在平稳解，序列非平稳。

故，针对(11)式，我们假设检验的原假设为$H_0:\delta =0$,备择假设为$H_1=\delta <0$。这里可通过普通最小二乘法下的t检验完成。

在零假设(序列非平稳)下，即使在大样本下统计量也是有偏误的，通常t检验无法使用。Fuller于1976年提出了这一情形下的t统计量服从的分布，即DF分布，DF分布临界值如下表所示：

ADF检验

在DF检验中，事实上假定了时间序列是由具有白噪声随机干扰项的一阶自回归过程生成的。但实际检验中，可能会出现以下3点情况：

①时间序列可能由更高阶的自回归过程生成，比如，时间序列存在滞后4阶自相关，即$Cov(X_t,X_{t-4})$显著不为0。

②随机干扰项并非是白噪声

比如，考虑模型：
$$X_t=\rho X_{t-1}+{\mu }_t\text{\hspace{1em}(12)}$$
若随机干扰项{${\mu }_t$}是m阶自回归过程:${\mu }_t={\beta }_1{\mu }_{t-1}+...+{\beta }_m{\mu }_{t-m}+{ϵ}_t$,在零假设$\rho =1$,${\mu }_t=X_t-X_{t-1}$的情况下，令$\delta =\rho -1$,则有：
$$X_t-X_{t-1}=\delta X_{t-1}+{\beta }_1(X_{t-1}-X_{t-2})+...+{\beta }_m(X_{t-m}-X_{t-m-1})+{ϵ}_t\text{\hspace{1em}(13)}$$
则差分形式为：
$$\Delta X_t=\delta X_{t-1}+\sum _{i=1}^m{\beta }_i\Delta X_{t-i}+{ϵ}_t\text{\hspace{1em}(14)}$$
③时间序列包含有明显的随时间变化的某种趋势，比如：线性趋势。

为了保证DF检验中随机干扰项的白噪声特性Fuller提出了增强的DF检验，即ADF检验。ADF检验是通过如下三个模型完成的(笔者不想打公式了！直接上图片)：

模型3中的t是时间变量，代表了时间序列随时间变化的某种趋势。这3个模型的原假设都是$H_0:\delta =0$,即存在1个单位根,模型1与另外两个模型的差别在于是否包含常数项和趋势项。

检验从模型3开始，然后模型2,最后是模型1。何时检验拒绝零假设，即原序列不存在单位根，为平稳序列，何时可停止检验。否则，就要继续检验。

一个简单的检验是同时估计出上述3个模型的适当形式，然后通过ADF临界值表检验是否拒绝原假设($H_0:\delta =0$).

R语言实现

• 模拟案例

按照定义模拟随机游走：

set.seed(1238)
n <- 1000
y <- vector(length = n)
y[1] = 0
for (i in 2:n){
  y[i] = y[i-1] + rnorm(1, 0, 1/n)
}
plot(y, type = 'l', xlab = 'time', main = '随机游走模拟')

图像：

ADF检验随机游走：

adf.test(y)

控制台输出：

	Augmented Dickey-Fuller Test

data:  y
Dickey-Fuller = -1.2817, Lag order = 9, p-value = 0.8824
alternative hypothesis: stationary

我们看到p值大于0.05的显著性水平，不能拒绝原假设，则序列是非平稳的。

向序列增加常数项和趋势，构造新序列new_y，并进行检验：

t = c(1:n)
new_y = y + 3 + 2*t
adf.test(new_y)

控制台输出：

	Augmented Dickey-Fuller Test

data:  new_y
Dickey-Fuller = -1.2817, Lag order = 9, p-value = 0.8824
alternative hypothesis: stationary

可以看到新序列的p值和随机游走序列的p值相同。

• 看源码

下面我们看一下部分源码，证明一下上面阐述的ADF理论部分内容。

adf.test源码1：

    table <- cbind(c(4.38, 4.15, 4.04, 3.99, 3.98, 3.96), c(3.95, 
        3.8, 3.73, 3.69, 3.68, 3.66), c(3.6, 3.5, 3.45, 3.43, 
        3.42, 3.41), c(3.24, 3.18, 3.15, 3.13, 3.13, 3.12), c(1.14, 
        1.19, 1.22, 1.23, 1.24, 1.25), c(0.8, 0.87, 0.9, 0.92, 
        0.93, 0.94), c(0.5, 0.58, 0.62, 0.64, 0.65, 0.66), c(0.15, 
        0.24, 0.28, 0.31, 0.32, 0.33))
    table <- -table
    tablen <- dim(table)[2]
    tableT <- c(25, 50, 100, 250, 500, 1e+05)
    tablep <- c(0.01, 0.025, 0.05, 0.1, 0.9, 0.95, 0.975, 0.99)

可以看到源码里有这么一个矩阵(table的类型是矩阵)，那这个6行8列的矩阵是啥呢？回顾一下DF检验，感觉这个表是不是和DF分布临界值表有点类似的意味？对呢，这可能就是Fuller于1996提出的临界值表，也可能是经过了后人改进的临界值表？这里我没有验证，欢迎大家告知[求教].JPG

adf.test源码2：

    k <- k + 1
    x <- as.vector(x, mode = "double")
    y <- diff(x)
    n <- length(y)
    z <- embed(y, k)
    yt <- z[, 1]
    xt1 <- x[k:n]
    tt <- k:n
    if (k > 1) {
        yt1 <- z[, 2:k]
        res <- lm(yt ~ xt1 + 1 + tt + yt1)
    }
    else res <- lm(yt ~ xt1 + 1 + tt)
    res.sum <- summary(res)

没错，这就是ADF检验中，模型3的回归方程啦。按照类比的思想，yt代表$\Delta X_t$，xt1表示$X_{t-1}$，1表示截距项(lm中的1代表包含截距项)，tt表示$t$， yt1表示包含$\Delta X_{t-1},...,\Delta X_{t-k}$的矩阵。

同时，我们看到if的判断条件为k>1,当k>1时，则执行包含$\Delta X_{t-1},...,\Delta X_{t-k}$矩阵的回归代码，否则执行不包含$\Delta X_{t-1},...,\Delta X_{t-k}$矩阵的回归代码，那么这里的k是啥意思？这个k其实是我们调用adf.test函数时需要自定义的参数，它代表着残差项{${\mu }_t$}有k阶自回归，当残差项不存在自回归时，我们设置k=0,存在一阶自回归则设置k=1，以此类推。那么if不是应该按照k>0来判断吗？NO!仔细看源码，可以看到第1行代码为k <- k+1，则if里的确应该用k>1来做判断。

---

自学结果，有错求教，谢谢