Nothing
#' VLSTAR Prediction
#'
#' One-step or multi-step ahead forecasts, with interval forecast, of a VLSTAR object.
#'
#' @param object An object of class \sQuote{\code{VLSTAR}} obtained through \command{VLSTAR()}
#' @param n.ahead An integer specifying the number of ahead predictions
#' @param conf.lev Confidence level of the interval forecast
#' @param st.new Vector of new data for the transition variable
#' @param M An integer with the number of errors sampled for the \code{Monte Carlo} method
#' @param B An integer with the number of errors sampled for the \code{bootstrap} method
#' @param st.num An integer with the index of dependent variable if \code{st.new} is \code{NULL}
#' and the transition variable is a lag of one of the dependent variables
#' @param method A character identifying which multi-step ahead method should be used among \code{naive}, \code{Monte Carlo} and \code{bootstrap}
#' @param newdata \code{data.frame} or \code{matrix} of new data for the exogenous variables
#' @param \dots further arguments to be passed to and from other methods
#' @return A \code{list} containing:
#'\item{forecasts}{\code{data.frame} of predictions for each dependent variable and the (1-\eqn{\alpha}) prediction intervals}
#'\item{y}{a matrix of values for y}
#' @aliases predict print.vlstarpred
#' @references Granger C.W.J. and Terasvirta T. (1993), Modelling Non-Linear Economic Relationships. \emph{Oxford University Press};
#'
#' Lundbergh S. and Terasvirta T. (2007), Forecasting with Smooth Transition Autoregressive Models. \emph{John Wiley and Sons};
#'
#' Terasvirta T. and Yang Y. (2014), Specification, Estimation and Evaluation of Vector Smooth Transition Autoregressive Models with Applications. \emph{CREATES Research Paper 2014-8}
#' @author Andrea Bucci and Eduardo Rossi
#' @seealso \code{\link{VLSTAR}} for log-likehood and nonlinear least squares estimation of the VLSTAR model.
#' @keywords VLSTAR
predict.VLSTAR <- function(object, ..., n.ahead = 1, conf.lev = 0.95, st.new = NULL, M = 5000, B = 1000,
st.num = NULL, newdata = NULL, method = c('naive', 'Monte Carlo', 'bootstrap')){
st_new <- st.new
method <- match.arg(method)
resid <- object$residuals
ncoly <- ncol(object$Data[[1]])
nrowy <- nrow(object$Data[[1]])
alpha = 1-conf.lev
k <- ncol(object$exo)
In <- diag(ncol(object$Data[[1]]))
datamat <- object$Data[[2]]
if(is.null(object$exo)){
if(object$constant == T){
constant <- matrix(rep(1, n.ahead), nrow = n.ahead, ncol = 1)
yy <- datamat[nrow(datamat),c(2:(ncoly*object$p+1))]
}else{
yy <- datamat[nrow(datamat),c(1:(ncoly*object$p))]
}
}else{
if(is.null(newdata) | length(newdata) != k){
stop('Please, provide valid new data for the exogenous variables')
}
if(object$constant == TRUE){
constant <- matrix(rep(1, n.ahead), nrow = n.ahead, ncol = 1)
yy <- datamat[nrow(datamat),c(2:(ncoly*object$p+1))]
}else{
yy <- datamat[nrow(datamat),c(1:(ncoly*object$p))]
}
newdata <- as.matrix(newdata)
}
BB <- object$B
BB1 <- object$B[,seq(1, ncol(BB))[c(rep(TRUE, ncoly), rep(FALSE, ncoly))]]
BB2 <- object$B[,seq(1, ncol(BB))[c(rep(FALSE, ncoly), rep(TRUE, ncoly))]]
Gtilde <- t(object$Gtilde[[length(object$Gtilde)]])
pred <- matrix(NA, ncol = ncoly, nrow = n.ahead)
Z <- c(constant[1, ], yy, newdata)
pred[1,] <- as.matrix((Gtilde%*%t(BB)%*%Z))
tmp <- pred[1, ]
yy <- c(tmp)
if(is.null(st.new)){
if(is.null(st.num)){
stop('Please, provide valid index for the transition variable')
}
st_new <-c(as.matrix(object$st[length(object$st)]), pred[1,st.num])
}
if(!is.null(st.new)){
if(length(st.new) != n.ahead){
stop('The length of the new data for the transition variable should be equal to the number of steps ahead')
}
}
lower <- matrix(NA, ncol = ncoly, nrow = n.ahead)
upper <- matrix(NA, ncol = ncoly, nrow = n.ahead)
for(j in 1:ncoly){
std_err <- sqrt((ncoly-object$p-1)/(nrowy-object$p-ncoly)*qf(1-alpha/2, df1 = ncoly, df2 = (nrowy-object$p-ncoly)))*sqrt((1+t(Z)%*%MASS::ginv(t(object$Data[[2]])%*%object$Data[[2]])%*%Z)*object$Omega[j,j])
lower[1,j] <- pred[1,j]-std_err
upper[1,j] <- pred[1,j]+std_err
}
if(n.ahead >1){
if(method == 'naive'){
for(i in 2 : n.ahead){
Z <- c(constant[i, ], yy, newdata)
GT <- matrix(NA, ncol = 1, nrow = ncoly)
for(j in 1:ncoly){
GT[j,] <- (1+exp(-object$Gammac[j,1]*(st_new[i]-object$Gammac[j,2])))^(-1)
}
Gt <- diag(GT[,1])
Gtilde <- t(cbind(In, Gt))
pred[i,] <- as.matrix((t(Gtilde)%*%t(BB)%*%Z))
tmp <- pred[i, ]
yy <- c(tmp)
if(is.null(st.new)){
st_new <- c(as.vector(st_new),tmp[st.num])
}
for(j in 1:ncoly){
std_err <- sqrt((ncoly-object$p-1)/(nrowy-object$p-ncoly)*qf(1-alpha/2, df1 = ncoly, df2 = (nrowy-object$p-ncoly)))*sqrt((1+t(Z)%*%MASS::ginv(t(object$Data[[2]])%*%object$Data[[2]])%*%Z)*object$Omega[j,j])
lower[i,j] <- pred[1,j]-std_err
upper[i,j] <- pred[1,j]+std_err
}
}
}else if(method == 'Monte Carlo'){
for(i in 2 : n.ahead){
epsilon_mc <- matrix(NA, nrow = M, ncol = ncoly)
for(k in 1:ncoly){
epsilon_mc[,k] <- rnorm(M,mean=0,sd=sd(resid[,k])*3)
}
nonlinear_MC <- matrix(NA, ncol = M, nrow = ncoly)
for(m in 1 : M){
ymc <- yy + epsilon_mc[m, ]
Zmc <- c(constant[i, ], ymc, newdata)
GT <- matrix(NA, ncol = 1, nrow = ncoly)
for(j in 1:ncoly){
GT[j,] <- (1+exp(-object$Gammac[j,1]*(st_new[i]-object$Gammac[j,2])))^(-1)
}
Gt <- diag(GT[,1])
Gtilde <- t(cbind(In, Gt))
nonlinear_MC[,m] <- t(Gtilde)%*%t(BB)%*%as.matrix(Zmc)
}
nonlinearMC <- t(nonlinear_MC)
Z <- c(constant[i, ], yy, newdata)
pred[i,] <- t(t(BB1)%*%Z + rowMeans(nonlinear_MC))
tmp <- pred[i, ]
yy <- c(tmp)
if(is.null(st.new)){
st_new <- c(as.vector(st_new),tmp[st.num])
}
for(j in 1:ncoly){
lower[i,j] <- t(t(BB1)%*%Z)[j] + quantile(nonlinearMC[,j], (alpha/2))
upper[i,j] <-t(t(BB1)%*%Z)[j] + quantile(nonlinearMC[,j], (1-alpha/2))
}
}
}else if(method == 'bootstrap'){
for(i in 2 : n.ahead){
epsilon_bo <- matrix(NA, nrow = B, ncol = ncoly)
nonlinear_BO <- matrix(NA, ncol = B, nrow = ncoly)
for(k in 1:ncoly){
for(b in 1 : B){
epsilon_bo[b,k] <- sample(resid[,k],1)
ybo <- yy + epsilon_bo[b, ]
Zbo <- c(constant[i, ], ybo, newdata)
GT <- matrix(NA, ncol = 1, nrow = ncoly)
for(j in 1:ncoly){
GT[j,] <- (1+exp(-object$Gammac[j,1]*(st_new[i]-object$Gammac[j,2])))^(-1)
}
Gt <- diag(GT[,1])
Gtilde <- t(cbind(In, Gt))
nonlinear_BO[,b] <- t(Gtilde)%*%t(BB)%*%Zbo
}}
nonlinearBO <- t(nonlinear_BO)
Z <- c(constant[i, ], yy, newdata)
pred[i,] <- t(t(BB1)%*%Z + rowMeans(nonlinear_BO))
tmp <- pred[i, ]
yy <- c(tmp)
if(is.null(st.new)){
st_new <- c(as.vector(st_new),tmp[st.num])
}
for(j in 1:ncoly){
lower[i,j] <- t(t(BB1)%*%Z)[j] + quantile(nonlinearBO[,j], (alpha/2))
upper[i,j] <-t(t(BB1)%*%Z)[j] + quantile(nonlinearBO[,j], (1-alpha/2))
}
}
}
}
colnames(pred) <- colnames(object$Bhat)
rownames(pred) <- paste("Step", 1:n.ahead)
forecasts <- list()
for(i in 1 : ncoly){
forecasts[[i]] <- cbind(pred[, i], lower[, i], upper[, i])
colnames(forecasts[[i]]) <- c("fcst", paste("lower ", (conf.lev*100), "%", sep = ''), paste("upper ", (conf.lev*100), "%", sep = ''))
}
names(forecasts) <- colnames(object$Bhat)
results <- list(forecasts = forecasts, y = object$Data[[1]])
class(results) <- "vlstarpred"
return(results)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.