R/autovar.R

In santoscs/pimfc: Prevendo Inflacao Com Modelo Fatores Comuns

#' Fit best VAR model to multivariate time series
#'
#'Returns best VAR 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 = freq)
    }
  }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 zoo
#' @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)))
}



#' 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.
#' @param h Length ahead of the forecast.
#' 
#' @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 stats
#' @export

comp.test <- function(y, fA, fB, h){
  yh <- y - glag(y, -h)
  fA <- fA - glag(y, -h)
  fB <- fB - glag(y, -h)
  fit <- stats::lm(yh ~ fA + fB) 
  return(lmtest::coeftest(fit, vcov. = sandwich::NeweyWest(fit)))
}

#' Generalized Lag 
#' 
#' Compute a lagged version of a vector like time series
#' 
#' @param x A vector or univariate time series
#' @param k	The number of lags
#' 
#' @return A vector object with the same length of x.
#' 
#' @import stats
#' 
#' @export
#' 
glag <- function(x, k = -1){
  h <- abs(k)
  if(k>0){
    xh <- c(x, rep(NA, h))
    ans <- embed(xh, h+1)[,1]
  }
    
  if(k<0){
    xh <- c(rep(NA, h), x) 
    ans <- embed(xh, h+1)[,(h+1)]
  }
  if(k==0){
    ans <- x
  }
  return(ans)
}