R/irf.r
In joergrieger/bvar: Estimation and forecasting of bayesian VAR models
Defines functions
irf.tvar
irf.msvar
irf.ftvar
irf.favar
irf.bvar
irf
compirf
Documented in
compirf
irf
irf.bvar
irf.favar
irf.ftvar
irf.msvar
irf.tvar
#' @export
#' @title Function to draw one single path for the impulse-response functions.
#' @param Alpha the estimate of the VAR-coefficients
#' @param Sigma the estimate of the Variance-covariance matrix
#' @param id_obj S3-object for identifying structural shocks
#' @param nolags Number of lags in the model
#' @param intercept Whether the model has an intercept or not
#' @param nhor the horizon of the impulse-response functions
compirf <- function(Alpha, Sigma,id_obj, nolags, intercept = TRUE, nhor){
# Preliminaries
K <- nrow(Sigma)
bigj <- matrix(0, K, K * nolags)
bigj[1:K,1:K] <- diag(K)
if(intercept == TRUE){
B <- Alpha[-c(1),] # remove first row if intercept
}
else{
B <- Alpha
}
PHI_Mat <- companionmatrix(B,nolags)
biga <- PHI_Mat
bigai <- biga
# identify the model
shock <- structural(id_obj, Alpha, Sigma)
impresp <- matrix(0,K,K * nhor)
impresp[1:K,1:K] <- shock
for(ii in 1:(nhor-1)){
impresp[,(ii * K + 1):(( ii + 1) * K)] <- (bigj %*% bigai %*% t(bigj) %*% shock)
bigai <- bigai %*% biga
}
imp <- array(0,dim=c(K,K,nhor))
jj <- 0
for(ii in 1:K){
for(ij in 1:nhor){
jj <- ii + K * (ij - 1)
imp[ii,,ij] <- impresp[,jj]
}
}
return(imp)
}
#' @export
#' @title Compute Impulse-Response Functions
#' @param obj an S3 object of class bvar
#' @param id_obj an S3 object with information about identifiaction of the model
#' @param nhor horizon of the impulse-response function
#' @param irfquantiles quantiles for the impulse-response functions
#' @param ncores number of cores used
#' @param bootrep bootstrap replication for generalized impulse-response functions
#' @param ... currently not used
#' @return returns an S3-object of the class bvirf, fvirf, msirf or tirf
#' @details The function irf computes impulse-response functions for estimatet linear and non-linear VAR models. For Threshold-VARs and Factor-augmented Threshold-VARs generalized impulse-response functions are computed. The parameter bootrep determines how many bootstrap replications are going to be computed.
#' @references Gary Koop, M. Hashem Pesaran and Simon M. Potter, 1996, Impulse response analysis in nonlinear multivariate model, Journal of Econometrics 74(1), 119-147
#' @rdname irf
irf <- function(obj,id_obj,nhor=12,ncores=1,irfquantiles = c(0.05,0.95),...) UseMethod("irf")
#' @importFrom stats quantile
#' @importFrom stats median
#' @export
#' @rdname irf
irf.bvar <- function(obj, id_obj, nhor = 12, ncores = 1, irfquantiles = c(0.05,0.95),...){
# Declare variables
intercept <- obj$general_info$intercept
Alphadraws <- obj$mcmc_draws$Alpha
Sigmadraws <- obj$mcmc_draws$Sigma
nolags <- obj$general_info$nolags
nreps <- ( obj$general_info$nreps - obj$general_info$burnin ) / obj$general_info$nthin
K <- dim(Sigmadraws[,,1])[1]
irfdraws <- array(0,dim=c(K,K,nhor,nreps))
irffinal <- array(0,dim=c(K,K,nhor,3))
#chck <- requireNamespace("doParallel") && requireNamespace("foreach") && requireNamespace("parallel")
if(ncores>1 && !requireNamespace("foreach",quietly=TRUE)){
stop("The foreach package cannot be loaded.")
}
if(ncores == 1){ # No parallelization
for(ii in 1:nreps){
Alpha <- Alphadraws[,,ii]
Sigma <- Sigmadraws[,,ii]
irf <- compirf(Alpha = Alpha, Sigma = Sigma, id_obj = id_obj, nolags = nolags, intercept = intercept, nhor = nhor)
irfdraws[,,,ii] <- irf
}
}
else{ # parallel version
# Register workers
cl <- parallel::makeCluster(ncores)
doParallel::registerDoParallel(cl)
`%dopar%` <- foreach::`%dopar%`
xtmp <- foreach::foreach(ii = 1:nreps) %dopar% {
Alpha <- Alphadraws[,,ii]
Sigma <- Sigmadraws[,,ii]
irf <- compirf(Alpha = Alpha, Sigma = Sigma, id_obj = id_obj, nolags = nolags, intercept = intercept, nhor = nhor)
}
# stop workers
parallel::stopCluster(cl)
for(ii in 1:nreps){
irfdraws[,,,ii] <- xtmp[[ii]]
}
}
irflower <- min(irfquantiles)
irfupper <- max(irfquantiles)
for(jj in 1:K){
for(kk in 1:K){
for(ll in 1:nhor){
irffinal[jj,kk,ll,1] <- stats::median(irfdraws[jj,kk,ll,])
irffinal[jj,kk,ll,2] <- stats::quantile(irfdraws[jj,kk,ll,],probs=irflower)
irffinal[jj,kk,ll,3] <- stats::quantile(irfdraws[jj,kk,ll,],probs=irfupper)
}
}
}
relist <- structure(list(irfdraws = irffinal,
irfhorizon = nhor,
varnames = obj$data_info$var_names,
no_variables = obj$data_info$no_variables),
class = "bvirf")
return(relist)
}
#' @export
#' @rdname irf
irf.favar <- function(obj, id_obj, nhor=12, ncores=1, irfquantiles = c(0.05,0.95),...){
# Preliminaries
intercept <- obj$general_info$intercept
Betadraws <- obj$mcmc_draws$Alpha
Sigmadraws <- obj$mcmc_draws$Sigma
Ldraws <- obj$mcmc_draws$L
nolags <- obj$general_info$nolags
nreps <- dim(Betadraws)[3]
k <- dim(Sigmadraws)[1]
dimXY <- dim(Ldraws)[1]
irf_small_draws <- array(0,dim=c(k,k,nhor,nreps))
irf_large_draws <- array(0,dim=c(k,dimXY,nhor,nreps))
if(ncores>1 && !requireNamespace("foreach",quietly=TRUE)){
stop("The foreach package cannot be loaded.")
}
if(ncores == 1){
for(ii in 1:nreps){
Alpha <- Betadraws[,,ii]
Sigma <- Sigmadraws[,,ii]
L <- Ldraws[,,ii]
irf <- compirf(Alpha = Alpha, Sigma = Sigma, id_obj = id_obj, nolags = nolags, intercept = intercept, nhor = nhor)
irf_small_draws[,,,ii] <- irf
for(jj in 1:k){
irf_large_draws[jj,,,ii] <- L[,] %*% irf[jj,,]
}
}
}
else{
# Get impulse-response for VAR-system
# Register workers
cl <- parallel::makeCluster(ncores)
doParallel::registerDoParallel(cl)
`%dopar%` <- foreach::`%dopar%`
xtmp <- foreach::foreach(ii = 1:nreps) %dopar% {
Alpha <- Betadraws[,,ii]
Sigma <- Sigmadraws[,,ii]
irf <- compirf(Alpha = Alpha, Sigma = Sigma, id_obj = id_obj, nolags = nolags, intercept = intercept, nhor = nhor)
} # End getting IRFs
# Transform it to larger system
for(ii in 1:nreps){
L <- Ldraws[,,ii]
irf_small_draws <- xtmp[[ii]]
for(jj in 1:k){
irf_large_draws[jj,,,ii] <- L[,] %*% irf_small_draws
}
}# end transformation
}# End loop over parallel version
# Store values
IrfSmallFinal <- array(0,dim=c(k,k,nhor,3))
IrfLargeFinal <- array(0,dim=c(k,dimXY,nhor,3))
irflower <- min(irfquantiles)
irfupper <- max(irfquantiles)
for(jj in 1:k){
for(kk in 1:k){
for(ll in 1:nhor){
IrfSmallFinal[jj,kk,ll,1] <- quantile(irf_small_draws[jj,kk,ll,],probs=0.5)
IrfSmallFinal[jj,kk,ll,2] <- quantile(irf_small_draws[jj,kk,ll,],probs=irflower)
IrfSmallFinal[jj,kk,ll,3] <- quantile(irf_small_draws[jj,kk,ll,],probs=irfupper)
}
}
}
for(jj in 1:k){
for(kk in 1:dimXY){
for(ll in 1:nhor){
IrfLargeFinal[jj,kk,ll,1] <- quantile(irf_large_draws[jj,kk,ll,],probs=0.5)
IrfLargeFinal[jj,kk,ll,2] <- quantile(irf_large_draws[jj,kk,ll,],probs=irflower)
IrfLargeFinal[jj,kk,ll,3] <- quantile(irf_large_draws[jj,kk,ll,],probs=irfupper)
}
}
}
# Returning values
relist <- structure(list(irf = IrfLargeFinal,irfhorizon = nhor,varnames = obj$data_info$varnames,
factor_varnames = obj$factordata_info$varnames),class = "fvirf")
return(relist)
}
#' @export
#' @rdname irf
irf.ftvar <- function(obj, id_obj, nhor=12, ncores = 1, irfquantiles = c(0.05,0.95), bootrep = 50,...){
nLength <- dim(obj$mcmc_draws$Alpha)[4]
Alphadraws <- obj$mcmc_draws$Alpha
Sigmadraws <- obj$mcmc_draws$Sigma
Ldraws <- obj$mcmc_draws$L
thVar <- obj$general_info$thVar
tardraws <- obj$mcmc_draws$tardraws
deldraws <- obj$mcmc_draws$deldraws
nolags <- obj$general_info$nolags
intercept <- obj$general_info$intercept
P <- dim(Sigmadraws[,,,1])[1]
dimxy <- dim(Ldraws)[1]
irf_small_draws <- array(0,dim=c(P,P,nhor,2,nLength))
irf_large_draws <- array(0,dim=c(P,dimxy,nhor,2,nLength))
# get data
factors <- obj$factordata_info$factors
data <- obj$data_info$data
fy <- cbind(data,factors)
if(ncores>1 && !requireNamespace("foreach",quietly = TRUE)){
stop("The foreach package cannot be loaded.")
}
if(ncores == 1){ # No parallelization
for(ii in 1:nLength){
print(ii)
Beta <- Alphadraws[,,,ii]
Sigma <- Sigmadraws[,,,ii]
tart <- tardraws[ii]
L <- Ldraws[,,,ii]
thDelay <- deldraws[ii]
xsplit <- splitVariables(y = fy, lags = nolags, thDelay = thDelay, thresh = thVar,
tart = tart, intercept = intercept)
for(jj in 1:P){
xx <- tirf(y = xsplit$ystar,ytest = xsplit$ytest,
beta1 = Beta[,,1],beta2 = Beta[,,2],
sigma1 = Sigma[,,1], sigma2 = Sigma[,,2],
tar = tart, thVar = thVar,thDelay = thDelay,
nolags = nolags, irfhor = nhor, intercept = intercept,
shockvar = jj,
bootrep = bootrep,
id_obj = id_obj)
irf_small_draws[jj,,,1,ii] <- xx$irf1
irf_small_draws[jj,,,2,ii] <- xx$irf2
irf_large_draws[jj,,,1,ii] <- L[,,1] %*% irf_small_draws[jj,,,1,ii]
irf_large_draws[jj,,,2,ii] <- L[,,2] %*% irf_small_draws[jj,,,2,ii]
}
}
}
else{
# Register clusters
# Warning: No check if number of cores are supported
cl <- parallel::makeCluster(ncores)
doParallel::registerDoParallel(cl)
irf_draws <- array(0,dim=c(P,P,nhor,2))
`%dopar%` <- foreach::`%dopar%`
xtmp <- foreach::foreach(ii = 1:nLength) %dopar% {
Beta <- Alphadraws[,,,ii]
Sigma <- Sigmadraws[,,,ii]
tart <- tardraws[ii]
L <- Ldraws[,,,ii]
thDelay <- deldraws[ii]
irftmp <- tirf1(y = fy,
Alpha = Beta, Sigma = Sigma,
tart = tart, thVar = thVar, thDelay = thDelay,
nolags = nolags, nhor = nhor, intercept = intercept,
bootrep, id_obj = id_obj, P)
} # End of parallel loop
# shut down workers
parallel::stopCluster(cl)
for(ii in 1:nLength){
irf_small_draws[,,,,ii] <- xtmp[[ii]]
L <- Ldraws[,,,ii]
for(jj in 1:P){
irf_large_draws[jj,,,1,ii] <- L[,,1] %*% irf_small_draws[jj,,,1,ii]
irf_large_draws[jj,,,2,ii] <- L[,,2] %*% irf_small_draws[jj,,,2,ii]
}
}# End projecting small system onto larger system
}# End getting impulse-response-functions
# Prepare return values
IrfSmallFinal <- array(0,dim=c(P,P,nhor,2,3))
IrfLargeFinal <- array(0,dim=c(P,dimxy,nhor,2,3))
irflower <- min(irfquantiles)
irfupper <- max(irfquantiles)
for(jj in 1:P){
for(kk in 1:P){
for(ll in 1:nhor){
# Regime 1
IrfSmallFinal[jj,kk,ll,1,1] <- median(irf_small_draws[jj,kk,ll,1,])
IrfSmallFinal[jj,kk,ll,1,2] <- quantile(irf_small_draws[jj,kk,ll,1,],probs=irflower)
IrfSmallFinal[jj,kk,ll,1,3] <- quantile(irf_small_draws[jj,kk,ll,1,],probs=irfupper)
# Regime 2
IrfSmallFinal[jj,kk,ll,2,1] <- median(irf_small_draws[jj,kk,ll,2,])
IrfSmallFinal[jj,kk,ll,2,2] <- quantile(irf_small_draws[jj,kk,ll,2,],probs=irflower)
IrfSmallFinal[jj,kk,ll,2,3] <- quantile(irf_small_draws[jj,kk,ll,2,],probs=irfupper)
}
}
}
for(jj in 1:P){
for(kk in 1:dimxy){
for(ll in 1:nhor){
# Regime 1
IrfLargeFinal[jj,kk,ll,1,1] <- mean(irf_large_draws[jj,kk,ll,1,])
IrfLargeFinal[jj,kk,ll,1,2] <- quantile(irf_large_draws[jj,kk,ll,1,],probs=irflower)
IrfLargeFinal[jj,kk,ll,1,3] <- quantile(irf_large_draws[jj,kk,ll,1,],probs=irfupper)
# Regime 2
IrfLargeFinal[jj,kk,ll,2,1] <- mean(irf_large_draws[jj,kk,ll,2,])
IrfLargeFinal[jj,kk,ll,2,2] <- quantile(irf_large_draws[jj,kk,ll,2,],probs=irflower)
IrfLargeFinal[jj,kk,ll,2,3] <- quantile(irf_large_draws[jj,kk,ll,2,],probs=irfupper)
}
}
}
relist <- structure(list(irf=IrfLargeFinal,irfhorizon=nhor,varnames=obj$data_info$varnames,factor_varnames = obj$factordata_info$var_names),class = "tvirf")
return(relist)
}
#' @export
#' @rdname irf
irf.msvar <- function(obj, id_obj, nhor = 12, ncores = 1, irfquantiles = c(0.05,0.95),...){
#
# Extract information Variables
#
intercept <- obj$general_info$intercept
nolags <- obj$general_info$nolags
nreps <- obj$general_info$nreps
burnin <- obj$general_info$burnin
nthin <- obj$general_info$nthin
noregimes <- obj$general_info$noregimes
Alphadraws <- obj$mcmc_draws$Alpha
Sigmadraws <- obj$mcmc_draws$Sigma
K <- dim(obj$mcmc_draws$Sigma)[2]
totdraws <- floor((nreps - burnin) / nthin)
irfdraws <- array(NA,dim=c(K,K,nhor,noregimes,totdraws))
irffinal <- array(NA, dim =c(K,K,nhor,3,noregimes))
# Write a function check_parallel
if(ncores>1 && !requireNamespace("foreach",quietly = TRUE)){
stop("The foreach package cannot be loaded.")
}
if(ncores == 1){ # No Parallelization
for(jj in 1:noregimes){ # Loop over regimes
for(ii in 1:totdraws){ # Loop over all draws
Alpha <- Alphadraws[,,jj,ii]
Sigma <- Sigmadraws[,,jj,ii]
irf <- compirf(Alpha = Alpha, Sigma = Sigma, id_obj = id_obj,nolags = nolags, intercept = intercept, nhor = nhor)
irfdraws[,,,jj,ii] <- irf
} # End loop over all draws
} # End loop over regimes
}
else{
for(jj in 1:noregimes){ # Loop over regimes, not parallelized
cl <- parallel::makeCluster(ncores)
doParallel::registerDoParallel(cl)
`%dopar%` <- foreach::`%dopar%`
xtmp <- foreach::foreach(ii = 1:totdraws) %dopar% {
Alpha <- Alphadraws[,,jj,ii]
Sigma <- Sigmadraws[,,jj,ii]
irf <- compirf(Alpha = Alpha, Sigma = Sigma, id_obj = id_obj, nolags = nolags, intercept = intercept, nhor = nhor)
} # End parallel loop
# stop workers
parallel::stopCluster(cl)
for(ii in 1:totdraws){
irfdraws[,,,jj,ii] <- xtmp[[ii]]
}
}
}
irflower <- min(irfquantiles)
irfupper <- max(irfquantiles)
for(ii in 1:noregimes){
for(jj in 1:K){
for(kk in 1:K){
for(ll in 1:nhor){
irffinal[jj,kk,ll,1,ii] <- median(irfdraws[jj,kk,ll,ii,])
irffinal[jj,kk,ll,2,ii] <- quantile(irfdraws[jj,kk,ll,ii,],probs=irflower)
irffinal[jj,kk,ll,3,ii] <- quantile(irfdraws[jj,kk,ll,ii,],probs=irfupper)
}
}
}
}
relist <- structure(list(irfdraws = irffinal,
irfhorizon = nhor,
varnames = obj$data_info$var_names,
no_variables = obj$data_info$no_variables,
noregimes = noregimes),
class = "msirf")
return(relist)
}
#' @export
#' @rdname irf
irf.tvar <- function(obj, id_obj, nhor = 12, ncores = 1, irfquantiles = c(0.05,0.95), bootrep=50,...){
nLength <- (obj$general_info$nreps - obj$general_info$burnin) / obj$general_info$nthin
Alphadraws <- obj$mcmc_draws$Alpha
Sigmadraws <- obj$mcmc_draws$Sigma
thVar <- obj$general_info$thVar
tardraws <- obj$mcmc_draws$tardraws
deldraws <- obj$mcmc_draws$deldraws
nolags <- obj$general_info$nolags
intercept <- obj$general_info$intercept
if(ncores>1 && !requireNamespace("foreach",quietly = TRUE)){
stop("The foreach package cannot be loaded.")
}
K <- dim(Sigmadraws[,,,1])[1]
irfdraws <- array(0,dim=c(K,K,nhor,2,nLength))
if(ncores == 1){
for(jj in 1:nLength){
print(jj)
Alpha <- Alphadraws[,,,jj]
Sigma <- Sigmadraws[,,,jj]
tart <- tardraws[jj]
thDelay <- deldraws[jj]
xsplit <- splitVariables(y = obj$data_info$data,lags = nolags, thDelay = thDelay, thresh = thVar, tart = tart, intercept = intercept)
for(ii in 1:K){
xx <- tirf(xsplit$ystar,xsplit$ytest,Alpha[,,1],Alpha[,,2],Sigma[,,1],Sigma[,,2],tart,thVar,thDelay,nolags,nhor,intercept,shockvar=ii,bootrep,id_obj)
irfdraws[ii,,,1,jj] <- xx$irf1
irfdraws[ii,,,2,jj] <- xx$irf2
}
}
}
else{ # Parallel version
# Register workers
cl <- parallel::makeCluster(ncores)
doParallel::registerDoParallel(cl)
irftmp <- array(0,dim=c(K,K,nhor,2))
`%dopar%` <- foreach::`%dopar%`
xtmp <- foreach::foreach(jj = 1:nLength) %dopar% {
Alpha <- Alphadraws[,,,jj]
Sigma <- Sigmadraws[,,,jj]
tart <- tardraws[jj]
thDelay <- deldraws[jj]
irftmp <- tirf1(y = obj$data_info$data,
Alpha = Alpha,Sigma = Sigma,
tart = tart,thVar = thVar, thDelay = thDelay,
nolags = nolags, nhor = nhor, intercept = intercept,
bootrep = bootrep, K = K, id_obj = id_obj)
}
# stop workers
parallel::stopCluster(cl)
for(jj in 1:nLength){
irfdraws[,,,,jj] <- xtmp[[jj]]
}
}
irffinal <- array(0,dim=c(K,K,nhor,2,3))
lowerquantile=min(irfquantiles)
upperquantile=max(irfquantiles)
for(ii in 1:K){
for(jj in 1:K){
for(kk in 1:nhor){
irffinal[ii,jj,kk,1,1] <- quantile(irfdraws[ii,jj,kk,1,],probs=0.5)
irffinal[ii,jj,kk,2,1] <- quantile(irfdraws[ii,jj,kk,2,],probs=0.5)
irffinal[ii,jj,kk,1,2] <- quantile(irfdraws[ii,jj,kk,1,],probs=lowerquantile)
irffinal[ii,jj,kk,2,2] <- quantile(irfdraws[ii,jj,kk,2,],probs=lowerquantile)
irffinal[ii,jj,kk,1,3] <- quantile(irfdraws[ii,jj,kk,1,],probs=upperquantile)
irffinal[ii,jj,kk,2,3] <- quantile(irfdraws[ii,jj,kk,2,],probs=upperquantile)
}
}
}
relist <- structure(list(irf=irffinal,irfhorizon=nhor,varnames=obj$general_info$var_names),class = "tvirf")
return(relist)
}
joergrieger/bvar documentation built on July 3, 2020, 5:34 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 irf.tvar irf.msvar irf.ftvar irf.favar irf.bvar irf compirf
Documented in compirf irf irf.bvar irf.favar irf.ftvar irf.msvar irf.tvar
#' @export
#' @title Function to draw one single path for the impulse-response functions.
#' @param Alpha the estimate of the VAR-coefficients
#' @param Sigma the estimate of the Variance-covariance matrix
#' @param id_obj S3-object for identifying structural shocks
#' @param nolags Number of lags in the model
#' @param intercept Whether the model has an intercept or not
#' @param nhor the horizon of the impulse-response functions
compirf <- function(Alpha, Sigma,id_obj, nolags, intercept = TRUE, nhor){
# Preliminaries
K <- nrow(Sigma)
bigj <- matrix(0, K, K * nolags)
bigj[1:K,1:K] <- diag(K)
if(intercept == TRUE){
B <- Alpha[-c(1),] # remove first row if intercept
}
else{
B <- Alpha
}
PHI_Mat <- companionmatrix(B,nolags)
biga <- PHI_Mat
bigai <- biga
# identify the model
shock <- structural(id_obj, Alpha, Sigma)
impresp <- matrix(0,K,K * nhor)
impresp[1:K,1:K] <- shock
for(ii in 1:(nhor-1)){
impresp[,(ii * K + 1):(( ii + 1) * K)] <- (bigj %*% bigai %*% t(bigj) %*% shock)
bigai <- bigai %*% biga
}
imp <- array(0,dim=c(K,K,nhor))
jj <- 0
for(ii in 1:K){
for(ij in 1:nhor){
jj <- ii + K * (ij - 1)
imp[ii,,ij] <- impresp[,jj]
}
}
return(imp)
}
#' @export
#' @title Compute Impulse-Response Functions
#' @param obj an S3 object of class bvar
#' @param id_obj an S3 object with information about identifiaction of the model
#' @param nhor horizon of the impulse-response function
#' @param irfquantiles quantiles for the impulse-response functions
#' @param ncores number of cores used
#' @param bootrep bootstrap replication for generalized impulse-response functions
#' @param ... currently not used
#' @return returns an S3-object of the class bvirf, fvirf, msirf or tirf
#' @details The function irf computes impulse-response functions for estimatet linear and non-linear VAR models. For Threshold-VARs and Factor-augmented Threshold-VARs generalized impulse-response functions are computed. The parameter bootrep determines how many bootstrap replications are going to be computed.
#' @references Gary Koop, M. Hashem Pesaran and Simon M. Potter, 1996, Impulse response analysis in nonlinear multivariate model, Journal of Econometrics 74(1), 119-147
#' @rdname irf
irf <- function(obj,id_obj,nhor=12,ncores=1,irfquantiles = c(0.05,0.95),...) UseMethod("irf")
#' @importFrom stats quantile
#' @importFrom stats median
#' @export
#' @rdname irf
irf.bvar <- function(obj, id_obj, nhor = 12, ncores = 1, irfquantiles = c(0.05,0.95),...){
# Declare variables
intercept <- obj$general_info$intercept
Alphadraws <- obj$mcmc_draws$Alpha
Sigmadraws <- obj$mcmc_draws$Sigma
nolags <- obj$general_info$nolags
nreps <- ( obj$general_info$nreps - obj$general_info$burnin ) / obj$general_info$nthin
K <- dim(Sigmadraws[,,1])[1]
irfdraws <- array(0,dim=c(K,K,nhor,nreps))
irffinal <- array(0,dim=c(K,K,nhor,3))
#chck <- requireNamespace("doParallel") && requireNamespace("foreach") && requireNamespace("parallel")
if(ncores>1 && !requireNamespace("foreach",quietly=TRUE)){
stop("The foreach package cannot be loaded.")
}
if(ncores == 1){ # No parallelization
for(ii in 1:nreps){
Alpha <- Alphadraws[,,ii]
Sigma <- Sigmadraws[,,ii]
irf <- compirf(Alpha = Alpha, Sigma = Sigma, id_obj = id_obj, nolags = nolags, intercept = intercept, nhor = nhor)
irfdraws[,,,ii] <- irf
}
}
else{ # parallel version
# Register workers
cl <- parallel::makeCluster(ncores)
doParallel::registerDoParallel(cl)
`%dopar%` <- foreach::`%dopar%`
xtmp <- foreach::foreach(ii = 1:nreps) %dopar% {
Alpha <- Alphadraws[,,ii]
Sigma <- Sigmadraws[,,ii]
irf <- compirf(Alpha = Alpha, Sigma = Sigma, id_obj = id_obj, nolags = nolags, intercept = intercept, nhor = nhor)
}
# stop workers
parallel::stopCluster(cl)
for(ii in 1:nreps){
irfdraws[,,,ii] <- xtmp[[ii]]
}
}
irflower <- min(irfquantiles)
irfupper <- max(irfquantiles)
for(jj in 1:K){
for(kk in 1:K){
for(ll in 1:nhor){
irffinal[jj,kk,ll,1] <- stats::median(irfdraws[jj,kk,ll,])
irffinal[jj,kk,ll,2] <- stats::quantile(irfdraws[jj,kk,ll,],probs=irflower)
irffinal[jj,kk,ll,3] <- stats::quantile(irfdraws[jj,kk,ll,],probs=irfupper)
}
}
}
relist <- structure(list(irfdraws = irffinal,
irfhorizon = nhor,
varnames = obj$data_info$var_names,
no_variables = obj$data_info$no_variables),
class = "bvirf")
return(relist)
}
#' @export
#' @rdname irf
irf.favar <- function(obj, id_obj, nhor=12, ncores=1, irfquantiles = c(0.05,0.95),...){
# Preliminaries
intercept <- obj$general_info$intercept
Betadraws <- obj$mcmc_draws$Alpha
Sigmadraws <- obj$mcmc_draws$Sigma
Ldraws <- obj$mcmc_draws$L
nolags <- obj$general_info$nolags
nreps <- dim(Betadraws)[3]
k <- dim(Sigmadraws)[1]
dimXY <- dim(Ldraws)[1]
irf_small_draws <- array(0,dim=c(k,k,nhor,nreps))
irf_large_draws <- array(0,dim=c(k,dimXY,nhor,nreps))
if(ncores>1 && !requireNamespace("foreach",quietly=TRUE)){
stop("The foreach package cannot be loaded.")
}
if(ncores == 1){
for(ii in 1:nreps){
Alpha <- Betadraws[,,ii]
Sigma <- Sigmadraws[,,ii]
L <- Ldraws[,,ii]
irf <- compirf(Alpha = Alpha, Sigma = Sigma, id_obj = id_obj, nolags = nolags, intercept = intercept, nhor = nhor)
irf_small_draws[,,,ii] <- irf
for(jj in 1:k){
irf_large_draws[jj,,,ii] <- L[,] %*% irf[jj,,]
}
}
}
else{
# Get impulse-response for VAR-system
# Register workers
cl <- parallel::makeCluster(ncores)
doParallel::registerDoParallel(cl)
`%dopar%` <- foreach::`%dopar%`
xtmp <- foreach::foreach(ii = 1:nreps) %dopar% {
Alpha <- Betadraws[,,ii]
Sigma <- Sigmadraws[,,ii]
irf <- compirf(Alpha = Alpha, Sigma = Sigma, id_obj = id_obj, nolags = nolags, intercept = intercept, nhor = nhor)
} # End getting IRFs
# Transform it to larger system
for(ii in 1:nreps){
L <- Ldraws[,,ii]
irf_small_draws <- xtmp[[ii]]
for(jj in 1:k){
irf_large_draws[jj,,,ii] <- L[,] %*% irf_small_draws
}
}# end transformation
}# End loop over parallel version
# Store values
IrfSmallFinal <- array(0,dim=c(k,k,nhor,3))
IrfLargeFinal <- array(0,dim=c(k,dimXY,nhor,3))
irflower <- min(irfquantiles)
irfupper <- max(irfquantiles)
for(jj in 1:k){
for(kk in 1:k){
for(ll in 1:nhor){
IrfSmallFinal[jj,kk,ll,1] <- quantile(irf_small_draws[jj,kk,ll,],probs=0.5)
IrfSmallFinal[jj,kk,ll,2] <- quantile(irf_small_draws[jj,kk,ll,],probs=irflower)
IrfSmallFinal[jj,kk,ll,3] <- quantile(irf_small_draws[jj,kk,ll,],probs=irfupper)
}
}
}
for(jj in 1:k){
for(kk in 1:dimXY){
for(ll in 1:nhor){
IrfLargeFinal[jj,kk,ll,1] <- quantile(irf_large_draws[jj,kk,ll,],probs=0.5)
IrfLargeFinal[jj,kk,ll,2] <- quantile(irf_large_draws[jj,kk,ll,],probs=irflower)
IrfLargeFinal[jj,kk,ll,3] <- quantile(irf_large_draws[jj,kk,ll,],probs=irfupper)
}
}
}
# Returning values
relist <- structure(list(irf = IrfLargeFinal,irfhorizon = nhor,varnames = obj$data_info$varnames,
factor_varnames = obj$factordata_info$varnames),class = "fvirf")
return(relist)
}
#' @export
#' @rdname irf
irf.ftvar <- function(obj, id_obj, nhor=12, ncores = 1, irfquantiles = c(0.05,0.95), bootrep = 50,...){
nLength <- dim(obj$mcmc_draws$Alpha)[4]
Alphadraws <- obj$mcmc_draws$Alpha
Sigmadraws <- obj$mcmc_draws$Sigma
Ldraws <- obj$mcmc_draws$L
thVar <- obj$general_info$thVar
tardraws <- obj$mcmc_draws$tardraws
deldraws <- obj$mcmc_draws$deldraws
nolags <- obj$general_info$nolags
intercept <- obj$general_info$intercept
P <- dim(Sigmadraws[,,,1])[1]
dimxy <- dim(Ldraws)[1]
irf_small_draws <- array(0,dim=c(P,P,nhor,2,nLength))
irf_large_draws <- array(0,dim=c(P,dimxy,nhor,2,nLength))
# get data
factors <- obj$factordata_info$factors
data <- obj$data_info$data
fy <- cbind(data,factors)
if(ncores>1 && !requireNamespace("foreach",quietly = TRUE)){
stop("The foreach package cannot be loaded.")
}
if(ncores == 1){ # No parallelization
for(ii in 1:nLength){
print(ii)
Beta <- Alphadraws[,,,ii]
Sigma <- Sigmadraws[,,,ii]
tart <- tardraws[ii]
L <- Ldraws[,,,ii]
thDelay <- deldraws[ii]
xsplit <- splitVariables(y = fy, lags = nolags, thDelay = thDelay, thresh = thVar,
tart = tart, intercept = intercept)
for(jj in 1:P){
xx <- tirf(y = xsplit$ystar,ytest = xsplit$ytest,
beta1 = Beta[,,1],beta2 = Beta[,,2],
sigma1 = Sigma[,,1], sigma2 = Sigma[,,2],
tar = tart, thVar = thVar,thDelay = thDelay,
nolags = nolags, irfhor = nhor, intercept = intercept,
shockvar = jj,
bootrep = bootrep,
id_obj = id_obj)
irf_small_draws[jj,,,1,ii] <- xx$irf1
irf_small_draws[jj,,,2,ii] <- xx$irf2
irf_large_draws[jj,,,1,ii] <- L[,,1] %*% irf_small_draws[jj,,,1,ii]
irf_large_draws[jj,,,2,ii] <- L[,,2] %*% irf_small_draws[jj,,,2,ii]
}
}
}
else{
# Register clusters
# Warning: No check if number of cores are supported
cl <- parallel::makeCluster(ncores)
doParallel::registerDoParallel(cl)
irf_draws <- array(0,dim=c(P,P,nhor,2))
`%dopar%` <- foreach::`%dopar%`
xtmp <- foreach::foreach(ii = 1:nLength) %dopar% {
Beta <- Alphadraws[,,,ii]
Sigma <- Sigmadraws[,,,ii]
tart <- tardraws[ii]
L <- Ldraws[,,,ii]
thDelay <- deldraws[ii]
irftmp <- tirf1(y = fy,
Alpha = Beta, Sigma = Sigma,
tart = tart, thVar = thVar, thDelay = thDelay,
nolags = nolags, nhor = nhor, intercept = intercept,
bootrep, id_obj = id_obj, P)
} # End of parallel loop
# shut down workers
parallel::stopCluster(cl)
for(ii in 1:nLength){
irf_small_draws[,,,,ii] <- xtmp[[ii]]
L <- Ldraws[,,,ii]
for(jj in 1:P){
irf_large_draws[jj,,,1,ii] <- L[,,1] %*% irf_small_draws[jj,,,1,ii]
irf_large_draws[jj,,,2,ii] <- L[,,2] %*% irf_small_draws[jj,,,2,ii]
}
}# End projecting small system onto larger system
}# End getting impulse-response-functions
# Prepare return values
IrfSmallFinal <- array(0,dim=c(P,P,nhor,2,3))
IrfLargeFinal <- array(0,dim=c(P,dimxy,nhor,2,3))
irflower <- min(irfquantiles)
irfupper <- max(irfquantiles)
for(jj in 1:P){
for(kk in 1:P){
for(ll in 1:nhor){
# Regime 1
IrfSmallFinal[jj,kk,ll,1,1] <- median(irf_small_draws[jj,kk,ll,1,])
IrfSmallFinal[jj,kk,ll,1,2] <- quantile(irf_small_draws[jj,kk,ll,1,],probs=irflower)
IrfSmallFinal[jj,kk,ll,1,3] <- quantile(irf_small_draws[jj,kk,ll,1,],probs=irfupper)
# Regime 2
IrfSmallFinal[jj,kk,ll,2,1] <- median(irf_small_draws[jj,kk,ll,2,])
IrfSmallFinal[jj,kk,ll,2,2] <- quantile(irf_small_draws[jj,kk,ll,2,],probs=irflower)
IrfSmallFinal[jj,kk,ll,2,3] <- quantile(irf_small_draws[jj,kk,ll,2,],probs=irfupper)
}
}
}
for(jj in 1:P){
for(kk in 1:dimxy){
for(ll in 1:nhor){
# Regime 1
IrfLargeFinal[jj,kk,ll,1,1] <- mean(irf_large_draws[jj,kk,ll,1,])
IrfLargeFinal[jj,kk,ll,1,2] <- quantile(irf_large_draws[jj,kk,ll,1,],probs=irflower)
IrfLargeFinal[jj,kk,ll,1,3] <- quantile(irf_large_draws[jj,kk,ll,1,],probs=irfupper)
# Regime 2
IrfLargeFinal[jj,kk,ll,2,1] <- mean(irf_large_draws[jj,kk,ll,2,])
IrfLargeFinal[jj,kk,ll,2,2] <- quantile(irf_large_draws[jj,kk,ll,2,],probs=irflower)
IrfLargeFinal[jj,kk,ll,2,3] <- quantile(irf_large_draws[jj,kk,ll,2,],probs=irfupper)
}
}
}
relist <- structure(list(irf=IrfLargeFinal,irfhorizon=nhor,varnames=obj$data_info$varnames,factor_varnames = obj$factordata_info$var_names),class = "tvirf")
return(relist)
}
#' @export
#' @rdname irf
irf.msvar <- function(obj, id_obj, nhor = 12, ncores = 1, irfquantiles = c(0.05,0.95),...){
#
# Extract information Variables
#
intercept <- obj$general_info$intercept
nolags <- obj$general_info$nolags
nreps <- obj$general_info$nreps
burnin <- obj$general_info$burnin
nthin <- obj$general_info$nthin
noregimes <- obj$general_info$noregimes
Alphadraws <- obj$mcmc_draws$Alpha
Sigmadraws <- obj$mcmc_draws$Sigma
K <- dim(obj$mcmc_draws$Sigma)[2]
totdraws <- floor((nreps - burnin) / nthin)
irfdraws <- array(NA,dim=c(K,K,nhor,noregimes,totdraws))
irffinal <- array(NA, dim =c(K,K,nhor,3,noregimes))
# Write a function check_parallel
if(ncores>1 && !requireNamespace("foreach",quietly = TRUE)){
stop("The foreach package cannot be loaded.")
}
if(ncores == 1){ # No Parallelization
for(jj in 1:noregimes){ # Loop over regimes
for(ii in 1:totdraws){ # Loop over all draws
Alpha <- Alphadraws[,,jj,ii]
Sigma <- Sigmadraws[,,jj,ii]
irf <- compirf(Alpha = Alpha, Sigma = Sigma, id_obj = id_obj,nolags = nolags, intercept = intercept, nhor = nhor)
irfdraws[,,,jj,ii] <- irf
} # End loop over all draws
} # End loop over regimes
}
else{
for(jj in 1:noregimes){ # Loop over regimes, not parallelized
cl <- parallel::makeCluster(ncores)
doParallel::registerDoParallel(cl)
`%dopar%` <- foreach::`%dopar%`
xtmp <- foreach::foreach(ii = 1:totdraws) %dopar% {
Alpha <- Alphadraws[,,jj,ii]
Sigma <- Sigmadraws[,,jj,ii]
irf <- compirf(Alpha = Alpha, Sigma = Sigma, id_obj = id_obj, nolags = nolags, intercept = intercept, nhor = nhor)
} # End parallel loop
# stop workers
parallel::stopCluster(cl)
for(ii in 1:totdraws){
irfdraws[,,,jj,ii] <- xtmp[[ii]]
}
}
}
irflower <- min(irfquantiles)
irfupper <- max(irfquantiles)
for(ii in 1:noregimes){
for(jj in 1:K){
for(kk in 1:K){
for(ll in 1:nhor){
irffinal[jj,kk,ll,1,ii] <- median(irfdraws[jj,kk,ll,ii,])
irffinal[jj,kk,ll,2,ii] <- quantile(irfdraws[jj,kk,ll,ii,],probs=irflower)
irffinal[jj,kk,ll,3,ii] <- quantile(irfdraws[jj,kk,ll,ii,],probs=irfupper)
}
}
}
}
relist <- structure(list(irfdraws = irffinal,
irfhorizon = nhor,
varnames = obj$data_info$var_names,
no_variables = obj$data_info$no_variables,
noregimes = noregimes),
class = "msirf")
return(relist)
}
#' @export
#' @rdname irf
irf.tvar <- function(obj, id_obj, nhor = 12, ncores = 1, irfquantiles = c(0.05,0.95), bootrep=50,...){
nLength <- (obj$general_info$nreps - obj$general_info$burnin) / obj$general_info$nthin
Alphadraws <- obj$mcmc_draws$Alpha
Sigmadraws <- obj$mcmc_draws$Sigma
thVar <- obj$general_info$thVar
tardraws <- obj$mcmc_draws$tardraws
deldraws <- obj$mcmc_draws$deldraws
nolags <- obj$general_info$nolags
intercept <- obj$general_info$intercept
if(ncores>1 && !requireNamespace("foreach",quietly = TRUE)){
stop("The foreach package cannot be loaded.")
}
K <- dim(Sigmadraws[,,,1])[1]
irfdraws <- array(0,dim=c(K,K,nhor,2,nLength))
if(ncores == 1){
for(jj in 1:nLength){
print(jj)
Alpha <- Alphadraws[,,,jj]
Sigma <- Sigmadraws[,,,jj]
tart <- tardraws[jj]
thDelay <- deldraws[jj]
xsplit <- splitVariables(y = obj$data_info$data,lags = nolags, thDelay = thDelay, thresh = thVar, tart = tart, intercept = intercept)
for(ii in 1:K){
xx <- tirf(xsplit$ystar,xsplit$ytest,Alpha[,,1],Alpha[,,2],Sigma[,,1],Sigma[,,2],tart,thVar,thDelay,nolags,nhor,intercept,shockvar=ii,bootrep,id_obj)
irfdraws[ii,,,1,jj] <- xx$irf1
irfdraws[ii,,,2,jj] <- xx$irf2
}
}
}
else{ # Parallel version
# Register workers
cl <- parallel::makeCluster(ncores)
doParallel::registerDoParallel(cl)
irftmp <- array(0,dim=c(K,K,nhor,2))
`%dopar%` <- foreach::`%dopar%`
xtmp <- foreach::foreach(jj = 1:nLength) %dopar% {
Alpha <- Alphadraws[,,,jj]
Sigma <- Sigmadraws[,,,jj]
tart <- tardraws[jj]
thDelay <- deldraws[jj]
irftmp <- tirf1(y = obj$data_info$data,
Alpha = Alpha,Sigma = Sigma,
tart = tart,thVar = thVar, thDelay = thDelay,
nolags = nolags, nhor = nhor, intercept = intercept,
bootrep = bootrep, K = K, id_obj = id_obj)
}
# stop workers
parallel::stopCluster(cl)
for(jj in 1:nLength){
irfdraws[,,,,jj] <- xtmp[[jj]]
}
}
irffinal <- array(0,dim=c(K,K,nhor,2,3))
lowerquantile=min(irfquantiles)
upperquantile=max(irfquantiles)
for(ii in 1:K){
for(jj in 1:K){
for(kk in 1:nhor){
irffinal[ii,jj,kk,1,1] <- quantile(irfdraws[ii,jj,kk,1,],probs=0.5)
irffinal[ii,jj,kk,2,1] <- quantile(irfdraws[ii,jj,kk,2,],probs=0.5)
irffinal[ii,jj,kk,1,2] <- quantile(irfdraws[ii,jj,kk,1,],probs=lowerquantile)
irffinal[ii,jj,kk,2,2] <- quantile(irfdraws[ii,jj,kk,2,],probs=lowerquantile)
irffinal[ii,jj,kk,1,3] <- quantile(irfdraws[ii,jj,kk,1,],probs=upperquantile)
irffinal[ii,jj,kk,2,3] <- quantile(irfdraws[ii,jj,kk,2,],probs=upperquantile)
}
}
}
relist <- structure(list(irf=irffinal,irfhorizon=nhor,varnames=obj$general_info$var_names),class = "tvirf")
return(relist)
}
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.