R/TVECM.gen.R

In tsDyn: Nonlinear Time Series Models with Regime Switching

as.matrix.ts <- function(x, ...)  {
    # A function implemented by Diethelm Wuertz
    ans = as.matrix.default(unclass(x))
    attr(ans, "tsp")<-NULL
     rownames(ans)<-NULL # colnames(ans)<-NULL
    ans
}

#'Simulation and bootstrap a VECM or bivariate TVECM
#'
#'Estimate or bootstraps a multivariate Threshold VAR
#'
#'This function offers the possibility to generate series following a
#'VECM/TVECM from two approaches: bootstrap or simulation. \code{VECM.sim} is
#'just a wrapper for \code{\link{TVECM.sim}}.
#'
#'When the argument \code{matrix} is given, on can only simulate a VECM
#'(\code{nthresh}=0) or TVECM (\code{nthresh}=1 or 2). One can have a
#'specification with constant (\code{"const"}), \code{"trend"}, \code{"both"}
#'or \code{"none"} (see argument \code{include}). Order for the parameters is
#'ECT/include/lags for VECM and ECT1/include1/lags1/ECT2/include2/lags2 for
#'TVECM. To be sure that once is using it correctly, setting \code{show.parMat
#'= TRUE} will show the matrix of parameters together with their values and
#'names.
#'
#'The argument \code{beta} is the cointegrating value on the right side of the
#'long-run relationship, and hence the function use the vector (1,-beta). The
#'\code{innov} argument specifies the innovations. It should be given as a
#'matrix of dim nxk, (here \var{n} does not include the starting values!), by
#'default it uses a multivariate normal distribution, with covariance matrix
#'specified by \code{varcov}.
#'
#'The starting values (of dim lags x k) can be given through argument
#'\code{starting}. The user should take care for their choice, since it is not
#'sure that the simulated values will cross the threshold even once. Notice
#'that only one cointegrating value is allowed. User interested in simulating a
#'VECM with more cointegrating values should do use the VAR representation and
#'use \code{\link{TVAR.sim}}.
#'
#'The second possibility is to bootstrap series. This is done on a object
#'generated by \code{\link{TVECM}} (or \code{\link{VECM}}). A simple residual
#'bootstrap is done, or one can simulate a series with the same parameter
#'matrix and with normal distributed residuals (with variance pre-specified),
#'corresponding to Monte-carlo simulations.
#'
#'One can alternatively give only the series, and then the function will call
#'internally \code{\link{TVECM}}.
#'
#'@param B Simulation: Matrix of coefficients to simulate
#'@param Thresh,nthresh,lag,include Simulation: parameters for the VECM/TVECM to simulate. 
#'See \code{\link{TVECM}} for their description.  
#'@param object Object computed by function \code{\link{TVECM}}
#'or linear \code{\link{VECM}}
#'@param beta The cointegrating value
#'@param n Simulation: Number of observations to simulate. 
#'@param starting Simulation: Starting values (same length as lag = 1)
#'@param innov Simulation: time series of innovations/residuals. 
#'@param boot.scheme Bootstrap: which resampling scheme to use for the residuals. See \code{\link{resample_vec}}. 
#'@template param_show.parMat
#'@template param_returnStarting
#'@param seed Bootstrap: seed used in the resampling
#'@param \dots additional arguments for the unexported \code{TVECM.gen}.  
#'@return A matrix with the simulated/bootstrapped series.
#'@author Matthieu Stigler
#'@seealso \code{\link{VECM}} or \code{\link{TVECM}} to estimate the VECM or TVECM.  
#'Similar \code{\link{TVAR.sim}} and \code{\link{TVAR.boot}} for \code{\link{TVAR}}, 
#'\code{\link{VAR.sim}} and \code{\link{VAR.boot}} for VAR models estimated with \code{\link{lineVar}} models. 
#'@keywords ts
#'@export
#'@examples
#'
#'
#'###reproduce example in Enders (2004, 2 edition) p. 350, 
#' # (similar example in Enders (2010, 3 edition) 301-302). 
#'
#'if(require(mnormt)){
#'#see that the full "VAR" coefficient matrix is:
#'  A <- matrix(c(-0.2, 0.2, 0.2, -0.2), byrow=TRUE, ncol=2)
#'
#'# but this is not the input of VECM.sim. You should decompose into the a and b matrix:
#'  a<-matrix(c(-0.2, 0.2), ncol=1)
#'  b<-matrix(c(1,-1), nrow=1)
#'
#'# so that:
#'  a%*%b
#'
#'# The a matrix is the input under argument B, while the b matrix is under argument beta: 
#' # (the other zeros in B are for the not-specified lags)
#'  innov<-rmnorm(100, varcov=diag(2))
#'  Bvecm <- rbind(c(-0.2, 0,0), c(0.2, 0,0))
#'  vecm1 <- VECM.sim(B=Bvecm, beta=1,n=100, lag=1,include="none", innov=innov)
#'  ECT <- vecm1[,1]-vecm1[,2]
#'
#'#add an intercept as in panel B
#'  Bvecm2 <- rbind(c(-0.2, 0.1,0,0), c(0.2,0.4, 0,0))
#'  vecm2 <- VECM.sim(B=Bvecm2,  n=100,beta=1, lag=1,include="const", innov=innov)
#'
#'  par(mfrow=c(2,1))
#'  plot(vecm1[,1], type="l", main="Panel a: no drift or intercept", ylab="", xlab="")
#'  lines(vecm1[,2], lty=2)
#'  plot(vecm2[,1], type="l", main="Panel b: drift terms (0.1)", ylab="", xlab="")
#'  lines(vecm2[,2], lty=2)
#'}
#'##Bootstrap a TVAR with 1 threshold (two regimes)
#'data(zeroyld)
#'TVECMobject <- TVECM(zeroyld, nthresh=1, lag=1, ngridBeta=20, ngridTh=20, plot=FALSE, trace = FALSE)
#'TVECM.boot(TVECMobject)
#'
#'##Check the bootstrap: do we get original series, when not resampling residuals?
#' TVECM.boot.check <- TVECM.boot(TVECMobject, boot.scheme = "check")
#' all.equal(as.data.frame(TVECM.boot.check), zeroyld)
#'
#' @name TVECM.sim
NULL # for roxygen, do not add 

VECM.gen <- function(B, beta, n=200, lag=1, 
                      include = c("const", "trend","none", "both"), 
                      starting=NULL, innov,
                      returnStarting = FALSE, 
                      show.parMat=FALSE, 
                      round_digits = 4){
  TVECM.gen(B=B, beta=beta, n=n, lag=lag, 
            include = include, 
            nthresh=0,
            starting=starting, innov = innov,
            returnStarting = returnStarting, show.parMat=show.parMat, round_digits = round_digits)
    
}
  
TVECM.gen <- function(B, beta, n=200, lag=1, 
                      include = c("const", "trend","none", "both"), 
                      nthresh=1, Thresh, 
                      starting=NULL, innov,
                      returnStarting = FALSE, 
                      show.parMat=FALSE, 
                      round_digits = 4){
  

  p <- lag
  include <- match.arg(include)
  
  ###check correct arguments
  if(!nthresh%in%c(0,1,2))
    stop("Arg nthresh should  be either 0, 1 or 2")
  if(missing(beta))
    stop("please provide arg beta (cointegrating value)")  
  
  
  ##include term
  ninc <- switch(include, "none"=0, "const"=1, "trend"=1, "both"=2)
  incVal <- switch(include, "none"=NULL, "const"="const", "trend"="trend", "both"=c("const","trend"))
  
  ### 
  k <- nrow(B)
  if(is.vector(beta)){
    if(length(beta)==k-1) beta <- c(1, -beta)
    tBETA<-matrix(beta, nrow=1)
    r <- 1
  } else {
    if(nrow(beta)!=k ) stop("beta should have k rows and r cols")
    r <- ncol(beta)
    tBETA <- t(beta)
  }
  esp <- p*k+r+ninc #number of lags +ecm
  
  ## naming of variables:
  pa <- switch(as.character(nthresh), "0"="", "1"=c("_low", "_upr"),"2"=c("_low", "_mid","_upr"))
  lags <- as.vector(outer("L{x", 1:k,  paste, sep=""))
  lags2 <- paste(rep(lags, times=p),"}{", rep(1:p,each=p),"}",sep="")
  
  if(esp*(nthresh+1)!=ncol(B)){
    colnames_Matrix_input <- as.vector(outer(c(rep("ECT", r), incVal, lags2), pa, paste, sep=""))
    cat("Matrix B badly specified: should have ", esp*(nthresh+1), "columns, but has", ncol(B), "\n")
    print(matrix(NA, nrow=k, ncol=esp*(nthresh+1), dimnames=list(paste("Equ x", 1:k, sep=""), colnames_Matrix_input)))
    stop()
  }
  rownames(B) <- paste("Equ x", 1:k, ":",sep="")
  

  if(!is.null(starting)){
    if(!is.matrix(starting)) stop("Provide 'starting' as matrix")
    if(!all(dim(starting) == c(p+1, k))) stop("Bad specification of starting values. Should have nrow = lag +1 and ncol = number of variables")
  }
  Bmat <- B
  
  ##### put coefficients vector in right form according to arg include (arg both need no modif)
  if(include!="both"){
    aa1 <- r+switch(include, "none"=1:2, "const"=2, "trend"=1, "both"=NULL)
    aa <- sort(rep(aa1, each=nthresh+1)+ (0:nthresh)*(p*k+max(aa1)))
    Bmat <- myInsertCol(Bmat, c=aa, 0)
  }
  nparBmat <- p * k + 2 + r
  
  
  ##############################
  ###Reconstitution boot/simul
  ##############################
  
  #initial values
  Yb <- matrix(0, nrow= n+p +1, ncol=k)
  if(!is.null(starting))  Yb[1:(p+1),] <- starting # needs 2 starting, as will take diff!
  
  trend <- c(rep(NA, p+1), seq_len(n))
  
  
  if(nthresh==0){
    for(i in (p+2):(n + p + 1)){
      ECT <- Bmat[,1:r]%*%tBETA%*%matrix(Yb[i-1,], ncol=1)
      Yb[i,] <- rowSums(cbind(Yb[i-1,],
                              Bmat[,r+1], 
                              Bmat[,r+2]*trend[i], 
                              ECT,
                              Bmat[,-c(1:(r+2))]%*%matrix(t(Yb[i-c(1:p),]-Yb[i-c(2:(p+1)),]), ncol=1),
                              innov[i-p-1,]))
    }
  } else if(nthresh==1){
    BD <- Bmat[,seq_len(nparBmat)]
    BU <- Bmat[,-seq_len(nparBmat)]
    
    for(i in (p+2):(n + p + 1)){
      ECT <- Yb[i-1,, drop = FALSE] %*% t(tBETA)
      if(round(ECT,round_digits)<=Thresh){
        B_here <-  BD
      }  else{
        B_here <-  BU
      }
      Yb[i,] <- rowSums(cbind(Yb[i-1,],
                              B_here[,1] * as.vector(ECT),
                              B_here[,2], 
                              B_here[,3]*trend[i],
                              B_here[,-c(1,2,3)] %*% matrix(t(Yb[i-c(1:p),]-Yb[i-c(2:(p+1)),]), ncol=1),
                              innov[i-p-1,]))
    }
  } else if(nthresh==2){
    BD <- Bmat[,seq_len(nparBmat)]
    BM <- Bmat[,seq_len(nparBmat)+nparBmat]
    BU <- Bmat[,seq_len(nparBmat)+2*nparBmat]
    for(i in (p+2):(n + p + 1)){
      ECT <- Yb[i-1,, drop = FALSE] %*% t(tBETA)
      if(round(ECT,round_digits)<=Thresh[1]){ 
        B_here <-  BD
      } else if(round(ECT,round_digits)>Thresh[2]) {
        B_here <-  BU
      } else{
        B_here <-  BM
      }
      Yb[i,] <- rowSums(cbind(Yb[i-1,],
                              B_here[,1] * as.vector(ECT),
                              B_here[,2], 
                              B_here[,3]*trend[i],
                              B_here[,-c(1,2,3)] %*% matrix(t(Yb[i-c(1:p),]-Yb[i-c(2:(p+1)),]), ncol=1),
                              innov[i-p-1,]))
    }
  }
  
  if(show.parMat){
    colnames_Matrix_system <- as.vector(outer(c(rep("ECT",r), "Const", "Trend", lags2), pa, paste, sep=""))
    colnames(Bmat)<- colnames_Matrix_system
    print(Bmat)
  }
  # res <- round(Yb, round_digits)
  if(!returnStarting) Yb <- Yb[-seq_len(p+1), ] 
  return(Yb)
}

#' @export
#' @rdname TVECM.sim
TVECM.sim <- function(B, n=200, lag=1, include = c("const", "trend","none", "both"),  
                      beta, 
                      nthresh = 1, Thresh,
                      starting=NULL, 
                      innov=rmnorm(n, varcov= diag(1,nrow(B))), 
                      show.parMat=FALSE, returnStarting=FALSE, ...){
  
  TVECM.gen(B=B, n=n, lag=lag, include = include,  
            beta = beta,
            nthresh = nthresh, Thresh = Thresh,
            starting=starting, innov=innov, 
            show.parMat=show.parMat, returnStarting=returnStarting, ...)
}

#' @export
#' @rdname TVECM.sim
VECM.sim <- function(B, n=200, lag=1, include = c("const", "trend","none", "both"),  
                     beta, 
                     starting=NULL, 
                     innov=rmnorm(n, varcov= diag(1,nrow(B))), 
                     show.parMat=FALSE, returnStarting=FALSE, ...){
  
  TVECM.gen(B=B, n=n, lag=lag, include = include,  
            beta = beta,
            nthresh = 0, 
            starting=starting, innov=innov, 
            show.parMat=show.parMat, returnStarting=returnStarting, ...)
}

#' @export
#' @rdname TVECM.sim
VECM.boot <- function(object, boot.scheme = c("resample", "resample_block", "wild1", "wild2", "check"),
                       seed = NULL, ...){
  TVECM.boot(object, boot.scheme = boot.scheme, seed = seed, ...)
}
  

#' @export
#' @rdname TVECM.sim
TVECM.boot <- function(object, boot.scheme = c("resample", "resample_block", "wild1", "wild2", "check"),
                       seed = NULL, ...){
  
  ## check bad cases
  if(!is.null(object$exogen) && object$exogen) stop("Cannot handle exo variables")
  if(object$model.specific$LRinclude != "none") stop("Cannot handle cases with LRinclude != 'none' ")
  
  
  
  boot.scheme <- match.arg(boot.scheme)
  
  B <- coef(object)
  if(inherits(object, "TVECM")) {
    B <-  do.call("cbind", B)
  }
  beta <-  object$model.specific$coint
  nthresh <- object$model.specific$nthresh
  t <- object$t
  k <- object$k
  lags <- object$lag
  if(lags==0) stop("Cannot handle case with lag = 0")
  startLag <- if(lags==0) 0 else 1
  include <- object$include
  # num_exogen <- object$num_exogen 
  
  YX <- object$model
  yorig <- YX[,1:k]
  starts <- as.matrix(yorig[seq_len(lags+1),, drop=FALSE])
  
  ## residuals, boot them
  resids <- residuals(object) ##[-seq_len(lags + 1),,drop = FALSE]
  if(!is.null(seed)) set.seed(seed) 
  innov <- resample_vec(resids, boot.scheme = boot.scheme, seed=seed)
  
  ## Exogen
  # if(num_exogen>0){
  #   nYX <- ncol(YX)
  #   exogen <- YX[-c(startLag:lags), -(num_exogen-1):0+nYX]
  #   #     starts <- cbind(starts, matrix(0, nrow=nrow(starts), ncol=num_exogen))
  # } else {
  #   exogen <- NULL
  # }  
  
  res <- TVECM.gen(B=B, n=t, lag=lags, include = include,  
                   beta = beta,
                   nthresh = nthresh,
                   Thresh = getTh(object), 
                   starting=starts, innov=innov, 
                   # exogen=exogen, 
                   returnStarting=TRUE,
                   ...)
  colnames(res) <- colnames(yorig)
  res
}

TVECM.boot.check <- function(object, ...) {
  dat_orig <-  as.data.frame(object$model[, 1:object$k])
  rownames(dat_orig) <- seq_len(nrow(dat_orig))
  boot <-  TVECM.boot(object, boot.scheme = "check", ...)
  all.equal(dat_orig, as.data.frame(boot))
}

if(FALSE){
  library(tsDyn)
  TVECM.boot.check <-  tsDyn:::TVECM.boot.check
  library(devtools)
  load_all()
  # environment(TVECM.gen)<-environment(star)
  # environment(TVECM.boot)<-environment(star)
  
  ##Simulation of a TVAR with 1 threshold
  B_0th <- rbind(c(-0.2, 0,0), c(0.2, 0,0))
  B_1th <- rbind(c(-0.2, 0.1,0.2, -0.1, 0.3, 0.4), c(0.2, 0.2,0.1, -0.3, 0.1, 0.15))
  B = B_1th
  beta <- 0.4
  nthresh = 1
  Thresh = 2.2
  lag = 1
  seed <-  NULL
  starting= matrix(c(2,1.5, 2.1, 1.9), nrow = 2)
  include="none"
  n= 10
  boot.scheme = "resample"
  innov <-  matrix(rnorm(n*2), ncol = 2)
  show.parMat = TRUE
  
  ## sim nthresh = 0
  a <- TVECM.gen(B=B_0th, nthresh=0, beta=1, lag=1,include="none", starting= starting, 
                 innov = innov, n = 10)
  
  ## sim nthresh = 1
  a <- TVECM.gen(B=B, nthresh=1, beta=1, lag=1,include="none", starting= starting, Thresh = 2.1,
                 innov = innov, n = 10)
  
  
  ### BOOT
  v_l1 <- VECM(zeroyld, lag=1, estim = "ML")
  v_l1_tr <- VECM(zeroyld, lag=1, estim = "ML", include = "trend")
  v_l1_r2 <- VECM(barry, lag=1, estim = "ML", r = 2)
  v_l1_lrInc <- VECM(zeroyld, lag=1, LRinclude = "const", estim = "ML")
  v_l1_lrInc$model.specific$coint
  v_l2 <- VECM(zeroyld, lag=2)
  
  TVECM.boot.check(object=v_l1)
  TVECM.boot.check(object = v_l2)
  TVECM.boot.check(object = v_l1_lrInc)
  TVECM.boot.check(object = v_l1_r2)
  TVECM.boot.check(object = v_l1_tr)
  
  
  tv_1th_l1 <- TVECM(zeroyld, lag=1, nthresh=1, plot=FALSE, trace=FALSE)
  tv_2th_l1 <- TVECM(zeroyld, lag=1, nthresh=2, plot=FALSE, trace=FALSE)
  tv_1th_l2 <- TVECM(zeroyld, lag=2, nthresh=1, plot=FALSE, trace=FALSE)
  tv_2th_l2 <- TVECM(zeroyld, lag=2, nthresh=2, plot=FALSE, trace=FALSE)
  
  TVECM.boot.check(tv_1th_l1)
  TVECM.boot.check(tv_2th_l1)
  TVECM.boot.check(tv_1th_l2)
  TVECM.boot.check(tv_2th_l2)
  
  TVECM.boot(object = tv_1th_l1, seed = 1)
  
  
  ECT<-a[,1]-a[,2]
  
  layout(matrix(1:2, ncol=1))
  plot(a[,1], type="l", xlab="", ylab="", ylim=range(a, ECT))
  lines(a[,2], col=2, type="l")
  
  plot(ECT, type="l")
  
  B<-rbind(c(0.2, 0.11928245, 1.00880447, -0.009974585, 0.3, -0.089316, 0.95425564, 0.02592617),
           c( -0.1, 0.25283578, 0.09182279,  0.914763741, 0.35,-0.0530613, 0.02248586, 0.94309347))
  sim <- TVECM.sim(B=B,beta=1, nthresh=1, n=500, Thresh=5, starting= matrix(c(5.2, 5.5, 5.1, 5.3), nrow = 2))
  
  #estimate the new serie
  TVECM(sim, lag=1)
  
  ##Bootstrap a TVAR with two threshold (three regimes)
  #data(zeroyld)
  TVECMobject <- TVECM(zeroyld, lag=1, nthresh=2, plot=FALSE, trace=FALSE)
  TVECM.boot(TVECMobject)
  
  #not working: (probably trend coefficients too small so digits errors)
  all(TVECM.boot(TVECMobject=lineVar(zeroyld, lag=1, model="VECM", include="trend"),type="check")==zeroyld)
  all(TVECM.sim(TVECMobject=lineVar(zeroyld, lag=1, model="VECM", include="both"),type="check")==zeroyld)
  
  #nthresh=1
  TVECMobject<-TVECM(zeroyld, nthresh=1, lag=1, ngridBeta=20, ngridTh=20, plot=FALSE)
  all(TVECM.sim(TVECMobject=TVECMobject,type="check")==zeroyld)
  
  all(TVECM.sim(TVECMobject=TVECM(zeroyld, nthresh=1, lag=2, ngridBeta=20, ngridTh=20, plot=FALSE),type="check")==zeroyld)
  all(TVECM.sim(TVECMobject=TVECM(zeroyld, nthresh=1, lag=1, ngridBeta=20, ngridTh=20, plot=FALSE, include="none"),type="check")==zeroyld)
  all(TVECM.sim(TVECMobject=TVECM(zeroyld, nthresh=1, lag=2, ngridBeta=20, ngridTh=20, plot=FALSE, include="none"),type="check")==zeroyld)
  
  #nthresh=2
  TVECMobject2<-TVECM(zeroyld, nthresh=2, lag=1, ngridBeta=20, ngridTh=20, plot=FALSE)
  all(TVECM.sim(TVECMobject=TVECMobject2,type="check")==zeroyld)
  all(TVECM.sim(TVECMobject=TVECM(zeroyld, nthresh=2, lag=2, ngridBeta=20, ngridTh=20, plot=FALSE),type="check")==zeroyld)
  
  all(TVECM.sim(TVECMobject=TVECM(zeroyld, nthresh=2, lag=1, ngridBeta=20, ngridTh=20, plot=FALSE, include="none"),type="check")==zeroyld) 
  #famous rounding problem...
  all(TVECM.sim(TVECMobject=TVECM(zeroyld, nthresh=2, lag=2, ngridBeta=20, ngridTh=20, plot=FALSE, include="none"),type="check")==zeroyld)
  
  ###TODO:
  #improve trend/both case
  #TVECM: possibility to give args!
  TVECM(zeroyld, nthresh=1, lag=2, ngridBeta=20, ngridTh=20, plot=FALSE, th1=list(exact=-1.4),include="none")
  TVECM(zeroyld, nthresh=1, lag=2, ngridBeta=20, ngridTh=20, plot=FALSE, th1=list(exact=-1.4),beta=list(exact=1),include="none")
  TVECM(zeroyld, nthresh=2, lag=2, ngridBeta=20, ngridTh=20, plot=FALSE, th1=list(exact=-1.4),th2=list(exact=0.5),include="none")
}