R/BigVARSupportFunctions.R

In BigVAR: Dimension Reduction Methods for Multivariate Time Series

### Support functions for BigVAR Package These are mostly utility functions that will not be seen by the user




#' Construct a VAR or VARX lag matrix
#' 
#' @param Y a \eqn{T \times k} matrix of endogenous (modeled) series
#' @param X a \eqn{T \times m} matrix of exogenous (unmodeled) series (default NULL)
#' @param p Endogenous Lag order
#' @param s exogenous lag order (default zero)
#' @param oos indicator as to whether the data should be constructed for out of sample prediction (i.e. last available entries of Y as final lags default FALSE)
#' @param contemp indicator as to whether to use contemporaneous exogenous predictors (for example, if exogenous series become available before exogenous default FALSE).
#' @return list with two entries:
#' \itemize{
#' \item{'Z'}{\eqn{kp+ms+1\times T-max(p,s)} VARX lag matrix}
#' \item{'Y'}{adjusted \eqn{k\times T-max(p,s)} endogenous series}
#' }
#' @details This function is not required unless you which to design your own cross validation routine. 
#' @references
#' See page 15 of Lutkepohl, 'A New Introduction to Multiple Time Series Analysis
#' @seealso \code{\link{MultVarSim}}
#' @examples
#' data(Y)
#' # construct VAR lag matrix with p=4
#' ZZ<-VARXLagCons(Y,X=NULL,p=4,s=0)
#' @export
VARXLagCons <- function(Y, X = NULL, p, s = 0, oos = FALSE, contemp = FALSE) {
    if (is.null(X)) {
        X <- matrix(0, nrow = nrow(Y))
    }
    if (nrow(Y) != nrow(X)) {
        stop("Y and X must have same dimensions")
    }
    if (s == 0 & !contemp) {
        m <- 0
    } else {
        m <- ncol(X)
    }

    if (p < 0 | m < 0) {
        stop("lag orders must be positive")
    }
    k <- ifelse(is.null(Y), 0, ncol(Y))
    XX <- VARXCons(Y, X, k, p, m, s, oos, contemp)
    Y <- t(Y[(max(c(p, s)) + 1):nrow(Y), ])
    return(list(Z = XX, Y = Y))
}

# VARX Construction for constructModel
VARXConsModel <- function(Y, p, VARX, tf) {
    k <- ncol(Y)
    if (length(VARX) != 0) {

        VARXI <- TRUE
        k1 <- VARX$k
        s <- VARX$s

        if(exists('contemp',where=VARX)){
            if(!is.logical(VARX$contemp)){
                stop("contemp must be logical")
            }
            if(VARX$contemp){
                contemp <- TRUE
                s1 <- 1
            }else{
                contemp <- FALSE
                s1 <- 0
            }
        }else{
            contemp <- FALSE
            s1 <- 0
        }

        m <- k - k1
        Y1 <- matrix(Y[, 1:k1, drop = F], ncol = k1)
        X <- matrix(Y[, (ncol(Y) - m + 1):ncol(Y), drop = F], ncol = m)

        if (!tf) {
            trainZ <- VARXCons(Y1, X, k1, p, m, s, contemp = contemp)

        } else {
            trainZ <- VARXCons(matrix(0, ncol = 1, nrow = nrow(X)), matrix(X, ncol = m), k = 0, p = 0, m = m, s = s, contemp = contemp,
                oos = FALSE)

        }

        trainZ <- trainZ[2:nrow(trainZ), , drop = F]

        trainY <- matrix(Y[(max(c(p, s)) + 1):nrow(Y), 1:k1, drop = F], ncol = k1)

    } else {
        # VAR setting
        VARXI <- FALSE
        contemp <- FALSE
        s1 <- 0

        Z1 <- VARXCons(Y, matrix(0, nrow = nrow(Y)), k, p, 0, 0)
        trainZ <- Z1[2:nrow(Z1), , drop = FALSE]
        trainY <- matrix(Y[(p + 1):nrow(Y), , drop = FALSE], ncol = k)

    }
    return(list(trainY = trainY, trainZ = trainZ, s1 = s1, contemp = contemp, VARXI = VARXI))
}


.check_is_matrix <- function(x) {
    !is.null(attr(x, "dim"))
}

.huber_loss <- function(r, delta) {
    l <- ifelse(abs(r) < delta, 1/2 * abs(r)^2, delta * (abs(r) - 1/2 * delta))
    return(l)
}

.calc.loss <- function(x, univ = FALSE, loss, delta) {
    if (loss == "L1") {
        l <- sum(abs(x))
    } else if (loss == "L2") {
        if (univ) {
            l <- x^2
        } else {
            l <- norm2(x)^2
        }

    } else if (loss == "Huber") {
        l <- .huber_loss(x, delta)

    }
    return(sum(l))
}
# mean benchmark
.evalMean <- function(Y, T1, T2, h = 1, loss = "L2", delta = 2.5) {

    ypredF <- NULL
    if (!is.matrix(Y)) {
        Y <- matrix(Y, ncol = 1)
    }

    MSFE <- c()

    k <- ncol(Y)

    for (u in (T1 - h + 2):T2) {

        if (h + u - 1 > T2) {
            break
        }

        trainY1 <- Y[1:(u - 1), ]

        if (k > 1) {

            ypred <- colMeans(trainY1)
            ypredF <- rbind(ypredF, ypred)
        } else {
            ypred <- mean(trainY1)
            ypredF <- c(ypredF, ypred)
        }

        uhat <- matrix(Y[u + h - 1, ] - ypred, ncol = k)

        MSFE <- c(MSFE, .calc.loss(uhat, univ = FALSE, loss, delta))

    }
    ypredF <- unname(ypredF)
    return(list(Mean = mean(na.omit(MSFE)), SD = sd(na.omit(MSFE))/sqrt(length(na.omit(MSFE))), preds = as.matrix(ypredF)))
}

# random walk benchmark
.evalRW <- function(Y, T1, T2, h = 1, loss = "L2", delta = 2.5) {

    if (!is.matrix(Y)) {
        Y <- matrix(Y, ncol = 1)
    }
    ypredF <- NULL
    MSFE <- c()

    k <- ncol(Y)

    for (u in (T1 - h + 2):T2) {

        if (h + u - 1 > T2) {
            break
        }

        trainY1 <- Y[u - 1, ]
        ypredF <- rbind(ypredF, trainY1)
        uhat <- matrix(Y[u + h - 1, ] - trainY1, ncol = k)

        MSFE <- c(MSFE, .calc.loss(uhat, univ = FALSE, loss, delta))

    }
    ypredF <- unname(ypredF)
    return(list(Mean = mean(na.omit(MSFE)), SD = sd(na.omit(MSFE))/sqrt(length(na.omit(MSFE))), preds = as.matrix(ypredF)))
}

# Construct Lambda Grid:
.LambdaGrid <- function(gran1, gran2, groups, Y, Z, group, p, k1, s, m, k, MN, alpha, C, intercept, tol, VARX = FALSE, separate_lambdas = FALSE,
    verbose = FALSE, gamma = 3, linear) {

    nseries <- ifelse(VARX == TRUE, k1, k)

    if (group == "Lag") {

        mat <- list()

        for (i in seq_len(length(groups))) {

            if (k > 1) {

                mat[[i]] <- norm2(Z[groups[[i]], ] %*% Y)

            } else {
                mat[[i]] <- norm2(t(Y) %*% Z[groups[[i]], ])
            }

        }

        lambdastart <- max(unlist(mat))

    }

    if (group == "Basic" | group == "BasicEN" | group == "Tapered") {
        if (!separate_lambdas) {
            if (group == "Basic" | group == "Tapered") {                
                lambdastart <- max(abs(t(Y) %*% t(Z)))
            } else {
                lambdastart <- max(abs(t(Y) %*% t(Z)))/max(c(alpha, 0.01))
            }
        } else {
            lambdastart <- c()
            for (i in 1:nseries) {

                if (group == "Basic" | group == "Tapered") {
                  lambdastart[i] <- max(abs(t(Y[, i, drop = F]) %*% t(Z)))
                } else {
                  lambdastart[i] <- max(abs(t(Y[, i, drop = F]) %*% t(Z)))/max(c(alpha, 0.01))
                }
            }

        }
    }


    if (group == "MCP" | group == "SCAD") {
        if (!separate_lambdas) {
            lambdastart <- max(abs(crossprod(t(Z), Y)))/nrow(Y)
        } else {
            lambdastart <- c()
            for (i in 1:nseries) {
                lambdastart[i] <- max(abs(crossprod(t(Z), Y[, i])))/nrow(Y)
            }
        }
    }


    if (group == "SparseLag") {

        mat <- list()

        if (alpha > 0) {

            for (i in seq_len(length(groups))) {

                if (k > 1) {

                  mat[[i]] <- norm2(Z[groups[[i]], ] %*% Y * (1/(alpha)))

                } else {

                  mat[[i]] <- norm2(t(Y) %*% Z[groups[[i]], ])

                }

            }

            lambdastart <- max(unlist(mat))

        } else {

            lambdastart <- max(t(Y) %*% t(Z))

        }


    }

    if (group == "OwnOther") {

        mat <- list()

        ZZ <- kronecker(t(Z), diag(nseries))

        for (i in seq_len(length(groups))) {

            mat[[i]] <- norm(as.vector(t(Y)) %*% ZZ[, groups[[i]]]/sqrt(length(groups[[i]])), "F")

        }

        lambdastart <- max(unlist(mat))

    }

    if (group == "SparseOO") {

        mat <- list()

        ZZ <- kronecker(t(Z), diag(nseries))

        if (alpha > 0) {
            for (i in seq_len(length(groups))) {

                mat[[i]] <- norm(1/(alpha) * as.vector(t(Y)) %*% ZZ[, groups[[i]]]/sqrt(length(groups[[i]])), "F")
            }

            lambdastart <- max(unlist(mat))

        } else {

            lambdastart <- max(t(Y) %*% t(Z))

        }

    }

    if (group == "EFX") {

        gmax <- c()

        for (i in 1:nseries) {

            gmax[i] <- norm2(Z %*% Y[, i])/sqrt(k * p)

        }

        lambdastart <- max(gmax)

    }

    if (group == "HLAGC" | group == "HLAGOO" | group == "HLAGELEM") {

        gmax <- c()

        for (i in 1:nseries) {

            gmax[i] <- norm2(Z %*% Y[, i])

        }

        if (!separate_lambdas) {
            lambdastart <- max(gmax)
        } else {
            lambdastart <- gmax
        }
    }

    if (VARX) {
        beta <- array(0, dim = c(k1, k1 * p + s * m + 1, 1))
    } else if (group == "Tapered") {

        beta <- array(0, dim = c(k, k * p + 1, 1))

    } else {
        beta <- array(0, dim = c(k, k * p + 1, 1))
    }

    if (!separate_lambdas) {
        lambdastart <- LGSearch(lambdastart, Y, Z, beta, group, k1, p, s, m, groups, k, MN, alpha, C, intercept, tol, VARX,
            gamma)

        if (!linear) {
            lambda <- exp(seq(from = log(lambdastart), to = log(lambdastart/gran1), length = gran2))
        } else {
            lambda <- seq(from = lambdastart, to = lambdastart/gran1, length = gran2)

        }

    } else {

        lambda <- matrix(NA, nrow = c(gran2), ncol = ncol(Y))

        for (i in seq_len(ncol(lambda))) {
            lambdastart[i] <- LGSearch(lambdastart[i], Y, Z, beta, group, k1, p, s, m, groups, k, MN, alpha, C, intercept,
                tol, VARX, gamma)
            lambdastart[i] <- ifelse(lambdastart[i] == 0, 1e-04, lambdastart[i])
            if (verbose & i%%20 == 0) {
                print(sprintf("determined lambda grid for series %s", i))
            }

            if (!linear) {
                lambda[, i] <- exp(seq(from = log(lambdastart[i]), to = log(lambdastart[i]/gran1), length = gran2))
            } else {
                lambda[, i] <- seq(from = lambdastart[i], to = lambdastart[i]/gran1, length = gran2)
            }

        }


    }


    return(lambda)

}

#' Converts a VAR coefficient matrix of order p to multiple companion form
#' 
#' @param B a \eqn{k \times kp} coefficient matrix
#' @param p Lag order
#' @param k Number of Series
#' @return Returns a \eqn{kp \times kp} coefficient matrix representing all coefficient matrices contained in Ai as a VAR(1).
#' @references See page 15 of Lutkepohl, 'A New Introduction to Multiple Time Series Analysis'
#' @seealso \code{\link{MultVarSim}}
#' @examples
#' k=3;p=6
#' B=matrix(0,nrow=k,ncol=p*k)
#' A1<- matrix(c(.4,-.02,.01,-.02,.3,.02,.01,.04,.3),ncol=3,nrow=3)
#' A2 <- matrix(c(.2,0,0,0,.3,0,0,0,.13),ncol=3,nrow=3)
#' B[,1:k]=A1
#' B[,(4*k+1):(5*k)]=A2
#' A <- VarptoVar1MC(B,p,k)
#' @export
VarptoVar1MC <- function(B, p, k) {

    Fp <- matrix(0, nrow = k * p, ncol = k * p)

    Fp[1:k, ] <- B

    Fp[-(1:k), 1:(k * (p - 1))] <- diag(k * (p - 1))
    # We require that the coefficient matrix generates a stationary VAR
    if (max(Mod(eigen(Fp)$values)) > 1) {
        warning("Coefficient Matrix is not stationary")
    }

    return(Fp)

}


#' Simulate a VAR
#' 
#' @param k Number of Series
#' @param A1 Either a \eqn{k \times k} coefficient matrix or a \eqn{kp \times kp} matrix created using \code{\link{VarptoVar1MC}}. 
#' @param p Maximum Lag Order
#' @param Sigma Residual Covariance Matrix of dimension \eqn{k\times k}
#' @param T Number of simulations
#' @return Returns a \eqn{T \times k} of realizations from a VAR.
#' @references Lutkepohl, 'A New Introduction to Multiple Time Series Analysis'
#' @seealso \code{\link{VarptoVar1MC}}
#' @examples
#' k=3;p=6
#' B=matrix(0,nrow=k,ncol=p*k)
#' A1<- matrix(c(.4,-.02,.01,-.02,.3,.02,.01,.04,.3),ncol=3,nrow=3)
#' A2 <- matrix(c(.2,0,0,0,.3,0,0,0,.13),ncol=3,nrow=3)
#' B[,1:k]=A1
#' B[,(4*k+1):(5*k)]=A2
#' A <- VarptoVar1MC(B,p,k)
#' Y <-MultVarSim(k,A,p,.1*diag(k),100)
#' @export
#' @importFrom MASS mvrnorm
MultVarSim <- function(k, A1, p, Sigma, T) {

    if (max(Mod(eigen(A1)$values)) > 1) {
        stop("Error: Generator Matrix is not stationary")
    }

    # add 500 observations for initialization purposes

    Y <- matrix(0, nrow = T + 500 + p, ncol = k)

    YY <- as.vector(Y)

    for (i in seq(from = (k * p + 1), to = (nrow(Y) * k - 1), by = k)) {

        u <- as.vector(c(mvrnorm(1, rep(0, k), Sigma), rep(0, k * p - k)))

        YY[(i + k):(i - k * p + 1 + k)] <- A1 %*% YY[(i):(i - k * p + 1)] + as.matrix(u, ncol = 1)

    }

    YY <- YY[YY != 0]

    Y <- matrix(YY, ncol = k, byrow = TRUE)

    Y <- Y[, c(ncol(Y):1)]

    Y <- Y[501:nrow(Y), ]

    return(Y)
}

# function to create subsets for lag group VARX-L
.groupfun <- function(p, k) {

    jjj <- list()
    jjj[[1]] <- 1:k

    if (p > 1) {

        for (i in 2:p) {
            jjj[[i]] <- jjj[[i - 1]] + k

        }

    }

    return(jjj)

}

# C++ groupings to account for 0 indexing
.groupfuncpp <- function(p, k) {
    jjj <- list()

    jjj[[1]] <- 0:(k - 1)

    if (p > 1) {
        for (i in 2:p) {

            jjj[[i]] <- jjj[[i - 1]] + k
        }
    }
    return(jjj)
}

# subsetting complement of groups in rcpp
.groupfuncomp <- function(p, k) {


    ownoth <- .groupfuncpp(p, k)

    kk2 <- list()

    pmax <- max(unlist(ownoth))

    to <- 0:(pmax)

    for (i in seq_len(length(ownoth))) {

        kk2[[i]] <- to[is.na(pmatch(to, ownoth[[i]]))]
    }

    return(kk2)

}

# Group indexing for own/other VARX-L
.lfunction2 <- function(p, k) {

    kk <- list()

    kk[[1]] <- 1:(k^2)

    if (p > 1) {
        for (i in 2:p) {
            kk[[i]] <- 1:(k^2) + tail(kk[[i - 1]], 1)

        }
    }
    return(kk)
}

.lfunction2cpp <- function(p, k) {

    kk <- list()

    kk[[1]] <- 0:(k^2 - 1)

    if (p > 1)

    {

        for (i in 2:p) for (i in 2:p) {

            kk[[i]] <- 0:(k^2 - 1) + tail(kk[[i - 1]], 1) + 1


        }

    }

    return(kk)


}

.lfunction3 <- function(p, k) {

    kk <- .lfunction2(p, k)

    oo <- list()

    pp <- list()

    for (i in seq_len(length(kk))) {
        j <- 0
        oo[[i]] <- kk[[i]][(seq(1, length(kk[[1]]), k + 1) + (j * k^2))]
        pp[[i]] <- kk[[i]][-(seq(1, length(kk[[1]]), k + 1) + (j * k^2))]
        j <- j + 1

    }

    ownoth <- c(oo, pp)
    return(ownoth)
}

.lfunction3cpp <- function(p, k) {

    kk <- .lfunction2cpp(p, k)
    oo <- list()
    pp <- list()

    for (i in seq_len(length(kk))) {
        j <- 0
        oo[[i]] <- kk[[i]][(seq(1, length(kk[[1]]), k + 1) + (j * k^2))]
        pp[[i]] <- kk[[i]][-(seq(1, length(kk[[1]]), k + 1) + (j * k^2))]
        j <- j + 1
    }

    ownoth <- c(oo, pp)

    return(ownoth)

}


.lfunctioncomp <- function(p, k) {

    ownoth <- .lfunction3cpp(p, k)
    kk2 <- list()
    pmax <- max(unlist(ownoth))
    to <- 0:(pmax)
    for (i in seq_len(length(ownoth))) {
        kk2[[i]] <- to[is.na(pmatch(to, ownoth[[i]]))]
    }

    return(kk2)
}


# This function should work for arbitrary groups
.lfunction <- function(groups, p) {

    H <- as.vector(do.call("cbind", groups))
    kk <- list()
    kk[[1]] <- H
    if (p > 1) {
        for (i in 2:p) {
            kk[[i]] <- as.vector(do.call("cbind", groups)) + tail(kk[[i - 1]], 1)
        }
    }
    return(kk)

}

# Indexing for HLAG
.vsubs <- function(p, k) {
    vi <- list()

    for (i in p:1) {
        g <- max(k * i - k, 0)
        vi[[i]] <- g:(k * (p) - 1)
    }
    return(vi)

}

# indexing for HLAG OO
.oofun <- function(p, k) {

    kk <- .lfunction2(p, k)
    oo <- list()
    pp <- list()

    for (i in seq_len(length(kk))) {
        j <- 0
        oo[[i]] <- kk[[i]][(seq(1, length(kk[[1]]), k + 1) + (j * k^2))]
        pp[[i]] <- kk[[i]][-(seq(1, length(kk[[1]]), k + 1) + (j * k^2))]
        j <- j + 1
    }
    ownoth <- list()
    jj <- .lfunction3(p, k)

    for (i in seq_len(length(jj))) {

        if (i == 1)

        {

            ownoth[[i]] <- jj[[i]]
            oo[[1]] <- NULL

        }

        if (i == length(jj)) {
            ownoth[[i]] <- tail(jj, 1)
            pp[[1]] <- NULL
        }

        if (i != 1 & i%%2 != 0) {

            ownoth[[i]] <- head(oo, 1)
            oo[[1]] <- NULL
        }

        if (i != length(jj) & i%%2 == 0) {
            ownoth[[i]] <- head(pp, 1)
            pp[[1]] <- NULL
        }


    }

    return(rev(ownoth))

}

.oocumfun <- function(p, k) {

    kk <- rev(.oofun(p, k))
    oogroups <- list()
    oogroups[[1]] <- unlist(kk)
    for (i in 2:length(kk)) {
        oogroups[[i]] <- unlist(kk[-(1:(i - 1))])

    }

    return(oogroups)
}


# indexing function for Own/Other HLAG
.vecoovars <- function(p, k, k1) {

    vv <- list()

    vv[[1]] <- 1:(p * k)

    vv[[2]] <- vv[[1]][-k1]

    q1 <- 1

    if (p > 1) {
        for (i in 3:(2 * p)) {
            if (i%%2 != 0) {
                vv[[i]] <- (q1 * k + 1):(k * p)
                q1 <- q1 + 1
            } else {
                vv[[i]] <- vv[[i - 1]][-k1]
            }
        }
    }
    return(vv)

}

# indexing to start at zero for use within rcpp
.vecoovarscpp <- function(p, k, k1) {

    vv <- list()
    vv[[1]] <- 0:(p * k - 1)
    vv[[2]] <- vv[[1]][-(k1)]
    q1 <- 1

    if (p > 1) {
        for (i in 3:(2 * p)) {
            if (i%%2 != 0) {
                vv[[i]] <- (q1 * k):(k * p - 1)
                q1 <- q1 + 1
            } else {

                vv[[i]] <- vv[[i - 1]][-(k1)]

            }

        }

    }

    return(vv)

}


# VARX Lag Group function
groupfunVARX <- function(p, k, k1, s) {

    jj <- list()
    m <- k - k1
    jj <- .groupfuncpp(p, k1)
    kp <- k1 * p + m * s - 1
    jj2 <- list()
    startjj <- max(unlist(jj)) + 1
    for (i in seq(startjj, kp, by = 1)) {
        jj[[i]] <- i
    }
    jj[sapply(jj, is.null)] <- NULL

    return(jj)
}

groupfunVARXcomp <- function(p, k, k1, s) {
    ownoth <- groupfunVARX(p, k, k1, s)
    kk2 <- list()
    pmax <- max(unlist(ownoth))
    to <- 0:(pmax)
    for (i in seq_len(length(ownoth))) {
        kk2[[i]] <- to[is.na(pmatch(to, ownoth[[i]]))]
    }
    return(kk2)

}


diaggroupfunVARX <- function(p, k, k1, s) {
    m <- k - k1
    jj <- list()
    jj <- .lfunction3cpp(p, k1)
    kp <- k1 * (p * k1 + s * m) - 1
    jj2 <- list()
    startjj <- max(unlist(jj)) + 1

    for (i in seq(startjj, kp, by = k1)) {
        jj[[i]] <- i:(i + k1 - 1)

    }

    jj[sapply(jj, is.null)] <- NULL


    return(jj)

}

diaggroupfunVARXcomp <- function(p, k, k1, s) {

    ownoth <- diaggroupfunVARX(p, k, k1, s)
    kk2 <- list()
    pmax <- max(unlist(ownoth))
    to <- 0:(pmax)

    for (i in seq_len(length(ownoth))) {
        kk2[[i]] <- to[is.na(pmatch(to, ownoth[[i]]))]
    }

    return(kk2)
}


diaggroupfunVARXL <- function(p, k, k1) {
    jj <- list()
    jj <- .lfunction3cpp(p, k1)
    kp <- k1 * p * k - 1
    jj2 <- list()
    startjj <- max(unlist(jj)) + 1

    for (i in seq(startjj, kp, by = 1)) {
        jj[[i]] <- i

    }

    jj[sapply(jj, is.null)] <- NULL
    return(jj)

}

diaggroupfunVARXcompL <- function(p, k, k1) {

    ownoth <- diaggroupfunVARXL(p, k, k1)
    kk2 <- list()
    pmax <- max(unlist(ownoth))

    to <- 0:(pmax)

    for (i in seq_len(length(ownoth))) {

        kk2[[i]] <- to[is.na(pmatch(to, ownoth[[i]]))]

    }

    return(kk2)

}


# iterative procedure to find a tighter bound for lambda starting value via binary search
LGSearch <- function(gstart, Y, Z, BOLD, group, k1, p, s, m, gs, k, MN, alpha, C, intercept, tol, VARX, gamma) {
    s1 <- 0
    palpha <- NULL
    tk <- 1/max(Mod(eigen(Z %*% t(Z))$values))
    lambdah <- gstart
    lambdal <- 0
    activeset <- list(rep(rep(list(0), length(gs))))

    gran2 <- 1
    grps <- create_group_indexes(group, p, k, gran2 * length(alpha), VARX, k1, s)
    groups <- grps$groups
    compgroups <- grps$compgroups
    activeset <- grps$activeset
    starting_eigvals <- grps$starting_eigvals
    nseries <- nrow(BOLD)
    while (max(abs(lambdah - lambdal)) > 10 * tol) {
        lambda <- (lambdah + lambdal)/2
        dual <- FALSE
        separate_lambdas <- FALSE
        temp <- .BigVAR.fit(group, BOLD, Z, Y, lambda, tol, p, m, k1, k, s, s1, MN, C, intercept, separate_lambdas, dual,
            activeset, starting_eigvals, groups, compgroups, VARX, alpha, palpha, gamma)
        # remove intercept from consideration
        if (group == "Tapered") {
            param <- adrop(BOLD[, -1, , drop = F], drop = 3)

        } else {

            BOLD <- temp$beta
            param <- adrop(BOLD[, -1, , drop = F], drop = 3)
        }

        activeset <- temp$activeset
        q1a <- temp$q1a

        if (MN) {
            submat <- param[1:nseries, 1:nseries, drop = F]
            diag(submat) <- ifelse(C == 0, diag(param[1:nseries, 1:nseries, drop = F]), diag(param[1:nseries, 1:nseries,
                drop = F]) - C)
            param[1:nseries, 1:nseries] <- submat
            BOLD[, -1, ] <- param

        }
        if (max(abs(param)) < tol) {
            lambdah <- lambda
        } else {
            lambdal <- lambda
        }
        
    }
    lambdah
}

#' Evaluate forecasts from a VAR or VARX with lag orders selected by AIC/BIC
#' 
#' @param Y a \eqn{T \times k} multivariate time series 
#' @param X a \eqn{T \times m} multivariate time series of unmodeled exogenous variables
#' @param p maximum lag order for endogenous series
#' @param s maximum lag order for exogenous series
#' @param T1 start of forecast evaluation period.
#' @param T2 end of forecast evaluation period
#' @param IC specifies whether to select lag order according to 'AIC' or 'BIC'
#' @param h desired forecast horizon
#' @param loss loss function (default 'L2', one of 'L1','L2','Huber')
#' @param delta delta for Huber loss function (default 2.5)
#' @param iterated indicator as to whether to use iterated or direct multistep forecasts (if applicable, VAR context only)
#' @return Returns the one-step ahead MSFE as well as the forecasts over the evaluation period and lag order selected.
#' @details This function evaluates the one-step ahead forecasts of a VAR or VARX fit by least squares over an evaluation period.  At every point in time, lag orders for the endogenous and exogenous series are selected according to AIC or BIC.  This function is run automatically when \code{\link{cv.BigVAR}} is called unless \code{ic} is set to \code{FALSE} in \code{\link{constructModel}}.      
#' @references Neumaier, Arnold, and Tapio Schneider. 'Estimation of parameters and eigenmodes of multivariate autoregressive models.' ACM Transactions on Mathematical Software (TOMS) 27.1 (2001): 27-57.
#' @seealso \code{\link{VARXFit}},\code{\link{constructModel}}, \code{\link{cv.BigVAR}}
#' @examples
#' data(Y)
#'
#' # Evaluate the performance of a VAR with lags selected by BIC.
#' p <- 4
#' T1 <- floor(nrow(Y))/3
#' T2 <- floor(2*nrow(Y))/3
#' # Matrix of zeros for X
#' X <- matrix(0,nrow=nrow(Y),ncol=ncol(Y))
#' BICMSFE <- VARXForecastEval(Y,X,p,0,T1,T2,'BIC',1)
#' 
#' @export
VARXForecastEval <- function(Y, X, p, s, T1, T2, IC, h, iterated = FALSE, loss = "L2", delta = 2.5) {


    if (T1 > nrow(Y) | T2 > nrow(Y) | T2 < T1) {
        stop("Training dates exceed series length")
    }

    if (!IC %in% c("AIC", "BIC")) {

        stop("IC must either be AIC or BIC")

    }

    MSFE <- c()
    predF <- NULL
    pvec <- NULL
    svec <- NULL
    k <- ncol(Y)
    m <- ifelse(s != 0, ncol(X), 0)
    for (i in (T1 - h + 2):T2) {

        if (h + i - 1 > T2) {
            break
        }

        testY <- as.matrix(Y[1:(i - 1), ])
        testX <- as.matrix(X[1:(i - 1), ])

        if (!iterated) {
            hd <- h

        } else {

            hd <- 1

        }
        if (IC == "BIC") {
            popt <- ICX(testY, testX, k, p, s, m, "BIC", h = hd)
        }
        if (IC == "AIC") {

            popt <- ICX(testY, testX, k, p, s, m, "AIC", h = hd)
        }
        B1 <- popt$B

        if (popt$p == 0 & popt$s == 0) {

            eZ <- matrix(rep(1, 1), ncol = 1)

            pred <- B1 %*% eZ

        } else {


            C <- max(popt$p, popt$s)

            ## # possibly memory leak in VARX lag matrix construction in Eigen if maxlag is 1.  # to be on the safe
            ## side, we will perform it in R
            eZ <- VARXCons(as.matrix(Y[(i - C):(i), ]), as.matrix(X[(i - C):(i), ]), k, popt$p, m, popt$s)

            ## }

            pred <- B1 %*% eZ

            # iterated multistep forecasts (if VAR and horizon greater than 1)
            if (h > 1 & s == 0 & iterated) {

                pred <- predictMS(matrix(pred, nrow = 1), Y, h - 1, B1, C, FALSE)

            }

        }
        predF <- rbind(predF, t(pred))
        MSFEi <- .calc.loss(Y[i + h - 1, ] - pred, univ = FALSE, loss, delta)

        MSFE <- c(MSFE, MSFEi)
        svec <- c(svec, popt$s)
        pvec <- c(pvec, popt$p)
    }

    return(list(MSFE = MSFE, pred = as.matrix(predF), p = pvec, s = svec))


}

#' Fit a VAR or VARX model by least squares
#' 
#' @param Y a \eqn{t \times k} multivariate time series
#' @param p maximum lag order
#' @param IC Information criterion indicator, if set to \code{NULL}, it will fit a least squares VAR(X) of orders p and s.  Otherwise, if set to 'AIC' or 'BIC' it return the model with lag orders that minimize the given IC. 
#' @param VARX a list of VARX specifications (as in \code{\link{constructModel}} (or NULL )
#' @return Returns a list with four entries:
#' \itemize{
#' \item{'Bhat'}{Estimated \eqn{k\times kp+ms} coefficient matrix}
#' \item{'SigmaU}{Estimated \eqn{k\times k} residual covariance matrix}
#' \item{'phat'}{Selected lag order for VAR component}
#' \item{'shat'}{Selected lag order for VARX component}
#' \item{'Y'}{multivariate time series retained for prediction purposes}
#' \item{'Y'}{number of endogenous (modeled) time series}
#' }
#' @details This function uses a modified form of the least squares technique proposed by Neumaier and Schneider (2001).  It fits a least squares VAR or VARX via a QR decomposition that does not require explicit matrix inversion.  This results in improved computational performance as well as numerical stability over the conventional least squares approach. 
#' @references Neumaier, Arnold, and Tapio Schneider. 'Estimation of parameters and eigenmodes of multivariate autoregressive models.' ACM Transactions on Mathematical Software (TOMS) 27.1 (2001): 27-57.
#' @seealso \code{\link{constructModel}}, \code{\link{cv.BigVAR}},\code{\link{BigVAR.fit}}
#' @examples
#' data(Y)
#' # fit a VAR_3(3)
#' mod <- VARXFit(Y,3,NULL,NULL)
#' # fit a VAR_3 with p= 6 and lag selected according to AIC
#' modAIC <- VARXFit(Y,6,'AIC',NULL)
#' # Fit a VARX_{2,1} with p=6, s=4 and lags selected by BIC
#' modXBIC <- VARXFit(Y,6,'BIC',list(k=1,s=4))
#' 
#' @export
VARXFit <- function(Y, p, IC, VARX = NULL) {

    if (!is.null(VARX)) {

        if (is.list(VARX) & !(exists("k", where = VARX) & exists("s", where = VARX))) {

            stop("VARX Specifications entered incorrectly")

        }

    }
    if (is.list(VARX) & (length(VARX) != 0)) {

        k1 <- VARX$k
        s <- VARX$s
        Y1 <- matrix(Y[, 1:k1], ncol = k1)
        m <- ncol(Y) - k1
        X <- matrix(Y[, (k1 + 1):ncol(Y)], ncol = m)

        if (exists("contemp", where = VARX)) {
            contemp <- VARX$contemp
        } else {
            contemp <- FALSE
        }
        Z <- VARXCons(Y1, X, k1, p, m, s, contemp = contemp)
        offset <- max(p, s) + 1
        YT <- matrix(Y1[offset:nrow(Y), ], ncol = k1)
        X <- matrix(X[offset:nrow(X), ], ncol = m)
    } else {

        k <- ncol(Y)
        k1 <- k
        s <- 0
        m <- 0
        offset <- p + 1
        X <- matrix(0, nrow = nrow(Y))

        Z <- VARXCons(Y, X, k, p, m, s)
        YT <- matrix(Y[(offset):nrow(Y), ], ncol = ncol(Y))

    }
    if (is.null(IC)) {

        Res <- ARFitVARXR(cbind(t(Z), YT), k1, p, m, s)

        shat <- s
        phat <- p

    } else {
        if (!IC %in% c("AIC", "BIC")) {

            stop("IC must either be AIC,BIC, or set to NULL")

        }

        Res <- ICX(YT, X, k1, p, s, m, IC)

        shat <- Res$s
        phat <- Res$p

    }
    if (is.null(VARX)) {
        k <- ncol(Y)
    } else {
        k <- VARX$k
    }
    list(Bhat = Res$B, SigmaU = Res$SigmaU, phat = phat, shat = shat, Y = Y, k = k)

}

#' One-step ahead predictions for VARX models
#' 
#' @param VARXRes the results from \code{\link{VARXFit}}
#' @return Returns a vector consisting of the out-of-sample forecasts for the provided \code{\link{VARXFit}} model.
#' @seealso \code{\link{VARXFit}}
#' @examples
#' data(Y)
#' # fit a VAR_3(3)
#' mod <- VARXFit(Y,3,NULL,NULL)
#' pred <-PredictVARX(mod)
#' 
#' @export
PredictVARX <- function(VARXRes) {

    B <- VARXRes$Bhat
    Y <- VARXRes$Y
    k <- VARXRes$k
    m <- ncol(Y) - k

    if (k < ncol(Y)) {
        Z <- VARXCons(Y[, 1:k, drop = FALSE], Y[, (k + 1):ncol(Y), drop = FALSE], k, VARXRes$phat, m, VARXRes$shat, oos = TRUE)
    } else {
        Z <- VARXCons(Y[, 1:k, drop = FALSE], matrix(0, nrow = nrow(Y)), k, VARXRes$phat, m, VARXRes$shat, oos = TRUE)
    }

    return(as.numeric(tail(t(B %*% Z), 1)))


}

# Recursive multi-step predictions
predictMS <- function(pred, Y, n.ahead, B, p, MN = FALSE, predict_all = FALSE, n.ahead_full = n.ahead) {

    # Augment Y with predictions, create lag matrix (no intercept if MN)
    Y <- rbind(Y, pred)

    Z <- VARXCons(Y, matrix(0, nrow = nrow(Y), ncol = 1), ncol(Y), p, 0, 0, oos = TRUE)

    if (MN) {
        Z <- Z[2:nrow(Z), , drop = F]
    }
    Z <- Z[, ncol(Z), drop = F]
    pred <- matrix(B %*% Z, ncol = ncol(Y), nrow = 1)

    if (n.ahead == 1) {
        if (predict_all) {
            return(rbind(Y[((nrow(Y) - n.ahead_full) + 1):nrow(Y), ], pred))
        } else {
            return(pred)
        }
    }

    predictMS(pred, Y, n.ahead - 1, B, p, MN, predict_all, n.ahead_full = n.ahead)

}

# Multi-step VARX with new data.
predictMSX <- function(pred, Y, n.ahead, B, p, newxreg, X, m, s, cumulative, MN, contemp = FALSE) {

    Y <- rbind(Y, pred)
    X <- rbind(X, matrix(newxreg[cumulative, ], ncol = m))

    if (nrow(Y) != nrow(X)) {
        stop("error, dimension issue")
    }
    if (!contemp) {
        Z <- VARXCons(as.matrix(Y), X, ncol(Y), p, m, s, oos = TRUE)
    } else {
        Z <- VARXCons(as.matrix(Y), as.matrix(X), ncol(Y), p, m, s, oos = FALSE, contemp = TRUE)
    }
    Z <- Z[, ncol(Z), drop = F]
    if (MN) {

        Z <- as.matrix(Z[2:nrow(Z), drop = F])
        pred <- matrix(B[, 2:ncol(B), drop = F] %*% Z, ncol = ncol(Y), nrow = 1)

    } else {

        pred <- matrix(B %*% Z, ncol = ncol(Y), nrow = 1)
    }
    if (n.ahead == 1) {
        return(pred)
    }

    predictMSX(pred, Y, n.ahead - 1, B, p, newxreg, X, m, s, cumulative + 1, MN)

}




# Find optimal values in 2-d gridsearch
findind <- function(opt, lambda1, lambda2) {
    if (opt < length(lambda2)) {
        lambda1ind <- 1
    } else {
        lambda1ind <- ceiling(opt/length(lambda2))
    }
    if (lambda1ind == 1) {
        jind <- opt
    } else {
        jind <- opt - (length(lambda2)) * (lambda1ind - 1)
    }
    return(c(lambda1ind, jind))
}



# Bayesian VAR with MN prior
BVARLitterman <- function(Y, Z, p, tau, mu, H, iRW) {
    T <- nrow(Y)
    k <- ncol(Y)

    # prior covariance based on univariate AR models
    sigmas <- c()
    for (i in 1:k) {
        Z1 <- VARXCons(Y[, i, drop = F], matrix(0, nrow = nrow(Y), ncol = 1), 1, p, 0, 0)
        K <- cbind(t(Z1), Y[(p + 1):nrow(Y), i])
        sigmas[i] <- sqrt(ARFitVARXR(K, 1, p, 0, 0)$SigmaU)

    }
    MMO <- colMeans(Y)

    # create prior random walk dummy
    Yrw1 <- diag(sigmas * iRW)
    Yrw2 <- matrix(0, nrow = k * (p - 1), ncol = k)
    Yrw <- tau * (rbind(Yrw1, Yrw2))
    Zrw <- tau * cbind(kronecker(diag(1:p), diag(sigmas)), matrix(0, nrow = k * p, ncol = 1))


    # create dummy for intercept
    epsilon <- 1e-05
    Ycs <- 1e-05 * matrix(0, nrow = 1, ncol = k)
    Zcs <- epsilon * cbind(matrix(0, ncol = k * p, nrow = 1), 1)

    # dummy on the sums of coefficients
    Ylr <- mu * diag(MMO * iRW)
    Zlr1 <- kronecker(matrix(1, nrow = 1, ncol = p), diag(MMO) * iRW)
    Zlr <- mu * (cbind(Zlr1, matrix(0, nrow = k, ncol = 1)))


    # Dummy for residual covariance matrix
    Ycv <- diag(sigmas)
    Zcv <- matrix(0, nrow = k, ncol = k * p + 1)

    Yprior <- rbind(Yrw, Ylr, Ycv, Ycs)
    Zprior <- rbind(Zrw, Zlr, Zcv, Zcs)

    Tstar <- nrow(Yprior)
    Z <- t(Z)
    Z <- cbind(Z[, 2:ncol(Z)], 1)

    ZZinv <- solve(t(Zprior) %*% Zprior + t(Z) %*% Z)
    ZY <- t(Zprior) %*% Yprior + t(Z) %*% Y
    beta <- ZZinv %*% ZY

    return(t(beta))

}

# grid search for BGR
BGRGridSearch <- function(Y, Z, p, grid, RWIND) {
    preds <- list()
    for (i in seq_len(length(grid))) {
        pi <- grid[i]
        mu <- pi * 0.1  # used in BGR paper
        preds[[i]] <- BVARLitterman(Y, Z, p, pi, mu, -1, RWIND)

    }
    preds <- array(unlist(preds), dim = c(nrow(preds[[1]]), ncol(preds[[1]]), length(preds)))

    return(preds)

}

# process MN prior
MN_prior <- function(Y, Z, C1) {
    C <- matrix(0, nrow = ncol(Y), ncol = nrow(Z))
    diag(C) <- C1
    Y <- t(Y)
    Y <- Y - C %*% Z
    Y <- t(Y)
    return(list(Y = Y, C = C))
}

# pre-process data (subtract intercept, adjust for shrinking toward constants)
pre_process <- function(Y, Z, C1, MN, intercept) {
    k <- ncol(Y)
    if (MN) {
        YC <- MN_prior(Y, Z, C1)
        Y <- YC$Y
        C <- YC$C
    }

    Y <- t(Y)

    if (intercept) {
        YMean <- c(apply(Y, 1, mean))
        ZMean <- c(apply(Z, 1, mean))
        if (k > 1) {
            Y <- Y - YMean %*% t(c(rep(1, ncol(Y))))
        } else {
            Y <- Y - mean(Y)
        }
        Z <- Z - ZMean %*% t(c(rep(1, ncol(Z))))

    } else {
        YMean <- rep(0, nrow(Y))
        ZMean <- rep(0, nrow(Z))
    }
    Y <- t(Y)
    return(list(Y = Y, Z = Z, YMean = YMean, ZMean = ZMean, C = C))

}


# s needs to be s+ s1 create indices for group structures
create_group_indexes <- function(group, p, k, gran2, VARX = FALSE, k1 = NULL, s = NULL) {
    starting_eigvals <- NULL
    groups <- NULL
    compgroups <- NULL
    activeset <- NULL
    if (VARX) {
        if (group == "Lag") {

            groups <- groupfunVARX(p, k, k1, s)
            compgroups <- groupfunVARXcomp(p, k, k1, s)
            activeset <- rep(list(rep(rep(list(0), length(groups)))), gran2)

        } else if (group == "SparseLag") {

            groups <- groupfunVARX(p, k, k1, s)
            compgroups <- groupfunVARXcomp(p, k, k1, s)
            starting_eigvals <- list()

            for (i in 1:(p + s)) {

                starting_eigvals[[i]] <- matrix(runif(k1, -1, 1), ncol = 1)

            }

            activeset <- rep(list(rep(rep(list(0), length(groups)))), gran2)


        } else if (group == "OwnOther") {
            groups <- diaggroupfunVARX(p, k, k1, s)
            compgroups <- diaggroupfunVARXcomp(p, k, k1, s)
            activeset <- rep(list(rep(rep(list(0), length(groups)))), gran2)


        } else if (group == "SparseOO") {
            groups <- diaggroupfunVARX(p, k, k1, s)
            compgroups <- diaggroupfunVARXcomp(p, k, k1, s)
            activeset <- rep(list(rep(rep(list(0), length(groups)))), gran2)

        }


    } else {
        if (group == "Lag") {

            groups <- .groupfuncpp(p, k)


            compgroups <- .groupfuncomp(p, k)


            activeset <- rep(list(rep(rep(list(0), length(groups)))), gran2)

        } else if (group == "SparseLag") {

            groups <- .groupfuncpp(p, k)
            compgroups <- .groupfuncomp(p, k)

            starting_eigvals <- list()

            for (i in 1:p) {

                starting_eigvals[[i]] <- matrix(runif(k, -1, 1), ncol = 1)
            }


            activeset <- rep(list(rep(rep(list(0), length(groups)))), gran2)



        } else if (group == "OwnOther") {

            groups <- .lfunction3cpp(p, k)
            compgroups <- .lfunctioncomp(p, k)
            activeset <- rep(list(rep(rep(list(0), length(groups)))), gran2)

        } else if (group == "SparseOO") {

            groups <- .lfunction3cpp(p, k)
            compgroups <- .lfunctioncomp(p, k)
            activeset <- rep(list(rep(rep(list(0), length(groups)))), gran2)
            starting_eigvals <- list()

            for (i in 1:(2 * p)) {

                starting_eigvals[[i]] <- matrix(runif(length(groups[[i]]), -1, 1), ncol = 1)
            }

        }
    }
    return(list(groups = groups, compgroups = compgroups, starting_eigvals = starting_eigvals, activeset = activeset))
}


create_lambda_grid <- function(trainY, trainZ, groups, gran1, gran2, group, p, k1, s, m, k, MN, alpha, C, intercept, tol,
    VARX, separate_lambdas, dual, gamma, linear, verbose) {
    # Constructs penalty grid if both alpha and lambda are selected
    if (dual) {
        lambda <- matrix(0, nrow = gran2, ncol = length(alpha))

        for (i in 1:length(alpha)) {

            lambda[, i] <- .LambdaGrid(gran1, gran2, groups, trainY, trainZ, group, p, k1, s, m, k, MN, alpha[i], C, intercept,
                tol, VARX = VARX, linear = linear)

        }

    } else {
        # Penalty parameter grid for just lambda
        if (group != "BGR") {

            lambda <- .LambdaGrid(gran1, gran2, groups, trainY, trainZ, group, p, k1, s, m, k, MN, alpha, C, intercept, tol,
                VARX = VARX, separate_lambdas, verbose, linear = linear)
        } else {
            # special handling for BGR
            lambda <- seq(1, 5, length = gran2)
            lambda <- lambda * sqrt(k * p)

        }
    }


    return(lambda)

}

# runs at each iteration of penalty parameter selection and validation stage
refine_and_forecast <- function(betaArray, eZ, trainZ, trainY, testY, lambda, h, recursive, MN, RVAR, refit_fraction, separate_lambdas,
    C, inds, loss, delta, k, p, k1, s, oos = FALSE) {
    if (separate_lambdas) {
        MSFE_temp <- matrix(0, nrow = dim(betaArray)[3], ncol = ncol(trainY))
    } else {
        MSFE_temp <- c()
    }
    if (MN) {

        eZ <- eZ[2:nrow(eZ), , drop = FALSE]

    }

    nlambdas <- dim(betaArray)[3]
    preds <- matrix(0, nrow = ncol(trainY), ncol = nlambdas)
    for (i in 1:nlambdas) {

        beta <- adrop(betaArray[, , i, drop = F], drop = 3)
        if (RVAR) {

            beta_rls <- RelaxedLS(cbind(t(trainZ), trainY), beta)
            beta <- (1 - refit_fraction) * beta + refit_fraction * beta_rls

        }
        if (!is.null(inds)) {
            if (i == inds) {
                beta_return <- beta
            }
        } else {
            # if we don't need beta matrix (i.e validation stage, just return the last value)
            beta_return <- beta
        }

        if (MN) {
            beta <- beta[, 2:ncol(beta)]
        }

        if (h > 1 & recursive) {
            ptemp <- beta %*% eZ

            pred <- matrix(ptemp, nrow = 1)

            pred <- predictMS(pred, trainY, h - 1, beta, p, MN)
            pred <- matrix(pred, ncol = 1)
        } else {
            pred <- beta %*% eZ
        }
        if (separate_lambdas) {
            if (!oos) {
                for (j in seq_len(ncol(testY))) {
                  MSFE_temp[i, j] <- .calc.loss(testY[, j] - pred[j, ], univ = TRUE, loss, delta)
                }
            } else {
                MSFE_temp[i] <- .calc.loss(t(testY) - beta %*% eZ, univ = FALSE, loss, delta)
            }
        } else {
            MSFE_temp[i] <- .calc.loss(t(testY) - beta %*% eZ, univ = FALSE, loss, delta)
        }


        preds[, i] <- pred
    }

    ## # just return optimal prediction in validation stage
    return(list(MSFE = MSFE_temp, betaArray = beta_return, preds = preds))
}

# determines optimal lambda from a grid of values
find_optimal_lambda <- function(MSFE, lambda, palpha, alpha, ONESE, dual, separate_lambdas, tapered) {

    if (tapered) {

        indopt <- which.min(colMeans(MSFE))

        if (indopt < length(lambda)) {

            alphaind <- 1

            alphaopt <- palpha[1]

        } else {

            alphaopt <- palpha[floor(indopt/length(lambda))]

            alphaind <- floor(indopt/length(lambda))

        }

        if (alphaind == 1) {

            lambdaopt <- lambda[indopt]

        } else if (indopt%%length(lambda) == 0) {

            lambdaopt <- lambda[length(lambda)]

        } else {
            lambdaind <- indopt - length(lambda) * alphaind
            lambdaopt <- lambda[lambdaind]

        }

        palpha <- alphaopt
        optind <- indopt
    } else {
        # one standard error correction
        if (ONESE & !dual & !separate_lambdas) {

            MSFE2 <- MSFE
            G2 <- colMeans(na.omit(MSFE2))
            G3 <- sd(na.omit(MSFE2))/sqrt(nrow(na.omit(MSFE2)))
            optind <- min(which(G2 < (min(G2) + G3)))
            lambdaopt <- lambda[optind]
        } else {

            if (!tapered & !dual) {
                # in rare cases in which MSFE is equal, the smaller penalty parameter is chosen.  This prevents
                # extremely sparse solutions
                if (separate_lambdas) {
                  if (ONESE) {
                    MSFES <- t(apply(MSFE, 3, colMeans))
                    sds <- t(apply(MSFE, 3, function(x) sd(na.omit(x))/sqrt(nrow(na.omit(x)))))
                    lambdaopt <- c()
                    optinds <- c()
                    for (i in 1:nrow(MSFES)) {
                      optinds[i] <- min(which(MSFES[i, ] < sds[i] + min(MSFES[i, ])))
                      lambdaopt[i] <- lambda[optinds[i], i, drop = F]
                    }
                    optind = optinds
                  } else {
                    MSFES <- t(apply(MSFE, 3, colMeans))
                    optinds <- apply(MSFES, 1, which.min)
                    lambdaopt <- c()
                    for (i in 1:nrow(MSFES)) {
                      lambdaopt[i] <- lambda[optinds[i], i]
                    }
                    optind = optinds


                  }
                } else {

                  optind <- max(which(colMeans(na.omit(MSFE)) == min(colMeans(na.omit(MSFE)))))
                  lambdaopt <- lambda[optind]
                }
            } else if (dual) {
                if (!ONESE) {

                  optind <- max(which(colMeans(na.omit(MSFE)) == min(colMeans(na.omit(MSFE)))))
                  inds <- findind(optind, lambda[, 1], alpha)
                } else {
                  G2 <- colMeans(na.omit(MSFE))

                  G3 <- sd(na.omit(MSFE))/sqrt(nrow(na.omit(MSFE)))

                  optind <- min(which(G2 < (min(G2) + G3)))
                  inds <- findind(optind, lambda[, 1], alpha)

                }
                lambdaopt <- lambda[inds[1], inds[2]]
                lambda <- lambda[, inds[2]]
                alphaopt <- alpha[inds[2]]
                optind <- inds

            }
        }
        if (!dual) {
            alphaopt <- alpha
        }

    }
    if (!exists("alphaopt")) {
        alphaopt <- NULL
    }
    return(list(optind = optind, lambdaopt = lambdaopt, alphaopt = alphaopt))
}



create_sigma_u <- function(Y, Z, B, h) {
    trainY <- Y[(h):(nrow(Y) - 1), , drop = F]
    trainZ <- Z[, 1:(ncol(Z) - h)]
    Sigma = 1/nrow(trainY) * sqrt(diag(tcrossprod(t(Y) - B %*% rbind(1, Z))))
    cis <- c(qnorm(0.025), qnorm(0.975))
    lower <- cis[1] * Sigma
    upper <- cis[2] * Sigma
    return(cbind(lower, upper))

}

# add C back to coefficient matrix in MN setting
adjust_mn_var <- function(beta, C) {
    for (i in 1:(dim(beta)[3])) beta[, 2:dim(beta)[2], i] <- beta[, 2:dim(beta)[2], i] + C
    return(beta)
}