Augmented Dickey-Fuller Test / Enders procedure when the data generating process is unknown (e.g., inclusion of deterministic trend, or not)

proc_unitroots.R

####################################
#                                  #
# UNIT ROOT TESTS                  #
# IN THE ABSENCE OF                #
# KNOWLEDGE ON THE ACTUAL          #
# DATA GENERATING PROCESS          #
#                                  #
# Enders 1995,                     #
# Applied Econometric Time Series  #
# pp. 254 - 258 and                #
# pp. 221 - 224                    #
####################################

# builds on R package urca for
# ADF tests
library(urca)

# borrowed from package dynamac (first-differencing)
dshift <- function(x, lags = 1) {
  stopifnot(is.numeric(x))
  if (lags == 0) {
    return(x)
  }
  out <- x
  for (i in 1:lags) {
    out <- c(NA, diff(out))
  }
  out
}

# conducts Enders procedure; returns level of integration of x (0 = no unit root; 1 = unit root)
.adf_enders <- function(x, maxlags = 6, pval = .1) {
  x <- x[!is.na(x)]
  res <- list(trend = 0, ur = 0, n = length(x), p = pval)

  if (pval == .1) signcol <- 3
  if (pval == .05) signcol <- 2
  if (pval == .01) signcol <- 1
  if (!pval %in% c(.1, .05, .01)) stop("invalid p value specified; choose .1, .05, or .01")

  critical_z <- qnorm((1 - pval / 2))

  m <- (ur.df(x, type = "trend", lags = maxlags, selectlags = c("BIC")))
  # summary(m)

  # is gamma = 0?
  gamma_is_zero <- (m@teststat[1] > m@cval[1, signcol])
  if (gamma_is_zero == F) {
    res$ur <- 0
    res$trend <- as.numeric(abs(m@testreg$coefficients["tt", 3]) > critical_z)
    return(res)
  }

  phi3_is_zero <- m@teststat[3] < m@cval[3, signcol]

  if (phi3_is_zero == F) {
    gamma_is_zero_normdist <- m@teststat[1] < (-critical_z)

    if (gamma_is_zero_normdist == F) {
      res$ur <- 0
      res$trend <- as.numeric(abs(m@testreg$coefficients["tt", 3]) > critical_z)
      return(res)
    }
    if (gamma_is_zero_normdist == T) {
      res$ur <- 1
      res$trend <- as.numeric(abs(m@testreg$coefficients["tt", 3]) > critical_z)
      return(res)
    }
  }

  # estimate without trend / test for presence of drift
  m <- (ur.df(x, type = "drift", lags = maxlags, selectlags = c("BIC")))

  # is gamma = 0?
  gamma_is_zero <- (m@teststat[1] > m@cval[1, signcol])
  if (gamma_is_zero == F) {
    res$ur <- 0
    res$trend <- 0
    return(res)
  }

  phi1_is_zero <- m@teststat[2] < m@cval[2, signcol]

  if (phi1_is_zero == F) {
    gamma_is_zero_normdist <- m@teststat[1] < (-critical_z)

    if (gamma_is_zero_normdist == F) {
      res$ur <- 0
      res$trend <- 0
      return(res)
    }
    if (gamma_is_zero_normdist == T) {
      res$ur <- 1
      res$trend <- 0
      return(res)
    }
  }

  m <- (ur.df(x, type = "none", lags = maxlags, selectlags = c("BIC")))

  gamma_is_zero <- (m@teststat[1] > m@cval[1, signcol])
  if (gamma_is_zero == F) {
    res$ur <- 0
    res$trend <- 0
    return(res)
  } else {
    res$ur <- 1
    res$trend <- 0
    return(res)
  }

  return(NA)
}

# returns level of integration of time series is up to I(2) (0 = no unit root; 1 = unit root; 2 = unit root in first-differences)
adf_enders <- function(x, ...) {
  if ((length(unique(x)) < 2) | all(is.na(x))) {
    return(NULL)
  }

  difforder <- 0:2
  res <- do.call("rbind", lapply(difforder, function(.lag) .adf_enders(dshift(x, .lag), ...)))

  out <- res[1, ]
  out$order <- 0 + out$ur
  if ((res[1, ]$ur == 1 & res[2, ]$ur == 1)) {
    out <- res[2, ]
    out$order <- 2
  }
  if ((res[1, ]$ur == 1 & res[2, ]$ur == 1 & res[3, ]$ur ==
    1)) {
    out <- res[3, ]
    out$order <- 3
  }

  return(out)
}

########################
# Running instructions #
########################

# Check for presence of unit root
set.seed(1234)
adf_enders(runif(100)) # no unit root
adf_enders(1:100 + rnorm(100)) # trend-stationary