R/forecast.R

In MSBVAR: Markov-Switching, Bayesian, Vector Autoregression Models

# Generate forecast methods for the models in the MSBVAR package
#
# Patrick T. Brandt
#
# 20120621 : Updated to include MSBVAR forecast functions.

"forecast" <- 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(varobj$exog.coefs)),
                       N1, N2)
{
    if(inherits(varobj,"VAR")){
        return(forecast.VAR(varobj, nsteps, A0=A0,
                            shocks=shocks, exog.fut=exog.fut))
    }

    if(inherits(varobj, "BVAR")){
        return(forecast.VAR(varobj, nsteps, A0=A0,
                            shocks=shocks, exog.fut=exog.fut))
    }

    if(inherits(varobj, "BSVAR")){
        return(forecast.VAR(varobj, nsteps, A0=solve(varobj$A0.mode),
                       shocks=shocks, exog.fut=exog.fut))
    }

    if(inherits(varobj, "MSBVAR")){
        return(forecast.MSBVAR(x=varobj, k=nsteps, N1, N2))
    }
}

# This is the generic VAR forecasting function.  The other
"forecast.VAR" <-
function(varobj, nsteps, A0, shocks, exog.fut)
{
   # Set up the initial parameters for the VAR forecast function from
   #  VAR object
  y <- varobj$y
  intercept <- varobj$intercept
  ar.coefs <- varobj$ar.coefs
  exog.coefs <- varobj$exog.coefs
  m<-dim(ar.coefs)[1]
  p<-dim(ar.coefs)[3]
  capT<-nrow(y)
  yhat<-rbind(y,matrix(0,ncol=m,nrow=nsteps))

   # Compute the deterministic part of the forecasts (less the intercept!)
   if(is.na(sum(varobj$exog.coefs))==F)
     {
       deterministic.VAR <- as.matrix(exog.fut) %*% exog.coefs
     }
   else
     { deterministic.VAR <- matrix(0,nrow=nsteps,ncol=m)
     }

   # Now loop over the forecast horizon
   for(h in 1:nsteps)
     {  yhat[capT + h, ] <- (yhat[capT + h - 1,] %*% ar.coefs[,,1] +
                             intercept + deterministic.VAR[h,] + (shocks[h,]%*%A0))
       if (p>1) {for(i in 2:p)
       { yhat[capT + h, ] <- (yhat[capT + h, ] +
                              (yhat[capT + h - i, ] %*% ar.coefs[,,i]))

       }}
     }
  output <- ts(yhat, start = start(varobj$y), frequency = frequency(varobj$y), names = colnames(varobj$y))
  attr(output, "class") <- c("forecast.VAR", "mts", "ts")
  attr(output, "eqnames") <- attr(varobj, "eqnames")
  return(output)
}

"forecast.BVAR" <- function(varobj, nsteps, A0, shocks, exog.fut)
{
    output <- forecast.VAR(varobj, nsteps, A0, shocks, exog.fut)
    attr(output, "class") <- c("forecast.BVAR", "mts", "ts")
    attr(output, "eqnames") <- attr(varobj, "eqnames")
    return(output)
}

"forecast.BSVAR" <- function(varobj, nsteps, A0=solve(varobj$A0.mode), shocks, exog.fut)
{
    output <- forecast.VAR(varobj, nsteps, A0, shocks, exog.fut)
    attr(output, "class") <- c("forecast.BSVAR", "mts", "ts")
    attr(output, "eqnames") <- attr(varobj, "eqnames")
    return(output)
}

"uc.forecast" <- function(varobj, nsteps, burnin, gibbs,
                          exog=NULL)
{
    if(inherits(varobj, "VAR"))
    {
        stop("Not implemented for VAR models!\nUse a BVAR with a flat-flat prior if you want this case.\n")
##         varobj$H0 <- matrix(0, nrow(varobj$Bhat), nrow(varobj$Bhat))
##         varobj$S0 <- matrix(0, ncol(varobj$Bhat), ncol(varobj$Bhat))
##         output <- uc.forecast.VAR(varobj, nsteps, burnin, gibbs, exog)
##         attr(output, "class") <- c("forecast.VAR")
##         return(output)
    }

    if(inherits(varobj, "BVAR"))
    {
        output <- uc.forecast.VAR(varobj, nsteps, burnin, gibbs, exog)
        attr(output, "class") <- c("forecast.VAR")
        return(output)
    }

    if(inherits(varobj, "BSVAR"))
    {
      stop("Not yet implemented for BSVAR models!\n")
##       output <- uc.forecast.VAR(varobj, nsteps, burnin, gibbs,exog)
##       attr(output, "class") <- c("uc.forecast.VAR", "mts", "ts")
##       return(output)
    }
}

"uc.forecast.VAR" <- function(varobj, nsteps, burnin, gibbs, exog)
  { # Extract all the elements from the VAR object
    y <- varobj$y
    ar.coefs <- varobj$ar.coefs
    intercept <- varobj$intercept
    A0 <- t(chol(varobj$mean.S))
    X <- varobj$X         # rhs variables for the model
    Y <- varobj$Y         # lhs variables for the model
    H0 <- varobj$H0
    S0 <- varobj$S0
#    mu <- varobj$hyperp
    exog.coefs <- varobj$exog.coefs
    z <- varobj$z
#    lambda0 <- varobj$prior[1]
#    lambda1 <- varobj$prior[2]
#    lambda3 <- varobj$prior[3]
#    lambda4 <- varobj$prior[4]
#    lambda5 <- varobj$prior[5]
#    mu5 <- varobj$prior[6]
#    mu6 <- varobj$prior[7]
    nu <- varobj$prior[8]
    prior <- varobj$prior.type
    num.exog <- varobj$num.exog
    qm <- varobj$qm
    ncoef <- nrow(varobj$Bhat)

  # Get some constants we are going to need
    starttime<-date()           # Starting time for simulation

    p<-dim(ar.coefs)[3]         # Capture the number of lags
                                # from input ar coefficients

    m<-ncol(y);                 # Number of endogenous
                                # variables in the VAR
    k<-m*nsteps;              # k = mh, the maximal number
                              # of forecasts

    capT<-nrow(y)               # Number of observations we
                                # are going to use.

  # Make arrays to hold the Gibbs sampler results
  yforc<-array(0,c(gibbs,nsteps,m))

  # Do the Gibbs draws.....
  for(i in 1:(burnin+gibbs))
    { # Step (a): Compute a draw of the conditional forecasts
      # COMPUTE INNOVATIONS: These are the structural innovations

      # First draw the innovations
      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,
                           exog)

      # Store draws that are past the burnin in the array
      if(i>burnin)
        { for(j in 1:m)
            { yforc[(i-burnin),(1:nsteps),j]<-ytmp[((capT+1):(capT+nsteps)),j] }
        }


      # Step (b): Compute the mode of the posterior for the
      # forecast distribution.  This is the "extended"
      # dataset that includes the i'th Gibbs sample forecast

      # Set up the updated Y Matrix
      # this is just "ytmp" from above

      Y.update <- ytmp[(capT-p+1):nrow(ytmp),]

      # Set up the updated X -- this is hard because we need to get
      # the RHS lags correct.  We do this by padding the existing Y
      # and then building the lags.  This reuses the code for the lag
      # construction in the szbvar code

      # 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 in exogenous coefficients if there are any.
      if(is.null(exog)==F)
        {
          X.update<-cbind(X.update,exog);
        }

      # Now, stack the original Y data and the augmented data.
      Y.update <- rbind(Y, ytmp[(capT+1):nrow(ytmp),])

      # Set up crossproducts and inverses we need
      X.update <- rbind(X, X.update)
      XX.update <- crossprod(X.update)     # Cross product of RHS variables
      hstar.update <- H0 + XX.update       # Prior + Cross product

      # Updated Regression estimates, Beta | Sigma
      B.update<-solve((hstar.update),(crossprod(X.update,Y.update) + H0[,1:m]))

      # Posterior mean of Sigma | Beta
      S.update <- (S0 + crossprod(Y.update)
                   + H0[1:m,1:m] - t(B.update)%*%(hstar.update)%*%(B.update))/(capT+nsteps+nu-m-1)

      # Posterior variance of B: VBh = diag(Sh.*.(inv((H0 + x'x)))
      hstarinv <- solve(hstar.update)
      vcv.Bh <- kronecker(S.update,hstarinv)

      # Draw from the conditional posterior pdfs of the parameters

      # This is only valid for just-identified models.
      df <- capT - m*p - m - 1 + nsteps
      wisharts <- rwishart(1, df, diag(m))

      # Generate the draws from the Wishart and the Beta
      # Wishart draw
      Sigmat <- (chol(S.update))
      Sigma.Draw <- t(Sigmat)%*%(df*solve(matrix(wisharts,m,m)))%*%Sigmat
      sqrtwish <- t(chol(Sigma.Draw))
      # Covariance of beta
      bcoefs.covar <- t(chol(vcv.Bh))

      # Draw beta|Sigma
      aplus <- matrix(B.update, ncol=1) + bcoefs.covar%*%matrix(rnorm(nrow(bcoefs.covar)), ncol=1)
      aplus <- matrix(aplus, ncol=m)

      aplus.coefs<-t(aplus[1:(m*p),]);        # extract the ar coefficients
      dim(aplus.coefs)<-c(m,m,p)                    # push ar coefs into M x M x P array
      aplus.coefs<-aperm(aplus.coefs,c(2,1,3))      # reorder array


      intercept <- aplus[(m*p+1),]       # get drawn intercept....
      ar.coefs<-aplus.coefs            # AR coefs
      A0 <- sqrtwish

      if(num.exog!=0)
        {
          exog.coefs <- aplus[(m*p+2):nrow(aplus),]
        }


      # Print some intermediate results to capture progress....
      # and tell us that things are still running
      if (i%%500==0)
        { cat("Gibbs Iteration = ", i, "     \n");
          if(i<=burnin)
            { cat("(Still a burn-in draw.)\n");
            }

        }
      # Back to the top of the Gibbs loop....
    }
  endtime<-date()
  # Print time stamp so we know how long everything took.
  cat("Start time : ", starttime, "\n");
  cat("End time   : ", endtime, "\n");
    # Returns a list object
    output <- list(forecast=yforc)
#    attr(output, "class") <- c("forecast.VAR")
    return(output)
}

"hc.forecast" <- function(varobj, yconst, nsteps, burnin, gibbs, exog=NULL)
{
    if(inherits(varobj, "VAR"))
    {
        stop("Not yet implemented for VAR models!\nUse a BVAR model.")
##         output <- hc.forecast.VAR(varobj, yconst, nsteps, burnin,
##                                   gibbs, exog)
##         attr(output, "class") <- c("forecast.VAR")
##         return(output)
    }

    if(inherits(varobj, "BVAR"))
    {
        output <- hc.forecast.VAR(varobj, yconst, nsteps, burnin,
                                  gibbs, exog)
        attr(output, "class") <- c("forecast.VAR")
        return(output)
    }

    if(inherits(varobj, "BSVAR"))
    {
      stop("Not yet implemented for B-SVAR models!\n")
##         output <- hc.forecast.VAR(varobj, yconst, nsteps, burnin,
##                                   gibbs, exog)
##         attr(output, "class") <- c("hc.forecast.VAR", "mts", "ts")
##         return(output)
    }
}

"hc.forecast.VAR" <-
function(varobj, yconst, nsteps, burnin, gibbs, exog=NULL)
{
    # Extract all the elements from the VAR object that we will need
    y <- varobj$y
    ar.coefs <- varobj$ar.coefs
    intercept <- varobj$intercept
    exog.coefs <- varobj$exog.coefs
    A0 <- t(chol(varobj$mean.S))
    mu <- varobj$hyperp
    prior <- varobj$prior
    ncoef <- nrow(varobj$Bhat)
    X <- varobj$X         # rhs variables for the model
    Y <- varobj$Y         # lhs variables for the model
    H0 <- varobj$H0       # precision for the var coefs
    S0 <- varobj$S0       # precision for the Sigma
    nu <- varobj$prior[8] # df

    # Get some constants we are going to need from the inputs
    starttime<-date()           # Starting time for simulation
                                # of forecasts

    q<-nrow(as.matrix(yconst));   # Number of restrictions

    p<-dim(ar.coefs)[3]         # Capture the number of lags
                                # from input ar coefficients

    m<-ncol(y);                 # Number of endogenous
                                # variables in the VAR

    capT<-nrow(y)               # Number of observations we
                                # are going to use.

    k<-m*nsteps;                # k = mh, the maximal number

    q <- nrow(yconst)             # Number of constraints

    # Make arrays to hold the Gibbs sampler results
    yforc<-array(0,c(gibbs,nsteps,m))


    # Do the Gibbs draws.....
    for(i in 1:(burnin+gibbs))
    {
      # Step (a): Compute a draw of the conditional forecasts
      # COMPUTE INNOVATIONS: These are the structural innovations
      # Solve the constraint equation for the updated forecast errors

      # Generate the forecasts without shocks
      ytmp<-as.matrix(coef.forecast.VAR(y, intercept, ar.coefs, exog.coefs, m, p,
                                     capT, nsteps, A0=(A0)))

      # Get the impulse responses that correspond to the forecasted data
      M <- irf.VAR(varobj, nsteps, A0)$mhat

      # Construct the draw of the orthogonalized innovations that
      # satisfy the hard condition.

      # These are the q constrained innovations
      r <- (yconst - ytmp[(capT+1):(capT+nsteps),])
      r<-matrix(r[1:nsteps,1],ncol=1)

      # Build the matrix of the impulses that define the constraint
      R <- matrix(0, k, q)

      # Put the g'th column of the impulse into the constraint matrix,
      # such that R * epsilon = r
      for (g in 1:q)
      {
        if(g==1)
          { R[1:length(M[,1,1]), 1] <- M[,1,1] }
        else
          {
            R[,g] <- c(M[,1,g], R[,g-1])[1:k]
          }
      }

      # Solve the minimization problem for the mean and variance of
      # the constrained innovations.

      RRinv<-solve(crossprod(R))
      mean.epsilon <- R%*%RRinv%*%r;
      var.epsilon <- diag(1,nrow=k) - (R%*%RRinv%*%t(R));

      # Draw from the singular MVN pdf of the constrained innovations.

      epsilon.i <- matrix(rmultnorm(1, mean.epsilon, var.epsilon),
                          nrow=nsteps, ncol=m, byrow=T)

      # Add the innovations to the forecasts
      ytmp[(capT+1):(capT +nsteps),] <- ytmp[(capT+1):(capT +nsteps),]+epsilon.i%*%A0

      # Store forecasts that are past the burnin point
      if(i>burnin)
        { for(j in 1:m)
            { yforc[(i-burnin),(1:nsteps),j]<-ytmp[(capT+1):(capT+nsteps),j] }
        }


      # Step (b): Compute the mode of the posterior for the
      # conditional forecast distribution.  This is the "extended"
      # dataset that includes the i'th Gibbs sample forecast


      # Build the augmented LHS and RHS matrices
      # 1) Get the nsteps+p observations we need to build the lagged
      # endogenous variables for the augmented system.

      Y.update <- ytmp[(capT-p+1):nrow(ytmp),]

      # 2) Build the updated X -- this is hard because we need to get
      # the RHS lags correct.  We do this by padding the existing Y
      # 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)

      # 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
      if(is.null(exog)==F)
        {
          X.update<-cbind(X.update,exog);
        }

      # Now, stack the original Y data and the augmented data.
      Y.update <- rbind(Y, ytmp[(capT+1):nrow(ytmp),])

      # Set up crossproducts and inverses we need
      X.update <- rbind(X, X.update)
      XX.update <- crossprod(X.update)     # Cross product of RHS variables
      hstar.update <- H0 + XX.update       # Prior + Cross product

      # Updated Regression estimates, Beta | Sigma
      B.update<-solve((hstar.update),(crossprod(X.update,Y.update) + H0[,1:m]))

      # Posterior mean of Sigma | Beta
      S.update <- (S0 + crossprod(Y.update)
                   + H0[1:m,1:m] - t(B.update)%*%(hstar.update)%*%(B.update))/(capT+nsteps+nu-m-1)

      # Posterior variance of B: VBh = diag(Sh.*.(inv((H0 + x'x)))
      hstarinv <- solve(hstar.update)
      vcv.Bh <- kronecker(S.update,hstarinv)

      # Draw from the conditional posterior pdfs of the parameters

      # This is only valid for just-identified models.
      df <- capT - m*p - m - 1 + nsteps
      wisharts <- rwishart(1, df, diag(m))

      # Generate the draws from the Wishart and the Beta
      # Wishart draw
      Sigmat <- (chol(S.update))
      Sigma.Draw <- t(Sigmat)%*%(df*solve(matrix(wisharts,m,m)))%*%Sigmat
      sqrtwish <- t(chol(Sigma.Draw))
      # Covariance of beta
      bcoefs.covar <- t(chol(vcv.Bh))

      # Draw of beta|Sigma ~ MVN(B.update, S.Update .*. Hstarinv)
      aplus <- matrix(B.update, ncol=1) +
          bcoefs.covar%*%matrix(rnorm(nrow(bcoefs.covar)), ncol=1)

      # Reshape and extract the coefs
      aplus <- matrix(aplus, ncol=m)
      aplus.coefs<-t(aplus[1:(m*p),]);          # extract the ar coefficients
      dim(aplus.coefs)<-c(m,m,p)                # push ar coefs into M x M x P array
      aplus.coefs<-aperm(aplus.coefs,c(2,1,3))  # reorder array

      intercept <- aplus[(m*p+1),]      # get drawn intercept....
      ar.coefs<-aplus.coefs             # AR coefs
      A0 <- sqrtwish

#      exog.coefs <- aplus[(m*p+2):nrow(aplus),]

      # Need to add something here to deal with the exogenous
      # regressors!


      # Print some intermediate results to capture progress....
      # and tell us that things are still running
      if (i%%1000==0)
        { cat("Gibbs Iteration = ", i, "     \n");
          if(i<=burnin)
            { cat("(Still a burn-in draw.)\n");
            }
        }
      # Back to the top of the Gibbs loop....
    }
  endtime<-date()
  # Print time stamp so we know how long everything took.
  cat("Start time : ", starttime, "\n");
  cat("End time   : ", endtime, "\n");
    # Returns a list object
    output <- list(forecast=yforc, orig.y=y) #llf=ts(llf),hyperp=c(mu,prior)))
    attr(output, "class") <- c("forecast.VAR")
    return(output)
}


# Forecasting function for MSBVAR models.
#
# 20110113 : Initial version
# 20110628 : Clean out remaining bugs on DF corrections for
#            observations
# 20120207 : Updated to work with revisions to MSBVAR

###### MSBVARfcast #####
# Forecast for one draw -- this is an internal function that only
# returns the forecasts averaged over the regimes.
#
# Inputs :
# y = data modeled + forecast steps
# k = number of forecast steps
# h = number of regimesq
# ar.coefs = array of c(m,m,p,h) for the MSBVAR
# intercepts = array of c(m,h) for the intercepts
# Sigma = array of c(m.m,h) for the variances
# shock = array of c(k,m,h) for the shocks
# ss = matrix of c(k,h) for the regimes

MSBVARfcast <- function(y, k, h, ar.coefs, intercepts, shocks, ss)
{
    # Get constants
    dims <- dim(ar.coefs)
    m <- dims[1]
    p <- dims[3]
    capT <- nrow(y)

    # Set up the object to hold the forecasts
    # Zero out forecast periods so we can cumulate sums over time,
    # lags, and regimes to get the correct weighted forecast.

    forcs <- rbind(y, matrix(0, k, m))

    for(j in 1:k)     # Loop over forecast steps
    {
        for(i in 1:h) # Now over the regime weights, multiplying by
                      # the regime weights and aggregating
        {
            forcs[capT + j,] <- forcs[capT+j,] +
                                ss[(capT+j),i]*(forcs[capT + j - 1,] %*%
                                                ar.coefs[,,1,i] + intercepts[,i] + shocks[j,,i])
            # Now add in the other lags -- be sure they are from the
            # right state!
            if(p>1) { for(lg in 2:p)
                      { forcs[capT + j, ] <- (forcs[capT + j, ] +
                                              ss[(capT+j-lg),i]*(forcs[capT + j - lg, ] %*% ar.coefs[,,lg,i]))
                    }
                  }
        }
    }

    return(forcs[(capT+1):(capT+k),])
}

# Betai2coefs
# Converts the draw Betai into arrays of AR coefficents for each
# regime and the intercepts
Betai2coefs <- function(Betai, m, p, h)
{
    ar.coefsi <- array(0, c(m,m,p,h))
    for(i in 1:h)
    {
        tmp <- t(Betai[1:(m*p),,i])
        dim(tmp) <- c(m,m,p)
        ar.coefsi[,,,i] <- aperm(tmp, c(2,1,3))  # lags then regimes.
#        ar.coefs[,,,i] <- tmp
}
    intercepts <- Betai[(m*p+1),,]
    return(list(ar.coefsi=ar.coefsi, intercepts=intercepts))
}


updateinit <- function(Y, init.model)
{
    n<-nrow(Y);	 	# of observations in data set
    m<-ncol(Y)	# of variables in data set
    p <- init.model$p

    # Compute the number of coefficients
    ncoef<-(m*p)+1;  # AR coefficients plus intercept in each RF equation
    ndum<-m+1;                # of dummy observations
    capT<-n-p+ndum;   	      # # of observations used in mixed estimation (Was T)
    Ts<-n-p;  		      # # of actual data observations @

    # Declare the endoegnous variables as a matrix.
    dat<-as.matrix(Y);

    # Create data matrix including dummy observations
    # X: Tx(m*p+1)
    # Y: Txm, ndum=m+1 prior dummy obs (sums of coeffs and coint).
	# mean of the first p data values = initial conditions
    if(p==1)
      { datint <- as.vector(dat[1,]) }
    else
      { datint<-as.vector(apply(dat[1:p,],2,mean)) }

    # Y and X matrices with m+1 initial dummy observations
    X<-matrix(0, nrow=capT, ncol=ncoef)
    Y<-matrix(0, nrow=capT, ncol=m)
    const<-matrix(1, nrow=capT)
    const[1:m,]<-0;	         # no constant for the first m periods
    X[,ncoef]<-const;

    # Build the dummy observations we need
    for(i in 1:m)
      {
        Y[ndum,i]<-datint[i];
        Y[i,i]<-datint[i];
        for(j in 1:p)
          { X[ndum,m*(j-1)+i]<-datint[i];
            X[i,m*(j-1)+i]<-datint[i];
          }
      }

    # Note that constant is put last here when the lags are made
    for(i in 1:p)
      { X[(ndum+1):capT,(m*(i-1)+1):(m*i)]<-matrix(dat[(p+1-i):(n-i),],ncol=m)
      }

     # Put on the exogenous regressors and make the constant the
     # first exog regressor after the AR coefs
    ## if(is.null(z)==F)
    ##   {
    ##     pad.z <- matrix(0,nrow=capT,ncol=ncol(z))
    ##     pad.z[(ndum+1):capT,] <- matrix(z[(p+1):n,], ncol=ncol(z))
    ##     X<-cbind(X,pad.z);
    ##   }

    # Get the corresponding values of Y
    Y[(ndum+1):capT,]<-matrix(dat[(p+1):n,],ncol=m);

    # Weight dummy observations
    X[1:m,] <- init.model$prior[6]*X[1:m,];
    Y[1:m,] <- init.model$prior[6]*Y[1:m,];
    X[ndum,] <- init.model$prior[7]*X[ndum,];
    Y[ndum,] <- init.model$prior[7]*Y[ndum,];

    init.model$X <- X
    init.model$Y <- Y
    init.model$hstar <- (init.model$H0 + crossprod(X))

    return(init.model)
}


# This is an amalgam of uc.forecast and gibbs.msbvar.  In this
# version, the same steps as gibbs.msbvar are used, but with the
# updated data from the forecast step.

# x = posterior mode object -- can either be from gibbs.msbvar or
#     msbvar.  Need to handle identification steps or permute as
#     appropriate.  Should warn against permuting!
# k = number of forecast steps.
#
# -----------------------------------------------------------------#
################## Steps in MSBVAR forecasting #####################
# -----------------------------------------------------------------#
# A) Augment the dataset from the mode -- forecasting
# B) Update the input data
# C) For given parameters,
#    1) draw statespace
#    2) update / draw Q
#    3) update regression -- need to remake the input matrices as part
#       of this.
#    4) sample variances
#    5) sample regression
#    6) update errors
#    7) Impose identification on the regimes (if necessary)
# top of loop
# -----------------------------------------------------------------#
#
# Initial version assumes user is inputing an identified object from
# gibbs.msbvar().  Need to add a warning when they using an msbvar()
#     object and reconcile this.
#
forecast.MSBVAR <- function(x, k, N1=1000, N2=1000)
{
    # Get some constants / objects
    init.model <- x$init.model  # initial model with the prior
                                # matrices.

    # Get other constants from the input object
    h <- x$h
    m <- x$m
    p <- x$p

    # Set the sampler method for Q based on inputs
    if(attr(x, "Qsampler")=="Gibbs") Qsampler <- Q.drawGibbs
    if(attr(x, "Qsampler")=="MH") Qsampler <- Q.drawMH

    # Input data
    y <- init.model$y

    # Get sample size of original dataset and the effective sample
    TT <- capT <- nrow(init.model$y)

    # Settle how we are going to handle input sample size versus the
    # other sample sizes for the inputs / outputs.  Do this the same
    # way as the earlier unconditional forecasting code.  In this
    # code, we treated capT as the full original sample size minus the
    # lags to make it match the estimator. We then
    # manipulate from that.  We do this by always forecasting off of
    # the original data from capT+1 to capT+k,  This is the principal
    # in uc.forecast() and forecast.VAR and coef.forecast.VAR() [in
    # hidden.R].  The difference here is that tbe value of capT has to
    # be adjusted compared to what is in the other code.
    #
    # Principal in this is just working with data augmentation off of
    # the original data setup.

    Tp <- TT-p # effective sample

    alpha.prior <- x$alpha.prior

    # Get the modes for the posterior
    Q <- matrix(apply(x$Q.sample, 2, mean), h, h)
    Betai <- array(apply(x$Beta.sample, 2, mean), c(m*p+1, m, h))
    tmp <- Betai2coefs(Betai, m, p, h)
    ar.coefsi <- tmp$ar.coefsi
    intercepts <- tmp$intercepts

    Sigmai <- array(apply(matrix(apply(x$Sigma.sample, 2, mean),
                                 m*(m+1)/2, h),
                          2, xpnd), c(m,m,h))
    intercepts <- Betai[(m*p+1),,]

    # Extend the dataset / input data
    # Residuals and regression
    hreg <- hregime.reg2(h, m, p, mean.SS(x), init.model)
    e <- array(0, c(Tp+k, m, h))

    # Now take random draws from the error process for each regime to
    # pad out initial values for the forecast periods.
    Sigmachol <- aperm(array(apply(Sigmai, 3, chol), c(m,m,h)),
                       c(2,1,3))

    for(i in 1:h)
    {
        etmp <-  t(Sigmachol[,,i]%*%matrix(rnorm(k*m), m, k))
        e[,,i] <- rbind(hreg$e[,,i], etmp)
    }

    # State-space initialization -- flat convolution of the last
    # state!

    tmpSS <- mean.SS(x)
    sstmp <- rbind(mean.SS(x), matrix(rep(NA, k*h), k, h))

    for(i in 1:k)
    {
        tmp <- sstmp[Tp+i-1,]%*%Q
        sstmp[Tp+i,] <- tmp/sum(tmp)

    }

    transtmp <- count.transitions(sstmp)
    ss <- list(SS=sstmp, transitions=transtmp)

    # Make list objects for Beta and Sigma
    Betai <- list(Betai=Betai)
    Sigmai <- list(Sigmai=Sigmai)

    # Storage
    forecasts <- array(0, c(N2, k, m))
    states <- vector("list", N2)

    # Burnin loop + Final loop
    for (j in 1:(N1+N2))
    {

        # Update the residuals
        shocks <- array(rnorm(m*k*h), c(k, m, h))

        Sigmai.chol <- aperm(array(apply(Sigmai$Sigmai, 3, chol), c(m,m,h)),
                             c(2,1,3))

        for(i in 1:h)
         {
#             shocks[,,i] <- Sigmai.chol[,,i]%*%shocks[,,i]
             shocks[,,i] <- t(tcrossprod(Sigmai.chol[,,i], shocks[,,i]))
             e[(Tp+1):(Tp+k),,i] <- (shocks[,,i])
         }


        # Forecast

        fcast <- MSBVARfcast(y[(p+1):nrow(y),], k, h, ar.coefsi,
                             intercepts, shocks, ss$SS)

#        print(round(cbind(shocks[,,1], shocks[,,2], fcast), 2))

        # Update / Draw statespace

        oldtran <- matrix(0,h,h)

        while(sum(diag(oldtran)==0)>0)
        {
            ss <- SS.ffbs(e, (Tp+k), m, p, h, Sigmai$Sigmai, Q)
            oldtran <- ss$transitions
        }

        #print(ss$transitions)
        ## plot(ts(ss$SS), plot.type="s", col=1:h,
        ##      main=paste("Iteration :", j))

        #if(sum(diag(ss$transitions)==0)>0) stop("Oops.")

        # Draw Q
        Q <- Qsampler(Q, ss$transitions, prior=alpha.prior, h)

        # Update regression
        # Update the model object input with the new data -- this is
        # not very efficient right now, but it works

        # Need to re-initialize the init.model object with the latest
        # forecast data so it can be an input to the next sample calls
        #
        # This should be really done with some BVAR setup function
        # that is then used in all of the later models in the pkg --
        # something to do later for speed.

        init.model <- updateinit(Y=rbind(y,fcast), init.model)

        # Update regression steps
        hreg <- hregime.reg2(h, m, p, ss$SS, init.model)

        # Draw variances

#        cat("Iteration :", j, "\n")
#        print(hreg$Sigmak)

        Sigmai <- Sigma.draw(m, h, ss, hreg, Sigmai$Sigmai)

#        print(Sigmai)

        # Draw VAR regression coefficients
        Betai <- Beta.draw(m, p, h, Sigmai$Sigmai, hreg, init.model,
                           Betai$Betai)

        # Split out AR coefs and intercepts
        tmp <- Betai2coefs(Betai$Betai, m, p, h)
        ar.coefsi <- tmp$ar.coefsi
        intercepts <- tmp$intercepts

        # Update error estimate
        e <- residual.update(m, h, init.model, Betai$Betai, e)

        # Permute regimes if necessary / order regimes

        # Save results if past burnin=N1.
        # Store forecasts and states
        if(j > N1)
        {
            for(i in 1:m)
            {
                forecasts[(j-N1),,i] <- t(fcast[,i])
                states[[(j-N1)]] <-
                    as.bit.integer(as.integer(ss$SS[,1:(h-1)]))
            }
        }

        # Print out iteration information
        if(j%%1000==0) cat("Iteration : ", j, "\n")
    }

    # Define the output object and its attributes
    class(forecasts) <- c("forecast.MSBVAR")
    class(states) <- c("SS")

    output <- list(forecasts=forecasts, ss.sample=states, k=k, h=h)

    # Add classing and attributes here.  These will be use later for
    # the plotting and other summary functions.
    attr(output, "eqnames") <- attr(x, "eqnames")
    tmp <- tsp(x$init.model$y)
    attr(output, "start") <- tmp[1]
    attr(output, "end") <- tmp[2]
    attr(output, "freq") <- tmp[3]

    return(output)
}