coef2array <- function(B, N, lag) {
    stopifnot(all.equal(dim(B), c(N, N * lag), check.attributes=FALSE))
    aperm(structure(B, dim = c(N, N, lag)), c(3, 1, 2))
}

coef2companion <- function(B) {
    N <- nrow(B)
    rbind(B, cbind(diag(ncol(B)-N), matrix(0, ncol(B)-N, N)))
}

getEVs <- function(B)
    eigen(coef2companion(B))$values

checkStable <- function(B)
    max(abs(getEVs(B))) < 1

## following functions adapted from VAR.etp
## source code for VAR.etp package
## https://github.com/cran/VAR.etp/tree/master/R

VAR.ys2 <- function (x, b, p, e, type) {
    ## simulate y
    ## faster version of VAR.ys (not expanding matrix y but pre-allocating)
    n <- nrow(x)
    k <- nrow(b)
    b0 <- b[, 1:(k * p), drop = FALSE]
    if (type == "const")
        b1 <- b[, (k * p) + 1]
    if (type == "const+trend") {
        b1 <- b[, (k * p) + 1]
        b2 <- b[, (k * p) + 2]
    }
    y <- matrix(0, n, k)
    colnames(y) <- colnames(x)
    y[1:p,] <- x[1:p, , drop = FALSE] ## initial values are first few observations of actual series
    for (i in (p + 1):n) {
        index <- 1:k
        d1 <- 0
        for (j in 1:p) {
            d1 <- d1 + b0[, index] %*% y[i - j,]
            index <- index + k
        }
        d1 <- d1 + e[i - p, ]
        if (type == "const")
            d1 <- d1 + b1
        if (type == "const+trend")
            d1 <- d1 + b1 + b2 * i
        y[i, ] <- d1
    }
    return(y)
}

VAR.adjust2 <- function (b, bias, p, type) {
    require(VAR.etp)
    ## modified version of VAR.adjust -- throw error if NAs/infinite
    k <- nrow(b)
    bs1 <- b - bias
    delta <- VAR.etp:::VAR.modul(bs1, p)[1]
    if (delta < 1)
        bs2 <- bs1[, 1:(k * p), drop = FALSE]
    if (delta >= 1) {
        delta1 <- 1
        while (delta >= 1) {
            delta1 <- delta1 - 0.01
            bias <- delta1 * bias
            bs2 <- b[, 1:(k * p), drop = FALSE] - bias[, 1:(k *
                p), drop = FALSE]
            if (is.nan(sum(bs2)) | is.infinite(sum(bs2))) {
                print(bs2)
                stop("no bueno")
                ## bs2 <- b[, 1:(p * k), drop = FALSE]
                ## break
            }
            delta <- VAR.etp:::VAR.modul(bs2, p)[1]
        }
    }
    if (type == "const" | type == "const+trend")
        bs2 <- cbind(bs2, bs1[, (p * k + 1):ncol(b), drop = FALSE])
    return(bs2)
}

VAR.est2 <- function (x, p, type = "const", coef.only=FALSE) {
    ## much faster version than VAR.est
    n <- nrow(x)
    k <- ncol(x)
    y <- t(x[(p + 1):n, ])
    z <- matrix(0, k*p, n-p)
    xp <- t(x)
    for (i in (p+1):n)
        z[, i-p] <- xp[,(i-1):(i-p)]
    if (type == "const")
        z <- rbind(z, 1)
    if (type == "const+trend")
        z <- rbind(z, 1, (p + 1):n)
    b <- tcrossprod(y, z) %*% solve(tcrossprod(z))
    rownames(b) <- colnames(x)
    colnames(b) <- VAR.etp:::VAR.names(x, p, type)
    if (coef.only) {
        return(list(coef=b))
    } else {
        e <- y - b %*% z
        rownames(e) <- colnames(x)  # necessary for k=1
        sigu <- cov(t(e))  # tcrossprod(e)/((n - p) - ncol(b))
        zz <- tcrossprod(z)/(n - p)
        tem1 = (n - p)^(-1) * solve(zz) %x% sigu
        tem2 = sqrt(diag(tem1))
        tem3 = matrix(tem2, nrow = k, ncol = ncol(b))
        tmat = b/tem3
        return(list(coef = b, resid = t(e), sigu = sigu, zzmat = zz,
                    zmat = z, tratio = tmat, p=p))
    }
}

VAR.Boot2 <- function(x, p, nb = 500, type = "const", seed=1) {
    ## reverse engineer VAR.Boot
    set.seed(seed)
    require(VAR.etp)
    n <- nrow(x)
    k <- ncol(x)
    var1 <- VAR.est(x, p, type)
    b <- var1$coef  # OLS estimates
    mat <- matrix(0, nrow = k, ncol = ncol(b))
    for (i in 1:nb) {
        ## es <- VAR.etp:::resamp(e)
        es <- var1$resid[sample(n-p, n, replace=TRUE),,drop=FALSE]
        ## xs <- VAR.etp:::VAR.ys(x, b, p, es, type)
        xs <- VAR.ys2(x, b, p, es, type)
        bs <- VAR.est2(xs, p, type, coef.only=TRUE)$coef
        mat <- mat + bs/nb
    }
    bias <- mat - b
    bs <- VAR.adjust2(b, bias, p, type)
    colnames(bs) <- VAR.etp:::VAR.names(x, p, type)
    es <- VAR.etp:::VAR.resid(x, bs, var1$zmat, p)
    colnames(es) <- rownames(b)
    sigu <- cov(es)  # t(es) %*% es/((n - p) - ncol(b))
    return(list(coef = bs, resid = es, sigu = sigu, Bias = bias, p=p))
}

