R help - Nov 2012 - Question about VAR (Vector Autoregression) in differences.

Folks,
  I have been using the VAR {vars} program to find a fit for the following
bi-variate time series (subset):

bivariateTS<-structure(c(0.950415958293559, 0.96077848972081,
0.964348957109053,
0.967852884998915, 0.967773510751625, 0.970342843688257, 0.97613937178359, 
0.980118627997436, 0.987059493773907, 0.99536830931504, 1.00622672085718, 
1.01198013845981, 1.01866618122606, 1.02039666233818, 1.02400374633675, 
1.02493877778841, 1.02865185614599, 1.03337978432753, 1.03244791539977, 
1.03237166203917, 1.02853704067597, 1.0271704031946, 1.02588255626357, 
1.02300209859012, 1.02397946952311, 1.02052043198189, 1.01593679111077, 
1.01535467405129, 1.01368552421158, 1.00833115261092, 1.00495328099247, 
1.00334161454411, 1.00029163432818, 0.999268774926758, 0.998174371988049, 
0.999643729403106, 1.00255061235855, 1.0036328278581, 1.00343133578339, 
1.00020064935273, 0.478881778679413, 0.404641679179503, 0.622804024030457, 
0.474677494404522, 0.433851413612414, 0.544920292447026, 0.561093836992123, 
0.563908924197768, 0.395583533713216, 0.273458857040352, 0.284537535240248, 
0.277701982320208, 0.451280071722745, 0.401744885050984, 0.533760449322162, 
0.488858723259857, 0.426799781286085, 0.503196832449693, 0.526637755320189, 
0.599897302279671, 0.519986312138689, 0.460526054741527, 0.550734747489121, 
0.487644621564283, 0.524566370547446, 0.700614666580712, 0.620626909353966, 
0.661635856684872, 0.631149920530168, 0.668383815717266, 0.646705357249337, 
0.685555413824799, 0.650175601510345, 0.679512593649786, 0.67242973935153, 
0.599346651479619, 0.630852735978974, 0.657541680143728, 0.694629910826878, 
0.654572164624051), .Dim = c(40L, 2L), .Dimnames = list(NULL, 
    c("a", "b")))

I have used the "usual" techniques to show that both of the series are
I(1), (integrated of order 1) AND that there is no co-integration. In addition I
used the VARselect program to estimate the number of lags.

To estimate the relationship between these variables the standard recipe (as far
as I know) is to take first differences of the individual series and apply the
VAR program to the new bivariate difference series.

The VAR program seems to derive a reasonable set of coefficients.

My Question:
How do I simulate future paths for the original pair of (undifferenced) series?

Thanks for your time and help,
KW

________________________________

The rest of this note is just how I have done this up until now...

The method I have been using is to take a simulated path of the differences
(note that in my case there are 2 simulated paths) say
> x.d
[1]   9  90 900

and then add back the starting value (1 in this case)
> c(1, 1+cumsum(x.d))
[1]    1   10  100 1000 

A kind of re-integration of the series.

In a real case I would use the final value(s) of bivariateTS (1.0002006,
0.6545722 ) in place of the "1" in the cumsum above.

My problem:
The paths that I get make no sense.

Here is the function I use (an alteration of Professor Pfaff's
predict.varest) to simulate var paths. My apologies if it is bad form to use
alter someone else's code for other purposes.

# Modified VAR simulation function jigged up from the predict function of vars
package
simVARpath<-function(object, n.ahead, mult = 1) {
    K <- object$K
    p <- object$p
    obs <- object$obs
    type <- object$type
    data.all <- object$datamat
    ynames <- colnames(object$y)
    n.ahead <- as.integer(n.ahead)
    Z <- object$datamat[, -c(1:K)]
    B <- Bcoef(object)
    if (type == "const") {
        Zdet <- matrix(rep(1, n.ahead), nrow = n.ahead, ncol = 1)
        colnames(Zdet) <- "const"
    }
    else if (type == "trend") {
        trdstart <- nrow(Z) + 1 + p
        Zdet <- matrix(seq(trdstart, length = n.ahead), nrow = n.ahead,
            ncol = 1)
        colnames(Zdet) <- "trend"
    }
    else if (type == "both") {
        trdstart <- nrow(Z) + 1 + p
        Zdet <- matrix(c(rep(1, n.ahead), seq(trdstart, length = n.ahead)), 
            nrow = n.ahead, ncol = 2)
        colnames(Zdet) <- c("const", "trend")
    }
    else if (type == "none") {
        Zdet <- NULL
    }
    if (!is.null(eval(object$call$season))) {
        season <- eval(object$call$season)
        seas.names <- paste("sd", 1:(season - 1), sep =
"")
        cycle <- tail(data.all[, seas.names], season)
        seasonal <- as.matrix(cycle, nrow = season, ncol = season - 
            1)
        if (nrow(seasonal) >= n.ahead) {
            seasonal <- as.matrix(cycle[1:n.ahead, ], nrow = n.ahead, 
                ncol = season - 1)
        }
        else {
            while (nrow(seasonal) < n.ahead) {
                seasonal <- rbind(seasonal, cycle)
            }
            seasonal <- seasonal[1:n.ahead, ]
        }
        rownames(seasonal) <- seq(nrow(data.all) + 1, length = n.ahead)
        if (!is.null(Zdet)) {
            Zdet <- as.matrix(cbind(Zdet, seasonal))
        }
        else {
            Zdet <- as.matrix(seasonal)
        }
    }
    if (!is.null(eval(object$call$exogen))) {
        if (is.null(dumvar)) {
            stop("\nNo matrix for dumvar supplied, but object varest
contains exogenous variables.\n")
        }
        if (!all(colnames(dumvar) %in% colnames(data.all))) {
            stop("\nColumn names of dumvar do not coincide with
exogen.\n")
        }
        if (!identical(nrow(dumvar), n.ahead)) {
            stop("\nRow number of dumvar is unequal to n.ahead.\n")
        }
        if (!is.null(Zdet)) {
            Zdet <- as.matrix(cbind(Zdet, dumvar))
        }
        else {
            Zdet <- as.matrix(dumvar)
        }
    }
    Zy <- as.matrix(object$datamat[, 1:(K * (p + 1))])

    cov<-summary(object)$covres
    forecast <- matrix(NA, ncol = K, nrow = n.ahead)
    lasty <- c(Zy[nrow(Zy), ])
    for (i in 1:n.ahead) {
        lasty <- lasty[1:(K * p)]
        Z <- c(lasty, Zdet[i, ])
        rv<-mult*matrix(rmvsnorm(1, dim = dim(cov)[1], Omega = cov), ncol =
1)
        forecast[i, ] <- B %*% Z + rv
        temp <- forecast[i, ]
        lasty <- c(temp, lasty)
    }
    forecast
}



--


	[[alternative HTML version deleted]]