餡子付゛録゛

ソフトウェア開発ツールの便利な使い方を紹介。

拡張Dickey-Fuller検定を横断面方向に拡張する

R言語

同一のデータ生成プロセスを踏んでいる複数の系列 y1, y2, y3, … を、それぞれ個別に単位根検定にかけると、定常であるのに標本数が不足して単位根の棄却に失敗しがちになる一方、c(y1, y2, y3, …) と単純に接続して単位根検定した場合、接続面に断続が生じるため単位根を過って棄却する事になります。

代表的な単位根検定のDickey-Fuller検定と、それに系列相関やトレンド項の影響も考慮したADF 検定では、差分データを前期の水準データで回帰して、前期の水準データの係数の標準誤差をτ分布で評価します。具体的には単位根が疑われる系列 Y_tを、\epsilonを誤差項として、

Y_t = \beta Y_{t-1} + \epsilon

と書き、以下のように差分 \Delta Y_tを説明する式に変形し、

Y_t - Y_{t-1} = \beta Y_{t-1} - Y_{t-1} + \epsilon
\Delta Y_t = (\beta - 1) Y_{t-1} + \epsilon = \gamma Y_{t-1} + \epsilon

\gamma = 0であるかを検定し、棄却されれば単位根が無いとします。ここで異なる系列でも同一のデータ生成プロセスであれば、前期の水準データ Y_{t-1}の係数 \gammaは同一となります。つまり、系列ごとに差分データと前期の水準データの組み合わせを作成後にデータを統合すれば、複数の系列を同時に単位根検定する事ができます。

一般に時系列方向と横断面方向にデータの広がりのあるパネルデータの場合、横断面方向にデータの広がりがあるため、単一系列の場合と異なり見せかけの相関が生じにくく、株価などの不確実性の高い高頻度取引などではなく、観察対象が低頻度に理屈を立てて行動した結果でああるため、単位根の存在は気にしない事が多い*1わけですが、たまたま稀な需要があったのでurcaパッケージのur.dfを拡張してみました。

1. ソースコード

urcaを呼び出したあと、関数ur.dfを上書きします。ちょっと長いですがデータ統合する部分以外は、元と同じです。ですから、元と同じオプションが使えます。元のが使いたくなった場合は、環境を指定するか、追加定義したグローバル領域のur.dfをrm(ur.df)して消してください。

library(urca)
ur.df <- function (..., type = c("none", "drift", "trend"), lags = 1, selectlags = c("Fixed", "AIC", "BIC"))
{
  lag <- as.integer(lags)
  if (lag < 0)
    stop("\nLags must be set to an non negative integer value.\n")
  lags <- lags + 1

  args <- list(...)
  len <- length(args)
  if(len <= 0)
    stop("\ny must be set.\n")
  tmp <- as.list(NULL)
  tmp[["n"]] <- 0

  for(i in 1:len){
    y <- args[[i]]
    y <- y[!is.na(y)]

    if (ncol(as.matrix(y)) > 1)
      stop("\ny is not a vector or univariate time series.\n")
    if (any(is.na(y)))
      stop("\nNAs in y.\n")

    y <- as.vector(y)
    z <- diff(y)
    n <- length(z)
    x <- embed(z, lags)
    z.diff <- x[, 1]
    z.lag.1 <- y[lags:n]

    tmp[["y"]] <- append(tmp[["y"]], y)
    tmp[["z"]] <- append(tmp[["z"]], z)
    tmp[["n"]] <- tmp[["n"]] + n
    tmp[["x"]] <- rbind(tmp[["x"]], x)
    tmp[["z.diff"]] <- append(tmp[["z.diff"]], z.diff)
    tmp[["z.lag.1"]] <- append(tmp[["z.lag.1"]], z.lag.1)
    tmp[["tt"]] <- append(tmp[["tt"]], lags:n)
  }

  y <- tmp[["y"]]
  z <- tmp[["z"]]
  n <- tmp[["n"]]
  x <- tmp[["x"]]
  z.diff <- tmp[["z.diff"]]
  z.lag.1 <- tmp[["z.lag.1"]]
  tt <- tmp[["tt"]]

  (n-length(z.lag.1)+1):n

  selectlags <- match.arg(selectlags)
  type <- match.arg(type)

  CALL <- match.call()
  DNAME <- deparse(substitute(y))
  x.name <- deparse(substitute(y))

  if (lags > 1) {
    if (selectlags != "Fixed") {
      critRes <- rep(NA, lags)
      for (i in 2:(lags)) {
        z.diff.lag = x[, 2:i]
        if (type == "none")
         result <- lm(z.diff ~ z.lag.1 - 1 + z.diff.lag)
        if (type == "drift")
         result <- lm(z.diff ~ z.lag.1 + 1 + z.diff.lag)
        if (type == "trend")
         result <- lm(z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
        critRes[i] <- AIC(result, k = switch(selectlags,
         AIC = 2, BIC = log(length(z.diff))))
      }
      lags <- which.min(critRes)
    }
    z.diff.lag = x[, 2:lags]
    if (type == "none") {
      result <- lm(z.diff ~ z.lag.1 - 1 + z.diff.lag)
      tau <- coef(summary(result))[1, 3]
      teststat <- as.matrix(tau)
      colnames(teststat) <- "tau1"
    }
    if (type == "drift") {
      result <- lm(z.diff ~ z.lag.1 + 1 + z.diff.lag)
      tau <- coef(summary(result))[2, 3]
      phi1.reg <- lm(z.diff ~ -1 + z.diff.lag)
      phi1 <- anova(phi1.reg, result)$F[2]
      teststat <- as.matrix(t(c(tau, phi1)))
      colnames(teststat) <- c("tau2", "phi1")
    }
    if (type == "trend") {
      result <- lm(z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
      tau <- coef(summary(result))[2, 3]
      phi2.reg <- lm(z.diff ~ -1 + z.diff.lag)
      phi3.reg <- lm(z.diff ~ z.diff.lag)
      phi2 <- anova(phi2.reg, result)$F[2]
      phi3 <- anova(phi3.reg, result)$F[2]
      teststat <- as.matrix(t(c(tau, phi2, phi3)))
      colnames(teststat) <- c("tau3", "phi2", "phi3")
    }
  }
  else {
    if (type == "none") {
      result <- lm(z.diff ~ z.lag.1 - 1)
      tau <- coef(summary(result))[1, 3]
      teststat <- as.matrix(tau)
      colnames(teststat) <- "tau1"
    }
    if (type == "drift") {
      result <- lm(z.diff ~ z.lag.1 + 1)
      phi1.reg <- lm(z.diff ~ -1)
      phi1 <- anova(phi1.reg, result)$F[2]
      tau <- coef(summary(result))[2, 3]
      teststat <- as.matrix(t(c(tau, phi1)))
      colnames(teststat) <- c("tau2", "phi1")
    }
    if (type == "trend") {
      result <- lm(z.diff ~ z.lag.1 + 1 + tt)
      phi2.reg <- lm(z.diff ~ -1)
      phi3.reg <- lm(z.diff ~ 1)
      phi2 <- anova(phi2.reg, result)$F[2]
      phi3 <- anova(phi3.reg, result)$F[2]
      tau <- coef(summary(result))[2, 3]
      teststat <- as.matrix(t(c(tau, phi2, phi3)))
      colnames(teststat) <- c("tau3", "phi2", "phi3")
    }
  }
  rownames(teststat) <- "statistic"
  testreg <- summary(result)
  res <- residuals(testreg)
  if (n < 25)
    rowselec <- 1
  if (25 <= n & n < 50)
    rowselec <- 2
  if (50 <= n & n < 100)
    rowselec <- 3
  if (100 <= n & n < 250)
    rowselec <- 4
  if (250 <= n & n < 500)
    rowselec <- 5
  if (n >= 500)
    rowselec <- 6
  if (type == "none") {
    cval.tau1 <- rbind(c(-2.66, -1.95, -1.6), c(-2.62, -1.95,
      -1.61), c(-2.6, -1.95, -1.61), c(-2.58, -1.95, -1.62),
      c(-2.58, -1.95, -1.62), c(-2.58, -1.95, -1.62))
    cvals <- t(cval.tau1[rowselec, ])
    testnames <- "tau1"
  }
  if (type == "drift") {
    cval.tau2 <- rbind(c(-3.75, -3, -2.63), c(-3.58, -2.93,
      -2.6), c(-3.51, -2.89, -2.58), c(-3.46, -2.88, -2.57),
      c(-3.44, -2.87, -2.57), c(-3.43, -2.86, -2.57))
    cval.phi1 <- rbind(c(7.88, 5.18, 4.12), c(7.06, 4.86,
      3.94), c(6.7, 4.71, 3.86), c(6.52, 4.63, 3.81), c(6.47,
      4.61, 3.79), c(6.43, 4.59, 3.78))
    cvals <- rbind(cval.tau2[rowselec, ], cval.phi1[rowselec,
      ])
    testnames <- c("tau2", "phi1")
  }
  if (type == "trend") {
    cval.tau3 <- rbind(c(-4.38, -3.6, -3.24), c(-4.15, -3.5,
      -3.18), c(-4.04, -3.45, -3.15), c(-3.99, -3.43, -3.13),
      c(-3.98, -3.42, -3.13), c(-3.96, -3.41, -3.12))
    cval.phi2 <- rbind(c(8.21, 5.68, 4.67), c(7.02, 5.13,
      4.31), c(6.5, 4.88, 4.16), c(6.22, 4.75, 4.07), c(6.15,
      4.71, 4.05), c(6.09, 4.68, 4.03))
    cval.phi3 <- rbind(c(10.61, 7.24, 5.91), c(9.31, 6.73,
      5.61), c(8.73, 6.49, 5.47), c(8.43, 6.49, 5.47),
      c(8.34, 6.3, 5.36), c(8.27, 6.25, 5.34))
    cvals <- rbind(cval.tau3[rowselec, ], cval.phi2[rowselec,
      ], cval.phi3[rowselec, ])
    testnames <- c("tau3", "phi2", "phi3")
  }
  colnames(cvals) <- c("1pct", "5pct", "10pct")
  rownames(cvals) <- testnames
  new("ur.df", y = y, model = type, cval = cvals, lags = lag,
    teststat = teststat, testreg = testreg, res = res, test.name = "Augmented Dickey-Fuller Test")
}

2. 使い方

多少はと言った程度ですが、個別に単位根検定をかけるよりは、系列をまとめて検定する方が正しくなります。

genS <- function(n=200, beta=1.0){
  x <- numeric(n)
  x[1] <- 0
  for(i in 2:n){
    x[i] <- beta*x[i-1] + rnorm(1)
  }
  x
}

set.seed(121) # 数字を固定しないと結果はマチマチ
y1 <- genS(30) # 非定常過程を生成
y2 <- genS(30)
y3 <- genS(30)
summary(ur.df(c(y1, y2, y3))) # 単位根を過って棄却する
summary(ur.df(y1, y2, y3)) # 単位根を棄却しない

set.seed(100) # 数字を固定しないと結果はマチマチ
y4 <- genS(30, 0.9) # 定常過程を生成
y5 <- genS(30, 0.9)
y6 <- genS(30, 0.9)
summary(ur.df(y4)) # 単位根を棄却しない
summary(ur.df(y5)) # 単位根を棄却しない
summary(ur.df(y6)) # 単位根を棄却しない
summary(ur.df(y4, y5, y6)) # 単位根を棄却する

*1: 同時性、不均一分散、系列相関はとても気にされる。