在 R 中使用 VECM 进行格兰杰因果检验

【问题标题】：Granger causality test using VECM in R在 R 中使用 VECM 进行格兰杰因果检验
【发布时间】：2014-07-14 15:46:19
【问题描述】：

我正在尝试使用 R 中的矢量纠错模型 (VECM) 计算 Granger 因果检验。我使用 tsDyn 包计算了 R 中的 VECM。由于我有 I(1) 和协整变量，因此假设 VECM 实现了 Granger 因果检验。但是我没有在 R 中找到任何可以对 VECM 执行 Granger Granger 因果检验的函数。我想问你，是否有人知道这样的功能。这是我的例子：

dols.est <- dynlm(ts_ln.API.real.1~ts_MR.var.nom.1+L(d(ts_MR.var.nom.1), -3:3)) # Estimate θ with DOLS

est.theta <- dols.est$coefficients[2]

int.mts <- ts.union(ts_ln.API.real.1, ts_MR.var.nom.1) # Create a multivariate time series

VEC.est <- VECM(int.mts, lag=1, r=1, include = c("both"), beta = est. theta)

任何帮助将不胜感激。 提前谢谢你！

【问题讨论】：

    标签： r causality


    【解决方案1】：

    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
}

    当然，我们可以在 causfinder::adfcstable 给出的表格中呈现整个相关的 ADF 统计数据（就像我们在 Procedia 的 CAD 论文中所做的那样）：

    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-因果关系等

    你也可以学习以下论文：

    1. “Causfinder：用于系统分析条件和部分格兰杰因果关系的 R 包”，国际科学与先进技术杂志，2014 年 10 月。http://www.ijsat.com/view.php?id=2014:October:Volume%204%20Issue%2010

    2. “土耳其经常账户赤字的决定因素：条件和部分格兰杰因果关系方法”（扩展版；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 包的 GENERALIZATION（您可以在 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：用于分析大脑功能整合的 R 包

    注意：

    在 R 中：

    可以进行CGC和PGC分析的包：FIAR和causfinder

    在 Matlab 中：

    可以进行CGC和PGC分析的包：

    GCCA（格兰杰因果连通性分析）（Anil SETH）2009：
    MVGC（多元格兰杰因果关系）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

    【讨论】：

      【解决方案2】：

      您可以使用户田-山本方法Toda-Yamamoto implementation in R。

      【讨论】：

        猜你喜欢
        • 1970-01-01
        • 2016-10-25
        • 2021-04-16
        • 2021-04-12
        • 1970-01-01
        • 2017-12-11
        • 2015-10-28
        • 2016-06-14
        • 1970-01-01
        相关资源
        最近更新 更多
        热门标签
        Java Python linux javascript Mysql C# Docker 算法 前端 SpringBoot Redis Vue spring 设计模式 .net core .net kubernetes c++ 数据库 数据结构 大数据 js 机器学习 微服务 Android Go 程序员 面试 JVM ASP.net core 云原生 人工智能 后端 PHP git CSS golang k8s Nginx Django mybatis 深度学习 多线程 React 架构 devops 爬虫 云计算 Spring Boot LeetCode