我fanchart从中获取了该功能vars并对其进行了修改，以添加您需要的一些功能（问题 1、3 和 4）。使用库中的功能已经可以解决问题 2 和 5。结果图如下所示。
fanchart这是我调用的被黑函数fanchart2。主要修改是：更改循环中的 y 轴标签（问题 1），将轴添加到循环中的每个图（问题 3），以及添加预测点（问题 4）。
#### hacked fanchart function ####
# base function taken from vars packaged
# hacked to allow more input as required for this problem
fanchart2 <- function(x,
colors = NULL,
cis = NULL,
names = NULL,
main = NULL,
ylab = NULL,
xlab = NULL,
col.y = NULL,
add.preds = NULL,
nc,
plot.type = c("multiple","single"),
mar = par("mar"),
oma = par("oma"),
lab.at=NULL,
lab.text=NULL,...) {
if (!(class(x) == "varprd")){
stop("\nPlease provide an object of class 'varprd',\ngenerated by predict-method for objects of class 'varest'.\n")
}
if (is.null(colors)){
colors <- gray(sqrt(seq(from = 0.05, to = 1, length = 9)))
}
if (is.null(cis)) {
cis <- seq(0.1, 0.9, by = 0.1)
}
else {
if ((min(cis) <= 0) || (max(cis) >= 1)){
stop("\nValues of confidence intervals must be in(0, 1).\n")
}
if (length(cis) > length(colors)){
stop("\nSize of 'colors' vector must be at least as long as\nsize of 'cis' vector\n")
}
}
n.regions <- length(cis)
n.ahead <- nrow(x$fcst[[1]])
K <- ncol(x$endog)
e.sample <- nrow(x$endog)
endog <- x$endog
fcst <- NULL
for (j in 1:n.regions) {
fcst[[j]] <- predict(x$model, n.ahead = n.ahead, ci = cis[j],
dumvar = x$exo.fcst)$fcst
}
xx <- seq(e.sample, length.out = n.ahead + 1)
xx <- c(xx, rev(xx))
op <- par(no.readonly = TRUE)
plot.type <- match.arg(plot.type)
ynames <- colnames(endog)
if (is.null(names)) {
names <- ynames
}
else {
names <- as.character(names)
if (!(all(names %in% ynames))) {
warning("\nInvalid variable name(s) supplied, using first variable.\n")
names <- ynames[1]
}
}
nv <- length(names)
ifelse(is.null(main), main <- paste("Fanchart for variable",
names), main <- rep(main, nv)[1:nv])
ifelse(is.null(ylab), ylab <- "", ylab <- ylab)
ifelse(is.null(xlab), xlab <- "", xlab <- xlab)
ifelse(is.null(col.y), col.y <- "black", col.y <- col.y)
if (plot.type == "single") {
if (nv > 1)
par(ask = TRUE)
par(mar = mar, oma = oma)
}
else if (plot.type == "multiple") {
if (missing(nc)) {
nc <- ifelse(nv > 4, 2, 1)
}
nr <- ceiling(nv/nc)
par(mfcol = c(nr, nc), mar = mar, oma = oma)
}
for (i in 1:nv) {
ymax <- max(c(fcst[[n.regions]][names[i]][[1]][, 3]),
endog[, names[i]])
ymin <- min(c(fcst[[n.regions]][names[i]][[1]][, 2]),
endog[, names[i]])
yy1 <- c(endog[e.sample, names[i]], fcst[[1]][names[i]][[1]][,
2], rev(c(endog[e.sample, names[i]], fcst[[1]][names[i]][[1]][,
3])))
plot.ts(c(endog[, names[i]], rep(NA, n.ahead)),
main = main[i],
ylim = c(ymin, ymax),
ylab = ylab[i],#### question 1 #### modivied ylab to depend on the loop counter
xlab = xlab,
col = col.y,
...)
polygon(xx, yy1, col = colors[1], border = colors[1])
if (n.regions > 1) {
for (l in 2:n.regions) {
yyu <- c(endog[e.sample, names[i]], fcst[[l]][names[i]][[1]][,
3], rev(c(endog[e.sample, names[i]], fcst[[l -
1]][names[i]][[1]][, 3])))
yyl <- c(endog[e.sample, names[i]], fcst[[l -
1]][names[i]][[1]][, 2], rev(c(endog[e.sample,
names[i]], fcst[[l]][names[i]][[1]][, 2])))
polygon(xx, yyu, col = colors[l], border = colors[l])
polygon(xx, yyl, col = colors[l], border = colors[l])
}
}
#### question 4 ####
# if a matrix of points at various times in the prediction is sent to the function
# they will be plotted here
# standard adjustments to color and pch are possible
# assumes a matrix of values is given with columns in order of the variables (price=col 1)
# and NA values are times without prediction and not plotted
if(is.null(add.preds)==F){
points(x=xx[2:(length(xx)/2)]
,y=add.preds[,i]
,col='gray48'
,pch=16)
}
#### question 3: adding axis to each plot inside the hacked function ####
# standard modifications can be done to this function
if(is.null(lab.at)==F){
axis(1,
at=lab.at,
labels=lab.text,
las=3,
line=1,
cex.axis=0.6)
}
}
on.exit(par(op))
}
将此作为用户定义的函数加载后，可以运行以下基于原始示例的代码来生成图形。我不确定将什么时间分配给预测，所以我选择了一些示例。
library(vars)
price <- as.vector(c(3755,3243,3109,2990,2949,3021,3104,2988,3014,2999,3090,3209,3039,2748,2671,2556,2554,2650,2627,2560))
people <- as.vector(c(4228,4966,4614,4752,4545,4851,4598,4597,4713,4672,4833,4790,4844,4995,5068,4918,4909,4807,5024,4898))
df <- cbind(price,people)
cointest <- ca.jo(df,
K = 5,
type = "eigen",
ecdet = "trend",
spec = "transitory")
vecm.level <- vec2var(cointest, r = 1)
vecm.pred <- predict(vecm.level, n.ahead = 6)
tmin <- as.Date("2016-08-01")
tmax <- as.Date("2018-03-01")
tlab <- seq(tmin, tmax, by="month")
time <- substr(tlab, 0, 7)
fanchart2(vecm.pred
,xaxt="n"
,ylab = c("Price (€)","Volume") #### question 1: see hacked function ###
,main = c("Price","People")
,add.preds = matrix(c(NA,2525,NA,NA,2500,NA,rep(NA,6)),byrow = F,ncol=2), #### question 4: see hacked function ####
,las=1 #### question 2: rotates text to horizontal for easier viewing ####
,colors=heat.colors(9) ##### question 5: capability in pre-hacked function ####
,lab.at = seq(1,length(tlab)+6,1) #### question 3: see hacked function, 6 from the n.ahead in vecm.pred ####
,lab.text = c( time, rep(NA,6)) #### question 3: see hacked function
)
