0首先，通过考虑所有滞后的普通样本的使用来执行 ADF 测试。例如：（causfinder::adfcs 和 causfinder::adfcstable）：
adfcs <- function (t, max = floor(12 * (length(t)/100)^(1/4)), type = c("c"))
# Augmented Dickey-Fuller function that takes into account the usage of common sample for all the lags
{
x <- ts(t)
x1d <- diff(x, differences = 1)
x1l <- lag(x, -1)
if (max == 0) {
x_names <- c("x1d", "x1l")
DLDlag <- ts.intersect(x1d, x1l)
DLDlag.df <- data.frame(DLDlag, obspts = c(time(DLDlag)))
}
else {
x_names <- c("x1d", "x1l", sapply(1:max, function(i) paste("x1d", i, "l", sep = "")))
}
if (max != 0) {
for (i in as.integer(1:max)) {
assign(x_names[i + 2], lag(x1d, -i))
}
DLDlag <- do.call(ts.intersect, sapply(x_names, as.symbol))
DLDlag.df <- data.frame(DLDlag, obspts = c(time(DLDlag)))
DifferenceLags <- as.vector(names(DLDlag.df), mode = "any")[3:(length(DLDlag.df) - 1)]
}
lmresults <- array(list())
SBCvalues <- array(list())
AICvalues <- array(list())
for (i in as.integer(0:max)) {
if (type == c("nc")) {
if (i == 0) {
lmresults[[max + 1]] <- lm(as.formula(paste("x1d ~x1l")), data = DLDlag.df)
SBCvalues[[max + 1]] <- BIC(lmresults[[max + 1]])
AICvalues[[max + 1]] <- AIC(lmresults[[max + 1]])
}
if (i > 0) {
lmresults[[i]] <- lm(as.formula(paste("x1d ~ x1l+", paste(DifferenceLags[1:i], collapse = "+"))), data = DLDlag.df)
SBCvalues[[i]] <- BIC(lmresults[[i]])
AICvalues[[i]] <- AIC(lmresults[[i]])
}
}
if (type == c("c")) {
if (i == 0) {
lmresults[[max + 1]] <- lm(as.formula(paste("x1d ~1+x1l")), data = DLDlag.df)
SBCvalues[[max + 1]] <- BIC(lmresults[[max + 1]])
AICvalues[[max + 1]] <- AIC(lmresults[[max + 1]])
}
if (i > 0) {
lmresults[[i]] <- lm(as.formula(paste("x1d ~ 1+x1l+", paste(DifferenceLags[1:i], collapse = "+"))), data = DLDlag.df)
SBCvalues[[i]] <- BIC(lmresults[[i]])
AICvalues[[i]] <- AIC(lmresults[[i]])
}
}
if (type == c("ct")) {
if (i == 0) {
lmresults[[max + 1]] <- lm(as.formula(paste("x1d ~ 1+x1l+seq_along(x1d)", collapse = "")), data = DLDlag.df)
SBCvalues[[max + 1]] <- BIC(lmresults[[max + 1]])
AICvalues[[max + 1]] <- AIC(lmresults[[max + 1]])
}
if (i > 0) {
lmresults[[i]] <- lm(as.formula(paste("x1d ~ 1+x1l+seq_along(x1d)+", paste(DifferenceLags[1:i], collapse = "+"))), data = DLDlag.df)
SBCvalues[[i]] <- BIC(lmresults[[i]])
AICvalues[[i]] <- AIC(lmresults[[i]])
}
}
}
out <- list()
out$optmins <- list(which.min(SBCvalues), which.min(AICvalues))
out$SBCAIC <- as.data.frame(cbind(SBCvalues, AICvalues))
typespecified <- type
if (which.min(SBCvalues) == max + 1) {
scs <- (max + 2) - (0 + 1)
out$adfcst <- unitrootTest(x[scs:length(x)], lags = 0,
type = typespecified)
}
else {
scs <- (max + 2) - (which.min(SBCvalues) + 1)
out$adfcst <- unitrootTest(x[scs:length(x)], lags = which.min(SBCvalues),
type = typespecified)
}
out
}
当然，可以在一个表格中呈现整个相关的 ADF 统计数据（就像我们在 Procedia 的 CAD 论文中所做的那样），该表格由 causfinder::adfcstable 给出：
adfcstable <- function (d, max = 5)
{
d <- as.data.frame(d)
LevelADFtable <- matrix(, nrow = dim(d)[[2]] * 3, ncol = 10)
FirstDiffADFtable <- matrix(, nrow = dim(d)[[2]] * 3, ncol = 9)
Result <- matrix(, nrow = dim(d)[[2]] * 3, ncol = 1)
ADFtable <- as.data.frame(cbind(LevelADFtable, FirstDiffADFtable, Result), stringsAsFactors = FALSE)
colnames(ADFtable) <- c("var", "type", "inc", "levelt", "Pc", "c", "Pt", "t", "prob", "omlo", "type", "inc", "1stDifft", "Pc", "c", "Pt", "t", "prob", "omlo", "intorder")
for (i in as.integer(1:dim(d)[[2]])) {
for (j in as.integer(1:3)) {
ADFtable[3 * (i - 1) + j, 1] <- colnames(d)[[i]]
}
ADFtable[3 * i - 2, 2] <- "dt"
ADFtable[3 * i - 2, 11] <- "dt"
ADFtable[3 * i - 1, 2] <- "d"
ADFtable[3 * i - 1, 11] <- "d"
ADFtable[3 * i, 2] <- "-"
ADFtable[3 * i, 11] <- "-"
ADFtable[3 * i - 2, 3] <- round(adfcs(d[, i], type = c("ct"))$adfcst@test$regression$coefficients[2, 1], digits = 3)
ADFtable[3 * i - 1, 3] <- round(adfcs(d[, i], type = c("c"))$adfcst@test$regression$coefficients[2, 1], digits = 3)
ADFtable[3 * i, 3] <- round(adfcs(d[, i], type = c("nc"))$adfcst@test$regression$coefficients[1, 1], digits = 3)
ADFtable[3 * i - 2, 12] <- round(adfcs(diff(d[, i], differences = 1), type = c("ct"))$adfcst@test$regression$coefficients[2, 1], digits = 3)
ADFtable[3 * i - 1, 12] <- round(adfcs(diff(d[, i], differences = 1), type = c("c"))$adfcst@test$regression$coefficients[2, 1], digits = 3)
ADFtable[3 * i, 12] <- round(adfcs(diff(d[, i], differences = 1), type = c("nc"))$adfcst@test$regression$coefficients[1, 1], digits = 3)
ADFtable[3 * i - 2, 4] <- round(adfcs(d[, i], type = c("ct"))$adfcst@test$statistic, digits = 3)
ADFtable[3 * i - 1, 4] <- round(adfcs(d[, i], type = c("c"))$adfcst@test$statistic, digits = 3)
ADFtable[3 * i, 4] <- round(adfcs(d[, i], type = c("nc"))$adfcst@test$statistic, digits = 3)
ADFtable[3 * i - 2, 13] <- round(adfcs(diff(d[, i], differences = 1), type = c("ct"))$adfcst@test$statistic, digits = 3)
ADFtable[3 * i - 1, 13] <- round(adfcs(diff(d[, i], differences = 1), type = c("c"))$adfcst@test$statistic, digits = 3)
ADFtable[3 * i, 13] <- round(adfcs(diff(d[, i], differences = 1), type = c("nc"))$adfcst@test$statistic, digits = 3)
ADFtable[3 * i - 2, 5] <- round(adfcs(d[, i], type = c("ct"))$adfcst@test$regression$coefficients[1, 4], digits = 3)
ADFtable[3 * i - 2, 7] <- round(adfcs(d[, i], type = c("ct"))$adfcst@test$regression$coefficients[3, 4], digits = 3)
ADFtable[3 * i - 1, 5] <- round(adfcs(d[, i], type = c("c"))$adfcst@test$regression$coefficients[1, 4], digits = 3)
ADFtable[3 * i - 1, 7] <- "X"
ADFtable[3 * i, 5] <- "X"
ADFtable[3 * i, 7] <- "X"
ADFtable[3 * i - 2, 14] <- round(adfcs(diff(d[, i], differences = 1), type = c("ct"))$adfcst@test$regression$coefficients[1, 4], digits = 3)
ADFtable[3 * i - 2, 16] <- round(adfcs(diff(d[, i], differences = 1), type = c("ct"))$adfcst@test$regression$coefficients[3, 4], digits = 3)
ADFtable[3 * i - 1, 14] <- round(adfcs(diff(d[, i], differences = 1), type = c("c"))$adfcst@test$regression$coefficients[1, 4], digits = 3)
ADFtable[3 * i - 1, 16] <- "X"
ADFtable[3 * i, 14] <- "X"
ADFtable[3 * i, 16] <- "X"
if (ADFtable[3 * i - 2, 5] < 0.05) {
ADFtable[3 * i - 2, 6] <- "s"
}
else {
ADFtable[3 * i - 2, 6] <- " "
}
if (ADFtable[3 * i - 2, 7] < 0.05) {
ADFtable[3 * i - 2, 8] <- "s"
}
else {
ADFtable[3 * i - 2, 8] <- " "
}
if (ADFtable[3 * i - 1, 5] < 0.05) {
ADFtable[3 * i - 1, 6] <- "s"
}
else {
ADFtable[3 * i - 1, 6] <- " "
}
ADFtable[3 * i - 1, 8] <- "X"
ADFtable[3 * i, 6] <- "X"
ADFtable[3 * i, 8] <- "X"
if (ADFtable[3 * i - 2, 14] < 0.05) {
ADFtable[3 * i - 2, 15] <- "s"
}
else {
ADFtable[3 * i - 2, 15] <- " "
}
if (ADFtable[3 * i - 2, 16] < 0.05) {
ADFtable[3 * i - 2, 17] <- "s"
}
else {
ADFtable[3 * i - 2, 17] <- " "
}
if (ADFtable[3 * i - 1, 14] < 0.05) {
ADFtable[3 * i - 1, 15] <- "s"
}
else {
ADFtable[3 * i - 1, 15] <- " "
}
ADFtable[3 * i - 1, 17] <- "X"
ADFtable[3 * i, 15] <- "X"
ADFtable[3 * i, 17] <- "X"
ADFtable[3 * i - 2, 9] <- round(adfcs(d[, i], type = c("ct"))$adfcst@test$p.value[[1]], digits = 3)
ADFtable[3 * i - 1, 9] <- round(adfcs(d[, i], type = c("c"))$adfcst@test$p.value[[1]], digits = 3)
ADFtable[3 * i, 9] <- round(adfcs(d[, i], type = c("nc"))$adfcst@test$p.value[[1]], digits = 3)
ADFtable[3 * i - 2, 18] <- round(adfcs(diff(d[, i], differences = 1), type = c("ct"))$adfcst@test$p.value[[1]], digits = 3)
ADFtable[3 * i - 1, 18] <- round(adfcs(diff(d[, i], differences = 1), type = c("c"))$adfcst@test$p.value[[1]], digits = 3)
ADFtable[3 * i, 18] <- round(adfcs(diff(d[, i], differences = 1), type = c("nc"))$adfcst@test$p.value[[1]], digits = 3)
ADFtable[3 * i - 2, 10] <- round(adfcs(d[, i], type = c("ct"))$adfcst@test$parameter, digits = 3)
ADFtable[3 * i - 1, 10] <- round(adfcs(d[, i], type = c("c"))$adfcst@test$parameter, digits = 3)
ADFtable[3 * i, 10] <- round(adfcs(d[, i], type = c("nc"))$adfcst@test$parameter, digits = 3)
ADFtable[3 * i - 2, 19] <- round(adfcs(diff(d[, i], differences = 1), type = c("ct"))$adfcst@test$parameter, digits = 3)
ADFtable[3 * i - 1, 19] <- round(adfcs(diff(d[, i], differences = 1), type = c("c"))$adfcst@test$parameter, digits = 3)
ADFtable[3 * i, 19] <- round(adfcs(diff(d[, i], differences = 1), type = c("nc"))$adfcst@test$parameter, digits = 3)
if (sum(as.numeric(c(ADFtable[3 * i - 2, 9] < 0.05 && ADFtable[3 * i - 2, 3] < 0, ADFtable[3 * i - 1, 9] < 0.05 && ADFtable[3 * i - 1, 3] < 0, ADFtable[3 * i, 9] < 0.05 && ADFtable[3 * i, 3] < 0))) > 1) {
ADFtable[3 * i - 1, 20] <- "I(0)"
}
else {
if (sum(as.numeric(c(ADFtable[3 * i - 2, 18] < 0.05 && ADFtable[3 * i - 2, 12] < 0, ADFtable[3 * i - 1, 18] < 0.05 && ADFtable[3 * i - 1, 12] < 0, ADFtable[3 * i, 18] < 0.05 && ADFtable[3 * i, 12] < 0))) > 1) {
ADFtable[3 * i - 1, 20] <- "I(1)"
}
else {
ADFtable[3 * i - 1, 20] <- "variableoi"
}
}
ADFtable[3 * i - 2, 20] <- ""
ADFtable[3 * i, 20] <- ""
}
ADFtable
}
1即使对于 VECM 模型（例如，我们的变量是 I(1) 和协整变量），我们也会根据时间序列水平上的 VAR 模型的信息标准来选择滞后数。相同职责不同能力的功能：
vars::VARselect # or
FIAR::ARorder # or
causfinder::ARorderG # or
causfinder::VARomlop (the last package is not free)
将在级别的变量上运行（没有差异）。
2要检查协整，请使用
ca.jo(..,K=cointegrationLength)
ca.jo 中的参数 K 控制 VECM 模型的滞后数。将在 VARomlop（或其他）中发现的滞后数作为参数 K 传递。使用 ca.jo 确定协整等级。ecdet 选项“none”表示协整方程中没有截距，“const”表示协整方程中的常数项，“trend”表示协整方程中的趋势变量。
长期关系的规范化是根据 mydata 中的第一列指定的，可以通过更改提交的 mydata 的列来更改（如果需要），例如 mydata[,"X2","X1","X3", “X4”]。
K 是 VAR 在级别上的滞后数，因此 K - 1 是 VECM 表示中的滞后数。例如，
summary(ca.jo(mydata, ecdet="none", type="eigen", K=29))
根据特征/迹检验的结果，判断协整是否存在，协整秩为r。
3如果在系统中的序列中检测到上述协整，则通过考虑协整等级来拟合矢量误差校正模型（VECM）。即，我们在上述步骤中使用 ca.jo 的协整向量来拟合 VECM 模型。ca.jo 的结果和协整向量的数量被传递给 cajorls。cajorls 有参数 r（协整秩）。
通过使用命令 cajorls() 估计受限 VECM 来生成归一化协整向量。例如，
cajorls(...,K=lagLength)
cajorls(ca.jo(mydata, ecdet="none", type="eigen", K=29),r=1)
纠错项可以只包含在 VECM 的每个方程中一次。它要么滞后 1，要么滞后 p，其中 p 是 VECM 的滞后阶；VECM 的相应表示被称为长期和暂时的；它仍然是相同的模型，只是表示不同；我们选择我们喜欢的那个。
4要将 VECM 转换为 VAR：
vars::vec2var
分析：
http ://www.r-bloggers.com/cointegration-r-irish-mortgage-debt-and-property-prices/
也可以回答您的问题。
如果您有一个 2 变量 VAR 系统，请保持经典 G 因果关系。如果你有一个“>2”变量的 VAR 系统，你必须去高级 G 因果关系：条件 G 因果关系、部分 G 因果关系、谐波 G 因果关系、规范 G 因果关系、全局 G 因果关系等。
你也可以研究以下论文：
“Causfinder：用于条件和部分 Granger 因果关系的系统分析的 R 包”，国际科学与先进技术杂志，2014 年 10 月。http://www.ijsat.com/view.php?id=2014:October:Volume % 204%20问题%2010
“土耳其经常账户赤字的决定因素：有条件和部分格兰杰因果关系法”（扩展版；33 页）
https://www.academia.edu/17698799/Determinants_of_Current_Account_Deficit_in_Turkey_The_Conditional_and_Partial_Granger_Causality_Approach_Extended_
“土耳其经常账户赤字的决定因素：条件和部分格兰杰因果关系方法”（9 页），Procedia 经济与金融，卷。2015 年 2 月 26 日，第 92-100 页
https://www.academia.edu/17057780/Determinants_of_Current_Account_Deficit_in_Turkey_The_Conditional_and_Partial_Granger_Causality_Approach
causfinder 是 FIAR 包的通用化（您可以在 CRAN 档案中找到 FIAR。FIAR 0.3、0.4 和 0.5。FIAR 完全免费）。
https://cran.r-project.org/src/contrib/Archive/FIAR
我强烈建议 FIAR 0.3，因为 0.3 版本显然比更高版本更具可扩展性。甚至你不需要分析 0.4 和 0.5 版本。因此，我在 FIAR 0.3 上构建了 causfinder。
在 FIAR 中，您会一一找到 CGC。在 causfinder 中，它系统地一次性提供所有 CGC。在 6 变量系统中，有 6*5=30 个 CGC 和 30 个 PGC。这 30+30=60 个 CGC 和 PGC 在 FIAR（60 个命令）中一一计算。在 causfinder 中，这 30+30 次 GC 只用 2 个命令计算。在 5 变量系统中，有 5*4=20 个 CGC 和 20 个 PGC。这 20+20=40 个 CGC 和 PGC 在 FIAR（40 个命令）中一一计算。在 causfinder 中，这 20+20 个 GC 只用 2 个命令计算。
causfinder（超过 FIAR）提供的是极快的速度/节奏、简单性、可视化和易于分析；没有其他的。
如果你想学习 CGC 或 PGC，你也可以通过 FIAR 来完成： 统计软件杂志 (JSS)：https://www.jstatsoft.org/article/view/v044i13
FIAR: An R Package for Analyzing Functional Integration in大脑
注意：
在 R 中：
可以进行CGC和PGC分析的包：FIAR和causfinder
在 Matlab 中：
可以进行 CGC 和 PGC 分析的包：
GCCA (Granger Causal Connectivity Analysis) (Anil SETH) 2009:
MVGC (Multivariate Granger Causality) 2014: GCCA
GrangerCausalityGUI 的新版本（冯建峰团队在 2008-2013 年期间与一些论文一起开发的成果工作。
2011 年，Roelstraete 和 Rosseel 处理高级 G 因果分析的 FIAR R 包发现了 GCCA 中的一个错误！
据我所知，在其他统计/计量经济学软件中，没有可以执行 CGC 和 PGC 的包/功能。
在 Matlab 中编程肯定比在 R 中更难。因此，我在 R 中编写了 causfinder（在我体验了 Gretl 和 Eviews 中的编码之后）。（我们认为因此我们 R！）
5在您获得 VAR（来自 VECM）后，其中由协整引起的限制被加载到 VAR 的系数上；（如果您在步骤 0 中有 ">2"-变量系统，请执行以下操作。如果没有，则 R 中已经有经典的 G 因果包；使用它们）
FIAR::condGranger
FIAR::partGranger
或者
causfinder::conditionalGgFp
causfinder::partialGgFp
如果你也想要引导，那么：
causfinder::conditionalGblup
causfinder::partialGblup