R/fevd.r
In joergrieger/bvar: Estimation and forecasting of bayesian VAR models
Defines functions
fevd.msvar
fevd.tvar
fevd.bvar
fevd
fevd1
gfevd
Documented in
fevd
fevd.bvar
fevd.msvar
fevd.tvar
#' @title internal function to calculate a generalized fevd
#' @param Alpha VAR-coefficients
#' @param Sigma Variance-Covariance matrix
#' @param h fevd-horizon
#' @param intercept whether the model had an intercept or not
#' @param nolags the number of lags in the model
#' @noRd
gfevd <- function(Alpha,Sigma,h,intercept,nolags){
nvar <- dim(Sigma)[1]
# if the model has an intercept delete the first row
if(intercept){
Alpha <- Alpha[-c(1),]
}
varcoef <- array(0,dim=c(nvar,nvar,nolags))
for(ii in 1:nolags){
varcoef[,,ii] <- Alpha[((ii - 1) * nvar + 1):(ii * nvar),]
}
# transform into vma form
vmacoef <- var2vma(varcoef,h,nolags,nvar)
# Calculate the generalized fevd
deltatemp <- array(0,dim=c(nvar,nvar))
for(i in 1:nvar){
ei <- array(0,dim=c(nvar))
ei[i] <- 1
tsum2 <- 0
for(hh in 1:h){
tmp <- t(ei) %*% vmacoef[,,hh] %*% Sigma[,] %*% t(vmacoef[,,hh]) %*% ei
tsum2 <- tsum2 + tmp
}
for(j in 1:nvar){
ej <- array(0, dim=c(nvar))
ej[j] <- 1
tsum <- 0
for(hh in 1:h){
tmp <- (ei %*% vmacoef[,,hh] %*% ej)^2
tsum <- tsum + tmp
}
deltatemp[i,j] <- Sigma[i,i]^(-1)*tsum/tsum2
}
}
return(deltatemp)
}
#' @title internal function to compute a forecast error variance decomposition
#' @param Alpha VAR-coefficients
#' @param Variance-Covariance matrix
#' @param h fevd-horizon
#' @param intercept whether the model had an intercept or not
#' @param nolags number of lags in the VAR model
#' @noRd
fevd1 <- function(Alpha,Sigma,h,intercept,nolags,id_obj){
nvar <- dim(Sigma)[1]
J <- array(0,dim=c(nvar,nvar*nolags))
J[1:nvar,1:nvar] <- diag(1,nvar)
# Coefficients in compannion form
if(intercept){
Alpha <- Alpha[-c(1),]
}
A <- companionmatrix(Alpha,nolags)
A0 <- diag(1,dim(A)[1])
# Start computing fevd
idSigma <- structural(id_obj, Alpha, Sigma)
TH1 <- J %*% A0 %*% t(J)
TH <- t(TH1 %*% idSigma)
TH2 <- (TH * TH)
TH3 <- TH2
if(h > 1){
B = A
for(i in 2:h){
TH <- t(J %*% B %*% t(J) %*% idSigma)
TH2 <- TH * TH
TH3 <- TH3 + TH2
B = B %*% A
}
}
return(TH3)
}
#' @export
#' @title Forecast error variance decomposition
#' @param obj S3 object generated by bvar, tvar or msvar
#' @param h forecast horizon
#' @param id_obj S3 object containing information about identifying the VAR model
#' @param type Type of forecast error variance decomposition. If type = "general" the generalized forecast error variance decomposition is used. For any other values the standard method. While the generalized fevd doesn't use an indefication method to identify structural shocks in the VAR model, this is not true for the standard method. Hence for the generalized fevd there is no need to supply the function via id_obj. If id_obj is not null, there will be a warning message but will still proceed.
#' @param normalize whether the fevd should be normalized or not
#' @param ... currently not used
#' @rdname fevd
fevd <- function(obj,h=12,id_obj=NULL,type="general",normalize = TRUE,...) UseMethod("fevd")
#' @export
#' @rdname fevd
fevd.bvar <- function(obj,h=12,id_obj=NULL,type="general",normalize = TRUE,...){
if(type == "general" && !is.null(id_obj)){
message("Ignoring identification method for generalized forecast error variance decomposition.\n")
}
Alpha <- obj$mcmc_draws$Alpha
Sigma <- obj$mcmc_draws$Sigma
nreps <- dim(Alpha)[3]
nVar <- dim(Alpha)[2]
intercept <- obj$general_info$intercept
nolags <- obj$general_info$nolags
fevd_return <- array(0,dim=c(nVar,nVar,1,nreps))
if(type == "general"){
for(ii in 1:nreps){
current_Alpha <- Alpha[,,ii]
current_Sigma <- Sigma[,,ii]
fevd_return[,,1,ii] <- gfevd(Alpha = current_Alpha,Sigma = current_Sigma, h = h, nolags = nolags, intercept = intercept)#deltatemp
}
}
else if(type == "standard"){
for(ii in 1:nreps){
current_Alpha <- Alpha[,,ii]
current_Sigma <- Sigma[,,ii]
fevd_return[,,1,ii] <- fevd1(Alpha = current_Alpha, Sigma = current_Sigma,h=h,nolags=nolags,intercept=intercept,id_obj=id_obj)
}
}
else{
stop("The type of the forecast error variance decomposition can be 'general' or 'standard'")
}
fevd_ret <- 0
for(ii in 1:nreps){
fevd_ret <- fevd_ret + fevd_return[,,1,ii]
}
fevd_ret <- fevd_ret / nreps
if(normalize){
t1 <- rowSums(fevd_ret)
for(ii in 1:nVar){
fevd_ret[ii,] <- fevd_ret[ii,] / t1[ii]
}
}
retlist = structure(list(fevd = fevd_ret,noregimes = 1),class="bvfevd")
return(retlist)
}
#' @export
#' @rdname fevd
fevd.tvar <- function(obj,h=12,id_obj=NULL,type="general",normalize = TRUE,...){
Alpha <- obj$mcmc_draws$Alpha
Sigma <- obj$mcmc_draws$Sigma
nreps <- dim(Alpha)[3]
nvar <- dim(Alpha)[2]
intercept <- obj$general_info$intercept
nolags <- obj$general_info$nolags
fevd_return <- array(0,dim=c(nvar,nvar,2,nreps))
if(type == "general"){
for(ii in 1:nreps){
current_Alpha <- Alpha[,,1,ii]
current_Sigma <- Sigma[,,1,ii]
fevd_return[,,1,ii] <- gfevd(Alpha = current_Alpha,Sigma = current_Sigma, h = h,
nolags = nolags, intercept = intercept)
current_Alpha <- Alpha[,,2,ii]
current_Sigma <- Sigma[,,2,ii]
fevd_return[,,2,ii] <- gfevd(Alpha = current_Alpha,Sigma = current_Sigma, h = h,
nolags = nolags, intercept = intercept)
}
}
else if(type == "standard"){
for(ii in 1:nreps){
current_Alpha <- Alpha[,,1,ii]
current_Sigma <- Sigma[,,1,ii]
fevd_return[,,1,ii] <- fevd1(Alpha = current_Alpha,Sigma = current_Sigma, h = h,
nolags = nolags, intercept = intercept)
current_Alpha <- Alpha[,,2,ii]
current_Sigma <- Sigma[,,2,ii]
fevd_return[,,2,ii] <- fevd1(Alpha = current_Alpha,Sigma = current_Sigma, h = h,
nolags = nolags, intercept = intercept)
}
}
else{
stop("Currently only generalized Forecast Error Variance Decomposition is implemented")
}
fevd_ret <- array(0,dim=c(nvar,nvar,2))
for(ii in 1:nreps){
fevd_ret[,,1] <- fevd_ret[,,1] + fevd_return[,,1,ii]
fevd_ret[,,2] <- fevd_ret[,,2] + fevd_return[,,2,ii]
}
fevd_ret <- fevd_ret / nreps
if(normalize){
t1 <- rowSums(fevd_ret[,,1])
t2 <- rowSums(fevd_ret[,,2])
for(ii in 1:nvar){
fevd_ret[ii,,1] <- fevd_ret[ii,,1] / t1[ii]
fevd_ret[ii,,2] <- fevd_ret[ii,,2] / t2[ii]
}
}
retlist = structure(list(fevd = fevd_ret,noregimes = 2),class="bvfevd")
return(retlist)
}
#' @export
#' @rdname fevd
fevd.msvar <- function(obj,h=12,id_obj=NULL,type="general",normalize = TRUE,...){
Alpha <- obj$mcmc_draws$Alpha
Sigma <- obj$mcmc_draws$Sigma
nreps <- dim(Alpha)[3]
nvar <- dim(Alpha)[2]
intercept <- obj$general_info$intercept
nolags <- obj$general_info$nolags
noregimes <- obj$general_info$noregimes
fevd_return <- array(0,dim=c(nvar,nvar,noregimes,nreps))
if(type == "general"){
for(ii in 1:nreps){
for(jj in 1:noregimes){
current_Alpha <- Alpha[,,jj,ii]
current_Sigma <- Sigma[,,jj,ii]
fevd_return[,,jj,ii] <- gfevd(Alpha = current_Alpha,Sigma = current_Sigma, h = h,
nolags = nolags, intercept = intercept)
}
}
}
else if(type == "standard"){
for(ii in 1:nreps){
for(jj in 1:noregimes){
current_Alpha <- Alpha[,,jj,ii]
current_Sigma <- Sigma[,,jj,ii]
fevd_return[,,jj,ii] <- fevd1(Alpha = current_Alpha,Sigma = current_Sigma, h = h,
nolags = nolags, intercept = intercept)
}
}
}
else{
stop("Currently only generalized Forecast Error Variance Decomposition is implemented")
}
fevd_ret <- array(0,dim=c(nvar,nvar,noregimes))
for(ii in 1:nreps){
for(jj in 1:noregimes){
fevd_ret[,,jj] <- fevd_ret[,,jj] + fevd_return[,,jj,ii]
}
}
fevd_ret <- fevd_ret / nreps
if(normalize){
for(jj in 1:noregimes){
t1 <- rowSums(fevd_ret[,,jj])
for(ii in 1:nvar){
fevd_ret[ii,,jj] <- fevd_ret[ii,,jj] / t1[ii]
}
}
}
retlist = structure(list(fevd = fevd_ret,noregimes = noregimes),class="bvfevd")
return(retlist)
}
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 fevd.msvar fevd.tvar fevd.bvar fevd fevd1 gfevd
Documented in fevd fevd.bvar fevd.msvar fevd.tvar
#' @title internal function to calculate a generalized fevd
#' @param Alpha VAR-coefficients
#' @param Sigma Variance-Covariance matrix
#' @param h fevd-horizon
#' @param intercept whether the model had an intercept or not
#' @param nolags the number of lags in the model
#' @noRd
gfevd <- function(Alpha,Sigma,h,intercept,nolags){
nvar <- dim(Sigma)[1]
# if the model has an intercept delete the first row
if(intercept){
Alpha <- Alpha[-c(1),]
}
varcoef <- array(0,dim=c(nvar,nvar,nolags))
for(ii in 1:nolags){
varcoef[,,ii] <- Alpha[((ii - 1) * nvar + 1):(ii * nvar),]
}
# transform into vma form
vmacoef <- var2vma(varcoef,h,nolags,nvar)
# Calculate the generalized fevd
deltatemp <- array(0,dim=c(nvar,nvar))
for(i in 1:nvar){
ei <- array(0,dim=c(nvar))
ei[i] <- 1
tsum2 <- 0
for(hh in 1:h){
tmp <- t(ei) %*% vmacoef[,,hh] %*% Sigma[,] %*% t(vmacoef[,,hh]) %*% ei
tsum2 <- tsum2 + tmp
}
for(j in 1:nvar){
ej <- array(0, dim=c(nvar))
ej[j] <- 1
tsum <- 0
for(hh in 1:h){
tmp <- (ei %*% vmacoef[,,hh] %*% ej)^2
tsum <- tsum + tmp
}
deltatemp[i,j] <- Sigma[i,i]^(-1)*tsum/tsum2
}
}
return(deltatemp)
}
#' @title internal function to compute a forecast error variance decomposition
#' @param Alpha VAR-coefficients
#' @param Variance-Covariance matrix
#' @param h fevd-horizon
#' @param intercept whether the model had an intercept or not
#' @param nolags number of lags in the VAR model
#' @noRd
fevd1 <- function(Alpha,Sigma,h,intercept,nolags,id_obj){
nvar <- dim(Sigma)[1]
J <- array(0,dim=c(nvar,nvar*nolags))
J[1:nvar,1:nvar] <- diag(1,nvar)
# Coefficients in compannion form
if(intercept){
Alpha <- Alpha[-c(1),]
}
A <- companionmatrix(Alpha,nolags)
A0 <- diag(1,dim(A)[1])
# Start computing fevd
idSigma <- structural(id_obj, Alpha, Sigma)
TH1 <- J %*% A0 %*% t(J)
TH <- t(TH1 %*% idSigma)
TH2 <- (TH * TH)
TH3 <- TH2
if(h > 1){
B = A
for(i in 2:h){
TH <- t(J %*% B %*% t(J) %*% idSigma)
TH2 <- TH * TH
TH3 <- TH3 + TH2
B = B %*% A
}
}
return(TH3)
}
#' @export
#' @title Forecast error variance decomposition
#' @param obj S3 object generated by bvar, tvar or msvar
#' @param h forecast horizon
#' @param id_obj S3 object containing information about identifying the VAR model
#' @param type Type of forecast error variance decomposition. If type = "general" the generalized forecast error variance decomposition is used. For any other values the standard method. While the generalized fevd doesn't use an indefication method to identify structural shocks in the VAR model, this is not true for the standard method. Hence for the generalized fevd there is no need to supply the function via id_obj. If id_obj is not null, there will be a warning message but will still proceed.
#' @param normalize whether the fevd should be normalized or not
#' @param ... currently not used
#' @rdname fevd
fevd <- function(obj,h=12,id_obj=NULL,type="general",normalize = TRUE,...) UseMethod("fevd")
#' @export
#' @rdname fevd
fevd.bvar <- function(obj,h=12,id_obj=NULL,type="general",normalize = TRUE,...){
if(type == "general" && !is.null(id_obj)){
message("Ignoring identification method for generalized forecast error variance decomposition.\n")
}
Alpha <- obj$mcmc_draws$Alpha
Sigma <- obj$mcmc_draws$Sigma
nreps <- dim(Alpha)[3]
nVar <- dim(Alpha)[2]
intercept <- obj$general_info$intercept
nolags <- obj$general_info$nolags
fevd_return <- array(0,dim=c(nVar,nVar,1,nreps))
if(type == "general"){
for(ii in 1:nreps){
current_Alpha <- Alpha[,,ii]
current_Sigma <- Sigma[,,ii]
fevd_return[,,1,ii] <- gfevd(Alpha = current_Alpha,Sigma = current_Sigma, h = h, nolags = nolags, intercept = intercept)#deltatemp
}
}
else if(type == "standard"){
for(ii in 1:nreps){
current_Alpha <- Alpha[,,ii]
current_Sigma <- Sigma[,,ii]
fevd_return[,,1,ii] <- fevd1(Alpha = current_Alpha, Sigma = current_Sigma,h=h,nolags=nolags,intercept=intercept,id_obj=id_obj)
}
}
else{
stop("The type of the forecast error variance decomposition can be 'general' or 'standard'")
}
fevd_ret <- 0
for(ii in 1:nreps){
fevd_ret <- fevd_ret + fevd_return[,,1,ii]
}
fevd_ret <- fevd_ret / nreps
if(normalize){
t1 <- rowSums(fevd_ret)
for(ii in 1:nVar){
fevd_ret[ii,] <- fevd_ret[ii,] / t1[ii]
}
}
retlist = structure(list(fevd = fevd_ret,noregimes = 1),class="bvfevd")
return(retlist)
}
#' @export
#' @rdname fevd
fevd.tvar <- function(obj,h=12,id_obj=NULL,type="general",normalize = TRUE,...){
Alpha <- obj$mcmc_draws$Alpha
Sigma <- obj$mcmc_draws$Sigma
nreps <- dim(Alpha)[3]
nvar <- dim(Alpha)[2]
intercept <- obj$general_info$intercept
nolags <- obj$general_info$nolags
fevd_return <- array(0,dim=c(nvar,nvar,2,nreps))
if(type == "general"){
for(ii in 1:nreps){
current_Alpha <- Alpha[,,1,ii]
current_Sigma <- Sigma[,,1,ii]
fevd_return[,,1,ii] <- gfevd(Alpha = current_Alpha,Sigma = current_Sigma, h = h,
nolags = nolags, intercept = intercept)
current_Alpha <- Alpha[,,2,ii]
current_Sigma <- Sigma[,,2,ii]
fevd_return[,,2,ii] <- gfevd(Alpha = current_Alpha,Sigma = current_Sigma, h = h,
nolags = nolags, intercept = intercept)
}
}
else if(type == "standard"){
for(ii in 1:nreps){
current_Alpha <- Alpha[,,1,ii]
current_Sigma <- Sigma[,,1,ii]
fevd_return[,,1,ii] <- fevd1(Alpha = current_Alpha,Sigma = current_Sigma, h = h,
nolags = nolags, intercept = intercept)
current_Alpha <- Alpha[,,2,ii]
current_Sigma <- Sigma[,,2,ii]
fevd_return[,,2,ii] <- fevd1(Alpha = current_Alpha,Sigma = current_Sigma, h = h,
nolags = nolags, intercept = intercept)
}
}
else{
stop("Currently only generalized Forecast Error Variance Decomposition is implemented")
}
fevd_ret <- array(0,dim=c(nvar,nvar,2))
for(ii in 1:nreps){
fevd_ret[,,1] <- fevd_ret[,,1] + fevd_return[,,1,ii]
fevd_ret[,,2] <- fevd_ret[,,2] + fevd_return[,,2,ii]
}
fevd_ret <- fevd_ret / nreps
if(normalize){
t1 <- rowSums(fevd_ret[,,1])
t2 <- rowSums(fevd_ret[,,2])
for(ii in 1:nvar){
fevd_ret[ii,,1] <- fevd_ret[ii,,1] / t1[ii]
fevd_ret[ii,,2] <- fevd_ret[ii,,2] / t2[ii]
}
}
retlist = structure(list(fevd = fevd_ret,noregimes = 2),class="bvfevd")
return(retlist)
}
#' @export
#' @rdname fevd
fevd.msvar <- function(obj,h=12,id_obj=NULL,type="general",normalize = TRUE,...){
Alpha <- obj$mcmc_draws$Alpha
Sigma <- obj$mcmc_draws$Sigma
nreps <- dim(Alpha)[3]
nvar <- dim(Alpha)[2]
intercept <- obj$general_info$intercept
nolags <- obj$general_info$nolags
noregimes <- obj$general_info$noregimes
fevd_return <- array(0,dim=c(nvar,nvar,noregimes,nreps))
if(type == "general"){
for(ii in 1:nreps){
for(jj in 1:noregimes){
current_Alpha <- Alpha[,,jj,ii]
current_Sigma <- Sigma[,,jj,ii]
fevd_return[,,jj,ii] <- gfevd(Alpha = current_Alpha,Sigma = current_Sigma, h = h,
nolags = nolags, intercept = intercept)
}
}
}
else if(type == "standard"){
for(ii in 1:nreps){
for(jj in 1:noregimes){
current_Alpha <- Alpha[,,jj,ii]
current_Sigma <- Sigma[,,jj,ii]
fevd_return[,,jj,ii] <- fevd1(Alpha = current_Alpha,Sigma = current_Sigma, h = h,
nolags = nolags, intercept = intercept)
}
}
}
else{
stop("Currently only generalized Forecast Error Variance Decomposition is implemented")
}
fevd_ret <- array(0,dim=c(nvar,nvar,noregimes))
for(ii in 1:nreps){
for(jj in 1:noregimes){
fevd_ret[,,jj] <- fevd_ret[,,jj] + fevd_return[,,jj,ii]
}
}
fevd_ret <- fevd_ret / nreps
if(normalize){
for(jj in 1:noregimes){
t1 <- rowSums(fevd_ret[,,jj])
for(ii in 1:nvar){
fevd_ret[ii,,jj] <- fevd_ret[ii,,jj] / t1[ii]
}
}
}
retlist = structure(list(fevd = fevd_ret,noregimes = noregimes),class="bvfevd")
return(retlist)
}
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.