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]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
