Nothing
# Generate the impulse responses -- requires a posterior sample of the
# A0 to be drawn already. Could probably put an option / overload in here that
# if there are no A0 supplied that it draws them.
"mc.irf" <- function(varobj, nsteps, draws=1000, A0.posterior=NULL, sign.list=rep(1,ncol(varobj$Y)))
{
if(inherits(varobj, "VAR")){
return(mc.irf.VAR(varobj=varobj, nsteps=nsteps, draws=draws))
}
if(inherits(varobj, "BVAR")){
return(mc.irf.BVAR(varobj=varobj, nsteps=nsteps, draws=draws))
}
if(inherits(varobj, "BSVAR")){
return(mc.irf.BSVAR(varobj=varobj, nsteps=nsteps,
A0.posterior=A0.posterior, sign.list=sign.list))
}
if(inherits(varobj, "MSBVAR")){
return(mc.irf.MSBVAR(varobj=varobj, nsteps=nsteps,
draws=length(varobj$ss.sample)))
}
}
"mc.irf.VAR" <- function(varobj, nsteps, draws)
{ output <- .Call("mc.irf.var.cpp", varobj, as.integer(nsteps),
as.integer(draws))
attr(output, "class") <- c("mc.irf", "mc.irf.VAR")
attr(output, "eqnames") <- attr(varobj, "eqnames")
return(output)
}
"mc.irf.BVAR" <- function(varobj, nsteps, draws)
{
output <- mc.irf.VAR(varobj, nsteps, draws)
attr(output, "class") <- c("mc.irf", "mc.irf.BVAR")
attr(output, "eqnames") <- attr(varobj, "eqnames")
return(output)
}
"mc.irf.BSVAR" <- function(varobj, nsteps, A0.posterior, sign.list)
{ m<-dim(varobj$ar.coefs)[1] # Capture the number of variablesbcoefs <- varobj$Bhat
# Check length of sign.list
if(length(sign.list)!=m) {
stop("sign.list has wrong number of elements. Must be m")
}
p<-dim(varobj$ar.coefs)[3] # Capture the number of lags
ncoef <- dim(varobj$B.posterior)[1]
n0 <- varobj$n0
n0cum <- c(0,cumsum(n0))
N2 <- A0.posterior$N2
# Get the covar for the coefficients
XXinv <- chol(solve(varobj$Hpinv.posterior[[1]]))
# storage for the impulses and the sampled coefficients.
impulse <- matrix(0,nrow=N2, ncol=(m^2*nsteps))
output <- .Call("mc.irf.bsvar.cpp", A0.posterior$A0.posterior,
as.integer(nsteps), as.integer(N2), as.integer(m),
as.integer(p), as.integer(ncoef), as.integer(n0),
as.integer(n0cum), XXinv, varobj$Ui, varobj$P.posterior,
sign.list)
attr(output, "class") <- c("mc.irf", "mc.irf.BSVAR")
attr(output, "eqnames") <- attr(varobj, "eqnames")
return(output)
}
######################################
# Helper functions for mc.irf.MSBVAR
######################################
# Get an A0 -- function to get the reduced form VAR A(0) to capture
# the contemp effects.
getA0 <- function(z, m, h)
{
tmp <- apply(matrix(z, ncol=h), 2, xpnd)
tmp <- array(tmp, c(m,m,h))
tmp <- array(apply(tmp, 3, chol), c(m,m,h))
aperm(tmp, c(2,1,3))
}
"irf.var.from.beta" <-
function(A0,bvec,nsteps)
{ m <- ncol(A0)
p <- (length(bvec))/m^2;
bmtx <- matrix(bvec,ncol=m)
ar.coefs<-t(bmtx) # extract the ar coefficients
dim(ar.coefs)<-c(m,m,p) # push ar coefs into M x M x P array
ar.coefs<-aperm(ar.coefs,c(2,1,3)) # reorder array so columns
# are for eqn
impulses <- irf.VAR(list(ar.coefs=ar.coefs),nsteps,A0=A0)$mhat
impulses <- aperm(impulses, c(3,1,2)) # flips around the responses
# to stack them as series
# for each response, so the
# first nstep elements are
# the responses to the first
# shock, the second are the
# responses of the first
# variable to the next
# shock, etc.
dim(impulses)<-c((m^2)*nsteps,1)
return(impulses)
}
"mc.irf.MSBVAR" <- function(varobj, nsteps, draws)
{
# Get constants
m <- varobj$m
p <- varobj$p
h <- varobj$h
N2 <- draws
# Get the A0 or structural shocks
Sigmavec <- varobj$Sigma.sample
cat("Running setup tasks to sort regime parameters.\n")
# Set them up in regime specific arrays
A0s <- array(sapply(1:N2, function(i) {getA0(Sigmavec[i,], m, h)}),
c(m, m, h, N2))
# Now set up the AR coefficients vector. Need to pluck out the
# intercepts from the posterior and assemble the AR coefs in a vector
# for each regime
# Intercept indices
mpplus1 <- m*p + 1
nc <- m*mpplus1
ii <- seq(mpplus1, by=mpplus1, length=m) # Intercept indices
# Split Beta by regimes
Beta <- array(varobj$Beta.sample, c(N2, nc, h))
# Remove the intercepts -- just the AR part
AR <- Beta[,-ii,]
rm(Beta)
cat("Computing long run regime probabilities.\n")
# Get summaries of the long run ergodic probabiities
lrQ <- sapply(1:N2,
function(i) {
steady.Q(matrix(varobj$Q.sample[i,],h,h))})
# Storage for the impulse responses: draws, steps, response-shock
# combinations, regimes array
tmp <- array(0, c(N2, nsteps, m^2, h))
outputavg <- array(0, c(N2, nsteps, m^2))
# Loop over the draws and regimes to compute the IRF for each
# regime and in equilibrium
# Computing the IRF for each regime
cat("Beginning to compute the IRFs for each regime, and averaged over the regime probabilities.\n")
for(i in 1:N2)
{
for(j in 1:h)
{
tmp[i,,,j] <- irf.var.from.beta(A0s[,,j,i], AR[i,,j], nsteps)
# Now take long run averages and store them
outputavg[i,,] <- tmp[i,,,j] + tmp[i,,,j]*lrQ[j,i]
}
if(i%%1000==0) cat("Finished ", 100*i/N2, " percent\n", sep="")
}
output <- list(shortrun=tmp, longrun=outputavg)
attr(output, "class") <- c("mc.irf", "mc.irf.MSBVAR")
attr(output, "eqnames") <- attr(varobj, "eqnames")
return(output)
}
######################################################################
# Main classes plotting function for impulse responses computed with
# mc.irf and it classes cousins.
######################################################################
"plot.mc.irf" <- function(x, method=c("Sims-Zha2"), component=1,
probs=c(0.16,0.84),
varnames=attr(x, "eqnames"),
regimelabels=NULL, ask=TRUE, ...)
{
# VAR
if(inherits(x, "mc.irf.VAR"))
{ tmp <- plot.mc.irf.VAR(x, method=method, component=component,
probs=probs, varnames=varnames,...) }
# BVAR
if(inherits(x, "mc.irf.BVAR"))
{ tmp <- plot.mc.irf.BVAR(x, method=method, component=component,
probs=probs, varnames=varnames,...) }
# BSVAR
if(inherits(x, "mc.irf.BSVAR"))
{ tmp <- plot.mc.irf.BSVAR(x, method=method, component=component,
probs=probs, varnames=varnames,...) }
# MSBVAR
if(inherits(x, "mc.irf.MSBVAR"))
{ tmp <- plot.mc.irf.MSBVAR(x, method=method,
component=component, probs=probs,
varnames=varnames,
regimelabels=regimelabels, ask,...) }
# Get the eigen-fractions, if relevant and available.
return(invisible(tmp))
}
# Use a single function to compute the IRF across each kind of VAR.
# Then the computation of the IRF can be separated from the plotting
# functions for more user control. This is hidden from the user, but
# called inside of the various functions needed to draw the IRF plots
# with error bands.
"compute.plot.mc.irf" <- function(x, method, component, probs)
{
mc.impulse <- x
m <- sqrt(dim(mc.impulse)[3])
nsteps <- dim(mc.impulse)[2]
draws <- dim(mc.impulse)[1]
# Storage
irf.ci <- array(0,c(nsteps,length(probs)+1,m^2))
# Compute the IRF confidence intervals for each element of MxM
# array
# There are multiple methods for doing this.
# Monte Carlo / Bootstrap percentile method
if (method=="Percentile")
{
eigen.sum <- 0
for(i in 1:m^2)
{
irf.bands <- t(apply(mc.impulse[,,i], 2, quantile, probs))
irf.mean <- apply(mc.impulse[,,i], 2, mean)
irf.ci[,,i] <- cbind(irf.bands, irf.mean)
}
}
if (method=="Normal Approximation")
{
eigen.sum <- 0
for (i in 1:m^2)
{
irf.mean <- apply(mc.impulse[,,i], 2, mean)
irf.var <- apply(mc.impulse[,,i], 2, var)
irf.bands <- irf.mean + matrix(rep(qnorm(probs), each=nsteps), nrow=nsteps)*irf.var
irf.ci[,,i] <- cbind(irf.bands, irf.mean)
}
}
# Sims and Zha symmetric eigen decomposition (assumes normality
# approximation) with no accounting for the correlation across the
# responses. Method does account for the correlation over time.
if (method=="Sims-Zha1")
{
eigen.sum <- matrix(0, m^2, nsteps)
for(i in 1:m^2)
{
decomp <- eigen(var(mc.impulse[,,i]), symmetric=T)
W <- decomp$vectors
lambda <- decomp$values
irf.mean <- apply(mc.impulse[,,i], 2, mean)
irf.bands <- irf.mean + W[,component]*matrix(rep(qnorm(probs), each=nsteps), nrow=nsteps)*sqrt(lambda[component])
irf.ci[,,i] <- cbind(irf.bands, irf.mean)
eigen.sum[i,] <- 100*lambda/sum(lambda)
}
}
# Sims and Zha asymmetric eigen decomposition (no normality
# assumption) with no accounting for the correlation across the
# responses. Method does account for correlation over time.
if (method=="Sims-Zha2")
{
eigen.sum <- matrix(0, m^2, nsteps)
for(i in 1:m^2)
{
decomp <- eigen(var(mc.impulse[,,i]), symmetric=T)
W <- decomp$vectors
lambda <- decomp$values
gammak <- mc.impulse[,,i]*(W[component,])
gammak.quantiles <- t(apply(gammak, 2, quantile, probs=probs))
irf.mean <- apply(mc.impulse[,,i], 2, mean)
irf.bands <- irf.mean + gammak.quantiles
irf.ci[,,i] <- cbind(irf.bands, irf.mean)
eigen.sum[i,] <- 100*lambda/sum(lambda)
}
}
# Sims and Zha asymmetric eigen decomposition with no normality
# assumption and an accounting of the temporal and cross response
# correlations.
if (method=="Sims-Zha3")
{
eigen.sum <- matrix(0, m^2, m^2*nsteps)
# Stack all responses and compute one eigen decomposition.
stacked.irf <- array(mc.impulse, c(draws, m^2*nsteps))
decomp <- eigen(var(stacked.irf), symmetric=T)
W <- decomp$vectors
lambda <- decomp$values
gammak <- stacked.irf*W[component,]
gammak.quantiles <- apply(gammak, 2, quantile, probs)
irf.mean <- matrix(apply(stacked.irf, 2, mean),
nrow=length(probs),
ncol=dim(stacked.irf)[2],
byrow=T)
irf.bands <- irf.mean + gammak.quantiles
# Reshape these....
irf.ci <- array((rbind(irf.bands,irf.mean[1,])),
c(length(probs)+1, nsteps, m^2))
irf.ci <- aperm(irf.ci, c(2, 1, 3))
eigen.sum <- 100*lambda/sum(lambda)
}
return(list(responses=irf.ci, eigenvector.fractions=eigen.sum,
m=m, nsteps=nsteps))
}
"plot.mc.irf.VAR" <- function(x, method=method, component=component,
probs=probs, varnames=varnames,...)
{
tmp <- compute.plot.mc.irf(x, method, component, probs)
# Get out the main components we need
m <- tmp$m
nsteps <- tmp$nsteps
irf.ci <- tmp$responses
eigen.sum <- tmp$eigenvector.fractions
# Compute the bounds for the plots
minmax <- matrix(0, nrow=m, ncol=2)
within.plots <- apply(irf.ci, 3, range)
tmp <- (c(1:m^2)%%m)
tmp[tmp==0] <- m
indices <- sort(tmp, index.return=T)$ix
dim(indices) <- c(m, m)
for(i in 1:m){minmax[i,] <- range(within.plots[,indices[,i]])}
# Now loop over each m columns to find the minmax for each column
# responses in the MAR plot.
j <- 1
# Plot the results
par(mfcol=c(m,m), mai=c(0.25,0.25,0.15,0.25),
omi=c(0.15,0.75,1,0.15))
for(i in 1:m^2)
{
lims <- ifelse((i-m)%%m==0, m, (i-m)%%m)
ts.plot(irf.ci[,,i],
gpars=list(xlab="",ylab="",ylim=minmax[lims,]), ...)
abline(h=0)
if(i<=m){ mtext(varnames[i], side=2, line=3)}
if((i-1)%%m==0){
mtext(varnames[j], side=3, line=2)
j <- j+1
}
}
mtext("Response in", side=2, line=3, outer=T)
mtext("Shock to", side=3, line=3, outer=T)
# Put response names on the eigenvector fractions
if(method == "Sims-Zha1" | method == "Sims-Zha2")
{
if(is.null(varnames)==T) varnames <- paste("V", seq(1:m), sep = "")
shock.name <- rep(varnames, m)
response.name <- rep(varnames, each=m)
eigen.sum <- cbind(shock.name, response.name, as.data.frame(eigen.sum))
colnames(eigen.sum) <- c("Shock","Response", paste("Component", seq(1:nsteps)))
}
if(method == "Sims-Zha3")
{ names(eigen.sum) <- c(paste("Component", seq(1:m^2*nsteps))) }
# Return
return(list(responses=irf.ci, eigenvector.fractions=eigen.sum))
}
"plot.mc.irf.BVAR" <- function(x, method=method, component=component,
probs=probs, varnames=varnames, ...)
{
plot.mc.irf.VAR(x, method, component,
probs, varnames, ...)
}
# BSVAR model IRFs
"plot.mc.irf.BSVAR" <- function(x, method=method, component=component,
probs=probs, varnames=varnames, ...)
{
m <- sqrt(dim(x)[3])
nsteps <- dim(x)[2]
draws <- dim(x)[1]
irf.ci <- array(0, c(nsteps, length(probs) + 1, m^2))
if (method == "Percentile") {
eigen.sum <- 0
for (i in 1:m^2) {
irf.bands <- t(apply(x[, , i], 2, quantile,
probs))
irf.mean <- apply(x[, , i], 2, mean)
irf.ci[, , i] <- cbind(irf.bands, irf.mean)
}
}
if (method == "Normal Approximation") {
eigen.sum <- 0
for (i in 1:m^2) {
irf.mean <- apply(x[, , i], 2, mean)
irf.var <- apply(x[, , i], 2, var)
irf.bands <- irf.mean + matrix(rep(qnorm(probs),
each = nsteps), nrow = nsteps) * irf.var
irf.ci[, , i] <- cbind(irf.bands, irf.mean)
}
}
if (method == "Sims-Zha1") {
eigen.sum <- matrix(0, m^2, nsteps)
for (i in 1:m^2) {
decomp <- eigen(var(x[, , i]), symmetric = T)
W <- decomp$vectors
lambda <- decomp$values
irf.mean <- apply(x[, , i], 2, mean)
irf.bands <- irf.mean + W[, component] * matrix(rep(qnorm(probs),
each = nsteps), nrow = nsteps) * sqrt(lambda[component])
irf.ci[, , i] <- cbind(irf.bands, irf.mean)
eigen.sum[i, ] <- 100 * lambda/sum(lambda)
}
}
if (method == "Sims-Zha2") {
eigen.sum <- matrix(0, m^2, nsteps)
for (i in 1:m^2) {
decomp <- eigen(var(x[, , i]), symmetric = T)
W <- decomp$vectors
lambda <- decomp$values
gammak <- x[, , i] * (W[component, ])
gammak.quantiles <- t(apply(gammak, 2, quantile,
probs = probs))
irf.mean <- apply(x[, , i], 2, mean)
irf.bands <- irf.mean + gammak.quantiles
irf.ci[, , i] <- cbind(irf.bands, irf.mean)
eigen.sum[i, ] <- 100 * lambda/sum(lambda)
}
}
if (method == "Sims-Zha3") {
eigen.sum <- matrix(0, m^2, m^2 * nsteps)
stacked.irf <- array(x, c(draws, m^2 * nsteps))
decomp <- eigen(var(stacked.irf), symmetric = T)
W <- decomp$vectors
lambda <- decomp$values
gammak <- stacked.irf * W[component, ]
gammak.quantiles <- apply(gammak, 2, quantile, probs)
irf.mean <- matrix(apply(stacked.irf, 2, mean), nrow = length(probs),
ncol = dim(stacked.irf)[2], byrow = T)
irf.bands <- irf.mean + gammak.quantiles
irf.ci <- array((rbind(irf.bands, irf.mean[1, ])), c(length(probs) +
1, nsteps, m^2))
irf.ci <- aperm(irf.ci, c(2, 1, 3))
eigen.sum <- 100 * lambda/sum(lambda)
}
minmax <- matrix(0, nrow = m, ncol = 2)
within.plots <- apply(irf.ci, 3, range)
tmp <- (c(1:m^2)%%m)
tmp[tmp==0] <- m
indices <- sort(tmp, index.return=T)$ix
dim(indices) <- c(m, m)
for(i in 1:m)
{ minmax[i,] <- range(within.plots[,indices[,i]])
}
# Now loop over each m columns to find the minmax for each column
# responses in the MAR plot.
j <- 1
# Plot the results
par(mfcol=c(m,m),mai=c(0.15,0.2,0.15,0.15), omi=c(0.15,0.75,1,0.15))
for(i in 1:m^2)
{
lims <- ifelse((i-m)%%m==0, m, (i-m)%%m)
ts.plot(irf.ci[,,i],
gpars=list(xlab="",ylab="",ylim=minmax[lims,]))
abline(h=0)
if(i<=m)
{ mtext(varnames[i], side=2, line=3) }
if((i-1)%%m==0)
{ mtext(varnames[j], side=3, line=2)
j <- j+1
}
}
# Add row and column labels for graph
mtext("Response in", side = 2, line = 3, outer = T)
mtext("Shock to", side = 3, line = 3, outer = T)
# Return eigendecomposition pieces if necessary
if (method == "Sims-Zha1" | method == "Sims-Zha2") {
if (is.null(varnames) == T)
varnames <- paste("V", seq(1:m), sep = "")
shock.name <- rep(varnames, m)
response.name <- rep(varnames, each = m)
eigen.sum <- cbind(shock.name, response.name, as.data.frame(eigen.sum))
colnames(eigen.sum) <- c("Shock", "Response", paste("Component",
seq(1:nsteps)))
}
if (method == "Sims-Zha3") {
names(eigen.sum) <- c(paste("Component", seq(1:m^2 *
nsteps)))
}
return(list(responses = irf.ci, eigenvector.fractions = eigen.sum))
}
# MSBVAR IRF plot
"plot.mc.irf.MSBVAR" <- function(x, method=method,
component=component, probs=probs,
varnames=varnames,
regimelabels=regimelabels, ask=ask, ...)
{
# Subset out the irfs by the number of regimes and
# variables. Recall that the dimensions of the irf array in "x"
# are N2 x nsteps x (shock*response) x no. regimes.
# Start by plotting each of the shortrun plots
d <- dim(x$shortrun)
h <- d[4]
# Storage / list for the IRF output summaries
out <- vector(mode="list", length=(h+1))
# Make regime labels if they have not been provided
if(is.null(regimelabels)) regimelabels <- paste("Regime", 1:h)
# Main loop for the plot.
for (i in 1:h)
{
# plot each short run irf, per regime, and a label as
# such.
# Store the summary / CI / eigendecomposition if it is done.
out[[i]] <- plot.mc.irf.VAR(x=x$shortrun[,,,i], method=method,
component=component,
probs=probs, varnames=varnames, ...)
# Add regime label to plot
mtext(regimelabels[i], side=3, line=4.5, outer=T)
devAskNewPage(ask=ask)
}
# Now plot the long-run of regime probability averaged IRF, with
# error bands
tmp <- plot.mc.irf.VAR(x$longrun, method=method, component=component,
probs=probs, varnames=varnames)
mtext("Regime averaged IRF", side=3, line=4.5, outer=TRUE)
out[[h+1]] <- tmp
# Add regime names to output
names(out) <- c(regimelabels, "Regime averaged")
return(out)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.