有关于时序协整分析

你是直接用GDP的数据吗？结果应该是一致的啊。我截了1980-2011年的年度GDP数据，分析如下：

> gdp[1:3,]
      cgdp syear
1 454562.4  1980
2 489156.1  1981
3 532335.1  1982


> fit=ur.df(gdp$cgdp,type="trend",lags=0)
> summary(fit)

###############################################
# Augmented Dickey-Fuller Test Unit Root Test #
###############################################

Test regression trend


Call:
lm(formula = z.diff ~ z.lag.1 + 1 + tt)

Residuals:
     Min       1Q   Median       3Q      Max
-2446389  -201311    11266   286045  1523872

Coefficients:
              Estimate Std. Error t value Pr(>|t|)   
(Intercept)  1.017e+05  3.051e+05   0.333    0.741   
z.lag.1      1.739e-01  2.414e-02   7.202 7.73e-08 ***
tt          -1.618e+04  2.903e+04  -0.557    0.582   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 669200 on 28 degrees of freedom
Multiple R-squared: 0.8825,        Adjusted R-squared: 0.8741
F-statistic: 105.1 on 2 and 28 DF,  p-value: 9.6e-14


Value of test-statistic is: 7.2015 122.7887 105.1156

Critical values for test statistics:
      1pct  5pct 10pct
tau3 -4.15 -3.50 -3.18
phi2  7.02  5.13  4.31
phi3  9.31  6.73  5.61

df统计量是7.2，明显大于-3.18，因此序列非平稳。

> library(tseries)
> adf.test(gdp$cgdp,k=0)

        Augmented Dickey-Fuller Test

data:  gdp$cgdp
Dickey-Fuller = 7.2015, Lag order = 0, p-value = 0.99
alternative hypothesis: stationary

这是用ADF.test()做的结果，和ur.df()的统计量是一致的，为7.2，显然接受原假设，而原假设是非平稳，因此序列非平稳。

这里详细讲一下单位根检验。
单方程的单位根检验一共有三个形式：
1.dx_{t}=bx_{t-1}+e_{t}，这种形式称为简单随机过程，关键在于没有截距，从图形上看它大体是围绕0。
我们用R生成模拟数据的程序可以写为：
x=rnorm(200)
x=cumsum(x)
2.dx_{t}=a+bx_{t-1}+e_{t}，这种形式称为带截距项的随机过程，从图形上看它是线性增长的。
我们用R生成模拟数据的程序可以写为：
x=rnorm(200)
x=cumsum(x)

3.dx_{t}=a+bx_{t-1}+t+e_{t}，这种形式称为带截距项和趋势项的随机过程，从图形上看它是加速增长的。
我们用R生成模拟数据的程序可以写为：
x=rnorm(200)
z=seq(0,2,,200)
x=cumsum(x+z)
1.DF的模型选择
这里提示我们可以根据时序数据的图形判断大体用什么形式来分析单位根。如果是随机游走的形式，则用第一种；如果有线性趋势，则用第二种；如果有加速增加趋势，则用第三种。大部分情况下，使用第二种线性趋势居多，因此很多程序默认的就是第二种形式。
下面的这段程序是用三个模型的模拟数据生成的图形：
x=rnorm(100)
y=cumsum(x+1)
z=seq(1,3,,100)
z=cumsum(x+z)
x=cumsum(x)
par(mfrow=c(2,2))
plot(1:100,x,type="l",main="简单随机游走")
plot(1:100,y,type="l",main="带截距的随机游走")
plot(1:100,z,type="l",main="带趋势的随机游走")
matplot(1:100,cbind(x,y,z),type="l",lty=rep(1,3),col=1:3,main="组合图")

通过比较，你会发现三个模型的特点。当然现实的序列没有这么完美，但即使模型使用有误，判断结果大多数情况下也是大体一致的。

2.ADF的滞后阶数：
首先我们要记住，单位根检验实际针对序列的一个AR模型，通常情况下，我们使用PACF图来确定应该滞后的阶数，当然也可以运用AIC准则来确定最合理的滞后阶数。假设现在我们得到的最优结果是滞后3阶，那么对于ur.df（）而言，本身就滞后了1阶，因此我们lags=2。

3.ADF检验的临界值：
ADF检验的统计量并不服从t分布，而是左偏的，而且三个模型的左偏程度是逐渐递增，我们不能使用传统的t统计量的临界值来判断。事实上，临界值与模型形式有关外，还与样本数有关，但主要与临界值有关。因此我们的t值实际是要与临界值比，如果大于临界值则非平稳，小于临界值则平稳。临界值可以通过蒙特卡罗模拟的方法得到。

再回到gdp的例子，事实上gdp是一个I(2)的时序数据，即要差分2次后才平稳。如果你经常阅读学术期刊的话，特别是宏观的做VECM的论文，就会知道这个是我国gdp数据的一个显著特征。下面给出分析代码：
> ngdp=gdp[-1,]
> ngdp$d=diff(gdp[,1])
> nngdp=ngdp[-1,]
> nngdp$d2=diff(ngdp[,3])
> summary(ur.df(nngdp$d,type="drift",lags=0))

###############################################
# Augmented Dickey-Fuller Test Unit Root Test #
###############################################

Test regression drift


Call:
lm(formula = z.diff ~ z.lag.1 + 1)

Residuals:
     Min       1Q   Median       3Q      Max
-2627881  -229811  -119239   230767  3133183

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.642e+05  2.181e+05   0.753    0.458
z.lag.1     5.936e-02  1.045e-01   0.568    0.575

Residual standard error: 886900 on 27 degrees of freedom
Multiple R-squared: 0.01181,        Adjusted R-squared: -0.02479
F-statistic: 0.3226 on 1 and 27 DF,  p-value: 0.5747


Value of test-statistic is: 0.568 1.2712

Critical values for test statistics:
      1pct  5pct 10pct
tau2 -3.58 -2.93 -2.60
phi1  7.06  4.86  3.94

> summary(ur.df(nngdp$d2,type="drift",lags=0))

###############################################
# Augmented Dickey-Fuller Test Unit Root Test #
###############################################

Test regression drift


Call:
lm(formula = z.diff ~ z.lag.1 + 1)

Residuals:
     Min       1Q   Median       3Q      Max
-2474313  -320867  -149410   146345  2841977

Coefficients:
              Estimate Std. Error t value Pr(>|t|)   
(Intercept)  2.776e+05  1.687e+05   1.646    0.111   
z.lag.1     -1.154e+00  1.931e-01  -5.974 2.27e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 881900 on 27 degrees of freedom
Multiple R-squared: 0.5693,        Adjusted R-squared: 0.5533
F-statistic: 35.69 on 1 and 27 DF,  p-value: 2.266e-06


Value of test-statistic is: -5.974 17.8687

Critical values for test statistics:
      1pct  5pct 10pct
tau2 -3.58 -2.93 -2.60
phi1  7.06  4.86  3.94

parazhu 发表于 2013-2-21 22:56
再回到gdp的例子，事实上gdp是一个I(2)的时序数据，即要差分2次后才平稳。如果你经常阅读学术期刊的话，特别 ...

十分感谢你详尽的回答。。另外，我后来用对数变换的名义GDP来做adf.test()（没做ur.df()），发现log(GDP)是一阶单整的，你觉得这有没有可能？。。。

还有一个问题：如何看ur.df()的输出结果？举个例子：

Value of test-statistic is: 7.2015 122.7887 105.1156

Critical values for test statistics:
      1pct  5pct 10pct
tau3 -4.15 -3.50 -3.18
phi2  7.02  5.13  4.31
phi3  9.31  6.73  5.61

以上是你的代码输出结果。。希望你还记得。。那么是不是对于任何情况，只将“Value of test-statistic is
”最左侧的数值与“tau3”的其中一个数值（取决于所选的显著性水平）比较即可？大于临界值则非平稳，小于临界值则平稳；如果是的话，能麻烦你给我讲一下上面的这些数值是啥意思吗？没时间的话就不用麻烦你了哈！。。

最后一个问题：你认为名义GDP要不要经过平减？。。即先平减，再取对数，再对这个新序列进行各种检验。。虽然我现在只取对数貌似也可以了。。

最后很感谢你的热心帮助！。。

parazhu 发表于 2013-2-21 22:56
再回到gdp的例子，事实上gdp是一个I(2)的时序数据，即要差分2次后才平稳。如果你经常阅读学术期刊的话，特别 ...

对了，我本来做ur.df()时，错把那个p-value当做判别准则了。。后来根据你的回复才知道是要比较统计值和临界值。。。具体怎么比较统计值和临界值需要你的解惑。。。顺带说一下那个ur.df()的p-value是干啥的吧。。它和adf.test()里的p-value不一样的？。。因为用adf.test()来检验平稳性直接看p只就可以了。。比ur.df()貌似方便一点。。

对GDP进行对数转化之后的计算结果，用ur.df()得到结果：
> gdp[1:3,]
      cgdp syear    lngdp
1 454562.4  1980 13.02709
2 489156.1  1981 13.10044
3 532335.1  1982 13.18503
> ngdp$d=diff(gdp[,3])
> ngdp[1:3,]
      cgdp syear    lngdp          d
2 489156.1  1981 13.10044 0.07334648
3 532335.1  1982 13.18503 0.08459150
4 596265.2  1983 13.29844 0.11341229
> nngdp=ngdp[-1,]
> nngdp$d2=diff(ngdp[,4])
> nngdp[1:3,]
      cgdp syear    lngdp         d         d2
3 532335.1  1982 13.18503 0.0845915 0.01124501
4 596265.2  1983 13.29844 0.1134123 0.02882079
5 720805.2  1984 13.48812 0.1896834 0.07627112

这里lngdp是对gdp进行对数转化，d对应的是对lngdp的一阶差分转化，d2是对d的一阶差分转化，即对lngdp的二阶差分转化。
然后我们使用ur.df()计算单位根：

> summary(ur.df(nngdp$lngdp,type="drift",lags=0))

###############################################
# Augmented Dickey-Fuller Test Unit Root Test #
###############################################

Test regression drift


Call:
lm(formula = z.diff ~ z.lag.1 + 1)

Residuals:
     Min       1Q   Median       3Q      Max
-0.08898 -0.04787  0.00261  0.03060  0.15187

Coefficients:
             Estimate Std. Error t value Pr(>|t|)  
(Intercept)  0.314492   0.133901   2.349   0.0264 *
z.lag.1     -0.010339   0.008635  -1.197   0.2416  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.06029 on 27 degrees of freedom
Multiple R-squared: 0.05042,        Adjusted R-squared: 0.01525
F-statistic: 1.434 on 1 and 27 DF,  p-value: 0.2416


Value of test-statistic is: -1.1973 96.2308

Critical values for test statistics:
      1pct  5pct 10pct
tau2 -3.58 -2.93 -2.60
phi1  7.06  4.86  3.94

统计量-1.1973通不过5%的显著性水平，因此存在单位根。

> summary(ur.df(nngdp$d,type="drift",lags=0))

###############################################
# Augmented Dickey-Fuller Test Unit Root Test #
###############################################

Test regression drift


Call:
lm(formula = z.diff ~ z.lag.1 + 1)

Residuals:
      Min        1Q    Median        3Q       Max
-0.082636 -0.042728  0.000812  0.046603  0.084284

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)  0.06412    0.02445   2.622   0.0142 *
z.lag.1     -0.40381    0.14930  -2.705   0.0117 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.04905 on 27 degrees of freedom
Multiple R-squared: 0.2132,        Adjusted R-squared: 0.184
F-statistic: 7.316 on 1 and 27 DF,  p-value: 0.01169


Value of test-statistic is: -2.7047 3.7031

Critical values for test statistics:
      1pct  5pct 10pct
tau2 -3.58 -2.93 -2.60
phi1  7.06  4.86  3.94

统计量-2.705 通不过5%的显著性水平，因此存在单位根。



> summary(ur.df(nngdp$d2,type="drift",lags=0))

###############################################
# Augmented Dickey-Fuller Test Unit Root Test #
###############################################

Test regression drift


Call:
lm(formula = z.diff ~ z.lag.1 + 1)

Residuals:
     Min       1Q   Median       3Q      Max
-0.10801 -0.03429  0.01051  0.03864  0.09040

Coefficients:
             Estimate Std. Error t value Pr(>|t|)   
(Intercept)  0.002414   0.010232   0.236    0.815   
z.lag.1     -0.896264   0.191347  -4.684 7.12e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.055 on 27 degrees of freedom
Multiple R-squared: 0.4483,        Adjusted R-squared: 0.4279
F-statistic: 21.94 on 1 and 27 DF,  p-value: 7.121e-05


Value of test-statistic is: -4.684 10.9706

Critical values for test statistics:
      1pct  5pct 10pct
tau2 -3.58 -2.93 -2.60
phi1  7.06  4.86  3.94

用adf.test()得到的结果：

> adf.test(nngdp$d)

        Augmented Dickey-Fuller Test

data:  nngdp$d
Dickey-Fuller = -3.3838, Lag order = 3, p-value = 0.07831
alternative hypothesis: stationary

> adf.test(nngdp$d2)

        Augmented Dickey-Fuller Test

data:  nngdp$d2
Dickey-Fuller = -3.7696, Lag order = 3, p-value = 0.03643
alternative hypothesis: stationary

这里有几点要说明：首先这个就是adf.test()的P值是可以直接用来判断的，因为这个是用插值得到的结果。我截取adf.test()的源代码做一个分析：

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)

    STAT <- res.sum$coefficients[2,1] / res.sum$coefficients[2,2]

    table <- cbind(c(4.38, 4.15, 4.04, 3.99, 3.98, 3.96),

                   c(3.95, 3.80, 3.73, 3.69, 3.68, 3.66),

                   c(3.60, 3.50, 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.80, 0.87, 0.90, 0.92, 0.93, 0.94),

                   c(0.50, 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, 100000)

    tablep <- c(0.01, 0.025, 0.05, 0.10, 0.90, 0.95, 0.975, 0.99)

    tableipl <- numeric(tablen)

    for(i in (1:tablen))

        tableipl <- approx(tableT, table[, i], n, rule=2)$y（这里是根据不同样本数对临界值进行插值调整）

        interpol <- approx(tableipl, tablep, STAT, rule=2)$y（这里是根据统计量的值在临界值中的位置，利用插值的方法计算出对应的概率值）

因此adf.test()的p值并不是我们一般OLS得到的P值结果。如需要详细了解，可阅读enders的《时间序列分析》。

2.对于adf.test()，需要注意的是：首先，它的模型假设就是一种，带趋势的模型。其次，它默认的滞后阶数是3阶。

既然是默认是三阶滞后，在使用adf.test()的时候要根据实际情况，对滞后阶数要调整。如果要使得ur.df()和adf.test()的结果一致，则GDP的这个例子，应该为：

> adf.test(nngdp$d,k=0)

        Augmented Dickey-Fuller Test

data:  nngdp$d
Dickey-Fuller = -2.7643, Lag order = 0, p-value = 0.2791
alternative hypothesis: stationary

> summary(ur.df(nngdp$d,type="trend",lags=0))

###############################################
# Augmented Dickey-Fuller Test Unit Root Test #
###############################################

Test regression trend


Call:
lm(formula = z.diff ~ z.lag.1 + 1 + tt)

Residuals:
      Min        1Q    Median        3Q       Max
-0.079723 -0.042838  0.001859  0.037072  0.083701

Coefficients:
              Estimate Std. Error t value Pr(>|t|)  
(Intercept)  0.0792738  0.0316573   2.504   0.0189 *
z.lag.1     -0.4201063  0.1519755  -2.764   0.0103 *
tt          -0.0008454  0.0011083  -0.763   0.4525  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.04944 on 26 degrees of freedom
Multiple R-squared: 0.2304,        Adjusted R-squared: 0.1712
F-statistic: 3.892 on 2 and 26 DF,  p-value: 0.03322


Value of test-statistic is: -2.7643 2.6244 3.8921

Critical values for test statistics:
      1pct  5pct 10pct
tau3 -4.15 -3.50 -3.18
phi2  7.02  5.13  4.31
phi3  9.31  6.73  5.61

以上是你的代码输出结果。。希望你还记得。那么是不是对于任何情况，只将“Value of test-statistic is

可能上面的回答没有直接针对你的问题，但上面的回答实际就是最重要的技术细节。
下面针对你的问题做个比较直接的回答：
1.Q:那么是不是对于任何情况，只将"Value of test-statistic is”最左侧的数值与“tau3”的其中一个数值（取决于所选的显著性水平）比较即可？大于临界值则非平稳，小于临界值则平稳?

A:这里三个临界值分别对应三个显著性水平0.01,0.05,0.1。在GDP的例子中，最右侧1% 的水平都通不过，当然无法通过更高的显著性水平了。因此不是比较最右侧的那个值，而是看你选取的显著性水平，再选择你的临界值。一般我们选取5%的显著性水平作为参考，SPSS就是默认5%。

2.Q:顺带说一下那个ur.df()的p-value是干啥的吧。。它和adf.test()里的p-value不一样的？。。因为用adf.test()来检验平稳性直接看p只就可以了。。比ur.df()貌似方便一点。
A: 在进行单位根检验的时候ur.df()的p-value不用看，这个pvalue就是回归系数t检验对应的p值，因此无法用于诊断是否平稳。以GDP为例，在ur.df()计算得到的p值0.0103，就是以下方法得到的。
> 2*pt(-2.764,26)
[1] 0.01035324
既然没有用，那么ur.df()为什么会出现呢？这是因为你使用了summary()这一个函数，其基本的输出格式是固定的，这里就是OLS的输出结构，所以虽然P值无用，这里也输出了。

3.Q:你认为名义GDP要不要经过平减？。。即先平减，再取对数，再对这个新序列进行各种检验。。虽然我现在只取对数貌似也可以了。
A:这里需要根据研究的需要进行。你看大部分期刊的文章，有平减的也有不平减的，牵涉到与其他解释变量以及被解释变量的一致性。只能说说一般情况，以做协整分析为例，因为这个是I(2)的，无论取对数还是不取对数，那么肯定要经过一次差分才能转化为I(1)，才能做VECM。而差分之后的意义就变为GDP增量了。

最后，建议可以读一下enders的计量经济学《时间序列分析》，它的一章非平稳序列对DF检验做了比较详细说明。

parazhu 发表于 2013-2-22 09:49
既然是默认是三阶滞后，在使用adf.test()的时候要根据实际情况，对滞后阶数要调整。如果要使得ur.df()和adf ...

谢谢你的回复！。。那么这样说的话，ur.df()函数更科学啊！因为它比adf.test()函数多了一个type参数选择。。。那我还是用ur.df()好了。。另外有一个问题，就是如果使用ur.df()的话，赤池信息准则怎么看？我没看到AIC啊？。。怎么根据ur.df()的输出结果选择最优的lags参数？。。。