R -应用adf.分组测试

Question

我有一个数据帧bbm与变量ticker, variable和value。我想通过adf应用增广迪基-富勒检验。测试函数按代码和变量分组。R应该在初始data.frame中添加一个新的列，其中包含相应的p值。

我试着

x <- with(bbm, tapply(value, list(ticker, variable), adf.test$p.value))
cbind(bbm, x)

生成Error in adf.test$p.value : object of type 'closure' is not subsettable。

Then I try

x <- with(bbm, tapply(value, list(ticker, variable), as.list(adf.test)$p.value))
cbind(bbm, x)

这会产生一个结果，但在新列中不是我想要的结果。即使我将代码中的p值更改为method，它仍然会产生一些奇数。

然后我尝试使用ddply:

bbm<-ddply(bbm, .(ticker, variable), mutate, df=adf.test(value)$p.value)

生成Error: wrong embedding Dimension。

我该如何解决这个问题?有什么建议吗?

下面是df:

的示例

            ticker                    variable   value
1  1002Z AV Equity        BS_CUSTOMER_DEPOSITS 29898.0
2  1002Z AV Equity        BS_CUSTOMER_DEPOSITS 31302.0
3  1002Z AV Equity        BS_CUSTOMER_DEPOSITS 29127.0
4  1002Z AV Equity        BS_CUSTOMER_DEPOSITS 24056.0
5  1002Z AV Equity        BS_CUSTOMER_DEPOSITS 22080.0
6  1002Z AV Equity        BS_CUSTOMER_DEPOSITS 22585.0
7  1002Z AV Equity        BS_CUSTOMER_DEPOSITS 22674.0
8  1002Z AV Equity        BS_CUSTOMER_DEPOSITS 21733.0
9  1002Z AV Equity        BS_CUSTOMER_DEPOSITS 22016.0
10 1002Z AV Equity        BS_CUSTOMER_DEPOSITS 21999.0
11 1002Z AV Equity        BS_CUSTOMER_DEPOSITS 22013.0
12 1002Z AV Equity        BS_CUSTOMER_DEPOSITS 21135.0
13 1002Z AV Equity                 BS_TOT_LOAN 28476.0
14 1002Z AV Equity                 BS_TOT_LOAN 29446.0
15 1002Z AV Equity                 BS_TOT_LOAN 29273.0
16 1002Z AV Equity                 BS_TOT_LOAN 27579.0
17 1002Z AV Equity                 BS_TOT_LOAN 20769.0
18 1002Z AV Equity                 BS_TOT_LOAN 21370.0
19 1002Z AV Equity                 BS_TOT_LOAN 22306.0
20 1002Z AV Equity                 BS_TOT_LOAN 21013.0
21 1002Z AV Equity                 BS_TOT_LOAN 21810.0
22 1002Z AV Equity          BS_TIER1_CAP_RATIO     6.5
23 1002Z AV Equity          BS_TIER1_CAP_RATIO     6.2
24 1002Z AV Equity          BS_TIER1_CAP_RATIO     7.9
25 1002Z AV Equity          BS_TIER1_CAP_RATIO     9.2
26 1002Z AV Equity          BS_TIER1_CAP_RATIO     8.5
27 1002Z AV Equity          BS_TIER1_CAP_RATIO     6.6
28 1002Z AV Equity          BS_TIER1_CAP_RATIO     9.6
29 1002Z AV Equity BS_TOT_CAP_TO_RISK_BASE_CAP    11.5
30 1002Z AV Equity BS_TOT_CAP_TO_RISK_BASE_CAP    10.9

 > dput(head(select(bbm, ticker, variable, value), 30))
structure(list(ticker = c("1002Z AV Equity", "1002Z AV Equity", 
"1002Z AV Equity", "1002Z AV Equity", "1002Z AV Equity", "1002Z AV Equity", 
"1002Z AV Equity", "1002Z AV Equity", "1002Z AV Equity", "1002Z AV Equity", 
"1002Z AV Equity", "1002Z AV Equity", "1002Z AV Equity", "1002Z AV Equity", 
"1002Z AV Equity", "1002Z AV Equity", "1002Z AV Equity", "1002Z AV Equity", 
"1002Z AV Equity", "1002Z AV Equity", "1002Z AV Equity", "1002Z AV Equity", 
"1002Z AV Equity", "1002Z AV Equity", "1002Z AV Equity", "1002Z AV Equity", 
"1002Z AV Equity", "1002Z AV Equity", "1002Z AV Equity", "1002Z AV Equity"
), variable = structure(c(4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 8L, 8L, 8L, 8L, 
8L, 8L, 8L, 9L, 9L), .Label = c("PX_LAST", "PE_RATIO", "VOL_MEAN", 
"BS_CUSTOMER_DEPOSITS", "BS_TOT_LOAN", "*", "RN366", "BS_TIER1_CAP_RATIO", 
"BS_TOT_CAP_TO_RISK_BASE_CAP", "RETURN_COM_EQY", "BS_LEV_RATIO_TO_TANG_CAP",
"NPLS_TO_TOTAL_LOANS"), class = "factor"), value = c(29898, 31302, 
29127, 24056, 22080, 22585, 22674, 21733, 22016, 21999, 22013, 
21135, 28476, 29446, 29273, 27579, 20769, 21370, 22306, 21013, 
21810, 6.5, 6.2, 7.9, 9.2, 8.5, 6.6, 9.6, 11.5, 10.9)), .Names = c("ticker", 
"variable", "value"), row.names = c(NA, 30L), class = "data.frame")

哦，使用模拟dplyr函数也会产生与ddply相同的误差。

Answer 1

这是一个整洁的解决方案:

bbm %>% 
    group_by(ticker,variable) %>% 
    summarise(pval = ifelse(n() <= 3,NA, adf.test(value)$p.value))
# A tibble: 4 x 3
# Groups:   ticker [?]
  ticker          variable                       pval
  <chr>           <fct>                         <dbl>
1 1002Z AV Equity BS_CUSTOMER_DEPOSITS         0.01  
2 1002Z AV Equity BS_TOT_LOAN                  0.951 
3 1002Z AV Equity BS_TIER1_CAP_RATIO           0.0118
4 1002Z AV Equity BS_TOT_CAP_TO_RISK_BASE_CAP NA     
Warning message:
In adf.test(value) : p-value smaller than printed p-value

您可以使用基本R ifelse函数来检查每组中是否存在少于3个点(这会将pval设置为NA)，否则您可以运行adf.test

我试了一下，@erasmortg似乎是正确的。错误"嵌入"来自没有足够的数据点来实际运行adf.test函数。

这需要至少四个数据点:

> adf.test(rnorm(1))
Error in embed(y, k) : wrong embedding dimension
> adf.test(rnorm(2))
Error in embed(y, k) : wrong embedding dimension
> adf.test(rnorm(3))
Error in res.sum$coefficients[2, 1] : subscript out of bounds
> adf.test(rnorm(4))
    Augmented Dickey-Fuller Test
data:  rnorm(4)
Dickey-Fuller = NaN, Lag order = 1, p-value = NA
alternative hypothesis: stationary

Answer 2

问题似乎出在组太小而无法处理。处理这个问题的一个选项是创建一个自定义函数来捕获错误(使用tryCatch，并通过lapply()调用传递该函数，如下所示:

testx <- function (x) {
  return(tryCatch(adf.test(x), error=function(e) NULL))
}
g<- lapply(split(bbm, bbm$variable), function(x) testx(x$value))
str(g)
#List of 12
# $ PX_LAST                    : NULL
# $ PE_RATIO                   : NULL
# $ VOL_MEAN                   : NULL
# $ BS_CUSTOMER_DEPOSITS       :List of 6
# ..$ statistic  : Named num -4.86
#  .. ..- attr(*, "names")= chr "Dickey-Fuller"
#  ..$ parameter  : Named num 2
#  .. ..- attr(*, "names")= chr "Lag order"
#  ..$ alternative: chr "stationary"
#  ..$ p.value    : num 0.01
#  ..$ method     : chr "Augmented Dickey-Fuller Test"
#  ..$ data.name  : chr "x"
#  ..- attr(*, "class")= chr "htest"
# $ BS_TOT_LOAN                :List of 6
#  ..$ statistic  : Named num -0.784
#  .. ..- attr(*, "names")= chr "Dickey-Fuller"
#  ..$ parameter  : Named num 2
#  .. ..- attr(*, "names")= chr "Lag order"
#  ..$ alternative: chr "stationary"
#  ..$ p.value    : num 0.951
#  ..$ method     : chr "Augmented Dickey-Fuller Test"
#  ..$ data.name  : chr "x"
#  ..- attr(*, "class")= chr "htest"
# $ *                          : NULL
# $ RN366                      : NULL
# $ BS_TIER1_CAP_RATIO         :List of 6
#  ..$ statistic  : Named num -4.33
#  .. ..- attr(*, "names")= chr "Dickey-Fuller"
#  ..$ parameter  : Named num 1
#  .. ..- attr(*, "names")= chr "Lag order"
#  ..$ alternative: chr "stationary"
#  ..$ p.value    : num 0.0118
#  ..$ method     : chr "Augmented Dickey-Fuller Test"
#  ..$ data.name  : chr "x"
#  ..- attr(*, "class")= chr "htest"
# $ BS_TOT_CAP_TO_RISK_BASE_CAP: NULL
# $ RETURN_COM_EQY             : NULL
# $ BS_LEV_RATIO_TO_TANG_CAP   : NULL
# $ NPLS_TO_TOTAL_LOANS        : NULL

这将创建一个长度为12(每个因子一个)的列表对象g，其中，对于有效的adf。测试调用，元素被相关的特征填充，其余的通过NULL。

如果感兴趣的参数只是每组的p.value，则可以将前面的lapply包裹在sapply()周围以获得以下对象:

h<- sapply(lapply(split(bbm, bbm$variable), function(x) testx(x$value)), function(x) print(x$p.value))
str(h)
#List of 12
# $ PX_LAST                    : NULL
# $ PE_RATIO                   : NULL
# $ VOL_MEAN                   : NULL
# $ BS_CUSTOMER_DEPOSITS       : num 0.01
# $ BS_TOT_LOAN                : num 0.951
# $ *                          : NULL
# $ RN366                      : NULL
# $ BS_TIER1_CAP_RATIO         : num 0.0118
# $ BS_TOT_CAP_TO_RISK_BASE_CAP: NULL
# $ RETURN_COM_EQY             : NULL
# $ BS_LEV_RATIO_TO_TANG_CAP   : NULL
# $ NPLS_TO_TOTAL_LOANS        : NULL

根据评论，如果ticker和variable都需要分组，这将产生期望的结果:

g<- lapply(split(bbm, list(bbm$variable, bbm$ticker)), function(x) testx(x$value))
#to remove the NULL which are not needed:
g[g != "NULL"]