if (FALSE) {
## We set up the same model as in the example of mvar(), but
## specify one time-varying parameter, and try to recover it with
## tvmvar()
# a) Specify time-varying VAR model
p <- 6 # Six variables
type <- c("c", "c", "c", "c", "g", "g") # 4 categorical, 2 gaussians
level <- c(2, 2, 4, 4, 1, 1) # 2 categoricals with 2 categories, 2 with 5
max_level <- max(level)
n_timepoints <- 4000
lags <- c(1, 3, 9) # include lagged effects of order 1, 3, 9
n_lags <- length(lags)
# Specify thresholds
thresholds <- list()
thresholds[[1]] <- matrix(0, ncol=level[1], nrow=n_timepoints)
thresholds[[2]] <- matrix(0, ncol=level[2], nrow=n_timepoints)
thresholds[[3]] <- matrix(0, ncol=level[3], nrow=n_timepoints)
thresholds[[4]] <- matrix(0, ncol=level[4], nrow=n_timepoints)
thresholds[[5]] <- matrix(0, ncol=level[5], nrow=n_timepoints)
thresholds[[6]] <- matrix(0, ncol=level[6], nrow=n_timepoints)
# Specify standard deviations for the Gaussians
sds <- matrix(NA, ncol=p, nrow=n_timepoints)
sds[, 5:6] <- 1
# Create coefficient array
coefarray <- array(0, dim=c(p, p, max_level, max_level, n_lags, n_timepoints))
# a.1) interaction between continuous 5<-6, lag=3
coefarray[5, 6, 1, 1, 2, ] <- seq(0, .4, length = n_timepoints) # only time-varying parameter
# a.2) interaction between 1<-3, lag=1
m1 <- matrix(0, nrow=level[2], ncol=level[4])
m1[1,1:2] <- 1
m1[2,3:4] <- 1
coefarray[1, 3, 1:level[2], 1:level[4], 1, ] <- m1 # constant across time
# a.3) interaction between 1<-5, lag=9
coefarray[1, 5, 1:level[1], 1:level[5], 3, ] <- c(0, 1) # constant across time
# b) Sample
set.seed(1)
dlist <- tvmvarsampler(coefarray = coefarray,
lags = lags,
thresholds = thresholds,
sds = sds,
type = type,
level = level,
pbar = TRUE)
# c.1) Recover: stationary
set.seed(1)
mvar_obj <- mvar(data = dlist$data,
type = type,
level = level,
lambdaSel = "CV",
lags = c(1, 3, 9),
signInfo = FALSE)
# Did we recover the true parameters?
mvar_obj$wadj[5, 6, 2] # cross-lagged effect of 6 on 5 over lag lags[2] (lag 3)
mvar_obj$wadj[1, 3, 1] # cross-lagged effect of 3 on 1 over lag lags[1] (lag 1)
mvar_obj$wadj[1, 5, 3] # cross-lagged effect of 1 on 5 over lag lags[3] (lag 9)
# c.2) Recover: time-varying
set.seed(1)
mvar_obj <- tvmvar(data = dlist$data,
type = type,
level = level,
estpoints = seq(0, 1, length=10),
bandwidth = .15,
lambdaSel = "CV",
lags = c(1, 3, 9),
signInfo = FALSE)
# Did we recover the true parameters?
mvar_obj$wadj[5, 6, 2, ] # true sort of recovered
mvar_obj$wadj[1, 3, 1, ] # yes
mvar_obj$wadj[1, 5, 3, ] # yes
# Plotting parameter estimates over time
plot(mvar_obj$wadj[5, 6, 2, ],
type="l", ylim=c(-.2,.7),
lwd=2, ylab="Parameter value", xlab="Estimation points")
lines(mvar_obj$wadj[1, 3, 1, ], col="red", lwd=2)
lines(mvar_obj$wadj[1, 5, 3, ], col="blue", lwd=2)
legend("bottomright", c("5 <-- 6", "1 <-- 3", "1 <-- 5"),
lwd = c(2,2,2), col=c("black", "red", "blue"))
# d) Predict values / compute nodewise error
mvar_pred_w <- predict.mgm(object=mvar_obj,
data=dlist$data,
tvMethod = "weighted")
mvar_pred_cM <- predict.mgm(object=mvar_obj,
data=dlist$data,
tvMethod = "closestModel")
mvar_pred_w$errors
mvar_pred_cM$errors
# For more examples see https://github.com/jmbh/mgmDocumentation
}
Run the code above in your browser using DataLab