#' Function to obtain h-step ahead direct forecast based on estimated VAR, VARX or VARMA model
#' @param fit Fitted sparse VAR, VARX or VARMA model.
#' @param h Desired forecast horizon. Default is h=1.
#' @export
#' @return Vector of length k containing the h-step ahead forecasts for the k time series.
#' @examples
#' data(var.example)
#' VARfit <- sparseVAR(Y=scale(Y.var), selection = "cv") # sparse VAR
#' VARforecast <- directforecast(fit=VARfit, h=1)
directforecast <- function(fit, h=1){
model <- "none"
if ("bigtime.VAR" %in% class(fit)) model <- "VAR"
if ("bigtime.VARX" %in% class(fit)) model <- "VARX"
if ("bigtime.VARMA" %in% class(fit)) model <- "VARMA"
if(h<=0){
stop("Forecast horizon h must be a strictly positive integer.")
}
if(fit$h!=h){
stop("Specify the same forecast horizon h as argument in the sparseVAR, sparseVARX, or sparseVARMA function and the directforecast function")
}
if(!is.element(model, c("VAR", "VARX", "VARMA"))){
stop("The model needs to be either VAR, VARX or VARMA")
}
# Preliminaries
k <- ncol(fit$Y)
if(model=="VAR"){
# We need to catch the situation in which no selection was done
if (fit$selection == "none"){
fcst <- array(dim = c(1, ncol(fit$Y), length(fit$lambdas)))
for (i in 1:length(fit$lambdas)){
mod_tmp <- fit
mod_tmp$selection <- "tmp"
mod_tmp$Phihat <- fit$Phihat[, , i]
mod_tmp$phi0hat <- fit$phi0hat[, , i]
fcst[, , i] <- directforecast(mod_tmp, h = h)
}
return(fcst)
}
Y <- fit$Y
p <- fit$p
Phi <- fit$Phihat
phi0 <- fit$phi0hat
if(is.null(Phi) | is.null(phi0)){
stop("Please provide a fitted VAR model")
}
VARFIT <- HVARmodel(Y=Y, p=p, h=h)
if(k==1){
VARFIT$fullY <- as.matrix(VARFIT$fullY)
Zpred <- c(c(t(VARFIT$fullY[nrow(VARFIT$fullY):(nrow(VARFIT$fullY)-p+1),1:ncol(VARFIT$fullY)])))
Phi <- matrix(Phi, nrow=1)
}else{
Zpred <- c(c(t(VARFIT$fullY[nrow(VARFIT$fullY):(nrow(VARFIT$fullY)-p+1),1:ncol(VARFIT$fullY)])))
}
Ypred <- matrix(Zpred, nrow=VARFIT$p, byrow=T)
Yhat <- phi0
for(i.p in 1:p){
if(k==1){
Yhat <- Yhat + matrix(Phi[, ((i.p-1)*k+1):(i.p*k) ], nrow=k)%*%Ypred[i.p,]
}else{
Yhat <- Yhat + Phi[, ((i.p-1)*k+1):(i.p*k) ]%*%Ypred[i.p,]
}
}
}
if(model=="VARX"| model=="VARMA"){
if(model=="VARX"){
m <- ncol(fit$X)
Y <- fit$Y
U <- fit$X
p <- fit$p
q <- fit$s
Phi <- fit$Phihat
Theta <- fit$Bhat
phi0 <- fit$phi0hat
if (fit$selection == "none"){
dir_fcst <- array(dim = c(1, ncol(Y), dim(Phi)[3]))
for (i in 1:dim(Phi)[3]){
mod_tmp <- fit
mod_tmp$Phihat <- fit$Phihat[,,i]
mod_tmp$phi0hat <- fit$phi0hat[,,i]
mod_tmp$Bhat <- fit$Bhat[,,i]
mod_tmp$selection = "tmp"
dir_fcst[,,i] <- directforecast(mod_tmp, h = h)
}
return(dir_fcst)
}
if(is.null(Phi) | is.null(Theta) | is.null(phi0)){
stop("Please provide a fitted VARX model")
}
}
if(model=="VARMA"){
m <- ncol(fit$U)
Y <- fit$Y
U <- fit$U
p <- fit$VARMAp
q <- fit$VARMAq
Phi <- fit$Phihat
Theta <- fit$Thetahat
phi0 <- fit$phi0hat
if (fit$VARMAselection == "none"){
# If no selection was done, then fitted values are tensors
dir_fcst <- array(dim = c(1, ncol(Y), dim(Phi)[3]))
for (i in 1:length(fit$PhaseII_lambdaPhi)){
mod_tmp <- fit
mod_tmp$Phihat <- fit$Phihat[, , i]
mod_tmp$Thetahat <- fit$Thetahat[, , i]
mod_tmp$phi0hat <- fit$phi0hat[, , i]
mod_tmp$VARMAselection <- "tmp"
dir_fcst[, , i] <- directforecast(mod_tmp, h = h)
}
return(dir_fcst)
}
if(is.null(Phi) | is.null(Theta) | is.null(phi0)){
stop("Please provide a fitted VARMA model")
}
}
VARXFIT <- HVARXmodelFORECAST(Y=Y, X=U, p=p, s=q, h=h)
if(k==1){
VARXFIT$fullY <- as.matrix(VARXFIT$fullY)
Phi <- matrix(Phi, nrow=1)
Theta <- matrix(Theta, nrow=1)
}
Zpred <- c(c(t(VARXFIT$fullY[nrow(VARXFIT$fullY):(nrow(VARXFIT$fullY)-p+1), 1:ncol(VARXFIT$fullY)])))
Ypred <- matrix(Zpred, nrow=VARXFIT$p, byrow=T)
if(m==1){
VARXFIT$fullXRESP <- as.matrix(VARXFIT$fullXRESP)
}
Xpred <- c(c(t(VARXFIT$fullXRESP[nrow(VARXFIT$fullXRESP):(nrow(VARXFIT$fullXRESP)-q+1),1:ncol(VARXFIT$fullXRESP)])))
Upred <- matrix(Xpred,nrow=VARXFIT$s,byrow=T)
Yhat <- phi0
for(i.p in 1:p){
if(k==1){
Yhat <- Yhat + matrix(Phi[, ((i.p-1)*k+1):(i.p*k) ], nrow=k)%*%Ypred[i.p,]
}else{
Yhat <- Yhat + Phi[, ((i.p-1)*k+1):(i.p*k) ]%*%Ypred[i.p,]
}
}
for(i.q in 1:q){
if(m==1){
Yhat <- Yhat + matrix(Theta[, ((i.q-1)*m+1):(i.q*m) ], nrow=k)%*%Upred[i.q,]
}else{
Yhat <- Yhat + Theta[, ((i.q-1)*m+1):(i.q*m) ]%*%Upred[i.q,]
}
}
}
return("Yhat"=c(Yhat))
}
HVARXmodelFORECAST<-function(Y, X, p, s, h=1){
# Preliminaries
k <- ncol(Y)
kX <-ncol(X)
# Response and predictor matrices
m <- max(s, p)
DATAY <- embed(Y, dimension=m+h)
fullY <- DATAY[, 1:k]
fullZNEW <- as.matrix(as.matrix(DATAY[,-c(1:k)])[,(1:((p+h-1)*k))])
fullZ <- t(fullZNEW[,(ncol(fullZNEW)-k*p+1):ncol(fullZNEW)])
DATAX <- embed(X,dimension=m+h)
fullXRESP <- DATAX[,1:kX]
fullXNEW <- as.matrix(as.matrix(DATAX[,-c(1:kX)])[,(1:((s+h-1)*kX))])
fullX <- t(fullXNEW[,(ncol(fullXNEW)-kX*s+1):ncol(fullXNEW)])
out <- list("fullY"=fullY, "fullX"=fullX, "fullZ"=fullZ, "k"=k, "kX"=kX, "p"=p, "s"=s,
"fullXRESP"=fullXRESP)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.