svars: Data-Driven Identification of SVAR Models

Documented in chow.test

#' Chow Test for Structural Break
#'
#' The Chow test for structural change is implemented as sample-split and break-point test (see Luetkepohl and Kraetzig, 2004, p. 135). An estimated VAR model and the presupposed structural break need to be provided.
#'
#' @param x An object of class 'vars', 'vec2var', 'nlVar'. Estimated VAR object. Or an object of class 'chowpretest' from stability()
#' @param SB Integer, vector or date character. The structural break is specified either by an integer (number of observations in the pre-break period),
#'                    a vector of ts() frequencies if a ts object is used in the VAR or a date character. If a date character is provided, either a date vector containing the whole time line
#'                    in the corresponding format or common time parameters need to be provided
#' @param nboot Integer. Number of bootstrap iterations to calculate quantiles and p-values
#' @param start Character. Start of the time series (only if dateVector is empty)
#' @param end Character. End of the time series (only if dateVector is empty)
#' @param frequency Character. Frequency of the time series (only if dateVector is empty)
#' @param format Character. Date format (only if dateVector is empty)
#' @param dateVector Vector. Vector of time periods containing SB in corresponding format
#'
#' @return A list of class "chow" with elements
#' \item{lambda_bp}{Test statistic of the Chow test with break point}
#' \item{testcrit_bp}{Critical value of the test statistic lambda_bp}
#' \item{p.value_bp}{p-value of the test statistic lambda_bp}
#' \item{lambda_sp}{Test statistic of the Chow test with sample split}
#' \item{testcrit_sp}{Critical value of the test statistic lambda_sp}
#' \item{p.value_sp}{p-value of the test statistic lambda_sp}
#' \item{SB}{Structural break tested}
#' \item{SBcharacter}{Structural break tested as character}
#' \item{p}{Number of lags used}
#'
#' @references Luetkepohl, H., 2005. New introduction to multiple time series analysis, Springer-Verlag, Berlin.\cr
#'      Luetkepohl, H., Kraetzig, M., 2004. Applied time series econometrics, Cambridge University Press, Cambridge.
#'
#' @seealso \code{\link{stability}}
#'
#'@examples
#' \donttest{
#' # Testing for structural break in USA data
#' #' # data contains quartlery observations from 1965Q1 to 2008Q2
#' # assumed structural break in 1979Q3
#' # x = output gap
#' # pi = inflation
#' # i = interest rates
#' set.seed(23211)
#' v1 <- vars::VAR(USA, lag.max = 10, ic = "AIC" )
#' z1 <- chow.test(v1, SB = 59)
#' summary(z1)
#'
#' #Using stability() to find potential break point and sample split
#' x1 <- stability(v1, type = "mv-chow-test")
#' plot(x1)
#' z1.1 <- chow.test(x1)
#' summary(z1.1)
#' #Or using sample split as benchmark
#' x1$break_point <- FALSE
#' z1.1 <- chow.test(x1)
#' summary(z1.1)
#'
#' #Structural brake via Dates
#' #given that time series vector with dates is available
#' dateVector <- seq(as.Date("1965/1/1"), as.Date("2008/7/1"), "quarter")
#' z2 <- chow.test(v1, SB = "1979-07-01", format = "%Y-%m-%d", dateVector = dateVector)
#' summary(z2)
#'
#' # alternatively pass sequence arguments directly
#' z3 <- chow.test(v1, SB = "1979-07-01", format = "%Y-%m-%d",
#'                 start = "1965-01-01", end = "2008-07-01",
#'                 frequency = "quarter")
#' summary(z3)
#'
#' # or provide ts date format (For quarterly, monthly, weekly and daily frequencies only)
#' z4 <- chow.test(v1, SB = c(1979,3))
#' summary(z4)
#' }
#' @import stats
#' @importFrom utils combn
#' @importFrom vars VAR
#' @importFrom expm sqrtm
#' @importFrom methods is
#'
#' @export
#'

## Chow tests of system structural break ##
#-----------------------------------------#

# Y     : data
# SB    : structural break
# nboot : number of bootstrap iterations
# lags  : maximum lag order

chow.test <- function(x, SB, nboot = 500, start = NULL, end = NULL,
                      frequency = NULL, format = NULL, dateVector = NULL){

  if(is(x, 'chowpretest')){
    if(x$break_point == TRUE){
      SB <-  which.max(x$teststat_bp)
    }else{
      SB <-  which.max(x$teststat_sp)
    }
    x <- x$var
  }

  u <- Tob <- p <- k <- residY <- coef_x <- yOut <- type <- y <-  NULL
  get_var_objects(x)

  # Null hypothesis of no sample split is rejected for large values of lambda
  #Tob <- x$obs
  if(is.numeric(SB)){
    SBcharacter <- NULL
  }

  if(!is.numeric(SB)){
    SBcharacter <- SB
    SB <- getStructuralBreak(SB = SB, start = start, end = end,
                             frequency = frequency, format = format, dateVector = dateVector, Tob = Tob, p = p)
  }

  if(length(SB) != 1 & inherits(yOut, "ts")){
    SBts = SB
    SB = dim(window(yOut, end = SB))[1]
    if(frequency(y == 4)){
      SBcharacter = paste(SBts[1], " Q", SBts[2], sep = "")
    }else if(frequency(yOut == 12)){
      SBcharacter = paste(SBts[1], " M", SBts[2], sep = "")
    }else if(frequency(yOut == 52)){
      SBcharacter = paste(SBts[1], " W", SBts[2], sep = "")
    }else if(frequency(yOut == 365.25)){
      SBcharacter = paste(SBts[1], "-", SBts[2], "-", SBts[3], sep = "")
    }else{
      SBcharacter = NULL
    }

  }
  y <- t(y)

  Full <- y

  # splitting sample
  sample1 <- y[1:SB, ]
  sample2 <- y[(SB+1+p):nrow(y), ]

  # estimating VAR for pre and post SB and for full series
  VAR.model <- VAR(Full, p = p, type = type)
  VAR1.model <- VAR(sample1, p = p, type = type)
  VAR2.model <- VAR(sample2, p = p, type = type)

  l1 <- VAR1.model$obs
  l2 <- VAR2.model$obs
  ll <- VAR.model$obs

  # calculating four covariance matrices
  Sigma.1 <- (1/l1)*t(residuals(VAR1.model))%*%(residuals(VAR1.model))
  Sigma.2 <- (1/l2)*t(residuals(VAR2.model))%*%(residuals(VAR2.model))

  Sigma <- (1/l1)*t(residuals(VAR.model)[1:l1,])%*%(residuals(VAR.model)[1:l1,]) +
    (1/l2)*t(residuals(VAR.model)[(ll-l2+1):ll,])%*%(residuals(VAR.model)[(ll-l2+1):ll,])

  Sigma.1.2 <- (1/(l1 + l2))*(t(residuals(VAR.model)[1:l1,])%*%(residuals(VAR.model)[1:l1,]) +
                                t(residuals(VAR.model)[(ll-l2+1):ll,])%*%(residuals(VAR.model)[(ll-l2+1):ll,]))

  # calculating the test statistic
  # lambda_bp : test statistic for break point (tests for change in covariance in addition)
  # lambda_SP : test statistic for structural break (tests only for parameter change)
  lambda_bp <- (l1 + l2)*log(det(Sigma)) - l1*log(det(Sigma.1)) - l2*log(det(Sigma.2))
  lambda_sp <- (l1 + l2)*(log(det(Sigma.1.2)) - log(det((1/(l1 + l2))*(l1*Sigma.1 + l2*Sigma.2))))

if(!is.null(nboot)){
  lambda_bpB <- rep(NA, nboot)
  lambda_spB <- rep(NA, nboot)
  TB <- l1


  # bootstrapping the test statistic to obtain empirical distribution
  # obtaining VAR parameter
  #coef_x <- coef(x)

  A <- matrix(0, nrow =k, ncol = k* p)
  #yl <- t(YLagCr(t(y), p))
  yl <- t(YLagCr(y, p))

  for(i in 1:k){
    A[i,] <- coef_x[[i]][1:(k*p),1]
  }

  if(type == "const"){
    v <- rep(1, k)

    for(i in 1:k){
      v[i] <- coef_x[[i]][(k*p+1), 1]
    }

    A <- cbind(v, A)
    Z_t <- rbind(rep(1, ncol(yl)), yl)
  }else if (type == "trend"){
    trend <- rep(1, k)

    for(i in 1:k){
      trend[i] <- coef_x[[i]][(k*p+1), 1]
    }

    A <- cbind(trend, A)
    Z_t <- rbind(seq(1, ncol(yl)), yl)
  }else if(type == "both"){
    v <- rep(1, k)

    for(i in 1:k){
      v[i] <- coef_x[[i]][(k*p+1), 1]
    }

    trend <- rep(1, k)
    Z_t <- rbind(rep(1, ncol(yl)), seq(1, ncol(yl)), yl)
    for(i in 1:k){
      trend[i] <- coef_x[[i]][(k*p+2), 1]
    }

    A <- cbind(v, trend, A)
  }else{
    Z_t <- yl
  }

  residF <- residuals(VAR.model)

  # creating new data
  datFull <- list()

  for(i in 1:nboot){
    u <- residF
    my <- rnorm(n = k)
    my <- (my > 0) - (my < 0)
    et <- u * my
    Ystar <- t(A %*% Z_t + t(et))
    datFull[[i]] <- Ystar
  }

  for(i in 1:nboot){
    xx <- datFull[[i]]
    # splitting sample
    sample1 <- xx[1:SB, ]
    sample2 <- xx[(SB+1+p):nrow(xx), ]

    VARB.model <- VAR(xx, p = p)
    VARB1.model <- VAR(sample1, p = p)
    VARB2.model <- VAR(sample2, p = p)

    l1 <- VARB1.model$obs
    l2 <- VARB2.model$obs
    ll <- VARB.model$obs

    Sigma.1 <- (1/l1)*t(residuals(VARB1.model))%*%(residuals(VARB1.model))
    Sigma.2 <- (1/l2)*t(residuals(VARB2.model))%*%(residuals(VARB2.model))

    Sigma <- (1/l1)*t(residuals(VARB.model)[1:l1,])%*%(residuals(VARB.model)[1:l1,]) +
      (1/l2)*t(residuals(VARB.model)[(ll-l2+1):ll,])%*%(residuals(VARB.model)[(ll-l2+1):ll,])

    Sigma.1.2 <- (1/(l1 + l2))*(t(residuals(VARB.model)[1:l1,])%*%(residuals(VARB.model)[1:l1,]) +
                                  t(residuals(VARB.model)[(ll-l2+1):ll,])%*%(residuals(VARB.model)[(ll-l2+1):ll,]))

    lambda_bpB[i] <- (l1 + l2)*log(det(Sigma)) - l1*log(det(Sigma.1)) - l2*log(det(Sigma.2))
    lambda_spB[i] <- (l1 + l2)*(log(det(Sigma.1.2)) - log(det((1/(l1 + l2))*(l1*Sigma.1 + l2*Sigma.2))))
  }


  df_bp <- p*k^2 + k + (k*(k + 1))/2
  df_sp <- p*k^2 + k

  # obtaining critical values and p-values for both tests
  testcrit_bp <- quantile(lambda_bpB, probs = 0.95)
  EmpDist_bp <- ecdf(lambda_bpB)
  p.value_bp <- 1 - EmpDist_bp(lambda_bp)

  testcrit_sp <- quantile(lambda_spB, probs = 0.95)
  EmpDist_sp <- ecdf(lambda_spB)
  p.value_sp <- 1 - EmpDist_sp(lambda_sp)
}else{
  df_bp <- p*k^2 + k + (k*(k + 1))/2
  df_sp <- p*k^2 + k

  # obtaining critical values and p-values for both tests
  testcrit_bp <- qchisq(p = 0.95, df = df_bp)
  p.value_bp <- 1 - unname(pchisq(lambda_bp, df = df_bp))

  testcrit_sp <- qchisq(p = 0.95, df = df_sp)
  p.value_sp <- 1 - unname(pchisq(lambda_sp, df = df_sp))
}


  chowTest <- list(lambda_bp = unname(lambda_bp),      # test statistic for break point test
                   testcrit_bp = testcrit_bp,  # critical value for 95% quantile
                   p.value_bp = p.value_bp,    # p-value
                   lambda_sp = unname(lambda_sp),      # test statistic for sample split
                   testcrit_sp = testcrit_sp,  # critical value for 95% quantile
                   p.value_sp = p.value_sp,     # p-value
                   SB = SB,                     # Structural breakpoint
                   SBcharacter = SBcharacter,   # Structural Break as character
                   p = p
  )
  class(chowTest) <- "chow"

  return(chowTest)

}