Тест причинности Грейнджера с использованием VECM в R

Question

Тест причинности Грейнджера с использованием VECM в R

Я пытаюсь вычислить тест причинности Грейнджера, используя модель коррекции ошибок вектора (VECM) в R. Я рассчитал VECM в R, используя пакет tsDyn. Поскольку у меня есть I(1) и коинтегрированные переменные, предполагается, что VECM реализует тест причинности Грейнджера. Однако я не нашел ни одной функции в R, которая могла бы выполнить тест причинности Грейнджер Грейнджер для VECM. Я хотел бы спросить Вас, знает ли кто-то такую функцию. Вот мой пример:

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)

Любая помощь будет оценена. Заранее спасибо!

1

r causality

Источник

user3837489 14 июл '14 в 15:46

1 ответ

Другие вопросы по тегам r causality

user2787767 03 май '16 в 13:01 2016-05-03 13:01 · Answer 1 · 2016-05-03 13:01

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 в виде таблицы (как мы это делали в нашей статье по CAD в Methoia), которую дает 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)

Аргумент K в ca.jo контролирует количество лагов модели VECM. Передайте число лагов, найденных в VARomlop (или в других), в качестве аргумента K. Определите ранг коинтеграции, используя ca.jo. Параметр ecdet имеет значение "none" для отсутствия перехвата в уравнении коинтеграции, "const" для постоянного члена в уравнении коинтеграции и "trend" для переменной тренда в уравнении коинтеграции.

Нормализация долгосрочных отношений задается в соответствии с 1-м столбцом в mydata, который может быть изменен (при желании) путем изменения столбцов отправленных mydata с помощью, например, mydata[,"X2","X1","X3","Х4"].

K - количество лагов для уровней VAR, поэтому K - 1 - количество лагов в представлении VECM. например,

summary(ca.jo(mydata, ecdet="none", type="eigen", K=29))

Основываясь на результатах собственных / следовых испытаний, определите наличие коинтеграции и ранг коинтеграции r.

3 Если выше обнаружена коинтеграция среди рядов в системе, установите модель коррекции ошибок вектора (VECM), принимая во внимание ранг коинтеграции. то есть мы подгоняем модель VECM, используя векторы коинтеграции ca.jo на предыдущем шаге. Результат ca.jo и количество векторов коинтеграции передаются в cajorls. у cajorls есть аргумент r (ранг коинтеграции).

Нормализованный вектор коинтеграции получают путем оценки ограниченного VECM с помощью команды cajorls(). например,

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: пакет R для системного анализа условных и частичных причин Грейнджера", Международный журнал науки и передовых технологий, октябрь 2014 года. http://www.ijsat.com/view.php?id=2014:October:Volume%204%20Issue%2010
"Детерминанты дефицита текущего счета в Турции: подход условной и частичной причинности Грейнджера" (расширенный; 33 страницы) https://www.academia.edu/17698799/Determinants_of_Current_Account_Deficit_in_Turkey_The_Conditional_and_Partial_Granger_Causality_Approach_Extended_

"Детерминанты дефицита текущего счета в Турции: подход условной и частичной причинности Грейнджера" (9 страниц), "Methoia Economics and Finance", Vol. 26, 2015, стр.92-100 https://www.academia.edu/17057780/Determinants_of_Current_Account_Deficit_in_Turkey_The_Conditional_and_Partial_Granger_Causality_Approach

causfinder - это ОБОБЩЕНИЕ пакета FIAR (FIAR можно найти в архиве CRAN. 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.

В FIAR вы найдете CGC один за другим. В causfinder он дает ВСЕ CGCs системно и сразу. В системе с 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 командами.

Каузфиндер обеспечивает (по сравнению с 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 (итоговая работа группы Jianfeng FENG, которая была разработана в сопровождении некоторых документов до 2008-2013 гг.

В 2011 году пакет Roelstraete и Rosseel FIAR R, который выполняет расширенный анализ G-причинности, выявили ошибку в GCCA!

Насколько мне известно, в других статистических / эконометрических программах нет пакетов / функций, которые могли бы выполнять CGC и PGC.
Программирование в Matlab определенно сложнее, чем в R. Поэтому я написал causfinder в R (после того, как я попробовал кодирование в Gretl и Eviews.). (Мы думаем поэтому мы R!)

5 После того, как вы получили VAR (от VECM), в котором ограничение, вызванное коинтеграцией, загружается на коэффициенты VAR; (сделайте следующее, если у вас на шаге 0 есть система с>> 2-переменной. Если нет, в R уже есть классические пакеты G-причинности; используйте их)

FIAR::condGranger
FIAR::partGranger

или же

causfinder::conditionalGgFp
causfinder::partialGgFp

Если вы также хотите начать загрузку, то:

causfinder::conditionalGblup
causfinder::partialGblup

user6423956 28 апр '20 в 00:00 2020-04-28 00:00 · Answer 2 · 2020-04-28 00:00

Вы можете использовать Тод-Ямамото подход реализацию Тода-Ямамото в R.

0

Источник

user6423956 28 апр '20 в 00:00