# Feature selection for a simulated data set

    require(sparsevar)
    require(fsMTS)
    require(plot.matrix)
    require(svMisc)
    require(MTS)

    k <- 5
    max.lag <- 3
    N <- 512
    show.progress = F
    genM <- function(k, sparsity_lag = 0.2){
      nonzero <- round(sparsity_lag*k*k)
      m = diag(k)*runif(k*k, min = 0, max=0.1)
      size <- k
      samples <- nonzero
      vals <- sample.int(size ^ 2, samples)
      pair<-cbind(vals %/% size + 1, vals %% size+1)
      pair[pair>k]<-k
      m[pair] <- runif(nonzero, min = -0.2, max=0.2)
      return(m)
    }
    set.seed(9)
    pars <- list()
    for (i in 1:max.lag){
      pars[[paste0("l",i)]] <- genM(k)
    }
    res<-data.frame()
    for (i in pars){
      res <-rbind(res,i)
    }
    mreal<-as.matrix(res)
    plot(abs(mreal), digits=2, col=rev(heat.colors(10)), key=NULL,
         main="Real structure of dependencies",
         xlab="Variables")

    mts.sim <- simulateVAR(N=k, p=max.lag, nobs=N,fixedMat=pars)
    data <- mts.sim$series
    colnames(data)<-paste0("V",1:ncol(data))

    shortest <- matrix(rexp(k*k, rate = 0.2), nrow=k)
    shortest <- shortest-diag(k)*shortest
    colnames(shortest)<-colnames(data)
    rownames(shortest)<-colnames(data)
    plot(shortest, digits=2, col=rev(heat.colors(10)), key=NULL,
         main="Shortest distances between nodes")

    data <- scale(data)
    data[5,]<-NA

    mIndep<-fsMTS(data, max.lag=max.lag, method="ownlags")
    plot(mIndep, digits=2, col=rev(heat.colors(10)), key=NULL,
         main="Only own lags")

    mCCF<-fsMTS(data, max.lag=max.lag, method="CCF",show.progress=show.progress)
    plot(mCCF, digits=2, col=rev(heat.colors(10)), key=NULL,
         main="Cross-correlations")

    mDistance<-fsMTS(data, max.lag=max.lag, method="distance", shortest = shortest, step = 1)
    plot(mDistance, digits=2, col=rev(heat.colors(10)), key=NULL,
         main="Distance-based feature selection")

    mGLASSO.localized<-fsMTS(data, max.lag=max.lag,method="GLASSO", rho = 0.1,show.progress=show.progress, localized= TRUE)
    plot(mGLASSO.localized, digits=2, col=rev(heat.colors(10)), key=NULL,
         main="Graphical LASSO-based (localized) feature selection")

    mGLASSO.global<-fsMTS(data, max.lag=max.lag,method="GLASSO", rho = 0.1,show.progress=show.progress, localized = FALSE)
    plot(mGLASSO.global, digits=2, col=rev(heat.colors(10)), key=NULL,
         main="Graphical LASSO-based (global) feature selection")

    mLARS<-fsMTS(data, max.lag=max.lag,method="LARS",show.progress=show.progress)
    plot(mLARS, digits=2, col=rev(heat.colors(10)), key=NULL,
         main="Least angle  regression-based feature selection")

    mRF.localized<-fsMTS(data, max.lag=max.lag,method="RF",show.progress=show.progress, localized = TRUE)
    plot(mRF.localized, digits=2, col=rev(heat.colors(10)), key=NULL,
         main="Random forest-based (localized) feature selection")

    mRF.global<-fsMTS(data, max.lag=max.lag,method="RF",show.progress=show.progress, localized = FALSE)
    plot(mRF.global, digits=2, col=rev(heat.colors(10)), key=NULL,
         main="Random forest-based (global) feature selection")

    mMI.localized<-fsMTS(data, max.lag=max.lag,method="MI",show.progress=show.progress, localized= TRUE)
    plot(mMI.localized, digits=2, col=rev(heat.colors(10)), key=NULL,
         main="Mutual information-based (localized) feature selection")

    mMI.global<-fsMTS(data, max.lag=max.lag,method="MI",show.progress=show.progress, localized= FALSE)
    plot(mMI.global, digits=2, col=rev(heat.colors(10)), key=NULL,
         main="Mutual information-based (global) feature selection")

    mPSC<-fsMTS(data, max.lag=max.lag,method="PSC",show.progress=show.progress)
    plot(mPSC, digits=2, col=rev(heat.colors(10)), key=NULL,
         main="PSC-based (global) feature selection")

    mlist <- list(Independent = mIndep,
                  Distance = mDistance,
                  CCF = mCCF,
                  GLASSO.localized = mGLASSO.localized,
                  GLASSO.global=mGLASSO.global,
                  LARS = mLARS,
                  RF.localized = mRF.localized,
                  RF.global = mRF.global,
                  MI.localized = mMI.localized,
                  MI.global = mMI.global,
                  PSC=mPSC)

    th<-0.1
    mE1 <- fsEnsemble(mlist, threshold = th, method="ranking")
    plot(mE1, digits=2, col=rev(heat.colors(10)), key=NULL,
         main="Ensemble feature selection  using Ranking")

    mlist[["EnsembleRank"]] <- mE1

    mE2 <- fsEnsemble(mlist, threshold = th, method="majority")
    plot(mE2, digits=2, col=rev(heat.colors(10)), key=NULL,
         main="Ensemble feature selection  using Majority Voting")

    mlist[["EnsembleMajV"]] <- mE2

    mlist[["Real"]] <- ifelse(mreal!=0,1,0)
    msimilarity <- fsSimilarityMatrix(mlist, threshold=th, method="Kuncheva")
    plot(msimilarity, digits=2, col=rev(heat.colors(ncol(msimilarity))), key=NULL,
         main="Pairwise comparison of feature sets", cex.axis=0.7)

Use obtained features for model estimation and prediction

    dat <- stats::na.omit(data)
    model <- VAR(dat, p=max.lag, include.mean = F, fixed = mE2)

    ## AR coefficient matrix
    ## AR( 1 )-matrix
    ##      [,1]  [,2] [,3]     [,4]    [,5]
    ## [1,]    0 0.000    0 -0.00352  0.0000
    ## [2,]    0 0.235    0  0.00000  0.0000
    ## [3,]    0 0.000    0  0.00000  0.0000
    ## [4,]    0 0.000    0  0.00000 -0.0544
    ## [5,]    0 0.000    0  0.00000  0.0000
    ## standard error
    ##      [,1]  [,2] [,3]   [,4]   [,5]
    ## [1,]    0 0.000    0 0.0448 0.0000
    ## [2,]    0 0.043    0 0.0000 0.0000
    ## [3,]    0 0.000    0 0.0000 0.0000
    ## [4,]    0 0.000    0 0.0000 0.0446
    ## [5,]    0 0.000    0 0.0000 0.0000
    ## AR( 2 )-matrix
    ##      [,1] [,2] [,3]  [,4] [,5]
    ## [1,]    0    0    0 0.000    0
    ## [2,]    0    0    0 0.000    0
    ## [3,]    0    0    0 0.000    0
    ## [4,]    0    0    0 0.134    0
    ## [5,]    0    0    0 0.000    0
    ## standard error
    ##      [,1] [,2] [,3]   [,4] [,5]
    ## [1,]    0    0    0 0.0000    0
    ## [2,]    0    0    0 0.0000    0
    ## [3,]    0    0    0 0.0000    0
    ## [4,]    0    0    0 0.0444    0
    ## [5,]    0    0    0 0.0000    0
    ## AR( 3 )-matrix
    ##      [,1] [,2] [,3]    [,4] [,5]
    ## [1,]    0    0    0  0.0000    0
    ## [2,]    0    0    0  0.0000    0
    ## [3,]    0    0    0  0.0000    0
    ## [4,]    0    0    0  0.0000    0
    ## [5,]    0    0    0 -0.0136    0
    ## standard error
    ##      [,1] [,2] [,3]   [,4] [,5]
    ## [1,]    0    0    0 0.0000    0
    ## [2,]    0    0    0 0.0000    0
    ## [3,]    0    0    0 0.0000    0
    ## [4,]    0    0    0 0.0000    0
    ## [5,]    0    0    0 0.0442    0
    ##
    ## Residuals cov-mtx:
    ##                               resi
    ##      1.0023482 0.5013753 0.2252566 0.1214393 0.1070205
    ##      0.5013753 0.9368659 0.4514411 0.2163830 0.1840746
    ## resi 0.2252566 0.4514411 0.9987328 0.4199933 0.3350474
    ##      0.1214393 0.2163830 0.4199933 0.9649700 0.4220635
    ##      0.1070205 0.1840746 0.3350474 0.4220635 0.9908030
    ##
    ## det(SSE) =  0.3285515
    ## AIC =  -1.093492
    ## BIC =  -1.05204
    ## HQ  =  -1.077242

    print(model$coef)

    ##               [,1]     [,2] [,3]        [,4]        [,5]
    ##  [1,]  0.000000000 0.000000    0  0.00000000  0.00000000
    ##  [2,]  0.000000000 0.235388    0  0.00000000  0.00000000
    ##  [3,]  0.000000000 0.000000    0  0.00000000  0.00000000
    ##  [4,] -0.003523185 0.000000    0  0.00000000  0.00000000
    ##  [5,]  0.000000000 0.000000    0 -0.05440679  0.00000000
    ##  [6,]  0.000000000 0.000000    0  0.00000000  0.00000000
    ##  [7,]  0.000000000 0.000000    0  0.00000000  0.00000000
    ##  [8,]  0.000000000 0.000000    0  0.00000000  0.00000000
    ##  [9,]  0.000000000 0.000000    0  0.13366447  0.00000000
    ## [10,]  0.000000000 0.000000    0  0.00000000  0.00000000
    ## [11,]  0.000000000 0.000000    0  0.00000000  0.00000000
    ## [12,]  0.000000000 0.000000    0  0.00000000  0.00000000
    ## [13,]  0.000000000 0.000000    0  0.00000000  0.00000000
    ## [14,]  0.000000000 0.000000    0  0.00000000 -0.01360898
    ## [15,]  0.000000000 0.000000    0  0.00000000  0.00000000

    VARpred(model,1)

    ## orig  511
    ## Forecasts at origin:  511
    ##         V1         V2         V3         V4         V5
    ##  0.0042468 -0.1575328  0.0000000  0.1612303 -0.0007172
    ## Standard Errors of predictions:
    ##                 resi
    ## 1.0012 0.9679 0.9994 0.9823 0.9954
    ## Root mean square errors of predictions:
    ##                 resi
    ## 1.0158 0.9820 1.0139 0.9966 1.0099