Skip to content

Commit

Permalink
version 0.2.2
Browse files Browse the repository at this point in the history
  • Loading branch information
Patrick T. Brandt authored and gaborcsardi committed Jan 14, 2007
1 parent 0961ed6 commit aeaa943
Show file tree
Hide file tree
Showing 20 changed files with 558 additions and 119 deletions.
9 changes: 5 additions & 4 deletions DESCRIPTION
Original file line number Diff line number Diff line change
@@ -1,10 +1,10 @@
Package: MSBVAR
Version: 0.2.2
Date: 2007-01-14
Title: Bayesian Vector Autoregression Models
Version: 0.2.1
Date: 2006-09-19
Author: Patrick T. Brandt <pbrandt@utdallas.edu>
Maintainer: Patrick T. Brandt <pbrandt@utdallas.edu>
Depends: R (>= 2.3.0), KernSmooth, xtable, coda
Depends: R (>= 2.4.0), KernSmooth, xtable, coda, methods
Description: Provides methods for estimating frequentist and
Bayesian Vector Autoregression (VAR) models. Functions for reduced
form and structural VAR models are also available. Includes
Expand All @@ -14,6 +14,7 @@ Description: Provides methods for estimating frequentist and
for plotting forecasts and impulse responses, and generating draws
from Wishart and singular multivariate normal densities. Future
versions will include some models with Markov switching.
LazyLoad: yes
License: GPL version 2 or newer
URL: http://www.utdallas.edu/~pbrandt/
Packaged: Wed Sep 20 08:43:41 2006; brandt
Packaged: Thu Jan 18 12:06:20 2007; brandt
4 changes: 4 additions & 0 deletions INDEX
Original file line number Diff line number Diff line change
@@ -1,4 +1,6 @@
A02mcmc Converts A0 objects to coda MCMC objects
BCF2006data Subset of Data from Brandt, Colaresi, and
Freeman (2006)
IsraelPalestineConflict
Weekly Goldstein Scaled Israeli-Palestinian
Conflict Data, 1979-2003
Expand Down Expand Up @@ -46,6 +48,8 @@ rmse Root mean squared error of a Monte Carlo / MCMC
sample of forecasts
rmultnorm Multivariate Normal Random Number Generator
rwishart Random deviates from a Wishart distribution
summary Summary functions for VAR / BVAR / B-SVAR model
objects
szbsvar Structural Sims-Zha Bayesian VAR model
estimation
szbvar Reduced form Sims-Zha Bayesian VAR model
Expand Down
15 changes: 11 additions & 4 deletions NAMESPACE
Original file line number Diff line number Diff line change
Expand Up @@ -52,18 +52,25 @@ export(dfev,
plot.mc.irf.VAR,
plot.mc.irf.BVAR,
plot.mc.irf.BSVAR,
plot.forecast.VAR,
plot.forecast,
plot.forc.ecdf,
print.posterior.fit.BVAR,
print.posterior.fit.BSVAR
)
print.posterior.fit.BSVAR,
summary.VAR,
summary.BVAR,
summary.BSVAR
)

S3Method(print, posterior.fit.BVAR)
S3Method(print, posterior.fit.BSVAR)
S3Method(print, dfev)

S3Method(summary, dfev)

S3Method(summary, VAR)
S3Method(summary, BVAR)
S3Method(summary, BSVAR)

S3Method(plot, irf.VAR)
S3Method(plot, irf.BVAR)
S3Method(plot, irf.BSVAR)
Expand All @@ -72,5 +79,5 @@ S3Method(plot, mc.irf.VAR)
S3Method(plot, mc.irf.BVAR)
S3Method(plot, mc.irf.BSVAR)

S3Method(plot, forecast.VAR)


129 changes: 87 additions & 42 deletions R/SZ.prior.evaluation.R
Original file line number Diff line number Diff line change
@@ -1,45 +1,90 @@
"SZ.prior.evaluation" <- function(Y, p, lambda0, lambda1, lambda3,
lambda4, lambda5, mu5, mu6, z=NULL, nu=ncol(Y)+1, qm, prior=0,
nstep, y.future)
{ combos <- length(lambda0)*length(lambda1)*length(lambda3)*length(lambda4)
results <- matrix(0, combos ,11)
h <- 0

for(i in 1:length(lambda0))
{ for(j in 1:length(lambda1))
{ for(k in 1:length(lambda3))
{ for(l in 1:length(lambda4))
{ fit <- szbvar(Y, p, z,
lambda0[i],
lambda1[j],
lambda3[k],
lambda4[l],
lambda5=lambda5,
mu5=mu5,
mu6=mu6,
nu=nu,
qm=qm,
prior=prior,
posterior.fit=T)

forecast <- forecast.VAR(fit, nstep)
eval.forecasts <-
cf.forecasts(forecast[(nrow(forecast)-nstep+1):nrow(forecast),], y.future)

# Now gather up all the parameters and results
tmp <- c(lambda0[i],lambda1[j],lambda3[k],lambda4[l],
lambda5, mu5, mu6, eval.forecasts[1],
eval.forecasts[2], fit$marg.llf[1],
fit$marg.post[1])

h <- h+1
results[h,] <- tmp

}}}}

colnames(results) <- c("lambda0", "lambda1", "lambda3", "lambda4",
"lambda5", "mu5", "mu6", "RMSE", "MAE", "MargLLF",
"MargPosterior")
# New version of this -- optimized for speed and more values of the
# hyperparameters

SZ.prior.evaluation <- function (Y, p, lambda0, lambda1, lambda3, lambda4, lambda5,
mu5, mu6, z = NULL, nu = ncol(Y) + 1, qm, prior = 0, nstep,
y.future)
{
combos <- length(lambda0) * length(lambda1) * length(lambda3) *
length(lambda4) * length(lambda5) * length(mu5) * length(mu6)

results <- matrix(0, combos, 11)

results[,1:7] <- as.matrix(expand.grid(lambda0=lambda0, lambda1=lambda1,
lambda3=lambda3, lambda4=lambda4,
lambda5=lambda5, mu5=mu5, mu6=mu6))

for (i in 1:nrow(results))
{
fit <- szbvar(Y, p, z, lambda0=results[i,1],
lambda1=results[i,2],
lambda3=results[i,3], lambda4=results[i,4],
lambda5 = results[i,5],
mu5 = results[i,6], mu6 = results[i,7], nu =
ncol(Y)+1, qm = qm, prior = prior,
posterior.fit = T)

forecast <- forecast.VAR(fit, nstep)
eval.forecasts <- cf.forecasts(forecast[(nrow(forecast) -
nstep + 1):nrow(forecast), ], y.future)

tmp <- c(eval.forecasts[1], eval.forecasts[2], fit$marg.llf[1], fit$marg.post[1])

results[i, 8:11] <- tmp

if(i%%100==T) { cat("Finished model", i, "of", nrow(results),"\n") }

}
colnames(results) <- c("lambda0", "lambda1", "lambda3", "lambda4",
"lambda5", "mu5", "mu6", "RMSE", "MAE", "LLF", "logMDD")
return(results)
}
}



## "SZ.prior.evaluation" <- function(Y, p, lambda0, lambda1, lambda3,
## lambda4, lambda5, mu5, mu6, z=NULL, nu=ncol(Y)+1, qm, prior=0,
## nstep, y.future)
## { combos <- length(lambda0)*length(lambda1)*length(lambda3)*length(lambda4)
## results <- matrix(0, combos ,11)
## h <- 0

## for(i in 1:length(lambda0))
## { for(j in 1:length(lambda1))
## { for(k in 1:length(lambda3))
## { for(l in 1:length(lambda4))
## { fit <- szbvar(Y, p, z,
## lambda0[i],
## lambda1[j],
## lambda3[k],
## lambda4[l],
## lambda5=lambda5,
## mu5=mu5,
## mu6=mu6,
## nu=nu,
## qm=qm,
## prior=prior,
## posterior.fit=T)

## forecast <- forecast.VAR(fit, nstep)
## eval.forecasts <-
## cf.forecasts(forecast[(nrow(forecast)-nstep+1):nrow(forecast),], y.future)

## # Now gather up all the parameters and results
## tmp <- c(lambda0[i],lambda1[j],lambda3[k],lambda4[l],
## lambda5, mu5, mu6, eval.forecasts[1],
## eval.forecasts[2], fit$marg.llf[1],
## fit$marg.post[1])

## h <- h+1
## results[h,] <- tmp

## }}}}

## colnames(results) <- c("lambda0", "lambda1", "lambda3", "lambda4",
## "lambda5", "mu5", "mu6", "RMSE", "MAE", "MargLLF",
## "MargPosterior")

## return(results)
## }
4 changes: 2 additions & 2 deletions R/dfev.R
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@
if (class(varobj)==c("VAR") || class(varobj)==c("BVAR"))
{ return(dfev.VAR(varobj, A0=t(chol(varobj$mean.S)), k))}
if (class(varobj)==c("BSVAR"))
{ return(dfev.VAR(varobj, A0=solve(varobj$A0.mode), k)) }
{ return(dfev.VAR(varobj, A0=t(solve(varobj$A0.mode)), k)) }
}

"dfev.BVAR" <- function(varobj, A0=t(chol(varobj$mean.S)), k)
Expand All @@ -13,7 +13,7 @@
return(output)
}

"dfev.BSVAR" <- function(varobj, A0=solve(varobj$A0.mode), k)
"dfev.BSVAR" <- function(varobj, A0=t(solve(varobj$A0.mode)), k)
{
output <- dfev.VAR(varobj,A0,k)
attr(output, "class") <- c("dfev")
Expand Down
69 changes: 45 additions & 24 deletions R/forecast.R
Original file line number Diff line number Diff line change
Expand Up @@ -3,27 +3,27 @@
if(inherits(varobj,"VAR")){
return(forecast.VAR(varobj, nsteps, A0=t(chol(varobj$mean.S)),
shocks=matrix(0,nrow=nsteps,ncol=dim(varobj$ar.coefs)[1]),
exog.fut=matrix(0,nrow=nsteps,ncol=nrow(exog.coefs))))
exog.fut=matrix(0,nrow=nsteps,ncol=nrow(varobj$exog.coefs))))
}

if(inherits(varobj, "BVAR")){
return(forecast.VAR(varobj, nsteps, A0=t(chol(varobj$mean.S)),
shocks=matrix(0,nrow=nsteps,ncol=dim(varobj$ar.coefs)[1]),
exog.fut=matrix(0,nrow=nsteps,ncol=nrow(exog.coefs))))
exog.fut=matrix(0,nrow=nsteps,ncol=nrow(varobj$exog.coefs))))

}

if(inherits(varobj, "BSVAR")){
return(forecast.VAR(varobj, nsteps, A0=solve(varobj$A0.mode),
return(forecast.VAR(varobj, nsteps, A0=t(solve(varobj$A0.mode)),
shocks=matrix(0,nrow=nsteps,ncol=dim(varobj$ar.coefs)[1]),
exog.fut=matrix(0,nrow=nsteps,ncol=nrow(exog.coefs))))
exog.fut=matrix(0,nrow=nsteps,ncol=nrow(varobj$exog.coefs))))
}
}

"forecast.VAR" <-
function(varobj, nsteps, A0=t(chol(varobj$mean.S)),
shocks=matrix(0,nrow=nsteps,ncol=dim(varobj$ar.coefs)[1]),
exog.fut=matrix(0,nrow=nsteps,ncol=nrow(exog.coefs)))
exog.fut=matrix(0,nrow=nsteps,ncol=nrow(varobj$exog.coefs)))
{
# Set up the initial parameters for the VAR forecast function from
# VAR object
Expand Down Expand Up @@ -62,16 +62,16 @@ function(varobj, nsteps, A0=t(chol(varobj$mean.S)),

"forecast.BVAR" <- function(varobj, nsteps, A0=t(chol(varobj$mean.S)),
shocks=matrix(0,nrow=nsteps,ncol=dim(varobj$ar.coefs)[1]),
exog.fut=matrix(0,nrow=nsteps,ncol=nrow(exog.coefs)))
exog.fut=matrix(0,nrow=nsteps,ncol=nrow(varobj$exog.coefs)))
{
output <- forecast.VAR(varobj, nsteps, A0, shocks, exog.fut)
attr(output, "class") <- c("forecast.BVAR", "mts", "ts")
return(output)
}

"forecast.BSVAR" <- function(varobj, nsteps, A0=solve(varobj$A0.mode),
"forecast.BSVAR" <- function(varobj, nsteps, A0=t(solve(varobj$A0.mode)),
shocks=matrix(0,nrow=nsteps,ncol=dim(varobj$ar.coefs)[1]),
exog.fut=matrix(0,nrow=nsteps,ncol=nrow(exog.coefs)))
exog.fut=matrix(0,nrow=nsteps,ncol=nrow(varobj$exog.coefs)))
{
output <- forecast.VAR(varobj, nsteps, A0, shocks, exog.fut)
attr(output, "class") <- c("forecast.BSVAR", "mts", "ts")
Expand Down Expand Up @@ -162,7 +162,7 @@ function(varobj, nsteps, A0=t(chol(varobj$mean.S)),
epsilon.i <- matrix(rnorm(nsteps*m),nrow=nsteps,ncol=m)

# Then construct a forecast using the innovations
ytmp <- forecast.VAR(varobj, nsteps, A0=A0, shocks=epsilon.i)
ytmp <- forecast.VAR(varobj, nsteps, A0=A0, shocks=epsilon.i, ...)

# Store draws that are past the burnin in the array
if(i>burnin)
Expand Down Expand Up @@ -507,7 +507,7 @@ function(varobj, yconst, nsteps, burnin, gibbs, exog=NULL)
return(output)
}

"plot.forecast.VAR" <-
"plot.forecast" <-
function(x,y=NULL,varnames=NULL,
start=c(0,1),
freq=1, probs=c(0.05,0.95),
Expand Down Expand Up @@ -536,34 +536,55 @@ function(x,y=NULL,varnames=NULL,
# }


par(las=1, mar=c(1,2,2.5,1))
par(las=1, mar=c(3,4,2.5,1.5), cex=1)
for(i in 1:m)
{ forc1.ci <- ts(fcast1.summary[,,i], start=start)
{ forc1.ci <- ts(fcast1.summary[,,i], start=start, freq=freq)
ylim <- c(floor(min(c(forc1.ci, compare.level[i]))),
ceiling(max(c(forc1.ci, compare.level[i]))))

if(is.null(fcasts2)==FALSE)
{ forc2.ci <- ts(fcast2.summary[,,i], start=start) }
{ forc2.ci <- ts(fcast2.summary[,,i], start=start, freq=freq)
ylim <- c(floor(min(c(forc1.ci,forc2.ci,compare.level[i]))),
ceiling(max(c(forc1.ci,forc2.ci,compare.level[i]))))

# if(is.null(fcasts3)==FALSE)
# { forc3.ci <- ts(fcast3.summary[,,i], start=start) }

ylim <- c(floor(min(c(forc1.ci,forc2.ci,compare.level[i]))),
ceiling(max(c(forc1.ci,forc2.ci,compare.level[i]))))
ts.plot(forc1.ci, forc2.ci,
gpars=list(lty=c(1,1,1,2,2,2),
ylim=ylim, xlab="",axes=FALSE, ... ))
ylim=ylim, xlab="", ylab=varnames[i], axes=FALSE, ... ))
axis(2,c(floor(min(c(forc1.ci,forc2.ci))),
ceiling(max(c(forc1.ci,forc2.ci)))))
mtext(varnames[i],side=3,line=1)

axis(1)
# mtext(varnames[i], side=2, line=1)

box();
if(i==1) { mtext(ylab, side=2, line=3, at=c(1.5*mean(ylim))) }
# if(i==1) { mtext(ylab, side=2, line=3, at=c(1.5*mean(ylim))) }
abline(h=0)

# put in the comparison level if one is provided
if (is.null(compare.level)==FALSE)
{ abline(h=compare.level[i], lty=c(2)) }

}
# par(oldpar)
}
} else {
ts.plot(forc1.ci,
gpars=list(lty=c(1,1,1),
ylim=ylim, xlab="", ylab=varnames[i], axes=FALSE, ... ))
axis(2,c(floor(min(c(forc1.ci))),
ceiling(max(c(forc1.ci)))))
axis(1)
# mtext(varnames[i], side=2, line=1)

box();
# if(i==1) { mtext(ylab, side=2, line=3, at=c(1.5*mean(ylim))) }
abline(h=0)

# put in the comparison level if one is provided
if (is.null(compare.level)==FALSE)
{ abline(h=compare.level[i], lty=c(2)) }

}

}
}

# if(is.null(fcasts3)==FALSE)
# { forc3.ci <- ts(fcast3.summary[,,i], start=start) }

0 comments on commit aeaa943

Please sign in to comment.