#' Fit best VAR model to multivariate time series
#'
#'Returns best ARIMA model according to either SC, HQ, AIC or FPE value.
#'The function conducts a search over possible model within the order
#'constraints provided.
#'
#'
#' @param y A multivariate time series
#' @param max.p Determines the highest lag order for lag length selection
#' according to the choosen ic.
#' @param ic Information criterion to be used in model selection.
#' @param seasonal If FALSE, restricts search to models without
#' seasonal dummies.
#'
#' @return A list with class attribute 'autovar' holding the following elements:
#'
#' \item{\code{fit}}{modelo estimado da class \code{varest}.}
#' \item{\code{d}}{ordem de diferenciacao das series.}
#' \item{\code{y}}{series das variaveis originais.}
#' \item{\code{x}}{series das variaveis diferenciadas \code{d} vezes.}
#' \item{\code{ic}}{criterios de informacao do modelo selecionado.}
#'
#' @import vars forecast stats
#' @export
#'
auto.var <-function(y, max.p=6, ic=c("SC", "HQ", "AIC", "FPE"), seasonal=TRUE){
# verifica NA e adiverte
if(anyNA(y))
warning("observacoes com NA foram excluidas")
# informacoes para estimacao
ic <- match.arg(ic)
ic <- paste(ic, "(n)", sep="")
k <- length(y[1,])
n <- length(y[,1])
freq <- frequency(y)
d <- vector()
x <- y
# primeira diferenca para estacionaridade
for(i in 1:k){
d[i] <- forecast::ndiffs(y[,i])
if(d[i]>0)
x[,i] <- stats::ts.union(y[,i]*NA, diff(y[,i], differences = d[i]))[,2]
else
x[,i] <- y[,i]
}
names(d) <- colnames(x) <- colnames(y)
# selecao e estimacao do modelo
x <- na.omit(x)
sele1 <- vars::VARselect(x, lag.max = max.p, type="const")
best1.ic <- apply(sele1$criteria, 1, min)
if(seasonal==TRUE){
sele2 <- vars::VARselect(x, lag.max = max.p, type="const", season = freq)
best2.ic <- apply(sele2$criteria, 1, min)
if(best1.ic[ic]<best2.ic[ic]){
best.ic <- best1.ic
p <- sele1$selection[ic]
# estima o modelo
model <- vars::VAR(x, p=p, type="const")
}else{
best.ic <- best2.ic
p <- sele2$selection[ic]
# estima o modelo
model <- vars::VAR(x, p=p, type="const", season = 12)
}
}else{
best.ic <- best1.ic
p <- sele1$selection[ic]
# estima o modelo
model <- vars::VAR(x, p=p, type="const")
}
result <- list(fit=model, d=d, y=y, x=x, ic=best.ic)
class(result) <- "autovar"
return(result)
}
#' Predict method for objects of class autovar
#'
#' Forecating a VAR object of class 'autovar' with confidence bands
#'
#' @param object An object of class 'autovar'; generated by auto.var().
#' @param n.ahead An integer specifying the number of forecast steps.
#' @param ci The forecast confidence interval.
#' @param ... Currently not used.
#'
#' @return A list with class attribute 'varprd' holding the following
#' elements:
#'
#' \item{\code{fcst}}{A list of matrices per variable containing
#' the forecasted values with lower and upper bounds
#' as well as the confidence interval.}
#' \item{model}{The estimated VAR object for auto.var().}
#'
#' @import stats
#'
#' @export
#'
predict.autovar <- function(object, ..., n.ahead=12, ci=0.95){
ny <- length(object$y[,1])
pr <- predict(object$fit, n.ahead = n.ahead, ci = ci)[[1]]
nomes <- names(pr)
ans <- vector("list")
for(i in 1:length(nomes)){
d <- object$d[i]
if(d>0){
fcst <- diffinv(pr[[i]][,'fcst'], xi = object$y[(ny-d+1):ny,i], differences=d)
lower <- diffinv(pr[[i]][,'lower'], xi = object$y[(ny-1+1):ny,i], differences=1)
upper <- diffinv(pr[[i]][,'upper'], xi = object$y[(ny-1+1):ny,i], differences=1)
}else{
fcst <- pr[[i]][,'fcst']
lower <- pr[[i]][,'lower']
upper <- pr[[i]][,'upper']
}
ans[[nomes[i]]]<- suppressWarnings(cbind(fcst, lower, upper))
if(d==2)
warning(paste("intervalos de predicao para", nomes[i], "considerando apenas primeira diferenca"))
}
result <- list(fcst=ans, model=object)
return(structure(result,class="autovarprd"))
}
#' Forecasting time series com VAR model
#'
#' Returns forecasts and other information for VAR models
#'
#' @param object An object of class "autovar".
#' @param h Number of periods for forecasting.
#' @param level Confidence level for prediction intervals.
#' @param fan If TRUE, level is set to seq(51,99,by=3). This is
#' suitable for fan plots.
#' @param ... Other arguments.
#'
#' @return An object of class "mforecast". An object of class
#' "mforecast" is a list containing at least the following elements:
#'
#' \item{model}{A list containing information about the fitted
#' model}
#' \item{method}{The name of the forecasting method as a character
#' string}
#' \item{mean}{Point forecasts as a time series}
#' \item{lower}{Lower limits for prediction intervals}
#' \item{upper}{Upper limits for prediction intervals}
#' \item{level}{The confidence values associated with the
#' prediction intervals}
#' \item{x}{The original time series (either \code{object} itself
#' or the time series used to create the model stored as \code{object}).}
#' \item{residuals}{Residuals from the fitted model.
#' That is x minus fitted values.}
#' \item{fitted}{Fitted values (one-step forecasts)}
#'
#' @import stats
#' @export
#'
forecast.autovar <- function(object, h=10, level=c(80,95), fan=FALSE, ...)
{
model <- object
object <- model$fit
out <- list(model=object,level=level,x=object$y)
# Get residuals and fitted values and fix the times
out$x <- model$y
out$residuals <- out$fitted <- ts(matrix(NA,nrow=nrow(out$x),ncol=ncol(out$x)))
tsp(out$residuals) <- tsp(out$fitted) <- tsp(out$x)
vres <- residuals(object)
vfits <- fitted(object)
out$residuals[ (nrow(out$residuals)-nrow(vres)+1):nrow(out$residuals), ] <- vres
out$fitted[ (nrow(out$fitted)-nrow(vfits)+1):nrow(out$fitted), ] <- vfits
# Add forecasts with prediction intervals
out$mean <- out$lower <- out$upper <- vector("list",object$K)
for(i in 1:(length(level)))
{
pr <- predict.autovar(model, n.ahead=h, ci=level[i]/100)
for(j in 1:object$K)
{
if(i==1)
out$mean[[j]] <- pr$fcst[[j]][,"fcst"]
out$lower[[j]] <- cbind(out$lower[[j]],pr$fcst[[j]][,"lower"])
out$upper[[j]] <- cbind(out$upper[[j]],pr$fcst[[j]][,"upper"])
}
}
names(out$mean) <- names(out$lower) <- names(out$upper) <- names(pr$fcst)
tspx <- tsp(object$y)
for(j in 1:object$K)
out$mean[[j]] <- ts(out$mean[[j]], frequency=tspx[3], start=tspx[2]+1/tspx[3])
out$method <- paste(paste("VAR(",object$p,")", sep=""), paste(names(pr$fcst), "d =", model$d, collapse = ", "), collapse = ", ")
return(structure(out,class="mforecast"))
}
#' Transform the forecast of an mforecast object for a list
#'
#' Transforma as previsoes forneccida por forecast.autovar para uma
#' lista contendo as previsoes e o intervalo de previsao
#'
#' @param x An object of class "mforecast".
#'
#' @return A list of matrices containing the forecasted values with
#' lower and upper bounds
#'
f2list <- function(x){
n <- length(x$mean)
y <- vector("list", n)
lo <- paste("Lo ",x$level, sep="")
hi <- paste("Hi ",x$level, sep="")
for(i in 1:n){
y[[i]] <- cbind(x$mean[[i]], x$lower[[i]][,1], x$upper[[i]][,1],
x$lower[[i]][,2], x$upper[[i]][,2])
colnames(y[[i]]) <- c('Point Forecast', lo[1], hi[1], lo[2], hi[2])
}
names(y) <- names(x$mean)
return(y)
}
#' Previsao fora da amostra com VAR
#'
#' outsample.var e usada para gerar previsoes fora da amostra
#' usando modelo VAR estimado automaticamente com a funcao
#' auto.var
#'
#' @inheritParams auto.var
#' @param h horizonte de previsao
#' @param k numero de previsoes fora da amostra
#'
#' @return uma lista com
#' fcast(ts): previsao pontual e intervalos de predicao usando arima
#' outdate(num): vetor contendo periodo fora da amostra
#' fit(Arima): ultimo modelo arima estimado
#'
#' @import stats
#'
#' @export
#'
outsample.var <- function(y, max.p = 6, ic = c("SC", "HQ", "AIC", "FPE"), seasonal = TRUE, h, k){
#input:
# x(ts): serie a ser prevista
# h(num): horizonte de previsao
# k(num): numero de previsoes fora da amostra
#output:
# fcast(ts): previsao pontual e intervalos de predicao usando VAR
# outdate(num): vetor contendo periodo fora da amostra
# fit(varest): ultimo modelo VAR estimado
x <- y
Tn <- dim(x)[1]
timeline <- time(x)
freq <- frequency(x)
fcast <- vector("list", dim(x)[2])
outdate <- timeline[(Tn-k+1):Tn]
for(i in 1:k){
# restringe os dados
x.ajuste <- window(x, end=timeline[Tn-k-h+i])
# estima o modelo
fit <- auto.var(x.ajuste, max.p = max.p, ic = ic, seasonal = seasonal)
# get prediction for ith value
aux <- forecast.autovar(fit, h = h)
aux <- f2list(aux)
for(j in 1:fit$fit$K){
df <- data.frame(aux[[j]])
fcast[[j]] <- rbind(fcast[[j]], df[h,])
}
}
fcast <- lapply(fcast, function(x){ts(x, start=outdate[1], frequency = freq)})
names(fcast) <- colnames(x)
return(list(fcast=fcast, outdate=outdate, fit=fit))
}
#' Test for Forecast Encompassing
#'
#' The Test for Forecast Encompassing compares the forecast of two forecast methods.
#'
#' The null hypothesis is that the method A forecast encompasses
#' of method B, i.e., method B have the same information
#' of method A. The alternative hypothesis is that method B
#' have additional information.
#'
#' @param fA Forecast from method A.
#' @param fB Forecast from method B.
#' @param y The observed time series.
#'
#' @return An object of class "coeftest" which is
#' essentially a coefficient matrix with columns
#' containing the estimates, associated standard
#' errors, test statistics and p values.
#'
#'
#' @import lmtest sandwich
#' @export
enc.test <- function(y, fA, fB){
requireNamespace("zoo")
fit <- dyn::dyn$lm(y ~ offset(fA) + I(fB-fA))
return(lmtest::coeftest(fit, vcov. = sandwich::NeweyWest(fit)))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.