R/impulseResponse.R
In svazzole/sparseVAR: Sparse VAR/VECM Models Estimation
Defines functions
getIRF
bootstrapOLS
jackknife
errorBandsIRF
checkImpulseZero
impulseResponse
Documented in
checkImpulseZero
errorBandsIRF
impulseResponse
#' @title Impulse Response Function
#'
#' @description A function to estimate the Impulse Response Function of a given VAR.
#'
#' @usage impulseResponse(v, len = 20)
#'
#' @param v the data in the for of a VAR
#' @param len length of the impulse response function
#'
#' @return \code{irf} a 3d array containing the impulse response function.
#'
#' @export
impulseResponse <- function(v, len = 20) {
# Check if v is a VAR object
if (!checkIsVar(v)) {
stop("Input v must be a VAR object")
}
# Numerical problems in the estimated variance covariance
e <- eigen(v$sigma)$values
if (!is.null(e[e <= 0])) {
P <- t(chol(v$sigma, pivot = TRUE))
} else {
P <- t(chol(v$sigma))
}
bigA <- companionVAR(v)
out <- getIRF(v, bigA, len = len, P)
out$cholP <- P # Add Choleski factorization to the output
return(out)
}
#' @title Check Impulse Zero
#'
#' @description A function to find which entries of the impulse response function
#' are zero.
#'
#' @usage checkImpulseZero(irf)
#'
#' @param irf irf output from impulseResponse function
#'
#' @return a matrix containing the indices of the impulse response function that
#' are 0.
#'
#' @export
checkImpulseZero <- function(irf) {
nx <- dim(irf)[1]
ny <- dim(irf)[2]
nz <- dim(irf)[3]
logicalIrf <- matrix(0, nx, ny)
for (z in 1:nz) {
logicalIrf <- logicalIrf + abs(irf[, , z])
}
logicalIrf <- logicalIrf == 0
return(which(logicalIrf == TRUE, arr.ind = TRUE))
}
#' @title Error bands for IRF
#'
#' @description A function to estimate the confidence intervals for irf and oirf.
#'
#' @usage errorBandsIRF(v, irf, alpha, M, resampling, ...)
#'
#' @param v a var object as from fitVAR or simulateVAR
#' @param irf irf output from impulseResponse function
#' @param alpha level of confidence (default \code{alpha = 0.01})
#' @param M number of bootstrapped series (default \code{M = 100})
#' @param resampling type of resampling: \code{"bootstrap"} or \code{"jackknife"}
#' @param ... some options for the estimation: \code{verbose = TRUE} or \code{FALSE},
#' \code{mode = "fast"} or \code{"slow"}, \code{threshold = TRUE} or \code{FALSE}.
#'
#' @return a matrix containing the indices of the impulse response function that
#' are 0.
#'
#' @export
errorBandsIRF <- function(v, irf, alpha = 0.01, M = 100, resampling = "bootstrap", ...) {
opt <- list(...)
verbose <- ifelse(!is.null(opt$verbose), opt$verbose, TRUE)
mode <- ifelse(!is.null(opt$mode), opt$mode, "fast")
threshold <- ifelse(!is.null(opt$threshold), opt$threshold, FALSE)
thresholdType <- ifelse(!is.null(opt$thresholdType), opt$thresholdType, "soft")
if (resampling == "bootstrap") {
lambda <- v$lambda
p <- length(v$A)
nr <- ncol(v$series)
nc <- ncol(v$A[[1]])
len <- dim(irf$irf)[3]
irfs <- array(data = rep(0, len * nc^2 * M), dim = c(nc, nc, len + 1, M))
oirfs <- array(data = rep(0, len * nc^2 * M), dim = c(nc, nc, len + 1, M))
if (verbose == TRUE) {
cat("Step 1 of 2: bootstrapping series and re-estimating VAR...\n")
pb <- utils::txtProgressBar(min = 0, max = M, style = 3)
}
for (k in 1:M) {
# create Xs and Ys (temp variables)
o <- bootstrappedVAR(v)
if (mode == "fast") {
if (v$penalty == "ENET") {
# fit ENET to a specific value of lambda
fit <- varENET(o, p, lambda, opt = list(method = v$method, penalty = v$penalty))
Avector <- stats::coef(fit, s = lambda)
A <- matrix(Avector[2:length(Avector)], nrow = nc, ncol = nc * p, byrow = TRUE)
} else if (v$penalty == "SCAD") {
fit <- varSCAD(o, p, lambda, opt = list(method = v$method, penalty = v$penalty))
Avector <- fit$beta[2:nrow(fit$beta), 1]
A <- matrix(Avector, nrow = nc, ncol = nc * p, byrow = TRUE)
} else {
fit <- varMCP(o, p, lambda, opt = list(method = v$method, penalty = v$penalty))
Avector <- fit$beta[2:nrow(fit$beta), 1]
A <- matrix(Avector, nrow = nc, ncol = nc * p, byrow = TRUE)
}
if (threshold == TRUE) {
applyThreshold(A, nr, nc, p, type = thresholdType)
}
M <- cbind(diag(x = 1, nrow = (nc * (p - 1)), ncol = (nc * (p - 1))), matrix(0, nrow = (nc * (p - 1)), ncol = nc))
bigA <- rbind(A, M)
} else {
# fit ENET on a series of lambdas
if (threshold == TRUE) {
fit <- fitVAR(o, p, penalty = v$penalty, method = v$method, threshold = TRUE)
} else {
fit <- fitVAR(o, p, penalty = v$penalty, method = v$method)
}
bigA <- companionVAR(fit)
}
tmpRes <- getIRF(v, bigA, len, irf$cholP)
irfs[, , , k] <- tmpRes$irf
oirfs[, , , k] <- tmpRes$oirf
if (verbose == TRUE) {
utils::setTxtProgressBar(pb, k)
}
}
if (verbose == TRUE) {
close(pb)
cat("Step 2 of 2: computing quantiles...\n")
pb <- utils::txtProgressBar(min = 0, max = (nc * nc), style = 3)
}
irfUB <- array(data = rep(0, len * nc^2), dim = c(nc, nc, len))
irfLB <- array(data = rep(0, len * nc^2), dim = c(nc, nc, len))
oirfUB <- array(data = rep(0, len * nc^2), dim = c(nc, nc, len))
oirfLB <- array(data = rep(0, len * nc^2), dim = c(nc, nc, len))
irfQUB <- array(data = rep(0, len * nc^2), dim = c(nc, nc, len))
irfQLB <- irfQUB
oirfQUB <- irfQUB
oirfQLB <- irfQUB
a <- alpha / 2
qLB <- stats::qnorm(a)
qUB <- stats::qnorm((1 - a))
for (i in 1:nc) {
for (j in 1:nc) {
for (k in 1:len) {
irfQUB[i, j, k] <- stats::quantile(irfs[i, j, k, ], probs = (1 - a), na.rm = TRUE)
oirfQUB[i, j, k] <- stats::quantile(oirfs[i, j, k, ], probs = (1 - a), na.rm = TRUE)
irfQLB[i, j, k] <- stats::quantile(irfs[i, j, k, ], probs = a, na.rm = TRUE)
oirfQLB[i, j, k] <- stats::quantile(oirfs[i, j, k, ], probs = a, na.rm = TRUE)
irfUB[i, j, k] <- qUB * stats::sd(irfs[i, j, k, ])
oirfUB[i, j, k] <- qUB * stats::sd(oirfs[i, j, k, ])
irfLB[i, j, k] <- qLB * stats::sd(irfs[i, j, k, ])
oirfLB[i, j, k] <- qLB * stats::sd(oirfs[i, j, k, ])
}
if (verbose == TRUE) {
utils::setTxtProgressBar(pb, (i - 1) * nc + j)
}
}
}
if (verbose == TRUE) {
close(pb)
}
output <- list()
output$irfUB <- irfUB
output$oirfUB <- oirfUB
output$irfLB <- irfLB
output$oirfLB <- oirfLB
output$irfQUB <- irfQUB
output$oirfQUB <- oirfQUB
output$irfQLB <- irfQLB
output$oirfQLB <- oirfQLB
attr(output, "class") <- "irfBands"
attr(output, "resampling") <- "bootstrap"
return(output)
} else if (resampling == "jackknife") {
output <- jackknife(v, irf, alpha = alpha, ...)
return(output)
} else if (resampling == "bootstrapOLS") {
output <- bootstrapOLS(v, irf, alpha = alpha, ...)
} else {
stop("Unknown resampling method. Possible values are \"bootstrap\" or \"jackknife\"")
}
}
jackknife <- function(v, irf, mode = "fast", alpha, ...) {
lambda <- v$lambda
p <- length(v$A)
nc <- ncol(v$A[[1]])
len <- dim(irf$irf)[3]
nr <- nrow(v$series)
opt <- list(...)
verbose <- ifelse(!is.null(opt$verbose), opt$verbose, TRUE)
mode <- ifelse(!is.null(opt$mode), opt$mode, "fast")
threshold <- ifelse(!is.null(opt$threshold), opt$threshold, FALSE)
thresholdType <- ifelse(!is.null(opt$thresholdType), opt$thresholdType, "soft")
irfs <- array(data = rep(0, len * nc^2 * nr), dim = c(nc, nc, len + 1, nr))
oirfs <- array(data = rep(0, len * nc^2 * nr), dim = c(nc, nc, len + 1, nr))
if (verbose == TRUE) {
cat("Step 1 of 2: jack knifing series and re-estimating VAR...\n")
pb <- utils::txtProgressBar(min = 0, max = nr, style = 3)
}
for (k in 1:nr) {
# create Xs and Ys (temp variables)
data <- v$series[-k, ]
trDt <- transformData(data, p, opt = list(method = v$method, penalty = v$penalty))
trDt$X <- trDt$X
trDt$y <- trDt$y
# data <- v$series
# trDt <- transformData(data, p, opt = list(method = v$method, penalty = v$penalty))
# trDt$X <- trDt$X[-k, ]
# trDt$y <- trDt$y[-k]
if (mode == "fast") {
if (v$penalty == "ENET") {
# fit ENET to a specific value of lambda
fit <- varENET(data, p, lambda, opt = list(method = v$method, penalty = v$penalty))
Avector <- stats::coef(fit, s = lambda)
A <- matrix(Avector[2:length(Avector)], nrow = nc, ncol = nc * p, byrow = TRUE)
} else if (v$penalty == "SCAD") {
fit <- varSCAD(data, p, lambda, opt = list(method = v$method, penalty = v$penalty))
Avector <- fit$beta[2:nrow(fit$beta), 1]
A <- matrix(Avector, nrow = nc, ncol = nc * p, byrow = TRUE)
} else {
fit <- varMCP(data, p, lambda, opt = list(method = v$method, penalty = v$penalty))
Avector <- fit$beta[2:nrow(fit$beta), 1]
A <- matrix(Avector, nrow = nc, ncol = nc * p, byrow = TRUE)
}
if (threshold == TRUE) {
applyThreshold(A, nr, nc, p, type = thresholdType)
}
M <- cbind(diag(x = 1, nrow = (nc * (p - 1)), ncol = (nc * (p - 1))), matrix(0, nrow = (nc * (p - 1)), ncol = nc))
bigA <- rbind(A, M)
} else {
# fit ENET on a series of lambdas
if (threshold == TRUE) {
fit <- fitVAR(data, p, penalty = v$penalty, method = v$method, threshold = TRUE)
} else {
fit <- fitVAR(data, p, penalty = v$penalty, method = v$method)
}
bigA <- companionVAR(fit)
}
tmpRes <- getIRF(v, bigA, len, irf$cholP)
irfs[, , , k] <- tmpRes$irf
oirfs[, , , k] <- tmpRes$oirf
if (verbose == TRUE) {
utils::setTxtProgressBar(pb, k)
}
}
if (verbose == TRUE) {
close(pb)
cat("Step 2 of 2: computing quantiles...\n")
pb <- utils::txtProgressBar(min = 0, max = (nc * nc), style = 3)
}
irfUB <- array(data = rep(0, len * nc^2), dim = c(nc, nc, len))
irfLB <- array(data = rep(0, len * nc^2), dim = c(nc, nc, len))
oirfUB <- array(data = rep(0, len * nc^2), dim = c(nc, nc, len))
oirfLB <- array(data = rep(0, len * nc^2), dim = c(nc, nc, len))
a <- alpha / 2
qLB <- stats::qnorm(a)
qUB <- stats::qnorm((1 - a))
for (i in 1:nc) {
for (j in 1:nc) {
for (k in 1:len) {
irfUB[i, j, k] <- base::mean(irfs[i, j, k, ]) + qUB * stats::sd(irfs[i, j, k, ])
oirfUB[i, j, k] <- base::mean(oirfs[i, j, k, ]) + qUB * stats::sd(oirfs[i, j, k, ])
irfLB[i, j, k] <- base::mean(irfs[i, j, k, ]) + qLB * stats::sd(irfs[i, j, k, ])
oirfLB[i, j, k] <- base::mean(oirfs[i, j, k, ]) + qLB * stats::sd(oirfs[i, j, k, ])
}
if (verbose == TRUE) {
utils::setTxtProgressBar(pb, (i - 1) * nc + j)
}
}
}
if (verbose == TRUE) {
close(pb)
}
output <- list()
output$irfUB <- irfUB
output$oirfUB <- oirfUB
output$irfLB <- irfLB
output$oirfLB <- oirfLB
attr(output, "class") <- "irfBands"
attr(output, "resampling") <- "jackknife"
return(output)
}
bootstrapOLS <- function(v, irf, mode = "fast", alpha, ...) {
lambda <- v$lambda
p <- length(v$A)
nc <- ncol(v$A[[1]])
len <- dim(irf$irf)[3]
nr <- nrow(v$series)
opt <- list(...)
verbose <- ifelse(!is.null(opt$verbose), opt$verbose, TRUE)
mode <- ifelse(!is.null(opt$mode), opt$mode, "fast")
threshold <- ifelse(!is.null(opt$threshold), opt$threshold, FALSE)
thresholdType <- ifelse(!is.null(opt$thresholdType), opt$thresholdType, "soft")
irfs <- array(data = rep(0, len * nc^2 * nr), dim = c(nc, nc, len + 1, nr))
oirfs <- array(data = rep(0, len * nc^2 * nr), dim = c(nc, nc, len + 1, nr))
if (verbose == TRUE) {
cat("Step 1 of 2: bootstrappingOLS series and re-estimating VAR...\n")
pb <- utils::txtProgressBar(min = 0, max = nr, style = 3)
}
for (k in 1:nr) {
# create Xs and Ys (temp variables)
o <- bootstrappedVAR(v)
N <- ncol(v$A[[1]])
nobs <- nrow(v$series)
bigA <- companionVAR(v)
trDt <- transformData(o, p = p, opt = list(method = v$method, scale = FALSE, center = TRUE))
nonZeroEntries <- as.matrix(bigA != 0)
## Create matrix R
t <- as.vector(nonZeroEntries)
n <- sum(t != 0)
ix <- which(t != 0)
j <- 1:n
R <- matrix(0, ncol = n, nrow = length(t))
for (zz in 1:n) {
R[ix[zz], j[zz]] <- 1
}
X <- as.matrix(trDt$X)
y <- as.vector(t(o[-(1:p), ]))
# Metodo A MANO
s <- corpcor::invcov.shrink(v$residuals, verbose = FALSE)
G <- t(o[-nobs, ]) %*% o[-nobs, ] / nobs
V <- solve(t(R) %*% (kronecker(G, s) %*% R))
VV <- nonZeroEntries
VV[nonZeroEntries] <- diag(V)
G1 <- solve(t(R) %*% (kronecker(t(o[-nobs, ]) %*% o[-nobs, ], s)) %*% R)
G2 <- t(R) %*% (kronecker(t(o[-nobs, ]), s))
g <- G1 %*% G2 # [ , (N+1):(length(y) + N)]
ga <- g %*% y
b1 <- vector(length = N * N)
b1 <- R %*% ga
A <- matrix(b1, ncol = N, byrow = F)
M <- cbind(diag(x = 1, nrow = (nc * (p - 1)), ncol = (nc * (p - 1))), matrix(0, nrow = (nc * (p - 1)), ncol = nc))
bigA <- rbind(A, M)
tmpRes <- getIRF(v, bigA, len, irf$cholP)
irfs[, , , k] <- tmpRes$irf
oirfs[, , , k] <- tmpRes$oirf
if (verbose == TRUE) {
utils::setTxtProgressBar(pb, k)
}
}
if (verbose == TRUE) {
close(pb)
cat("Step 2 of 2: computing quantiles...\n")
pb <- utils::txtProgressBar(min = 0, max = (nc * nc), style = 3)
}
irfUB <- array(data = rep(0, len * nc^2), dim = c(nc, nc, len))
irfLB <- array(data = rep(0, len * nc^2), dim = c(nc, nc, len))
oirfUB <- array(data = rep(0, len * nc^2), dim = c(nc, nc, len))
oirfLB <- array(data = rep(0, len * nc^2), dim = c(nc, nc, len))
a <- alpha / 2
qLB <- stats::qnorm(a)
qUB <- stats::qnorm((1 - a))
for (i in 1:nc) {
for (j in 1:nc) {
for (k in 1:len) {
irfUB[i, j, k] <- base::mean(irfs[i, j, k, ]) + qUB * stats::sd(irfs[i, j, k, ])
oirfUB[i, j, k] <- base::mean(oirfs[i, j, k, ]) + qUB * stats::sd(oirfs[i, j, k, ])
irfLB[i, j, k] <- base::mean(irfs[i, j, k, ]) + qLB * stats::sd(irfs[i, j, k, ])
oirfLB[i, j, k] <- base::mean(oirfs[i, j, k, ]) + qLB * stats::sd(oirfs[i, j, k, ])
}
if (verbose == TRUE) {
utils::setTxtProgressBar(pb, (i - 1) * nc + j)
}
}
}
if (verbose == TRUE) {
close(pb)
}
output <- list()
output$irfUB <- irfUB
output$oirfUB <- oirfUB
output$irfLB <- irfLB
output$oirfLB <- oirfLB
attr(output, "class") <- "irfBands"
attr(output, "resampling") <- "jackknife"
return(output)
}
getIRF <- function(v, bigA, len = 20, P) {
nr <- nrow(v$A[[1]])
irf <- array(data = rep(0, len * nr^2), dim = c(nr, nr, len + 1))
oirf <- array(data = rep(0, len * nr^2), dim = c(nr, nr, len + 1))
Atmp <- diag(nrow = nrow(bigA), ncol = ncol(bigA))
irf[, , 1] <- Atmp[1:nr, 1:nr]
oirf[, , 1] <- Atmp[1:nr, 1:nr] %*% P
for (k in 1:len) {
Atmp <- Atmp %*% bigA
irf[, , (k + 1)] <- as.matrix(Atmp[1:nr, 1:nr])
oirf[, , (k + 1)] <- as.matrix(Atmp[1:nr, 1:nr] %*% P)
}
## TODO: add cumulative response functions
out <- list()
out$irf <- irf
out$oirf <- oirf
attr(out, "class") <- "irf"
return(out)
}
svazzole/sparseVAR documentation built on April 19, 2021, 2:11 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions getIRF bootstrapOLS jackknife errorBandsIRF checkImpulseZero impulseResponse
Documented in checkImpulseZero errorBandsIRF impulseResponse
#' @title Impulse Response Function
#'
#' @description A function to estimate the Impulse Response Function of a given VAR.
#'
#' @usage impulseResponse(v, len = 20)
#'
#' @param v the data in the for of a VAR
#' @param len length of the impulse response function
#'
#' @return \code{irf} a 3d array containing the impulse response function.
#'
#' @export
impulseResponse <- function(v, len = 20) {
# Check if v is a VAR object
if (!checkIsVar(v)) {
stop("Input v must be a VAR object")
}
# Numerical problems in the estimated variance covariance
e <- eigen(v$sigma)$values
if (!is.null(e[e <= 0])) {
P <- t(chol(v$sigma, pivot = TRUE))
} else {
P <- t(chol(v$sigma))
}
bigA <- companionVAR(v)
out <- getIRF(v, bigA, len = len, P)
out$cholP <- P # Add Choleski factorization to the output
return(out)
}
#' @title Check Impulse Zero
#'
#' @description A function to find which entries of the impulse response function
#' are zero.
#'
#' @usage checkImpulseZero(irf)
#'
#' @param irf irf output from impulseResponse function
#'
#' @return a matrix containing the indices of the impulse response function that
#' are 0.
#'
#' @export
checkImpulseZero <- function(irf) {
nx <- dim(irf)[1]
ny <- dim(irf)[2]
nz <- dim(irf)[3]
logicalIrf <- matrix(0, nx, ny)
for (z in 1:nz) {
logicalIrf <- logicalIrf + abs(irf[, , z])
}
logicalIrf <- logicalIrf == 0
return(which(logicalIrf == TRUE, arr.ind = TRUE))
}
#' @title Error bands for IRF
#'
#' @description A function to estimate the confidence intervals for irf and oirf.
#'
#' @usage errorBandsIRF(v, irf, alpha, M, resampling, ...)
#'
#' @param v a var object as from fitVAR or simulateVAR
#' @param irf irf output from impulseResponse function
#' @param alpha level of confidence (default \code{alpha = 0.01})
#' @param M number of bootstrapped series (default \code{M = 100})
#' @param resampling type of resampling: \code{"bootstrap"} or \code{"jackknife"}
#' @param ... some options for the estimation: \code{verbose = TRUE} or \code{FALSE},
#' \code{mode = "fast"} or \code{"slow"}, \code{threshold = TRUE} or \code{FALSE}.
#'
#' @return a matrix containing the indices of the impulse response function that
#' are 0.
#'
#' @export
errorBandsIRF <- function(v, irf, alpha = 0.01, M = 100, resampling = "bootstrap", ...) {
opt <- list(...)
verbose <- ifelse(!is.null(opt$verbose), opt$verbose, TRUE)
mode <- ifelse(!is.null(opt$mode), opt$mode, "fast")
threshold <- ifelse(!is.null(opt$threshold), opt$threshold, FALSE)
thresholdType <- ifelse(!is.null(opt$thresholdType), opt$thresholdType, "soft")
if (resampling == "bootstrap") {
lambda <- v$lambda
p <- length(v$A)
nr <- ncol(v$series)
nc <- ncol(v$A[[1]])
len <- dim(irf$irf)[3]
irfs <- array(data = rep(0, len * nc^2 * M), dim = c(nc, nc, len + 1, M))
oirfs <- array(data = rep(0, len * nc^2 * M), dim = c(nc, nc, len + 1, M))
if (verbose == TRUE) {
cat("Step 1 of 2: bootstrapping series and re-estimating VAR...\n")
pb <- utils::txtProgressBar(min = 0, max = M, style = 3)
}
for (k in 1:M) {
# create Xs and Ys (temp variables)
o <- bootstrappedVAR(v)
if (mode == "fast") {
if (v$penalty == "ENET") {
# fit ENET to a specific value of lambda
fit <- varENET(o, p, lambda, opt = list(method = v$method, penalty = v$penalty))
Avector <- stats::coef(fit, s = lambda)
A <- matrix(Avector[2:length(Avector)], nrow = nc, ncol = nc * p, byrow = TRUE)
} else if (v$penalty == "SCAD") {
fit <- varSCAD(o, p, lambda, opt = list(method = v$method, penalty = v$penalty))
Avector <- fit$beta[2:nrow(fit$beta), 1]
A <- matrix(Avector, nrow = nc, ncol = nc * p, byrow = TRUE)
} else {
fit <- varMCP(o, p, lambda, opt = list(method = v$method, penalty = v$penalty))
Avector <- fit$beta[2:nrow(fit$beta), 1]
A <- matrix(Avector, nrow = nc, ncol = nc * p, byrow = TRUE)
}
if (threshold == TRUE) {
applyThreshold(A, nr, nc, p, type = thresholdType)
}
M <- cbind(diag(x = 1, nrow = (nc * (p - 1)), ncol = (nc * (p - 1))), matrix(0, nrow = (nc * (p - 1)), ncol = nc))
bigA <- rbind(A, M)
} else {
# fit ENET on a series of lambdas
if (threshold == TRUE) {
fit <- fitVAR(o, p, penalty = v$penalty, method = v$method, threshold = TRUE)
} else {
fit <- fitVAR(o, p, penalty = v$penalty, method = v$method)
}
bigA <- companionVAR(fit)
}
tmpRes <- getIRF(v, bigA, len, irf$cholP)
irfs[, , , k] <- tmpRes$irf
oirfs[, , , k] <- tmpRes$oirf
if (verbose == TRUE) {
utils::setTxtProgressBar(pb, k)
}
}
if (verbose == TRUE) {
close(pb)
cat("Step 2 of 2: computing quantiles...\n")
pb <- utils::txtProgressBar(min = 0, max = (nc * nc), style = 3)
}
irfUB <- array(data = rep(0, len * nc^2), dim = c(nc, nc, len))
irfLB <- array(data = rep(0, len * nc^2), dim = c(nc, nc, len))
oirfUB <- array(data = rep(0, len * nc^2), dim = c(nc, nc, len))
oirfLB <- array(data = rep(0, len * nc^2), dim = c(nc, nc, len))
irfQUB <- array(data = rep(0, len * nc^2), dim = c(nc, nc, len))
irfQLB <- irfQUB
oirfQUB <- irfQUB
oirfQLB <- irfQUB
a <- alpha / 2
qLB <- stats::qnorm(a)
qUB <- stats::qnorm((1 - a))
for (i in 1:nc) {
for (j in 1:nc) {
for (k in 1:len) {
irfQUB[i, j, k] <- stats::quantile(irfs[i, j, k, ], probs = (1 - a), na.rm = TRUE)
oirfQUB[i, j, k] <- stats::quantile(oirfs[i, j, k, ], probs = (1 - a), na.rm = TRUE)
irfQLB[i, j, k] <- stats::quantile(irfs[i, j, k, ], probs = a, na.rm = TRUE)
oirfQLB[i, j, k] <- stats::quantile(oirfs[i, j, k, ], probs = a, na.rm = TRUE)
irfUB[i, j, k] <- qUB * stats::sd(irfs[i, j, k, ])
oirfUB[i, j, k] <- qUB * stats::sd(oirfs[i, j, k, ])
irfLB[i, j, k] <- qLB * stats::sd(irfs[i, j, k, ])
oirfLB[i, j, k] <- qLB * stats::sd(oirfs[i, j, k, ])
}
if (verbose == TRUE) {
utils::setTxtProgressBar(pb, (i - 1) * nc + j)
}
}
}
if (verbose == TRUE) {
close(pb)
}
output <- list()
output$irfUB <- irfUB
output$oirfUB <- oirfUB
output$irfLB <- irfLB
output$oirfLB <- oirfLB
output$irfQUB <- irfQUB
output$oirfQUB <- oirfQUB
output$irfQLB <- irfQLB
output$oirfQLB <- oirfQLB
attr(output, "class") <- "irfBands"
attr(output, "resampling") <- "bootstrap"
return(output)
} else if (resampling == "jackknife") {
output <- jackknife(v, irf, alpha = alpha, ...)
return(output)
} else if (resampling == "bootstrapOLS") {
output <- bootstrapOLS(v, irf, alpha = alpha, ...)
} else {
stop("Unknown resampling method. Possible values are \"bootstrap\" or \"jackknife\"")
}
}
jackknife <- function(v, irf, mode = "fast", alpha, ...) {
lambda <- v$lambda
p <- length(v$A)
nc <- ncol(v$A[[1]])
len <- dim(irf$irf)[3]
nr <- nrow(v$series)
opt <- list(...)
verbose <- ifelse(!is.null(opt$verbose), opt$verbose, TRUE)
mode <- ifelse(!is.null(opt$mode), opt$mode, "fast")
threshold <- ifelse(!is.null(opt$threshold), opt$threshold, FALSE)
thresholdType <- ifelse(!is.null(opt$thresholdType), opt$thresholdType, "soft")
irfs <- array(data = rep(0, len * nc^2 * nr), dim = c(nc, nc, len + 1, nr))
oirfs <- array(data = rep(0, len * nc^2 * nr), dim = c(nc, nc, len + 1, nr))
if (verbose == TRUE) {
cat("Step 1 of 2: jack knifing series and re-estimating VAR...\n")
pb <- utils::txtProgressBar(min = 0, max = nr, style = 3)
}
for (k in 1:nr) {
# create Xs and Ys (temp variables)
data <- v$series[-k, ]
trDt <- transformData(data, p, opt = list(method = v$method, penalty = v$penalty))
trDt$X <- trDt$X
trDt$y <- trDt$y
# data <- v$series
# trDt <- transformData(data, p, opt = list(method = v$method, penalty = v$penalty))
# trDt$X <- trDt$X[-k, ]
# trDt$y <- trDt$y[-k]
if (mode == "fast") {
if (v$penalty == "ENET") {
# fit ENET to a specific value of lambda
fit <- varENET(data, p, lambda, opt = list(method = v$method, penalty = v$penalty))
Avector <- stats::coef(fit, s = lambda)
A <- matrix(Avector[2:length(Avector)], nrow = nc, ncol = nc * p, byrow = TRUE)
} else if (v$penalty == "SCAD") {
fit <- varSCAD(data, p, lambda, opt = list(method = v$method, penalty = v$penalty))
Avector <- fit$beta[2:nrow(fit$beta), 1]
A <- matrix(Avector, nrow = nc, ncol = nc * p, byrow = TRUE)
} else {
fit <- varMCP(data, p, lambda, opt = list(method = v$method, penalty = v$penalty))
Avector <- fit$beta[2:nrow(fit$beta), 1]
A <- matrix(Avector, nrow = nc, ncol = nc * p, byrow = TRUE)
}
if (threshold == TRUE) {
applyThreshold(A, nr, nc, p, type = thresholdType)
}
M <- cbind(diag(x = 1, nrow = (nc * (p - 1)), ncol = (nc * (p - 1))), matrix(0, nrow = (nc * (p - 1)), ncol = nc))
bigA <- rbind(A, M)
} else {
# fit ENET on a series of lambdas
if (threshold == TRUE) {
fit <- fitVAR(data, p, penalty = v$penalty, method = v$method, threshold = TRUE)
} else {
fit <- fitVAR(data, p, penalty = v$penalty, method = v$method)
}
bigA <- companionVAR(fit)
}
tmpRes <- getIRF(v, bigA, len, irf$cholP)
irfs[, , , k] <- tmpRes$irf
oirfs[, , , k] <- tmpRes$oirf
if (verbose == TRUE) {
utils::setTxtProgressBar(pb, k)
}
}
if (verbose == TRUE) {
close(pb)
cat("Step 2 of 2: computing quantiles...\n")
pb <- utils::txtProgressBar(min = 0, max = (nc * nc), style = 3)
}
irfUB <- array(data = rep(0, len * nc^2), dim = c(nc, nc, len))
irfLB <- array(data = rep(0, len * nc^2), dim = c(nc, nc, len))
oirfUB <- array(data = rep(0, len * nc^2), dim = c(nc, nc, len))
oirfLB <- array(data = rep(0, len * nc^2), dim = c(nc, nc, len))
a <- alpha / 2
qLB <- stats::qnorm(a)
qUB <- stats::qnorm((1 - a))
for (i in 1:nc) {
for (j in 1:nc) {
for (k in 1:len) {
irfUB[i, j, k] <- base::mean(irfs[i, j, k, ]) + qUB * stats::sd(irfs[i, j, k, ])
oirfUB[i, j, k] <- base::mean(oirfs[i, j, k, ]) + qUB * stats::sd(oirfs[i, j, k, ])
irfLB[i, j, k] <- base::mean(irfs[i, j, k, ]) + qLB * stats::sd(irfs[i, j, k, ])
oirfLB[i, j, k] <- base::mean(oirfs[i, j, k, ]) + qLB * stats::sd(oirfs[i, j, k, ])
}
if (verbose == TRUE) {
utils::setTxtProgressBar(pb, (i - 1) * nc + j)
}
}
}
if (verbose == TRUE) {
close(pb)
}
output <- list()
output$irfUB <- irfUB
output$oirfUB <- oirfUB
output$irfLB <- irfLB
output$oirfLB <- oirfLB
attr(output, "class") <- "irfBands"
attr(output, "resampling") <- "jackknife"
return(output)
}
bootstrapOLS <- function(v, irf, mode = "fast", alpha, ...) {
lambda <- v$lambda
p <- length(v$A)
nc <- ncol(v$A[[1]])
len <- dim(irf$irf)[3]
nr <- nrow(v$series)
opt <- list(...)
verbose <- ifelse(!is.null(opt$verbose), opt$verbose, TRUE)
mode <- ifelse(!is.null(opt$mode), opt$mode, "fast")
threshold <- ifelse(!is.null(opt$threshold), opt$threshold, FALSE)
thresholdType <- ifelse(!is.null(opt$thresholdType), opt$thresholdType, "soft")
irfs <- array(data = rep(0, len * nc^2 * nr), dim = c(nc, nc, len + 1, nr))
oirfs <- array(data = rep(0, len * nc^2 * nr), dim = c(nc, nc, len + 1, nr))
if (verbose == TRUE) {
cat("Step 1 of 2: bootstrappingOLS series and re-estimating VAR...\n")
pb <- utils::txtProgressBar(min = 0, max = nr, style = 3)
}
for (k in 1:nr) {
# create Xs and Ys (temp variables)
o <- bootstrappedVAR(v)
N <- ncol(v$A[[1]])
nobs <- nrow(v$series)
bigA <- companionVAR(v)
trDt <- transformData(o, p = p, opt = list(method = v$method, scale = FALSE, center = TRUE))
nonZeroEntries <- as.matrix(bigA != 0)
## Create matrix R
t <- as.vector(nonZeroEntries)
n <- sum(t != 0)
ix <- which(t != 0)
j <- 1:n
R <- matrix(0, ncol = n, nrow = length(t))
for (zz in 1:n) {
R[ix[zz], j[zz]] <- 1
}
X <- as.matrix(trDt$X)
y <- as.vector(t(o[-(1:p), ]))
# Metodo A MANO
s <- corpcor::invcov.shrink(v$residuals, verbose = FALSE)
G <- t(o[-nobs, ]) %*% o[-nobs, ] / nobs
V <- solve(t(R) %*% (kronecker(G, s) %*% R))
VV <- nonZeroEntries
VV[nonZeroEntries] <- diag(V)
G1 <- solve(t(R) %*% (kronecker(t(o[-nobs, ]) %*% o[-nobs, ], s)) %*% R)
G2 <- t(R) %*% (kronecker(t(o[-nobs, ]), s))
g <- G1 %*% G2 # [ , (N+1):(length(y) + N)]
ga <- g %*% y
b1 <- vector(length = N * N)
b1 <- R %*% ga
A <- matrix(b1, ncol = N, byrow = F)
M <- cbind(diag(x = 1, nrow = (nc * (p - 1)), ncol = (nc * (p - 1))), matrix(0, nrow = (nc * (p - 1)), ncol = nc))
bigA <- rbind(A, M)
tmpRes <- getIRF(v, bigA, len, irf$cholP)
irfs[, , , k] <- tmpRes$irf
oirfs[, , , k] <- tmpRes$oirf
if (verbose == TRUE) {
utils::setTxtProgressBar(pb, k)
}
}
if (verbose == TRUE) {
close(pb)
cat("Step 2 of 2: computing quantiles...\n")
pb <- utils::txtProgressBar(min = 0, max = (nc * nc), style = 3)
}
irfUB <- array(data = rep(0, len * nc^2), dim = c(nc, nc, len))
irfLB <- array(data = rep(0, len * nc^2), dim = c(nc, nc, len))
oirfUB <- array(data = rep(0, len * nc^2), dim = c(nc, nc, len))
oirfLB <- array(data = rep(0, len * nc^2), dim = c(nc, nc, len))
a <- alpha / 2
qLB <- stats::qnorm(a)
qUB <- stats::qnorm((1 - a))
for (i in 1:nc) {
for (j in 1:nc) {
for (k in 1:len) {
irfUB[i, j, k] <- base::mean(irfs[i, j, k, ]) + qUB * stats::sd(irfs[i, j, k, ])
oirfUB[i, j, k] <- base::mean(oirfs[i, j, k, ]) + qUB * stats::sd(oirfs[i, j, k, ])
irfLB[i, j, k] <- base::mean(irfs[i, j, k, ]) + qLB * stats::sd(irfs[i, j, k, ])
oirfLB[i, j, k] <- base::mean(oirfs[i, j, k, ]) + qLB * stats::sd(oirfs[i, j, k, ])
}
if (verbose == TRUE) {
utils::setTxtProgressBar(pb, (i - 1) * nc + j)
}
}
}
if (verbose == TRUE) {
close(pb)
}
output <- list()
output$irfUB <- irfUB
output$oirfUB <- oirfUB
output$irfLB <- irfLB
output$oirfLB <- oirfLB
attr(output, "class") <- "irfBands"
attr(output, "resampling") <- "jackknife"
return(output)
}
getIRF <- function(v, bigA, len = 20, P) {
nr <- nrow(v$A[[1]])
irf <- array(data = rep(0, len * nr^2), dim = c(nr, nr, len + 1))
oirf <- array(data = rep(0, len * nr^2), dim = c(nr, nr, len + 1))
Atmp <- diag(nrow = nrow(bigA), ncol = ncol(bigA))
irf[, , 1] <- Atmp[1:nr, 1:nr]
oirf[, , 1] <- Atmp[1:nr, 1:nr] %*% P
for (k in 1:len) {
Atmp <- Atmp %*% bigA
irf[, , (k + 1)] <- as.matrix(Atmp[1:nr, 1:nr])
oirf[, , (k + 1)] <- as.matrix(Atmp[1:nr, 1:nr] %*% P)
}
## TODO: add cumulative response functions
out <- list()
out$irf <- irf
out$oirf <- oirf
attr(out, "class") <- "irf"
return(out)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.