做数据序列的惩罚最大F检验

Question

数据序列 惩罚最大 F检验( the penalized maximal F test，PMFT)，帮个忙，十分感谢。

Answer 1

单位根过程是特征方程含有单位根的数据序列，如随机游走模型就是一个单位根过程，它的特征方程为，其根为。检验数据序列是否存在单位根的方法是DF检验。

1 随机游走过程的自相关系数
1.1 理论推导
若，其中，则与的相关系数绝对绝对值为1，这是一个很自然的推论，但在时间序列分析中却并非如此。

对于：

当时，序列为AR(1)平稳序列，其自相关系数为，可知与的（自）相关系数即，其绝对值小于1；

而当时，序列为随机游走过程，那么其自相关系数是否为呢？下面就来简单推导一下随机游走模型的自相关系数表达式。

假设其初始值为，使用累积法可得，

方差，

同理，

协方差

因此，自相关系数

这说明随机游走模型的自相关系数的数学期望并非是1，而是随而出现衰减趋势，这与接近于1的AR平稳过程难以区分，因此无法通过ACF图象来判断序列是否存在单位根。

另外，随机游走模型的方差会随t无限增大，不符合等方差的假设，使用OLS估计会得出有偏的结果，因此t检验也不能用于单位根存在的判别方法。

1.2 示例
下面生成两个序列：

两个序列的初始值均设定为0，并使用同一个{}序列。

set.seed(2231)
epsilon = rnorm(200)
y1 <- y2  <- NULL
y1[1] = y2[1] = 0
for(i in c(2:200)) {
  y1[i] = 0.9*y1[i-1] + epsilon[i]   
  y2[i] = y2[i-1] + epsilon[i]   
}
y1 = y1[51:200]
y2 = y2[51:200]

对比两个序列的走势：

par(mfrow = c(2,1), family = "mono",
    plt = c(0.1,0.9,0.12,0.9),
    mgp = c(2,0.5,0))
plot(y1, type = "l")
plot(y2, type = "l")

对比两个序列的ACF图象：

par(mfrow = c(2,1), family = "mono",
    plt = c(0.1,0.9,0.12,0.9),
    mgp = c(2,0.5,0))
acf(y1, ylab = "y1 ACF")
acf(y2, ylab = "y2 ACF")

从走势图来看，两个序列没有明显的区别；而从ACF图象来看，二者的区别在于随机游走过程的ACF衰减速度更慢，但在实际中很难判断“慢”到何种程度才是平稳过程向单位根过程过渡的临界状态，不过当ACF图象衰减很慢时能提示数据序列可能是单位根过程。

下面使用OLS回归拟合和之间的关系：

model.1 <- lm(y2[2:150] ~ y2[1:149])
model.2 <- lm(y2[2:150] ~ y2[1:149] - 1)
summary(model.1)
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.28865    0.16940   1.704   0.0905 .  
## y2[1:149]    0.94812    0.02521  37.611   <2e-16 ***
summary(model.2)
## Coefficients:
##           Estimate Std. Error t value Pr(>|t|)    
## y2[1:149]  0.98592    0.01205   81.83   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

-假设OLS估计是无偏的，若序列存在单位根，则模型的一次项系数应该不显著区别于1：
model.1包含截距。一次项系数拟合结果为0.94918，标准误为0.02645，，表明其在95%置信水平下（t统计量临界绝对值为1.96）显著区别于1，这与数据的实际生成过程相去甚远（实际生成过程的一次项系数就是1）；

• model.2不包含截距。，在90%置信水平下（t统计量临界绝对值为1.65）不显著区别于1，与数据生成过程比较贴近。

上述结果说明，t检验的估计是很不稳定的。正如1.1节的推导表明，OLS对单位根过程的估计本身就是有偏的，t检验不能直接运用于单位根的判别。
2 DF检验的基本思想
DF检验由Dickey和Fuller提出。它的基本思想是使用蒙特卡洛方法来模拟t统计量的分布，然后再将由实际数据计算的结果与之对比。

但是，DF检验并不是一个规范检验，对于不同数据生成过程具有不同的临界值，并且还与样本量有关，样本量越大，临界值越小。

上文中使用的序列样本量为150，若进行1000次蒙特卡洛试验，需要先生成1000个长度为150的{}序列，简便的方法是生成一个长度为1000*150的序列，再把它平均分为1000份。

set.seed(123)
epsilon0 <- rnorm(200000, 0, 1)
t1 <- t2 <- NULL
c <- NULL
for(i in c(1:1000)) {
  m = 200*i - 199; n = 200*i
  epsilon = epsilon0[m:n]
  y = cumsum(epsilon)[51:200]
  model.1 <- lm(y[2:150] ~ y[1:149])
  model.2 <- lm(y[2:150] ~ y[1:149] - 1)
  t1[i] = (summary(model.1)$coeff[2,1]-1)/summary(model.1)$coeff[2,2]
  t2[i] = (summary(model.2)$coeff[1,1]-1)/summary(model.2)$coeff[1,2]
}
 
par(mfrow = c(1,2), family = "mono",
    plt = c(0.15,0.9,0.12,0.9),
    mgp = c(1.8,0.5,0))
plot(density(t1))
plot(density(t2))

计算t统计量的临界值：

quantile(t1, probs = c(0.01, 0.05, 0.1))
##        1%        5%       10% 
## -3.685972 -2.846401 -2.535439
quantile(t2, probs = c(0.01, 0.05, 0.1))
##        1%        5%       10% 
## -2.501684 -1.848091 -1.594933

根据模拟结果显示：

• 对于带漂移项的随机游走过程（即带截距的模型），其计算出的t统计量有99%的几率大于-3.685972，95%的几率大于-2.846401，90%的几率大于-2.535439；

• 对于随机游走过程（即不带截距的模型），其计算出的t统计量有99%的几率大于-2.501684，95%的几率大于-1.848091，90%的几率大于-1.594933。
  在model.1中，t统计量为-2.0579，大于对应的90%临界值-2.535439，因此不能拒绝数据序列存在单位根的假设；在model.2中，t统计量为-1.16846，也大于对应的90%临界值-1.594933，同样不能拒绝数据序列存在单位根的假设。

3 DF检验
蒙特卡洛方法每次模拟的结果虽然并非完全一致，但在模型形式和样本量确定的情况下会相对稳定，因此在进行DF检验时并不需要每次都进行蒙特卡洛模拟，而是可以直接使用Dickey和Fuller等人的模拟结果。

3.1 理论基础
由

令，则

因此，检验是否显著异于1就等价为检验是否显著异于0。

从前面的介绍可以看出，模型中有无截距对t统计量的计算结果存在很大的影响。实际上，DF检验一共讨论了如下三种情况：

• 仅包含系数，即假设过程为随机游走过程；

• 包含截距，即假设过程为带漂移项的随机游走过程；

• 包含截距和时间变量，即假设过程包含时间趋势。

针对上述三种情况，不仅要确定数据序列是否存在单位根，即是否显著异于0，还要确定这三种情况哪种才是拟合数据生成过程的“最佳”形式。

对于前者，DF检验使用的是类似于t检验的方法，并把三种形式得到的t统计量分别记为、、。

对于后者，DF检验则使用的是类似于F检验的方法，统计量记为：

有约束无约束有约束无约束无约束无约束

• 有约束模型表示估计参数较少者，如式1相比于式2，其不需要估计参数，实际上相当于加了约束条件；

• 为模型的残差平方和，为模型自由度（样本量减去估计参数的数量）。

DF检验设计了三个统计量：

• 对式2应用，约束条件为；

• 对式3为应用，约束条件为；

• 对式3为应用，约束条件为。

当统计量大于对应的临界值时，表明约束条件不成立。

3.2 R中的函数
在R中，可以使用urca工具包中的ur.df()函数进行DF检验。该函数的语法结构如下：

ur.df(y, type = c("none", "drift", "trend"),
      lags = 1,
      selectlags = c("Fixed", "AIC", "BIC"))

• y：数据序列；

• type：有"none"、"drift"、"trend"三个取值，分别对应上述三种情况；

• lags：在DF检验中设置为0；

• 该函数输出结果会给出对应统计量的计算结果及其临界值。

如下分别构造三种模型形式：

library(urca)
model.01 <- ur.df(y2, type = "none", lags = 0)
model.02 <- ur.df(y2, type = "drift", lags = 0)
model.03 <- ur.df(y2, type = "trend", lags = 0)

以model.03为例，使用summary()函数输出的结果分为如下几个部分：

summary(model.03)

模型形式

## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression trend 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt)

• z.diff表示；

• z.lag.1表示；

• +1表示包含截距，若为-1则表示不含截距；

• tt表示时间趋势项。

残差分布

## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.71275 -0.59136  0.02147  0.75579  2.12529

参数估计结果及相关统计量

## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  0.454794   0.317119   1.434   0.1537  
## z.lag.1     -0.062278   0.030319  -2.054   0.0417 *
## tt          -0.001395   0.002249  -0.620   0.5360  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.984 on 146 degrees of freedom
## Multiple R-squared:  0.03056,    Adjusted R-squared:  0.01728 
## F-statistic: 2.301 on 2 and 146 DF,  p-value: 0.1038

和统计量计算结果及其临界值

## Value of test-statistic is: -2.0541 1.551 2.3011 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau3 -3.99 -3.43 -3.13
## phi2  6.22  4.75  4.07
## phi3  8.43  6.49  5.47

这部分是DF检验的重点，具体看下文“结果解析”。

3.3 结果解析
在的基础上加上项形成序列，则为带漂移项的随机游走过程，即，也即。

y3 = 2*c(51:200) + y2
model.11 <- ur.df(y3, type = "none", lags = 0)
model.12 <- ur.df(y3, type = "drift", lags = 0)
model.13 <- ur.df(y3, type = "trend", lags = 0)

先看model.13的结果：

summary(model.13)
## Value of test-statistic is: -2.0541 203.0029 2.3011 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau3 -3.99 -3.43 -3.13
## phi2  6.22  4.75  4.07
## phi3  8.43  6.49  5.47

• 的计算结果-2.0541大于10%显著水平（即90%置信水平）所有临界值-3.13，因此不能拒绝的假设；

• 计算结果为203.0029，大于1%显著水平临界值6.22，因此可以拒绝的假设；计算结果2.3011小于10%显著水平的临界值5.47，因此不能拒绝的假设。

• 综合可得，，。

再看model.12的结果：

summary(model.12)
## Value of test-statistic is: 0.5476 295.7266 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1  6.52  4.63  3.81

的计算结果0.5476大于10%显著水平的临界值-2.57，因此不能拒绝的假设；

的计算结果295.7266大于1%显著水平的临界值6.52，因此可以拒绝的假设；

综合可得，，。

最后看model.11的结果：

summary(model.11)
## Value of test-statistic is: 20.153 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau1 -2.58 -1.95 -1.62

的计算结果20.153大于10%显著水平下的临界值-1.62，因此不能拒绝的假设。
综上所述，可以判断该序列是一个带漂移项的随机游走过程，这与数据生成生成过程相符合。