Rの関数:PP.test {stats}

args(PP.test)

PP.test_commented <- function(x, lshort = TRUE) {
  # --- 1. 入力チェック ---
  # 単変量データであることを確認します。
  if (NCOL(x) > 1) {
    stop("x is not a vector or univariate time series")
  }
  DNAME <- deparse1(substitute(x)) # データ名を保存

  # --- 2. データの準備 (ラグデータの作成) ---
  # embed関数を使って、現在の値(yt)と1期前の値(yt1)のペアを作ります。
  # x = (x1, x2, ..., xn)
  # z[,1] = (x2, x3, ..., xn) -> yt
  # z[,2] = (x1, x2, ..., xn-1) -> yt1
  z <- embed(x, 2)
  yt <- z[, 1]
  yt1 <- z[, 2]
  n <- length(yt) # 有効サンプルサイズ (元の長さ - 1)

  # --- 3. 基礎となる回帰分析 (OLS) ---
  # 時間トレンド項 u を作成 (中心化した時間)
  u <- (1L:n) - n / 2

  # 線形回帰モデル: y_t = alpha + beta * t + rho * y_{t-1} + e_t
  # ここでの係数 rho (yt1の係数) が 1 かどうかを検定します。
  res <- lm(yt ~ 1 + u + yt1)

  if (res$rank < 3) {
    stop("singularities in regression")
  } # 共線性などで計算不能な場合

  # --- 4. ディッキー・フラー統計量の計算 ---
  # 回帰係数と標準誤差を取得
  cf <- coef(summary(res))

  # yt1の係数 (cf[3, 1]) に対するt統計量を計算
  # 帰無仮説は係数 = 1 なので、(係数 - 1) / 標準誤差
  # これが通常のt値ですが、PP検定ではこれを後で補正します。
  tstat <- (cf[3, 1] - 1) / cf[3, 2]

  # 残差 u の取得
  u <- residuals(res)
  # 残差の分散 (単純な分散)
  ssqru <- sum(u^2) / n

  # --- 5. Newey-West 推定量による長期分散の推定 ---
  # 打ち切りラグパラメータ l の決定
  # lshort=TRUE なら 4*(n/100)^0.25, FALSE なら 12*(n/100)^0.25
  l <- if (lshort) {
    trunc(4 * (n / 100)^0.25)
  } else {
    trunc(12 * (n / 100)^0.25)
  }

  # C言語の関数 .Call(C_pp_sum, u, l) を呼び出して、
  # 重み付き自己共分散の和を計算しています。
  # ssqrtl は「長期分散 (long-run variance)」の推定値です。
  ssqrtl <- ssqru + .Call(stats:::C_pp_sum, u, l)

  # --- 6. 統計量の補正 (Phillips-Perron Z統計量の導出) ---
  # ここからの計算は、通常のt統計量(tstat)を、誤差項の自己相関や
  # 不均一分散を考慮して修正するための数式です。

  n2 <- n^2

  # 補正項の計算に必要な要素 (回帰行列の行列式に関連する項など)
  trm1 <- n2 * (n2 - 1) * sum(yt1^2) / 12
  trm2 <- n * sum(yt1 * (1L:n))^2
  trm3 <- n * (n + 1) * sum(yt1 * (1L:n)) * sum(yt1)
  trm4 <- (n * (n + 1) * (2 * n + 1) * sum(yt1)^2) / 6
  Dx <- trm1 - trm2 + trm3 - trm4 # 行列式に相当する値

  # 統計量 STAT (Z_t) の計算
  # 第1項: tstat * (単純分散 / 長期分散)^(1/2) -> 分散の違いでスケーリング
  # 第2項: 自己相関によるバイアスを差し引く補正項
  STAT <- sqrt(ssqru) / sqrt(ssqrtl) * tstat - (n^3) / (4 * sqrt(3) *
    sqrt(Dx) * sqrt(ssqrtl)) * (ssqrtl - ssqru)

  # --- 7. P値の計算 ---
  # ディッキー・フラー分布の近似表
  # このテーブルは有限サンプルでのクリティカル値を補間するために使われます。
  table <- cbind(c(4.38, 4.15, 4.04, 3.99, 3.98, 3.96), c(
    3.95,
    3.8, 3.73, 3.69, 3.68, 3.66
  ), c(
    3.6, 3.5, 3.45, 3.43,
    3.42, 3.41
  ), c(3.24, 3.18, 3.15, 3.13, 3.13, 3.12), c(
    1.14,
    1.19, 1.22, 1.23, 1.24, 1.25
  ), c(
    0.8, 0.87, 0.9, 0.92,
    0.93, 0.94
  ), c(0.5, 0.58, 0.62, 0.64, 0.65, 0.66), c(
    0.15,
    0.24, 0.28, 0.31, 0.32, 0.33
  ))
  table <- -table # 統計量は負の値をとるため符号反転

  tablen <- dim(table)[2L]
  tableT <- c(25, 50, 100, 250, 500, 1e+05) # サンプルサイズの基準点
  tablep <- c(0.01, 0.025, 0.05, 0.1, 0.9, 0.95, 0.975, 0.99) # 確率の基準点

  # 現在のサンプルサイズ n に合わせてクリティカル値を補間
  tableipl <- numeric(tablen)
  for (i in (1L:tablen)) {
    tableipl[i] <- approx(tableT, table[
      ,
      i
    ], n, rule = 2)$y
  }

  # 計算された統計量 STAT に対応する P値を補間
  PVAL <- approx(tableipl, tablep, STAT, rule = 2)$y

  # --- 8. 結果の出力 ---
  PARAMETER <- l
  METHOD <- "Phillips-Perron Unit Root Test"
  names(STAT) <- "Dickey-Fuller"
  names(PARAMETER) <- "Truncation lag parameter"

  structure(list(
    statistic = STAT, parameter = PARAMETER, p.value = PVAL,
    method = METHOD, data.name = DNAME
  ), class = "htest")
}