## File: TVECM_HSTest.R
## Source code of the function TVCECM.HSTest in the version 0.7-40 of the
## tsDyn package.
##
## Reference: Berner Larsen (2012). A threshold cointegration analysis
##            of Norwegian interest rates, Master Thesis in Statistics,
##            University of Tromso. Norway. Available at:
##            https://munin.uit.no/handle/10037/4370
#########################################################################

library(foreach) # Provides Foreach Looping Construct 
library(tsDyn)   # Nonlinear Time Series Models with Regime Switching

TVECM.HStest <- function(data, lag = 1, ngridTh = 300, trim = 0.05, nboot = 100,
fixed.beta=NULL, intercept=TRUE, boot.type=c("FixedReg", "ResBoot"),
hpc = c("none", "foreach"), common=c("All", "ECT"), type=c("2Reg", "SymReg")) {

## Check args:
boot.type <- match.arg(boot.type)
hpc       <- match.arg(hpc)
common    <- match.arg(common)
type      <- match.arg(type)

dir=FALSE # internal value, was used to try different implementation of lmtest

### Organize Data
data <- as.matrix(data)
if(ncol(data)>2) {warning("Please no more than two equations")}
if(is.null(colnames(data))){colnames(data)<-paste("Var", c(1:2), sep="")}
T <- nrow(data)
p <- lag
y <- diff(data)[(p+1):(T-1),]

DeltaX <- embed(diff(data),p+1)[,-(1:2)]
if(intercept) DeltaX <- cbind(1, DeltaX)
x <- DeltaX
t <- nrow(y)

### Compute beta with VECM() and extract ect
if(is.null(fixed.beta)){
  ve <- VECM(data, lag=lag, include="const", estim="ML")
  }else{
  ve <- VECM(data, lag=lag, include="const", beta=fixed.beta, estim="2OLS")
}

ect <- ve$model[,grep("ECT", colnames(ve$model))]
w0  <- matrix(ect[!is.na(ect)], ncol=1)

### Set up of the grid
q <- if(type=="2Reg") sort(w0) else sort(abs(w0))
if(ngridTh>(1-2*trim)*t) {
  newNgridTh <- round((1-2*trim)*t-1)
  warning("ngridTh (", ngridTh,
  ") bigger than number of potential threshold values, set to ",newNgridTh,"\n")
  ngridTh <- newNgridTh
}

gamma2  <- q[round(seq(from=trim*t, to=(1-trim)*t,length.out=ngridTh))]
gamma2  <- unique(gamma2)
ngridTh <- length(gamma2)

###########
### LM Test
###########
lmtest02 <- function(y,x,w0,gammas,dir=dir){
  # y: y var, x: intercept and lags matrix, w0: ECT term, gammas: potential thresholds
  X <- cbind(w0,x) #X: ECT and intercept and lags
  if(dir){
    q         <- qr(X)
    res_restr <- qr.resid(q,y)
  } else{
  z0zz      <- X%*%solve(t(X)%*%X)
  res_restr <- lm.fit(X,y)$residuals # residuals from the linear VECM given b0
  }

  res_restr1 <- res_restr[,1]
  res_restr2 <- res_restr[,2]

  store   <- rep(NA, ngridTh)
  Ttrim   <- trim*t
  ngridTh <- min(t*(1-2*trim), length(gammas))

  for(j in 1:ngridTh){
    d1 <- if(type=="2Reg") ifelse(w0<=gammas[j],1,0) else
    ifelse(w0<= - gammas[j] | w0> gammas[j],1,0) # d1: dummy variable
  n1 <- sum(d1)
  if (min(c(n1,(t-n1)))>Ttrim){
    z1 <- if(common=="All") c(d1)*X else c(d1)*w0
    res_unrestr <- if(dir) qr.resid(q, z1) else z1-z0zz%*%(t(X)%*%z1)
  # z11: residuals from unrestricted model (with threhsold)
  zea <- res_restr1*res_unrestr
  zeb <- res_restr2*res_unrestr
  ze  <- cbind(zea,zeb)
#[(z11.*(e1*ones(1,length(z11(1,:))))),(z11.*(e2*ones(1,length(z11(1,:)))))];
  v     <- crossprod(ze)
  z11y  <- crossprod(res_unrestr,y)
  s     <- matrix(c(z11y), ncol=1) # vectorization of the parameter matrix z11y
  VV    <- crossprod(v)
  VVinv <- try(solve(VV), silent=TRUE)
  
  if(inherits(VVinv, "try-error")) VVinv <- ginv(VV)
     store[j] <- t(s)%*%VVinv%*%t(v)%*%s
  } # end of the if
  
} # end of the whole loop
return(store)
} # end of the function lmtest01

#### LM test for fixed regressor bootstrap: X'X^-1 evaluated only once
lmtest02_boot<-function(y,x,w0,gammas,dir=dir){
  X <- cbind(w0,x) # X: ECT and intercept and lags

  if(dir){
    res_restr <- qr.resid(q,y)
  } else{
    res_restr <- lm.fit(X,y)$residuals # residuals from the linear VECM given b0
  }

  res_restr1 <- res_restr[,1]
  res_restr2 <- res_restr[,2]
  store      <- rep(0, ngridTh)
  ngridTh    <- min(t*(1-2*trim), length(gammas))

  for(j in 1:ngridTh){
    d1 <- if(type=="2Reg") ifelse(w0<=gammas[j],1,0) else
    ifelse(w0<= - gammas[j] | w0> gammas[j],1,0) #d1: dummy variable
    n1 <- sum(d1)
    if (min(c(n1,(t-n1)))>Ttrim){
      z1<-if(common=="All") c(d1)*X else c(d1)*w0
      res_unrestr <- if(dir) qr.resid(q, z1) else z1-z0zz%*%(t(X)%*%z1)
      # z11: residuals from unrestricted model (with threhsold)
      zea <- res_restr1*res_unrestr
      zeb <- res_restr2*res_unrestr
      ze  <- cbind(zea,zeb)
      # [(z11.*(e1*ones(1,length(z11(1,:))))),(z11.*(e2*ones(1,length(z11(1,:)))))];
      v     <- crossprod(ze)
      z11y  <- crossprod(res_unrestr,y)
      s     <- matrix(c(z11y), ncol=1) # vectorization of the parameter matrix z11y
      VV    <- crossprod(v)
      VVinv <- try(solve(VV), silent=TRUE)
      if(inherits(VVinv, "try-error")) VVinv <- ginv(VV)
        store[j]<-t(s)%*%VVinv%*%t(v)%*%s
      } # end of the if
  } # end of the whole loop

  lm01 <- max(store, na.rm=TRUE)
  lm01
} # end of the function lmtest02
###

lm01     <- lmtest02(y,x,w0,gamma2, dir=dir)
teststat <- max(lm01, na.rm=TRUE)

##################################
### Bootstraps
##################################
if(nboot==0){
  CriticalValBoot <- NULL
  PvalBoot        <- NULL
  boots.reps      <- NULL
  if(hpc=="foreach") warning("hpc='foreach' used only when nboot>0\n")
   } else if (nboot>0){
  ##################################
  ### Fixed Regressor Bootstrap %

  ##################################
  if(boot.type=="FixedReg"){
    X     <- cbind(w0,x) #X: ECT and intercept and lags
    Ttrim <- trim*t
    if(dir){
      q <- qr(X)
      } else{
      z0zz <- X%*%solve(t(X)%*%X)
    }
  
    lmtest_withBoot<-function(e){
      yr <- rnorm(n=t,0,1)*e
      return(lmtest02_boot(yr,x,w0,gamma2,dir=dir))
     }

   boots.reps <- if(hpc=="none") replicate(nboot, lmtest_withBoot(e=residuals(ve)))
   else foreach(i=1:nboot, .export="lmtest_withBoot", .combine="c")%dopar%lmtest_withBoot(e=residuals(ve))

  ##################################
  ### Residual Bootstrap
  ##################################
  } else{
    lmtest_with_resBoot<-function(ve){
    # bootstrap it
    data.boot <- TVECM.sim(TVECMobject=ve, type="boot")
    # estimate VECM
    if(is.null(fixed.beta)){
      ve.boot <- VECM(data.boot, lag=lag, include="const", estim="ML")
    }else{
      ve.boot <- VECM(data.boot, lag=lag, include="const", beta=fixed.beta,
      estim="2OLS")
    }
  
  # extract w0, y and x
  ect.boot  <- ve.boot$model[,"ECT"]
  which.nas <- 1:(p+1)
  w0.boot   <- matrix(ect.boot[-which.nas], ncol=1)
  x.boot    <- ve.boot$model[-which.nas,-c(1:3)]
  y.boot    <- ve.boot$model[,c(1:2)]
  y.boot    <- diff(y.boot)[(p+1):(T-1),]
 
  # set-up grid
  w0.ord.boot <- sort(w0)
  gamma2.boot <- w0.ord.boot[round(seq(from=trim*T,
  to=(1-trim)*T,length.out=ngridTh))]
  gamma2.boot  <- unique(gamma2.boot)
  ngridTh.boot <- length(gamma2.boot)
  test.boot    <- lmtest02(y.boot,x.boot,w0.boot,gamma2.boot,dir=dir)
  return(max(test.boot, na.rm=TRUE))
}

boots.reps <- if(hpc=="none") replicate(nboot, lmtest_with_resBoot(ve)) else
foreach(i=1:nboot,.export="lmtest_with_resBoot", .combine="c")%dopar% lmtest_with_resBoot
} # end if boot= ResBoot
## result: compute p values and critical values:

PvalBoot        <- mean(ifelse(boots.reps>teststat,1,0))
CriticalValBoot <- quantile(boots.reps, probs= c(0.9, 0.95,0.99))
} # end if boot > 0

####### Return args
  args           <- list()
  args$nboot     <- nboot
  args$boot.type <- boot.type
  ret            <- list()

  ret$args     <- args
  ret$stat     <- teststat
  ret$values   <- lm01
  ret$ths      <- gamma2
  ret$maxTh    <- gamma2[which(lm01==ret$stat)]
  ret$PvalBoot <- PvalBoot

  ret$CriticalValBoot <- CriticalValBoot
  ret$allBoots        <- boots.reps
  class(ret)          <- "TVECMHanSeo02Test"
return(ret)
} # End of the whole function

### Print method
print.TVECMHanSeo02Test <- function(x,...){
  cat("## Test of linear versus threshold cointegration",
  " of Hansen and Seo (2002) ##\n\n", sep="")
  cat("Test Statistic:\t", x$stat)
  cat("\t(Maximized for threshold value:", x$maxTh, ")\n")
  if(x$args$nboot>0){
    boot.name<-switch(x$args$boot.type, "FixedReg"="Fixed regressor bootstrap",
    "ResBoot"="Residual Bootstrap")
    cat("P-Value:\t", x$PvalBoot, "\t\t(",boot.name, ")\n")
   }
}