version 0.3.2

cran · Oct 11, 2008 · 70c00d1 · 70c00d1
1 parent 66670b1
commit 70c00d1
Show file tree

Hide file tree

Showing 35 changed files with 140 additions and 220 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,10 +1,10 @@
 Package: MSBVAR
-Version: 0.3.1
-Date: 2007-10-23
+Version: 0.3.2
+Date: 2008-10-11
 Title: Markov-Switching Bayesian Vector Autoregression Models
-Author: Patrick T. Brandt <pbrandt@utdallas.edu>, Justin Appleby <appleby@utdallas.edu>
+Author: Patrick T. Brandt <pbrandt@utdallas.edu>
 Maintainer: Patrick T. Brandt <pbrandt@utdallas.edu>
-Depends: R (>= 2.4.0), KernSmooth, xtable, coda
+Depends: R (>= 2.7.0), KernSmooth, xtable, coda
 Description: Provides methods for estimating frequentist and 
   Bayesian Vector Autoregression (VAR) models.  Functions for reduced
   form and structural VAR models are also available. Includes 
@@ -17,4 +17,4 @@ Description: Provides methods for estimating frequentist and
 LazyLoad: yes
 License: GPL version 2 or newer
 URL: http://www.utdallas.edu/~pbrandt/
-Packaged: Tue Oct 23 20:33:27 2007; brandt
+Packaged: Fri Oct 10 16:44:53 2008; ptb
diff --git a/INDEX b/INDEX
@@ -50,8 +50,6 @@ 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

diff --git a/R/forecast.R b/R/forecast.R
@@ -126,14 +126,7 @@ function(varobj, nsteps, A0=t(chol(varobj$mean.S)),
     num.exog <- varobj$num.exog
     qm <- varobj$qm
     ncoef <- nrow(varobj$Bhat)
-#     # build the extended regressions if necessary
-#     if(is.na(sum(exog.fut))==F)
-#       {
-#         z.fut <- rbind(z,exog.fut)
-#       }
-#     else
-#       { z.fut <- NULL
-#       }
+
   # Get some constants we are going to need
     starttime<-date()           # Starting time for simulation
 
@@ -183,17 +176,19 @@ function(varobj, nsteps, A0=t(chol(varobj$mean.S)),
       # and then building the lags.  This reuses the code for the lag
       # construction in the szbvar code
 
-      X.update <- matrix(0, nsteps, ncoef)
-      X.update[,ncoef] <- matrix(1, nsteps, 1)
+      # Now build rhs -- start with an empty matrix
+      X.update <- matrix(0, nsteps, m*p+1)
+        # Add in the constant
+      X.update[, (m*p+1)] <- matrix(1, nsteps, 1)
+        # Put in the lagged y's
 
       # Note that constant is put last here when the lags are made
       for(j in 1:p)
         {
           X.update[1:nsteps,(m*(j-1)+1):(m*j)] <- matrix(Y.update[(p+1-j):(nsteps+p-j),],ncol=m)
         }
 
-      # Put on the exogenous regressors and make the constant the
-      # first exog regressor after the AR coefs
+      # Put in exogenous coefficients if there are any.
       if(is.null(exog)==F)
         {
           X.update<-cbind(X.update,exog);

diff --git a/R/mc.irf.R b/R/mc.irf.R
@@ -53,15 +53,21 @@
                             probs=c(0.16,0.84), varnames=NULL,...)
 {
     if(inherits(x, "mc.irf.VAR"))
-    { plot.mc.irf.VAR(x, method, component, probs, varnames) }
+    { plot.mc.irf.VAR(x, method=method, component=component,
+                            probs=probs, varnames=varnames) }
+
     if(inherits(x, "mc.irf.BVAR"))
-    { plot.mc.irf.BVAR(x, method, component, probs, varnames) }
+    { plot.mc.irf.BVAR(x, method=method, component=component,
+                       probs=probs, varnames=varnames) }
+
     if(inherits(x, "mc.irf.BSVAR"))
-    { plot.mc.irf.BSVAR(x, method, component, probs, varnames) }
+    { plot.mc.irf.BSVAR(x, method=method, component=component,
+                        probs=probs, varnames=varnames) }
+
 }
 
-"plot.mc.irf.VAR" <- function(x, method=c("Sims-Zha2"), component=1,
-                            probs=c(0.16,0.84), varnames=NULL,...)
+"plot.mc.irf.VAR" <- function(x, method, component,
+                              probs, varnames=NULL,...)
 { mc.impulse <- x
   m <- sqrt(dim(mc.impulse)[3])
   nsteps <- dim(mc.impulse)[2]
@@ -214,16 +220,17 @@
   invisible(list(responses=irf.ci, eigenvector.fractions=eigen.sum))
 }
 
-"plot.mc.irf.BVAR" <- function(x, method=c("Sims-Zha2"), component=1, probs=c(0.16,0.84),varnames=NULL, ...)
+"plot.mc.irf.BVAR" <- function(x, method, component,
+                               probs, varnames, ...)
 {
-    plot.mc.irf.VAR(x, method=c("Sims-Zha2"), component=1,
-                    probs=c(0.16,0.84), varnames=varnames, ...)
+    plot.mc.irf.VAR(x, method, component,
+                    probs, varnames=varnames, ...)
 }
 
 # This version uses a median, not a mean for the central tendency for
 # the percentile method.  It is based on the plot.mc.irf code.
 
-"plot.mc.irf.BSVAR" <- function(x, method = c("Sims-Zha2"), component = 1, probs = c(0.16, 0.84), varnames = NULL, ...)
+"plot.mc.irf.BSVAR" <- function(x, method, component, probs, varnames, ...)
 {
     m <- sqrt(dim(x)[3])
     nsteps <- dim(x)[2]
@@ -350,112 +357,3 @@
     invisible(list(responses = irf.ci, eigenvector.fractions = eigen.sum))
 }
 
-
-## "mc.irf.BSVAR.DEPRECATED" <- function(varobj, nsteps, A0.posterior)
-## { m<-dim(varobj$ar.coefs)[1]  # Capture the number of variablesbcoefs <- varobj$Bhat
-##   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))
-
-##   # Loop for the impulse responses
-##   for(i in 1:N2)
-##     {
-##       # Set up the A0 for this iteration
-##       A0 <- A0.get(A0.posterior$A0.posterior, i)
-##       A0inv <- solve(A0)
-
-##       bj <- a2b(A0, varobj$Ui)
-
-##       F.draw <- matrix(0, ncoef, m)
-
-##       for(j in 1:m)
-##         { btmp <- bj[(n0cum[j]+1):(n0cum[(j+1)])]
-##           F.draw[,j] <- varobj$P.posterior[[j]]%*%(btmp)
-##         }
-
-##       F.draw <- F.draw + XXinv%*%matrix(rnorm(m*ncoef), ncoef, m)
-##       B.draw <- F.draw%*%(A0inv)
-##       B.draw <- B.draw[1:(m*p),]
-##       dim(B.draw) <- c(m^2*p, 1)
-
-##       # Now compute the irfs
-##       impulse[i,] <- t(irf.var.from.beta(t(A0inv),
-##                                          B.draw,
-##                                          nsteps))
-##       if (i%%1000==0)
-##         { cat("Monte Carlo IRF Iteration = ", i, "\n"); }
-##     }
-
-##   # Put the results into an array of the impulses.
-##   impulse <- array(impulse,c(N2,nsteps,m^2))
-##   attr(impulse, "class") <- c("mc.irf.BSVAR")
-##   return(impulse)
-## }
-
-## "mc.irf.VAR.DEPRECATED" <- function(varobj, nsteps, draws)
-## { m<-dim(varobj$ar.coefs)[1]  # Capture the number of variablesbcoefs <- varobj$Bhat
-##   p<-dim(varobj$ar.coefs)[3]    # Capture the number of lags
-
-##   bcoefs <- t(varobj$Bhat[1:(m*p),])
-##   dim(bcoefs) <- c((m^2*p),1)
-##   capT <- nrow(varobj$Y)
-##   X <- varobj$X[,1:(m*p)]
-##   XXinv <- solve(crossprod(X))
-
-##   # Get the correct moment matrix for the BVAR object
-##   if(is.null(varobj$prior)==FALSE)
-##     { XXinv <- solve(varobj$hstar[1:(m*p),1:(m*p)]) }
-
-##   # Cholesky of the covariance
-##   Sigmat <- chol(varobj$mean.S)
-
-##   # Matrices to hold stuff
-##   impulse <- matrix(0,nrow=draws,ncol=(m^2*nsteps))
-
-##   # Do all the Wishart draws
-##   # DF from Box and Tiao Section 8.5.1 p. 460. and Sims and Zha 1998.
-##   df <- capT - m*p - m - 1
-##   wisharts <- rwishart(draws, df, diag(m))
-##   XXinv <- t(chol(XXinv))
-
-##   # Main loop -- uses antithetic acceleration to cut the number of
-##   # effective draws in half.
-
-##   for(i in 1:draws)
-##   {
-##       # Generate the draws from the Wishart and the Beta
-##       Sigma.Draw <- t(Sigmat)%*%(capT*solve(matrix(wisharts[,,i],m,m)))%*%Sigmat
-##       sqrtwish <- t(chol(Sigma.Draw))
-
-##       # Here we exploit the fact that since there is a Kronecker
-##       # product structure to the covariance mtx of beta.  This means
-##       # we can take the Cholesky of each component in the Kronecker
-##       # product of the VCV matrix of Beta and combine them with random
-##       # normal draws to get a draw from beta.  This is much much easier
-##       # than using a draw from a full MVN pdf
-
-##       bcoefs.covar <- kronecker(sqrtwish, XXinv)
-##       coef.u <- bcoefs.covar%*%matrix(rnorm(m^2*p), ncol=1)
-##       coef.mtx.odd <- bcoefs + coef.u
-##       coef.mtx.even <- bcoefs - coef.u
-
-##       # Now compute the irfs
-##       impulse[i,] <- t(irf.var.from.beta(sqrtwish, coef.mtx.odd, nsteps))
-##       if (i%%1000==0)
-##       { cat("Monte Carlo IRF Iteration = ", i, "\n"); }
-##   }
-
-##   # Put the results into an array of the impulses.
-##   output <- array(impulse,c(draws,nsteps,m^2))
-##   attr(output, "class") <- c("mc.irf.VAR")
-##   return(output)
-## }
-
diff --git a/R/sanity.R b/R/sanity.R
@@ -3,8 +3,9 @@
     if(args$lambda0<=0 || args$lambda0>1){
         stop("\n\n\t-- Overall prior tightness (lambda0) must be between 0 and 1\n")
     }
-    else if(args$lambda1<0 || args$lambda1>1){
-        stop("\n\n\t-- Prior tightness around AR(1) parameters (lambda1) must be between 0 and 1\n")
+    else if(args$lambda1<0){
+        stop("\n\n\t--
+        Prior tightness around AR(1) parameters (lambda1) must be greater than 0\n")
     }
     else if(args$lambda3<0){
         stop("\n\n\t-- Order of lag decay (lambda3) must be greater than 0\n")
@@ -13,7 +14,8 @@
         stop("\n\n\t-- Prior tightness around intercept (lambda4) must be greater than 0\n")
     }
     else if(args$lambda5<0){
-        stop("\n\n\t-- Prior tightness around exogenous parameters (lambda5) must be between 0 and 1\n")
+        stop("\n\n\t--
+        Prior tightness around exogenous parameters (lambda5) must be greater than 0\n")
     }
     else if(args$mu5<0){
         stop("\n\n\t-- Prior weight of sum of coefficients must be non-negative\n")
@@ -44,10 +46,10 @@
         m <- ncol(Y)
     }
 
-    # Check time series properties of input z
-    if(is.null(args$z)==FALSE && is.null(tsp(args$z))==TRUE){
-        stop("\n\n\t-- Input data [z] must be NULL, a ts() object, or an mts() object\n")
-    }
+##     # Check time series properties of input z
+##     if(is.null(args$z)==FALSE || is.ts(args$z)==FALSE){
+##         stop("\n\n\t-- Input data [z] must be NULL, a ts() object, or an mts() object\n")
+##     }
 
     # Check if num lags is possible given the amount of data supplied
     if(args$p<=0 || args$p>=T){
@@ -104,7 +106,11 @@
 
     m <- ncol(args[[1]])
 
-    # Check identification
+    # Check identification MATRIX
+
+    if(is.matrix(args$ident)==FALSE) {
+        stop("\n\n\t-- Check Identification: ident argument should be of the class matrix, not list or dataframe\n")
+    }
 
     if(length(which(args$ident==1))>(m*(m+1))/2){
         stop("\n\n\t-- Check Identification: no more than m*(m+1)/2 free parameters allowed\n")
@@ -126,7 +132,7 @@
     methodlist <- c("DistanceMLA", "DistanceMLAhat", "Euclidean", "PositiveDiagA", "PositiveDiagAinv")
 
     if(length(which(methodlist==args$normalization))==0){
-        warning("\n\n\t-- Specified normalization method does not exist: defalut set to 'DistanceMLA'\n")
+        warning("\n\n\t-- Specified normalization method does not exist: default set to 'DistanceMLA'\n")
         return("DistanceMLA")
     }
     else { return(1) }
@@ -151,11 +157,13 @@
     if(args$nsteps<=0) stop("\n\n\t-- 'nsteps' must be greater than 0\n")
 
     if(attr(args$varobj,"class")=="VAR" && !(args$nsteps>0)){
-        stop("\n\n\t-- For VAR models, argument 'nsteps' must be defined and non-negative\n")
+        stop("\n\n\t--
+        For VAR models, argument 'nsteps' must be defined and non-negative integer\n")
     }
 
     if(attr(args$varobj,"class")=="BVAR" && !(args$draws>0)){
-        stop("\n\n\t-- For BVAR models, argument 'draws' must be defined and non-negative\n")
+        stop("\n\n\t--
+        For BVAR models, argument 'draws' must be defined and non-negative integer\n")
     }
 
     if(attr(args$varobj,"class")=="BSVAR" && is.null(args$A0.posterior)){

diff --git a/R/szbsvar.R b/R/szbsvar.R
@@ -321,6 +321,13 @@ function(Y, p, z=NULL, lambda0, lambda1, lambda3, lambda4, lambda5,
       # compute the structural innovations
       structural.innovations <- Y1%*%A0.mode - X1%*%F.posterior
 
+      # reduced form exogenous coefficients
+      if(nexog==0){
+          exog.coefs <- NA
+      } else {
+          exog.coefs <- B.posterior[((m*p)+2):nrow(B.posterior),]
+      }
+
       # Now build an output list / object for the B-SVAR model
       output <- list(XX=XX,                               # data matrix moments with dummy obs
                      XY=XY,
@@ -341,11 +348,12 @@ function(Y, p, z=NULL, lambda0, lambda1, lambda3, lambda4, lambda5,
                      B.posterior=B.posterior,
                      ar.coefs=AR.coefs.posterior,
                      intercept=F.posterior[(m*p+1),],
-                     prior=c(lambda0,lambda1,lambda3,lambda4,lambda5,
-                       mu5,mu6),
+                     exog.coefs=exog.coefs,
+                     prior=c(lambda0,lambda1,lambda3,lambda4,lambda5,mu5,mu6),
                      df=capT,
                      n0=n0,
-                     ident=ident
+                     ident=ident,
+                     b=max.obj$par
                      )
       class(output) <- c("BSVAR")
       return(output)

diff --git a/R/zzz.R b/R/zzz.R
@@ -1,8 +1,8 @@
 .onLoad <- function(...) {
-   cat("##\n## MSBVAR Package v.0.3.1\n")
-   cat("## 2007-10-23\n")
-   cat("## Copyright (C) 2005-2007 Patrick T. Brandt\n")
-   cat("## Written by Patrick T. Brandt and Justin Appleby\n")
+   cat("##\n## MSBVAR Package v.0.3.2\n")
+   cat("## 2008-10-11\n")
+   cat("## Copyright (C) 2005-2008 Patrick T. Brandt\n")
+   cat("## Written by Patrick T. Brandt with help from Justin Appleby\n")
    cat("##\n## Support provided by the U.S. National Science Foundation\n")
    cat("## (Grants SES-0351179, SES-0351205, and SES-0540816)\n##\n")
    require(KernSmooth, quietly=TRUE)