![ちなみにcausality関数は実際にはこう動いてます（汗）
2015/8/6 34
function (x, cause = NULL, vcov. = NULL, boot = FALSE, boot.runs = 100)
{
if (!(class(x) == "varest")) {
stop("?nPlease provide an object of class 'varest', generated by 'var()'.?n")
}
K <- x$K
p <- x$p
obs <- x$obs
type <- x$type
obj.name <- deparse(substitute(x))
y <- x$y
y.names <- colnames(x$y)
if (is.null(cause)) {
cause <- y.names[1]
warning("?nArgument 'cause' has not been specified;?nusing first variable in 'x$y' (",
cause, ") as cause variable.?n")
}
else {
if (!all(cause %in% y.names))
stop("Argument cause does not match variables names.?n")
}
y1.names <- subset(y.names, subset = y.names %in% cause)
y2.names <- subset(y.names, subset = !(y.names %in% cause))
Z <- x$datamat[, -c(1:K)]
xMlm <- toMlm(x)
PI <- coef(xMlm)
PI.vec <- as.vector(PI)
R2 <- matrix(0, ncol = ncol(PI), nrow = nrow(PI))
g <- which(gsub("??.l[[:digit:]]", "", rownames(PI)) %in%
cause)
j <- which(colnames(PI) %in% cause)
R2[g, -j] <- 1
w <- which(as.vector(R2) != 0)
N <- length(w)
R <- matrix(0, ncol = ncol(PI) * nrow(PI), nrow = N)
for (i in 1:N) R[i, w[i]] <- 1
sigma.pi <- if (is.null(vcov.))
vcov(xMlm)
else if (is.function(vcov.))
vcov.(xMlm)
else vcov.
df1 <- p * length(y1.names) * length(y2.names)
df2 <- K * obs - length(PI)
STATISTIC <- t(R %*% PI.vec) %*% solve(R %*% sigma.pi %*%
t(R)) %*% R %*% PI.vec/N
if (boot) {
co.names <- Bcoef(x)
k <- which(gsub("??.l[[:digit:]]", "", colnames(co.names)) %in%
cause)
l <- which(rownames(co.names) %in% cause)
R2inv <- matrix(1, ncol = nrow(PI), nrow = ncol(PI))
R2inv[-l, k] <- 0
xres <- restrict(x, method = "man", resmat = R2inv)
pred <- sapply(xres$varresult, predict)
res <- residuals(xres)
if (is.null(vcov.)) {
Zmlm <- model.matrix(xMlm)
cross <- crossprod(Zmlm)
inside <- solve(R %*% sigma.pi %*% t(R))
boot.fun <- function(x = 1) {
Ynew <- pred + res * rnorm(n = obs, mean = 0,
sd = x)
PI.boot <- solve(cross, crossprod(Zmlm, Ynew))
PI.boot.vec <- as.vector(PI.boot)
t(R %*% PI.boot.vec) %*% inside %*% (R %*% PI.boot.vec)/N
}
}
else {
X <- x$datamat
if (x$type %in% c("const", "both"))
X <- X[, -grep("const", colnames(X))]
boot.fun <- function(x = 1) {
X[, 1:K] <- pred + res * rnorm(n = obs, sd = x,
mean = 1)
xMlm.boot <- update(xMlm, . ~ .)
sigma.pi.boot <- if (is.function(vcov.))
vcov.(xMlm.boot)
else {
vcov.
warning("vcov. should be function, not an object, when used with boot=TRUE")
}
PI.boot.vec <- as.vector(coef(xMlm.boot))
t(R %*% PI.boot.vec) %*% solve(R %*% sigma.pi.boot %*%
t(R)) %*% R %*% PI.boot.vec/N
}
}
res.rep <- replicate(boot.runs, boot.fun(x = 1))
pval <- mean(res.rep > as.numeric(STATISTIC))
}
names(STATISTIC) <- "F-Test"
if (!boot) {
PARAMETER1 <- df1
PARAMETER2 <- df2
names(PARAMETER1) <- "df1"
names(PARAMETER2) <- "df2"
PVAL <- 1 - pf(STATISTIC, PARAMETER1, PARAMETER2)
PARAM <- c(PARAMETER1, PARAMETER2)
}
else {
PARAMETER1 <- boot.runs
names(PARAMETER1) <- "boot.runs"
PVAL <- pval
PARAM <- PARAMETER1
}
METHOD <- paste("Granger causality H0:", paste(y1.names,
collapse = " "), "do not Granger-cause", paste(y2.names,
collapse = " "))
result1 <- list(statistic = STATISTIC, parameter = PARAM,
p.value = PVAL, method = METHOD, data.name = paste("VAR object",
obj.name))
class(result1) <- "htest"
sigma.u <- crossprod(resid(x))/(obs - ncol(Z))
colnames(sigma.u) <- y.names
rownames(sigma.u) <- y.names
select <- sigma.u[rownames(sigma.u) %in% y2.names, colnames(sigma.u) %in%
y1.names]
sig.vech <- sigma.u[lower.tri(sigma.u, diag = TRUE)]
index <- which(sig.vech %in% select)
N <- length(index)
Cmat <- matrix(0, nrow = N, ncol = length(sig.vech))
for (i in 1:N) {
Cmat[i, index[i]] <- 1
}
Dmat <- .duplicate(K)
Dinv <- ginv(Dmat)
lambda.w <- obs %*% t(sig.vech) %*% t(Cmat) %*% solve(2 *
Cmat %*% Dinv %*% kronecker(sigma.u, sigma.u) %*% t(Dinv) %*%
t(Cmat)) %*% Cmat %*% sig.vech
STATISTIC <- lambda.w
names(STATISTIC) <- "Chi-squared"
PARAMETER <- N
names(PARAMETER) <- "df"
PVAL <- 1 - pchisq(STATISTIC, PARAMETER)
METHOD <- paste("H0: No instantaneous causality between:",
paste(y1.names, collapse = " "), "and", paste(y2.names,
collapse = " "))
result2 <- list(statistic = STATISTIC, parameter = PARAMETER,
p.value = PVAL, method = METHOD, data.name = paste("VAR object",
obj.name))
class(result2) <- "htest"
result2
return(list(Granger = result1, Instant = result2))
}
<environment: namespace:vars>
causality関数のR内部実装（読めば色々理解できるかと思います）](https://image.slidesharecdn.com/causality20150806tjo-150806033703-lva1-app6891/85/Granger-2015-34-320.jpg)
Granger因果による 時系列データの因果推定（因果フェス2015）
- 1. Granger因果による 時系列データの因果推定 株式会社リクルートコミュニケーションズ 尾崎 隆 (Takashi J. OZAKI, Ph. D.) 2015/8/6 1
- 3. 今更ですが、自己紹介を… ? 実はGranger因果を使った論文を大昔書いたことがあります 2015/8/6 3 (Ozaki, PLoS One, 2011)
- 7. 下準備 2015/8/6 7 サンプルデータとして{RFinanceYJ}で株価時系列を取ってくる > library(RFinanceYJ) > st1<-quoteStockTsData('4751.t',since='2012-01-01',date.end='2013-07-31') > st2<-quoteStockTsData('2432.t',since='2012-01-01',date.end='2013-07-31') > st3<-quoteStockTsData('3632.t',since='2012-01-01',date.end='2013-07-31') > st4<-quoteStockTsData('4755.t',since='2012-01-01',date.end='2013-07-31') > st5<-quoteStockTsData('3793.t',since='2012-01-01',date.end='2013-07-31') > stock1<-ts(st1$adj_close) > stock2<-ts(st2$adj_close) > stock3<-ts(st3$adj_close) > stock4<-ts(st4$adj_close) > stock5<-ts(st5$adj_close)
- 10. 2015/8/6 10 ※面倒なので単変量時系列、即ち AR, MA, ARIMA過程の説明は飛ばします 自己回帰の発想は後で話します
- 11. VAR (Vector Auto Regressive)モデルとは ? VAR(p)過程 ? ベクトル自己回帰過程 ? ラグ次数pと各係数について ? 「見かけ上無相関なモデル」(SUR)で各説明変数が相互に共 通なので推定は容易、ただOLSを解くだけ（AICで次数選択） ? 多変量時系列モデリングの基本 ? AR過程を並べるだけでパラメータいっぱい ? だいたいの多変量時系列はこれで十分表現できる 2015/8/6 11
- 12. Rでは{vars}パッケージで実践できる ? とりあえず株価3系列のVARモデルを推定してみる ※ある程度スパイクの少ない3系列に限定しました 2015/8/6 12 > vstock<-cbind(stock2,stock3,stock4) ?株価3系列のマトリクスを組む > VARselect(vstock,lag.max=5) ?VAR次数pを推定 $selection AIC(n) HQ(n) SC(n) FPE(n) 1 1 1 1 ?VAR次数p = 1ということに $criteria 1 2 3 4 5 AIC(n) 2.245074e+01 2.247979e+01 2.249565e+01 2.252026e+01 2.251009e+01 HQ(n) 2.249960e+01 2.256531e+01 2.261782e+01 2.267908e+01 2.270557e+01 SC(n) 2.257395e+01 2.269543e+01 2.280370e+01 2.292071e+01 2.300296e+01 FPE(n) 5.626414e+09 5.792362e+09 5.885075e+09 6.031926e+09 5.971320e+09 > vstock.var<-VAR(vstock,p=1) ?VAR(){vars}関数で次数p=1を与えて係数推定
- 13. {vars}：多変量時系列をモデリング?予測?因果性の推定 2015/8/6 13 > summary(vstock.var) VAR Estimation Results: ========================= Endogenous variables: stock2, stock3, stock4 Deterministic variables: const Sample size: 389 Log Likelihood: -6010.006 Roots of the characteristic polynomial: 0.9926 0.9926 0.9615 Call: VAR(y = vstock, p = 1) Estimation results for equation stock2: ======================================= stock2 = stock2.l1 + stock3.l1 + stock4.l1 + const Estimate Std. Error t value Pr(>|t|) stock2.l1 0.96768 0.01167 82.953 < 2e-16 *** stock3.l1 -0.01560 0.01005 -1.552 0.121417 stock4.l1 -0.08362 0.02780 -3.008 0.002800 ** const 172.01591 48.83806 3.522 0.000479 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 80.51 on 385 degrees of freedom Multiple R-Squared: 0.953, Adjusted R-squared: 0.9526 F-statistic: 2601 on 3 and 385 DF, p-value: < 2.2e-16 （以下省略） VAR係数の概要はsummary()で
- 14. {vars}：多変量時系列をモデリング?予測?因果性の推定 2015/8/6 14 > plot(vstock.var) VARモデルの推定残差?自己相関?偏自己相関はplot.varest()で
- 15. {vars}：多変量時系列をモデリング?予測?因果性の推定 2015/8/6 15 > vstock.prd<-predict(vstock.var,n.ahead=100,ci=0.95) > plot(vstock.prd) VARモデルの予測はpredict.varest()で
- 17. Granger因果を考え出した人 2015/8/6 17 Clive W. J. Granger (1934-2009) 2003年ノーベル経済学賞 *Photo from Wikimedia Commons **実は授賞理由は共和分に関する理論的研究の成果です
- 18. Granger因果の基本的な発想 2015/8/6 18 時系列Aの未来値予測が、 時系列Bを予測式に入れることに よって改善する場合、 「Grangerの意味でB?Aなる因果関係にある」 と言う
- 32. Granger因果性検定はcausality()で 2015/8/6 32 > dstock<-diff(vstock,lag=1) ?Granger因果性検定なので1階階差を取る > VARselect(dstock,lag.max=5) > dstock.var<-VAR(dstock,p=1) > causality(dstock.var,cause=“stock2”) ?stock2からの因果性を見る $Granger Granger causality H0: stock2 do not Granger-cause stock3 stock4 data: VAR object dstock.var F-Test = 1.8703, df1 = 2, df2 = 1152, p-value = 0.1545 $Instant H0: No instantaneous causality between: stock2 and stock3 stock4 data: VAR object dstock.var Chi-squared = 124.7876, df = 2, p-value < 2.2e-16
- 33. 一応有意なGranger因果が得られる例を… 2015/8/6 33 > data(Canada) ?{vars}パッケージ同梱の”Canada”データより > var.2c <- VAR(Canada, p = 2, type = “const”) > causality(var.2c, cause = “e”)?”e”からの因果性を見る $Granger Granger causality H0: e do not Granger-cause prod rw U data: VAR object var.2c F-Test = 6.2768, df1 = 6, df2 = 292, p-value = 3.206e-06 $Instant H0: No instantaneous causality between: e and prod rw U data: VAR object var.2c Chi-squared = 26.0685, df = 3, p-value = 9.228e-06
- 34. ちなみにcausality関数は実際にはこう動いてます（汗） 2015/8/6 34 function (x, cause = NULL, vcov. = NULL, boot = FALSE, boot.runs = 100) { if (!(class(x) == "varest")) { stop("?nPlease provide an object of class 'varest', generated by 'var()'.?n") } K <- x$K p <- x$p obs <- x$obs type <- x$type obj.name <- deparse(substitute(x)) y <- x$y y.names <- colnames(x$y) if (is.null(cause)) { cause <- y.names[1] warning("?nArgument 'cause' has not been specified;?nusing first variable in 'x$y' (", cause, ") as cause variable.?n") } else { if (!all(cause %in% y.names)) stop("Argument cause does not match variables names.?n") } y1.names <- subset(y.names, subset = y.names %in% cause) y2.names <- subset(y.names, subset = !(y.names %in% cause)) Z <- x$datamat[, -c(1:K)] xMlm <- toMlm(x) PI <- coef(xMlm) PI.vec <- as.vector(PI) R2 <- matrix(0, ncol = ncol(PI), nrow = nrow(PI)) g <- which(gsub("??.l[[:digit:]]", "", rownames(PI)) %in% cause) j <- which(colnames(PI) %in% cause) R2[g, -j] <- 1 w <- which(as.vector(R2) != 0) N <- length(w) R <- matrix(0, ncol = ncol(PI) * nrow(PI), nrow = N) for (i in 1:N) R[i, w[i]] <- 1 sigma.pi <- if (is.null(vcov.)) vcov(xMlm) else if (is.function(vcov.)) vcov.(xMlm) else vcov. df1 <- p * length(y1.names) * length(y2.names) df2 <- K * obs - length(PI) STATISTIC <- t(R %*% PI.vec) %*% solve(R %*% sigma.pi %*% t(R)) %*% R %*% PI.vec/N if (boot) { co.names <- Bcoef(x) k <- which(gsub("??.l[[:digit:]]", "", colnames(co.names)) %in% cause) l <- which(rownames(co.names) %in% cause) R2inv <- matrix(1, ncol = nrow(PI), nrow = ncol(PI)) R2inv[-l, k] <- 0 xres <- restrict(x, method = "man", resmat = R2inv) pred <- sapply(xres$varresult, predict) res <- residuals(xres) if (is.null(vcov.)) { Zmlm <- model.matrix(xMlm) cross <- crossprod(Zmlm) inside <- solve(R %*% sigma.pi %*% t(R)) boot.fun <- function(x = 1) { Ynew <- pred + res * rnorm(n = obs, mean = 0, sd = x) PI.boot <- solve(cross, crossprod(Zmlm, Ynew)) PI.boot.vec <- as.vector(PI.boot) t(R %*% PI.boot.vec) %*% inside %*% (R %*% PI.boot.vec)/N } } else { X <- x$datamat if (x$type %in% c("const", "both")) X <- X[, -grep("const", colnames(X))] boot.fun <- function(x = 1) { X[, 1:K] <- pred + res * rnorm(n = obs, sd = x, mean = 1) xMlm.boot <- update(xMlm, . ~ .) sigma.pi.boot <- if (is.function(vcov.)) vcov.(xMlm.boot) else { vcov. warning("vcov. should be function, not an object, when used with boot=TRUE") } PI.boot.vec <- as.vector(coef(xMlm.boot)) t(R %*% PI.boot.vec) %*% solve(R %*% sigma.pi.boot %*% t(R)) %*% R %*% PI.boot.vec/N } } res.rep <- replicate(boot.runs, boot.fun(x = 1)) pval <- mean(res.rep > as.numeric(STATISTIC)) } names(STATISTIC) <- "F-Test" if (!boot) { PARAMETER1 <- df1 PARAMETER2 <- df2 names(PARAMETER1) <- "df1" names(PARAMETER2) <- "df2" PVAL <- 1 - pf(STATISTIC, PARAMETER1, PARAMETER2) PARAM <- c(PARAMETER1, PARAMETER2) } else { PARAMETER1 <- boot.runs names(PARAMETER1) <- "boot.runs" PVAL <- pval PARAM <- PARAMETER1 } METHOD <- paste("Granger causality H0:", paste(y1.names, collapse = " "), "do not Granger-cause", paste(y2.names, collapse = " ")) result1 <- list(statistic = STATISTIC, parameter = PARAM, p.value = PVAL, method = METHOD, data.name = paste("VAR object", obj.name)) class(result1) <- "htest" sigma.u <- crossprod(resid(x))/(obs - ncol(Z)) colnames(sigma.u) <- y.names rownames(sigma.u) <- y.names select <- sigma.u[rownames(sigma.u) %in% y2.names, colnames(sigma.u) %in% y1.names] sig.vech <- sigma.u[lower.tri(sigma.u, diag = TRUE)] index <- which(sig.vech %in% select) N <- length(index) Cmat <- matrix(0, nrow = N, ncol = length(sig.vech)) for (i in 1:N) { Cmat[i, index[i]] <- 1 } Dmat <- .duplicate(K) Dinv <- ginv(Dmat) lambda.w <- obs %*% t(sig.vech) %*% t(Cmat) %*% solve(2 * Cmat %*% Dinv %*% kronecker(sigma.u, sigma.u) %*% t(Dinv) %*% t(Cmat)) %*% Cmat %*% sig.vech STATISTIC <- lambda.w names(STATISTIC) <- "Chi-squared" PARAMETER <- N names(PARAMETER) <- "df" PVAL <- 1 - pchisq(STATISTIC, PARAMETER) METHOD <- paste("H0: No instantaneous causality between:", paste(y1.names, collapse = " "), "and", paste(y2.names, collapse = " ")) result2 <- list(statistic = STATISTIC, parameter = PARAMETER, p.value = PVAL, method = METHOD, data.name = paste("VAR object", obj.name)) class(result2) <- "htest" result2 return(list(Granger = result1, Instant = result2)) } <environment: namespace:vars> causality関数のR内部実装（読めば色々理解できるかと思います）
- 36. 長所：時系列データのみから仮説フリーに因果性を推論できる ? 元々は計量経済学上の要請 ? 従来は複雑な経済?ファイナンス理論的な仮説に基づいてAとB との間に因果関係があるかないかを論じてきた ? だが、それでは時系列データのみに着目して計量分析を行う 利点がない ? データドリブン＆仮説フリーではやれないの？ ? Grangerが打ち建てたのは「何の理論にも基づかない、仮説 フリーで予測を基準とする因果性」のアイデア 2015/8/6 36
- 37. 短所：真の意味での「因果」ではない ? Granger因果性は通常の因果性が存在する 必要条件であるが、十分条件ではない ? VARモデルに依存する以上、あくまでも「時間差」（ラグ）に着目した指標 ? 時系列上のラグが不明確なケースでは使えない ? 逆に言えば、ラグの存在する時系列同士ならGranger因果性検定さえ通過して しまえば何でも因果関係があることになりかねない ? 例えば、「稲妻が光る?雷鳴が轟く」 ? 確かに「稲妻が光」ってから「雷鳴は轟」く ? だが雷鳴が轟くのは稲妻が光ることそれ自体によるものではない ? あくまでも高電圧放電が空気中を伝わることで空気が瞬時に高温に熱せられ、 その空気の振動により雷鳴が発生する ? 真に因果関係を持つのは「高電圧放電による空気の過熱?雷鳴」であり、「稲 光」ではない 2015/8/6 37
- 38. 短所：単位根過程に対してGranger因果は使えない ? 単位根過程（ランダムウォークを含む非定常な時系列）を含むVARモ デルのOLS推定では、F検定が使えない（漸近的にF分布から外れる） ? 沖本本p.122 ? このためF検定に依拠するGranger因果性検定は、単位根過程を含 むVARモデルに対しては適用できない ? 一階階差VARモデルを用いるという解決策もあるが、その原理を鑑みるに本 質的にその感度は下がると見るべき（株価データでの実践例は実は非推奨 だったりする…） ? さらに時系列A、Bが共和分の関係にある場合、一階階差VARモデル 自体が定義できないため（共和分の定義による）、Granger因果性 検定はそもそも成立しない ? それらの場合はインパルス応答など別の手法を用いる必要がある 2015/8/6 38
- 41. 第三者効果？ 2015/8/6 41 A?Bなる因果関係がある際に、 A?C?Bなる因果関係があるとすると、 A?Bなる因果関係は実はないのではないか？ ?偏相関のアナロジーでキャンセルしたい
- 43. 最後に 2015/8/6 43 計量時系列分析の世界には因果推論の手法として、他にも ? インパルス応答 ? 分散分解 などの手法がありますが、これらはおいおい… おしまい