Multivariate break test

In this tutorial, I explain how to implement, in a flexible way, the algorithm of Bai, Lumsdaine, and Stock (1998).

Step 1: Lag variables.

This function takes as argument a matrix of time series observations and lags it by an order (q).

Code

compute_lags <- function(Y                                                   #time series matrix Y
                         , q)                                                #lag order q
{
  p <- dim(Y)[1]                                                             #get the dimensions
  n <- dim(Y)[2]

  myDates <- rownames(Y)[(q + 1) : p]                                        #optional: keep the rownames dates of the data frame with final matching
  Y <- data.matrix(Y)                                                        #matrix conversion

  YLAG <- matrix(data = NA, nrow <- (p - q), ncol <- (n * (q + 1)))          #create an empty matrix

  for(i in 0:q)
  {
    YLAG[ , (n * i + 1):(n * (i + 1))] <- Y[(q - i + 1):(p - i), ]
  }

  Y <- YLAG[,1:n]
  YLAG <- YLAG[,(n + 1):dim(YLAG)[2]]
  return(list(Y = Y, YLAG = YLAG, myDates = myDates))
}

Step 2: Define all matrices

This function takes as argument, a matrix Y (of multivariate time series), an order q for the VAR, an optional matrix X of contemporaneous covariates, it adds an optional trend and determines whether we test for a break at the mean level (intercept = TRUE) or for all the parameters (intercept = FALSE).

Code

matrix_conformation <- function(Y, q, X, trend, intercept)
{
  pInit <-  dim(Y)[1]                                                           #get the original number of observations

  lY <- compute_lags(Y, q)                                                      #get the list of contemporaneous and lags objects

  Y <- lY$Y                                                                     #get the matching dependent matrix
  YLAG <- lY$YLAG                                                               #get the lagged dependent variables matrix
  myDates <- lY$myDates                                                         #get the original matching rowname vector (of dates)
  
  p <- dim(Y)[1]                                                                #length of the matrix
  nEq <- dim(Y)[2]                                                              #number of equations of the vAR

  print(p)
  In <- diag(nEq)                                                               #identity matrix of the number of equations of the VAR

  Gt <- as.matrix(cbind(rep(1, p), YLAG))                                       #create a unique matrix transpose of G with one intercept and autoregressive terms
  
  n <- dim(Gt)[2]                                                               #incremental number of regressors

  if(!is.null(X))                                                               #if additional covariates matrix is passed as argument
  {
    if(pInit==dim(X)[1])                                                        #and if its size is equal to the original matrix Y
    {
      Gt <- cbind(Gt, data.matrix(X[(q + 1) : pInit, ]))                        #increment Gt by the contemporaneous covariate matrix X
      n <- dim(Gt)[2]                                                           #increment the total number of regressors
    }
    else
      print("The number of observations of X does not match the one of Y")
  }

  if(trend)                                                                     #check if we add a trend
  {
    Gt <- cbind(Gt, seq(1, p, by = 1))                                          #add trend
    n <- dim(Gt)[2]                                                             #increment the total number of regressors
  }

  if(intercept)                                                                 #if only the intercept is allowed to break
  {
    s <- t(data.matrix(c(rep(0, n))))                                           #create the selection vector
    s[1] <- 1                                                                   #vector only select the first element (intercept) to test the shift
    S <- kronecker(s, In)                                                       #create the selection matrix
  }

  if(!intercept)                                                                #full parameters structural change estimation
  {                                                                             #get the full dimension of the test
    r <- nEq * n
    S <- diag(r)                                                                #identity matrix of the size of the test is the selection matrix
  }

  G <- t(Gt)                                                                    #transpose Gt to get G as in BLS

  Yex <- data.matrix(c(t(Y)))                                                   #Expend and vectorize Y
  
  Gex <- kronecker(t(G), In)                                                    #Expend G
  return(list(Yex = Yex, Gex = Gex, p = p, G = G, S = S, myDates = myDates, Y = Y, nEq = nEq))
}

Step 3: Optimal VAR(q) order

We determine which is the optimal order of the VAR(q). For this we keep the option of using the Akaike Information Criterion (AIC) or Bayesian Information Criterion (BIC). We add the argument qMax which determines up until which lag we are going to compute the criteria.

Code

compute_aicbic <- function(Y, qMax, X, trend, intercept)           #compute the AIC and BIC criteria for lags from 1 to qMax
{
  library(stats)                                                   #load stats package
  AICBIC = matrix(data <- NA, nrow = 2, ncol = qMax)               #create empty matrix for the AIC / BIC criteria
  
  for(q in 1:qMax)
  {
    print(paste0("Testing lags number : ", q))
    
    lConfMatrix <- matrix_conformation(Y, q, X, trend, intercept) #create a list of conformed objects for the estimation
    Yex <- lConfMatrix$Yex
    Gex <- lConfMatrix$Gex
    
    mod <- lm(Yex~Gex)                                            #estimate the model with lm
    AICBIC[1,q] <- AIC(mod)                                       #get AIC
    AICBIC[2,q] <- BIC(mod)                                       #get BIC
  }
  rownames(AICBIC) <- c("AIC", "BIC")
  colnames(AICBIC) <- paste0("lags = ", 1:qMax)
  return(AICBIC)
}

Step 4: Compute VAR(q) parameters.

This function is based on the pre-computed conformed matrix. We keep the option to compute it by Ordinary Least Squares (OLS), Feasible Generalized Least Squares (FGLS) or Iterative Generalized Least Squares (IGLS). FGLS and IGLS are required if we add additional covariates that are not common to all equations as in a Seemingly Unrelated Regression (SUR). The arguments are a matrix of covariates Z, which includes original and breaking covariates. Yex is an expended vector of the dependent variables Y, nEq is the number of equations, p is the number of observations, estMode is a string which can switch between “OLS”, “FGLS” and “IGLS”. In the case of “IGLS”, iter is an integer which determines the number of iterations used.

Code

compute_beta <- function(Z, Yex, nEq, p, estMode, iter)                   #function which takes the regressor matrices Z and Yex, number of equations nEq
  #number of observations p, mode of estimation "estMode" and number of iterations
  #in the case of iterative feasible general least square
{
  if(estMode == "OLS")                                                    #if OLS
  {
    Beta <- solve(Z %*% t(Z), tol = 0) %*% Z %*% Yex                      #solve the system and get the vector of betas
  }

  if(estMode == "FGLS")                                                   #if FGLS
  {
    Beta <- solve(Z %*% t(Z), tol = 0) %*% Z %*% Yex                      #solve the system and get the vector of betas
    Sigma <- compute_sigma(Z, Yex, Beta, nEq)                             #get the covariance matrix of errors
 
    Omega <- kronecker(diag(p), Sigma)                                    #get Omega
    Beta <- solve(Z %*% solve(Omega, tol = 0) %*% t(Z), tol = 0) %*% Z %*% solve(Omega, tol = 0) %*% Yex        #solve the system and get the vector of betas
  }

  if(estMode == "IGLS")                                                   #if IGLS
  {
    Beta <- solve(Z %*% t(Z), tol = 0) %*% Z %*% Yex                      #solve the system and get the vector of betas
    
    for (i in 1:iter)
    {
      Sigma <- compute_sigma(Z, Yex, Beta, nEq)                           #get the covariance matrix of errors
      Omega <- kronecker(diag(p), Sigma)                                  #get Omega
      Beta <- solve(Z %*% solve(Omega, tol = 0) %*% t(Z), tol = 0) %*% Z %*% solve(Omega, tol = 0) %*% Yex      #solve the system and get the vector of betas
      # print(Sigma)                                                      #control for convergence / non divergence
    }
  }
  Sigma <- compute_sigma(Z, Yex, Beta, nEq)                               #final estimation of Sigma
  return(list(Beta = Beta, Sigma = Sigma))
}

Step 5: Compute Sigma

This function simply compute the errors, and their covariance matrix, generated by the VAR in OLS, FGLS or IGLS mode.

Code

compute_sigma <- function(Z, Yex, Beta, nEq)                    #compute the covariance matrix of errors as in BLS (1998)
{
  errors <- (Yex - t(Z) %*% Beta)                               #get the n*p vector of residuals
  Errors <- matrix(errors, ncol = nEq, byrow = T)               #reset as matrix to obtain Sigma
  Sigma <- cov(Errors)

  return(Sigma)
}

Step 6: Compute F-statistic

This function computes the F-statistic of the selected breaking parameters. ### Code

compute_fstat <- function(R, Beta, Z, p, Sigma)                   #compute the f-statistic as in BLS (1998)
{
  RBeta <- R %*% Beta                         #pre compute the RBeta matrix with the selected coefficients allowed to break

  if(!is.null(Sigma))                         #if the covariance matrix of error is passed as argument, compute F-stat
  {
    Omega <- kronecker(diag(p), Sigma)        #get Omega
    Fk <- p * t(RBeta) %*% solve(R %*% solve((Z %*% solve(Omega, tol = 0) %*% t(Z)) / p, tol = 0) %*% t(R), tol = 0) %*% RBeta
  }
  return(Fk)
}

Step 7: Critical values

Before stepping into the confidence interval computation, let’s write a function computing the alphas based on critical values of the limiting distribution V (see, Picard, 1985). This function takes a vector of critical values (e.g., 90%, 95% and 99%) and get the corresponding alpha levels.

Code

compute_vdistr_cv = function(ci = c(0,9, 0.95, 0,99))          #computes the critical values for a vector of confidence intervals proposed, ci
{
  u <- length(ci)                         #get the number of ci elements
  target <- 1 - (1 - ci) / 2              #redefine target for a two tail CI
  
  print(paste0("The vdistr targets are: ",target))
  x <- seq(-200, 200, 0.01)               #define the support sequence "x" for the CDF of V
  
  V <- (3 / 2) * exp(abs(x)) * pnorm( (-3 / 2) * abs(x)^0.5 )  - (1 / 2) * pnorm( (-1 / 2) * abs(x)^0.5 ) #compute V
  
  cumV <- cumsum(V)/sum(V)                #scale the CDF of V to reach one
  
  
  # dev.new()                                               
  # plot(x, cumsum(gamma)/sum(gamma), t = 'l')                    #optionally plot V (nice shape!)
  
  cv <- rep(NA, u)
  k <- 1
  
  print(target)
  for(i in 2:length(x))
  {
    if(cumV[i-1] < target[k] && cumV[i] >= target[k])
    {
      cv[k] <- x[i]
      k <- k + 1
      if(k>u)
        break
    }  
  }
  print(cv)
  return(cv)
}

Step 8: Confidence intervals

We compute the corresponding confidence intervals for a potential break date. This function uses the original G matrix of covariates, the selection matrix S, the covariance matrix Sigma and the alphas from the critical values computed above.

Code

compute_ci <- function(G, S, Sigma, R, Beta, cv, p)             #as per Bekeart Harvey Lumsdaine (2002) // recall that RBeta = S %*% deltaT
{
  u <- length(cv)                                               #get the number of critical values from the vdistr
  ciDelta <- rep(NA, u)                                         #create empty vector

  RBeta <- R %*% Beta                                           #pre compute the selection of (breaking) parameters to test

  tCi <- solve(t(RBeta) %*% S %*% kronecker((G %*% t(G)) / p, solve(Sigma, tol = 0)) %*% t(S) %*% RBeta, tol = 0)    #compute the conf interval factor

  for(i in 1:u)
  {
    ciDelta[i] <- cv[i] * tCi                                   #get the vector of critical values
  }
  # print(ciDelta)
  return(ciDelta)
}

Step 9: Graphical output

One step before launching the main function, get a nice graphical output for our results.

Code

compute_plot_stats <- function(myDates, Variables, fstat, CI, Y, meanShift)     #function to compute summary statistics and deliver plots
{
  library(ggplot2)
  library(reshape2)
  library(ggpubr)
  library(stargazer)
  
  myDates <- as.Date(myDates, format = "%d.%m.%Y")
  
  
  maxF <- which.max(fstat)                                         #get the index when the break occurs
  cis <- round(CI[maxF, ])                                         #get the confidence intervals around the break
  
                                                                   #hard code for three different confidence intervals
  startDate90 <- myDates[maxF - cis[1]]                                       
  endDate90 <- myDates[maxF + cis[1]] 
  
  startDate95 <- myDates[maxF - cis[2]]                                             
  endDate95 <- myDates[maxF + cis[2]]                            
                       
  startDate99 <- myDates[maxF - cis[3]]                                              
  endDate99 <- myDates[maxF + cis[3]]                            
  
  p <- dim(Y)[1]
  n <- dim(Y)[2]
  Y <- Y * 100                                                     #get values in percentage
  
  Y <- apply(Y, 2, as.double)
  Y <- data.frame(cbind(myDates, Y))
  
  colnames(Y)[1] <- "Date"
  colnames(Y)[2:(n + 1)] <- Variables
 
  Y[,1] <- myDates
 
  Y <- melt(Y, id.var = "Date")                                   #reshape to long format
  names(Y)[2] = "Variables"
  
  dev.new()
  g1 <- ggplot(Y, aes(x = Date, y = value, group = Variables, colour = Variables))
  g1 <- g1  + geom_line()
  
  g1 <- g1 + ggtitle("") + xlab("Date") + ylab("Intensity (%)")
  
  g1 <- g1 + coord_cartesian(ylim = c(0, max(Y$value, na.rm = T)))
  g1 <- g1 + scale_y_continuous(expand = c(0,0))                                                    #force the y axis to start at zero
  
  g1 <- g1 + scale_x_date(breaks = scales::pretty_breaks(n = 10))

  g1 <- g1 + theme_bw()
                                                                                                    #add shaded area for various ci
  g1 <- g1+annotate("rect", xmin = startDate90, xmax = endDate90, ymin = 0, ymax = Inf,
                   alpha = .6)
  g1 <- g1+annotate("rect", xmin = startDate95, xmax = endDate95, ymin = 0, ymax = Inf,
                    alpha = .4)
  g1 <- g1+annotate("rect", xmin = startDate99, xmax = endDate99, ymin = 0, ymax = Inf,
                    alpha = .2)
  
  d <- data.frame(date = myDates[maxF], event = "The event")
  
  g1 <- g1 + geom_vline(data = d, mapping = aes(xintercept = date), color = "black", size = 1)      #add the break line
  
  return(g1)
}

Step 10: Main entry point

Finally, the main function, wrapping up all others and running a loop for all potential k breaks considered This function takes all the aforementioned arguments, adding a trim parameter (in percent) to compute the burn-in and “burn-out” periods. In addition, I add a boolean parameter (posBreak), which if set to TRUE will discard all breaks arising from a decrease in intercept or full parameters.

Code

main <- function(Y                                      #Y is a matrix or vector which will be lagged by (q) to compute a VAR(q)
                 , X = NULL                             #X is a matrix of (contemporaneous) covariates
                 , trend = FALSE                        #trend is a boolean indicating whether a trend vector should be added to the VAR
                 , intercept = TRUE                     #intercept is a boolean indicating whether the test applies on the mean shift (TRUE) or all parameters (FALSE)
                 , ci = c(0.9, 0.95, 0.99)              #ci is the vector of confidence intervals (in growing order) to compute based on the CDF of a V distr.
                 , estMode = "OLS"                      #estMode can take values of "OLS", "FGLS", "IGLS"
                 , iter = 3                             #in the case of "IGLS", the number of iteration "iter" can be specified.
                 , aicbicMode = "AIC"                   #AicbicMode can be "AIC" or "BIC" depending on the minimum criterion to select
                 , qMax = 6                             #qMax is the number of lags (from 1 to qMax) tested to determine the AIC / BIC maximum.
                 , trim = 0.15                          #trim is the percentage parameter to start and end the sample analysis
                 , posBreak = FALSE                     #if we want the algorithm to only detect positive breaks
)
{
  myVars <- colnames(Y)                                                   #get the variable names

  qOpt <- compute_aicbic(Y, qMax, X, trend, intercept)                    #return the AIC and BIC criteria for lags from 1 to 6
  
  
  print(qOpt)                                                             #print the matrix of AIC and BIC for each lags

  q <- as.numeric(which.min(qOpt[aicbicMode,]))                           #choose the lag q according to the min AIC
  print(paste0("lag with the minimum ", aicbicMode, " is: ", q))

  lconf <- matrix_conformation(Y, q, X, trend, intercept)                 #create a list of conform objects for the estimation

  Yex <- lconf$Yex                                                        #get the conformed (expanded) Yex matrix (for the system, in vector form)
  Gex <- lconf$Gex                                                        #get the conformed G matrix of regressors for the system
  p <- lconf$p                                                            #final number of observations
  G <- lconf$G                                                            #original matrix of regressors
  S <- lconf$S                                                            #selection matrix
  Y <- lconf$Y                                                            #matching original dependent variables matrix
  nEq <- lconf$nEq                                                        #original number of equations / dependent variables
  myDates <- lconf$myDates                                                #matching dates

  print(paste0("The number of equations in the system is: ", nEq))

  fstat <- rep(NA, p)                                                     #create a vector of f_statistics for each k tested
  meanShift <- rep(NA, p)                                                 #create a vector with the evaluated size of the intercept difference
  CI <- matrix(data = NA, nrow = p, ncol = length(ci))                    #create a matrix of confidence intervals for each k tested

  cv <- compute_vdistr_cv(ci)                                             #compute critical values for the vector of confidence intervals proposed

  startInd <- round(trim * p)                                             #start index
  endInd <- round(p - trim * p)                                           #end index

  for(k in startInd:endInd)                                               #loop over the k with a trimming date / burn period
  {
    if(k%%10==0)
      print(paste0("The iteration is at the level k = ", k))              #get an idea of where we are in the loop every 10 iterations
    
    
    GexB <- Gex %*% t(S)
    GexB[1:((k - 1) * (nEq)),] <- 0                                       #force filling the GexB matrix with 0 before and original values after k

    Z <- t(cbind(Gex, GexB))                                              #bind the regressor and breaking regressor matrices together

    lbetaSigma <- compute_beta(Z, Yex, nEq, p, estMode, iter)             #compute the BetaSigma object list
    
    Beta <- lbetaSigma$Beta                                               #get the vector of betas
    Sigma <- lbetaSigma$Sigma                                             #get the covariance matrix of errors
    pBeta <- length(Beta)                                                 #get the length of the vector of betas

    #create a selection matrix to get only the betas of interest (breakings)
    if(intercept)                                                         #1 - case where only shift in intercept
    {
      R <- matrix(data = 0, nrow = nEq, ncol = pBeta)
      R[,(pBeta - nEq + 1):pBeta] <- diag(nEq)
    }

    if(!intercept)                                                        #2 - case where all parameters break
    {
      R <- matrix(data = 0, nrow = pBeta / 2 , ncol = pBeta)
      R[,(pBeta / 2 + 1):pBeta] <- diag(pBeta / 2)
    }

    fstat[k] <- compute_fstat(R, Beta, Z, p, Sigma)                      #compute the F-statistic for the current k
    CI[k, ] <- compute_ci(G, S, Sigma, R, Beta, cv, p)                   #compute the confidence interval for the current k

    meanShift[k] <- mean(R %*% Beta)                                       #get the mean intercept shift
  }

  if(posBreak)                                                           #if posBreak is TRUE, limit to positive break detection
    fstat[meanShift < 0] <- 0

  dev.new()                                                              
  plot(fstat)                                                            #plots the generated sequence of F-statistics

  g1 <- compute_plot_stats(myDates, myVars, fstat, CI, Y)
  
  breakInd <- which.max(fstat)
  breakDate <- myDates[breakInd]
  breakCi <- CI[breakInd, ]
  rownames(Y) <- myDates
  Gt <- t(G)
  rownames(Gt) <- myDates
  meanShift <- meanShift[breakInd]
  maxF = max(fstat, na.rm=T)
  trimDates = matrix(data=c(myDates[startInd], myDates[endInd]),nrow = 1, ncol = 2)
  colnames(trimDates) <- c("begin trim date", "end trim date")
  
  return(list(fstat = fstat                                               #return a "multibreak" class object
              , maxF = maxF
              , confInterval = breakCi
              , criticalValues = cv
              , breakDate = breakDate
              , Y = data.frame(Y)
              , G = data.frame(Gt)
              , breakInd = breakInd
              , meanShift = meanShift
              , aicbic = qOpt
              , g1 = g1
              , trimDates = trimDates))
}

Simulation

Now that everything is ready, let’s generate some clean data to test whether we are able to identify break dates. This extra function generates a matrix of p time series of n observations, based on the R rnorm generation. The argument intensity is a double value which is added to the generated data after the break that we can set to occur at any percent of the dataset (here at 35%)

Code

compute_simul <- function(n, p, intensity = 1, whenbreak = 0.35)
{
  set.seed(123)                                     #optional
  
  prebreakN <- round(n * whenbreak)
  postbreakN <- round(n * (1 - whenbreak))
  
  Y <- rbind(matrix(rnorm(prebreakN * p) + 5, nrow = prebreakN, ncol = p), matrix(rnorm(postbreakN * p) + 5 + intensity, nrow = postbreakN, ncol = p))
  startDate <- as.Date("22.01.2020", format = "%d.%m.%Y")                                            #date when I wrote this code
  
  simulDates <- 0:(n - 1) + startDate
  simulDates <- format(simulDates, format = "%d.%m.%Y")

  rownames(Y) <- simulDates
  colnames(Y) <- LETTERS[1:p]
  
  print(head(Y))
  return(Y)
}

Result with three simulated variables

Fig 1. Break detection for three simulated variables experiencing a common break at 35% of the sample, with an intensity parameter of four. The order of the VAR(q) is 1, after selection by the BIC, with no trend or additional covariates added. The break is easy to identify visually and is well captured by the algorithm. Max F-statistic: 159.34, and the confidence intervals for 90%, 95% and 99%, represented by the shaded areas are very concentrated around the break.

Fig 1. Break detection for three simulated variables experiencing a common break at 35% of the sample, with an intensity parameter of four. The order of the VAR(q) is 1, after selection by the BIC, with no trend or additional covariates added. The break is easy to identify visually and is well captured by the algorithm. Max F-statistic: 159.34, and the confidence intervals for 90%, 95% and 99%, represented by the shaded areas are very concentrated around the break.

Result with five simulated variables

Fig 2. Break detection for five simulated variables experiencing a common break at 35% of the sample, with an intensity parameter of 0.5. The order of the VAR(q) is 1, after selection by the BIC, with no trend or additional covariates added. The break is impossible to detect visually. However, the algorithm accurately captures it. Max F-statistic: 22.72 and the confidence intervals for 90%, 95% and 99%, represented by the shaded areas are much wider.

Fig 2. Break detection for five simulated variables experiencing a common break at 35% of the sample, with an intensity parameter of 0.5. The order of the VAR(q) is 1, after selection by the BIC, with no trend or additional covariates added. The break is impossible to detect visually. However, the algorithm accurately captures it. Max F-statistic: 22.72 and the confidence intervals for 90%, 95% and 99%, represented by the shaded areas are much wider.

To get the full working package, including data, you can clone or import the project from my github repository

Avatar
Loïc Maréchal
PhD candidate in Finance

My research interests include finance, statistical learning and natural science originated data applied to finance.

Related