#' Forecast future observations of BVAR model
#'
#' This function returns a forecast future observations of BVAR-SV-fatTail model.
#' @param Chain The fatBVARSV object from command BVAR.
#' @param t_pred The time prediction horizon.
#' @return The list of forecast future observations of BVAR model
#' @export
#' @examples
#' \dontrun{
#' forecast1 <- get_forecast(Chain1)
#' }
#'
get_forecast <- function(Chain, y0 = NULL, t_pred = 12, t_current, Nfsample = NULL, thin = 1){
K <- Chain$K
p <- Chain$p
if( is.null(y0)) {
# y <- Chain$y
y0 <- tail(Chain$y, p) # Get last p obs
}
# if( is.null(t_current)) {
# t_current <- which(Chain$y == tail(y0,p)[p,1])
# }
dist <- Chain$dist
SV <- Chain$prior$SV
ndraws <- nrow(Chain$mcmc)
mcmc <- Chain$mcmc[seq(from = 1, to = ndraws, by = thin),]
if (is.null(Nfsample) || (Nfsample == nrow(mcmc)) ){
Nfsample <- nrow(mcmc)
frange <- c(1:Nfsample)
} else {
frange <- sample(nrow(mcmc), Nfsample, replace = T)
}
y_pred <- array(NA, dim = c(t_pred, K, Nfsample))
ymean_pred <- array(NA, dim = c(t_pred, K, Nfsample)) # Predictive mean
yvar_pred <- array(NA, dim = c(t_pred, 0.5*K*(K+1), Nfsample)) # Predictive Chol Sigma
vol_pred <- array(NA, dim = c(t_pred, K, Nfsample)) # Predictive diag(Sigma)
B_mat <- get_post(mcmc, element = "B")
A_mat <- get_post(mcmc, element = "a")
Sigma_mat <- get_post(mcmc, element = "sigma") # No SV
Sigma_H <- sqrt(get_post(mcmc, element = "sigma_h")) # SV
Sigma_G <- sqrt(get_post(mcmc, element = "sigma_G")) # SV of gamma
H_mat <- get_post(mcmc, element = "h")
if (ncol(H_mat)>0) H_mat <- H_mat[,( (t_current - 1) * K +1):(t_current * K)]
Nu_mat <- get_post(mcmc, element = "nu")
Gamma_mat <- get_post(mcmc, element = "gamma")
count_id <- 0
seed_set <- sample.int(n = 50000, size = Nfsample, replace = F)
for (i in frange){
count_id <- count_id+1
b0 <- B_mat[i,]
a0 <- A_mat[i,]
if (SV) {
sigma_h <- Sigma_H[i,]
if (K == 1) {
h <- H_mat[i]
} else {
h <- H_mat[i,]
}
if (dist == "Gaussian") {
pred_tmp <- sim.VAR.Gaussian.SV(K = K, p = p, t_max = t_pred,
b0 = b0, a0 = a0, h = h, sigma_h = sigma_h,
y0 = y0,
seednum = seed_set[count_id], burn_in = 0)
}
if (dist == "Student") {
nu <- Nu_mat[i]
pred_tmp <- sim.VAR.Student.SV(K = K, p = p, t_max = t_pred,
b0 = b0, a0 = a0, h = h, sigma_h = sigma_h,
y0 = y0, nu = nu,
seednum = seed_set[count_id], burn_in = 0)
}
if (dist == "Skew.Student") {
nu <- Nu_mat[i]
gamma <- Gamma_mat[i,]
pred_tmp <- sim.VAR.Skew.Student.SV(K = K, p = p, t_max = t_pred,
b0 = b0, a0 = a0, h = h, sigma_h = sigma_h,
y0 = y0, nu = nu, gamma = gamma,
seednum = seed_set[count_id], burn_in = 0)
}
if (dist == "MT") {
nu <- Nu_mat[i,]
pred_tmp <- sim.VAR.MT.SV(K = K, p = p, t_max = t_pred,
b0 = b0, a0 = a0, h = h, sigma_h = sigma_h,
y0 = y0, nu = nu,
seednum = seed_set[count_id], burn_in = 0)
}
if (dist == "MST") {
nu <- Nu_mat[i,]
gamma <- Gamma_mat[i,]
pred_tmp <- sim.VAR.MST.SV(K = K, p = p, t_max = t_pred,
b0 = b0, a0 = a0, h = h, sigma_h = sigma_h,
y0 = y0, nu = nu, gamma = gamma,
seednum = seed_set[count_id], burn_in = 0)
}
if (dist == "OT") {
nu <- Nu_mat[i,]
pred_tmp <- sim.VAR.OT.SV(K = K, p = p, t_max = t_pred,
b0 = b0, a0 = a0, h = h, sigma_h = sigma_h,
y0 = y0, nu = nu,
seednum = seed_set[count_id], burn_in = 0)
}
if (dist == "OST") {
nu <- Nu_mat[i,]
gamma <- Gamma_mat[i,]
pred_tmp <- sim.VAR.OST.SV(K = K, p = p, t_max = t_pred,
b0 = b0, a0 = a0, h = h, sigma_h = sigma_h,
y0 = y0, nu = nu, gamma = gamma,
seednum = seed_set[count_id], burn_in = 0)
}
if (dist == "dynSkew.Student") {
nu <- Nu_mat[i]
gamma <- tail(Gamma_mat[i,], K)
sigma_G <- Sigma_G[i,]
pred_tmp <- sim.VAR.dynSkew.Student.SV(K = K, p = p, t_max = t_pred,
b0 = b0, a0 = a0, h = h,
y0 = y0, nu = nu, gamma = gamma, sigma_G = sigma_G,
seednum = seed_set[count_id], burn_in = 0)
}
} else {
sigma <- Sigma_mat[i,]
h <- log(sigma)
if (dist == "Gaussian") {
pred_tmp <- sim.VAR.Gaussian.novol(K = K, p = p, t_max = t_pred,
b0 = b0, a0 = a0, h = h,
y0 = y0,
seednum = seed_set[count_id], burn_in = 0)
}
if (dist == "Student") {
nu <- Nu_mat[i]
pred_tmp <- sim.VAR.Student.novol(K = K, p = p, t_max = t_pred,
b0 = b0, a0 = a0, h = h,
y0 = y0, nu = nu,
seednum = seed_set[count_id], burn_in = 0)
}
if (dist == "Skew.Student") {
nu <- Nu_mat[i]
gamma <- Gamma_mat[i,]
pred_tmp <- sim.VAR.Skew.Student.novol(K = K, p = p, t_max = t_pred,
b0 = b0, a0 = a0, h = h,
y0 = y0, nu = nu, gamma = gamma,
seednum = seed_set[count_id], burn_in = 0)
}
if (dist == "MT") {
nu <- Nu_mat[i,]
pred_tmp <- sim.VAR.MT.novol(K = K, p = p, t_max = t_pred,
b0 = b0, a0 = a0, h = h,
y0 = y0, nu = nu,
seednum = seed_set[count_id], burn_in = 0)
}
if (dist == "MST") {
nu <- Nu_mat[i,]
gamma <- Gamma_mat[i,]
pred_tmp <- sim.VAR.MST.novol(K = K, p = p, t_max = t_pred,
b0 = b0, a0 = a0, h = h,
y0 = y0, nu = nu, gamma = gamma,
seednum = seed_set[count_id], burn_in = 0)
}
if (dist == "OT") {
nu <- Nu_mat[i,]
pred_tmp <- sim.VAR.OT.novol(K = K, p = p, t_max = t_pred,
b0 = b0, a0 = a0, h = h,
y0 = y0, nu = nu,
seednum = seed_set[count_id], burn_in = 0)
}
if (dist == "OST") {
nu <- Nu_mat[i,]
gamma <- Gamma_mat[i,]
pred_tmp <- sim.VAR.OST.novol(K = K, p = p, t_max = t_pred,
b0 = b0, a0 = a0, h = h,
y0 = y0, nu = nu, gamma = gamma,
seednum = seed_set[count_id], burn_in = 0)
}
if (dist == "dynSkew.Student") {
nu <- Nu_mat[i]
gamma <- tail(Gamma_mat[i,], K)
sigma_G <- Sigma_G[i,]
pred_tmp <- sim.VAR.dynSkew.Student.novol(K = K, p = p, t_max = t_pred,
b0 = b0, a0 = a0, h = h,
y0 = y0, nu = nu, gamma = gamma, sigma_G = sigma_G,
seednum = seed_set[count_id], burn_in = 0)
}
}
y_pred[,,count_id] <- pred_tmp$y
ymean_pred[,,count_id] <- pred_tmp$y_mean
yvar_pred[,,count_id] <- pred_tmp$y_var
vol_pred[,,count_id] <- pred_tmp$volatility
}
return(list(y_pred = y_pred,
ymean_pred = ymean_pred,
yvar_pred = yvar_pred,
vol_pred = vol_pred))
}
#' Forecast density of future observations of BVAR model
#'
#' This function returns a forecast density of future observations of BVAR-SV-fatTail model.
#' @param Chain The fatBVARSV object from command BVAR.
#' @param y_current The current values of y.
#' @param y_obs_future The future observable values of y.
#' @param t_pred The time prediction horizon.
#' @return The forecast density of future observations of BVAR model
#' @export
#' @examples
#' \dontrun{
#' forecast_dens1 <- forecast_density(Chain1, y_current, y_obs_future, t_current)
#' }
#'
#' @export
forecast_density <- function(Chain, y_current = NULL, y_obs_future, t_current){
K <- Chain$K
p <- Chain$p
if (! ncol(y_obs_future) == K) { stop("ncol(y_obs_future) != K") }
if (is.null(y_current)) y_current <- Chain$y
#if (is.null(t_current)) t_current = which(Chain$y == tail(y_current,p)[p,1] )
t_pred = nrow(y_obs_future)
Nfsample = nrow(Chain$mcmc)
predictive_samples <- get_forecast(Chain = Chain, y0 = tail(y_current,p), t_pred = t_pred, t_current = t_current, thin = 1, Nfsample = Nfsample) # Nfsample
#############################################################
# Non cummulative
#############################################################
log_pred<- matrix(NA, nrow = t_pred, ncol = K)
bilog_pred<- matrix(NA, nrow = t_pred, ncol = K*(K-1)*0.5)
emp_CDF<- matrix(NA, nrow = t_pred, ncol = K)
for (i in c(1:K)){
for (j in c(1:t_pred)){
y_obs <- as.numeric(y_obs_future[j,i])
log_pred[j,i] <- log(mean(dnorm(y_obs, mean = predictive_samples$ymean_pred[j,i,], sd = sqrt(predictive_samples$vol_pred[j,i,]))))
emp_CDF[j,i] <- ecdf(x = predictive_samples$y_pred[j,i,])(y_obs)
}
}
if (K > 1){
Sigma2 <- array(0, dim = c(K, K, Nfsample))
Sigma <- matrix(0, nrow = K, ncol = K)
for (j in c(1:t_pred)){
for (sample_id in c(1:Nfsample)){
Sigma[lower.tri(Sigma, diag = T)] <- predictive_samples$yvar_pred[j,,sample_id]
Sigma2[,,sample_id] <- Sigma %*% t(Sigma)
}
for (i in c(1: (K-1))){
for (k in c((i+1):K)){
pred_bidens <- rep(NA, Nfsample)
for (sample_id in c(1:Nfsample)){
pred_bidens[sample_id] <- dmvn(X = c(y_obs_future[j,i], y_obs_future[j,k]),
mu = predictive_samples$ymean_pred[j, c(i,k),sample_id],
sigma = Sigma2[c(i,k),c(i,k), sample_id], log = FALSE, isChol = FALSE)
}
bilog_pred[j,(i-1)*K + k-i*(i+1)*0.5] <- log(mean(pred_bidens))
}
}
}
}
#############################################################
# Cummulative
#############################################################
clog_pred<- matrix(NA, nrow = t_pred, ncol = K)
cbilog_pred<- matrix(NA, nrow = t_pred, ncol = K*(K-1)*0.5)
cemp_CDF<- matrix(NA, nrow = t_pred, ncol = K)
cy_pred <- apply(predictive_samples$y_pred, MARGIN = c(2,3), FUN = cumsum)
cy_obs_future <- apply(y_obs_future, MARGIN = 2, FUN = cumsum)
# for (i in c(1:K)){
# for (j in c(1:t_pred)){
# y_obs <- as.numeric(cy_obs_future[j,i])
# predict_den <- stats::density(cy_pred[j,i,], adjust = 3)
# id <- which(predict_den$y == 0);
# if (length(id)>0) {
# for (fix in id){
# predict_den$y[fix] <- mean(predict_den$y[(fix-1):(fix+1)])
# }
# }
# appximate_density <- smooth.spline(x = predict_den$x, y = log(predict_den$y))
# pred_dens <- predict(appximate_density, y_obs)$y
# clog_pred[j,i] <- pred_dens
# cemp_CDF[j,i] <- ecdf(x = cy_pred[j,i,])(y_obs)
# }
# }
CRPS <- matrix(NA, nrow = t_pred, ncol = K)
qwCRPS_2t <- matrix(NA, nrow = t_pred, ncol = K) # Quantile weighted CRPS - 2 tails
qwCRPS_rt <- matrix(NA, nrow = t_pred, ncol = K) # Quantile weighted CRPS - right tail
qwCRPS_lt <- matrix(NA, nrow = t_pred, ncol = K) # Quantile weighted CRPS - left tail
alpha <- seq(0.05, 0.95, by = 0.05)
for (i in c(1:K)){
for (j in c(1:t_pred)){
y_obs <- as.numeric(y_obs_future[j,i])
samples1 <- predictive_samples$y_pred[j,i,]
samples2 <- sample(samples1, size = length(samples1), replace = T)
CRPS[j,i] <- - mean(abs(samples1 - y_obs)) + 0.5*mean(abs(samples1 - samples2))
Q_tau <- quantile(predictive_samples$y_pred[j,i,], probs = alpha)
QS <- (y_obs - Q_tau) * ( (y_obs < Q_tau) - alpha )
qwCRPS_2t[j,i] <- 2/length(alpha) * sum( (2 * alpha - 1)^2 * QS )
qwCRPS_rt[j,i] <- 2/length(alpha) * sum( alpha^2 * QS )
qwCRPS_lt[j,i] <- 2/length(alpha) * sum( (1 - alpha)^2 * QS )
}
}
return(list(log_pred = log_pred,
bilog_pred = bilog_pred,
emp_CDF = emp_CDF,
MSFE = (apply(predictive_samples$y_pred, MARGIN = c(1,2), FUN = mean) - y_obs_future)^2,
MAFE = abs(apply(predictive_samples$y_pred, MARGIN = c(1,2), FUN = mean) - y_obs_future),
# predictive_samples = predictive_samples,
predictive_quantile05 = apply(predictive_samples$y_pred, MARGIN = c(1,2), FUN = quantile, probs = 0.05),
predictive_quantile10 = apply(predictive_samples$y_pred, MARGIN = c(1,2), FUN = quantile, probs = 0.10),
predictive_quantile50 = apply(predictive_samples$y_pred, MARGIN = c(1,2), FUN = quantile, probs = 0.50),
predictive_quantile90 = apply(predictive_samples$y_pred, MARGIN = c(1,2), FUN = quantile, probs = 0.90),
predictive_quantile95 = apply(predictive_samples$y_pred, MARGIN = c(1,2), FUN = quantile, probs = 0.95),
t_current = t_current,
clog_pred = clog_pred,
cbilog_pred = cbilog_pred,
cemp_CDF = cemp_CDF,
cMSFE = (apply(cy_pred, MARGIN = c(1,2), FUN = mean) - cy_obs_future)^2,
cMAFE = abs(apply(cy_pred, MARGIN = c(1,2), FUN = mean) - cy_obs_future),
CRPS = CRPS,
qwCRPS_2t = qwCRPS_2t,
qwCRPS_rt = qwCRPS_rt,
qwCRPS_lt = qwCRPS_lt
))
}
#' Recursive forecast of BVAR model
#'
#' This function returns a recursive forecast of future observations of BVAR-SV-fatTail model.
#' @param y the whole data matrix T x k.
#' @param t_start the start point for forecast.
#' @param t_pred The time prediction horizon.
#' @param K The number of variables in y.
#' @param p The number of lags in BVAR model.
#' @return The recursive forecast density of future observations of BVAR model
#' @export
#' @examples
#' \dontrun{
#' recuresive_model1 <- recursive_forecast(y, t_start, t_pred, K, p, dist = "MST", SV = T)
#' }
#'
#' @export
recursive_forecast <- function(y, t_start = 100, t_pred = 12, K, p, dist = "MST", SV = T, numCores = NULL){
t_max = nrow(y)
max_rolling <- t_max - (t_start + t_pred) + 1
log_pred <- array(NA, dim = c(t_pred, K, max_rolling))
bilog_pred <- array(NA, dim = c(t_pred, K*(K-1)*0.5, max_rolling))
MSFE <- array(NA, dim = c(t_pred, K, max_rolling))
MAFE <- array(NA, dim = c(t_pred, K, max_rolling))
PIT <- array(NA, dim = c(t_pred, K, max_rolling))
if (is.null(numCores)) numCores <- parallel::detectCores()*0.5
out_recursive <- parallel::mclapply(1:max_rolling,
FUN = function(j) {
time_current <- j + t_start - 1
y_current <- y[c(1:time_current), ]
prior <- get_prior(tail(y_current, time_current - p), p = p, dist=dist, SV = SV)
inits <- get_init(prior, samples = 15000, burnin = 5000, thin = 1)
if (SV) {
Chain <- BVAR.SV(y = tail(y_current, time_current - p), K = K, p = p, dist = dist,
y0 = head(y_current, p), prior = prior, inits = inits)
} else {
Chain <- BVAR.novol(y = tail(y_current, time_current - p), K = K, p = p, dist = dist,
y0 = head(y_current, p), prior = prior, inits = inits)
}
y_obs_future <- y[c((time_current+1):(time_current+t_pred)), ]
forecast_err <- forecast_density(Chain = Chain, y_obs_future = y_obs_future)
list(time_id = j, forecast_err = forecast_err)
},
mc.cores = numCores)
for (forecast_element in out_recursive){
time_id = forecast_element$time_id
forecast_err = forecast_element$forecast_err
log_pred[,,time_id] <- forecast_err$log_pred
bilog_pred[,,time_id] <- forecast_err$bilog_pred
MSFE[,,time_id] <- forecast_err$MSFE
MAFE[,,time_id] <- forecast_err$MAFE
PIT[,,time_id] <- forecast_err$emp_CDF
}
mlog_pred <- apply(log_pred, MARGIN = c(1,2), FUN = mean)
mbLP <- apply(bilog_pred, MARGIN = c(1,2), FUN = mean)
mMSFE <- apply(MSFE, MARGIN = c(1,2), FUN = mean)
mMAFE <- apply(MAFE, MARGIN = c(1,2), FUN = mean)
return( list( mLP = mlog_pred,
mbLP = mbLP,
mMSFE = mMSFE,
mMAFE = mMAFE,
log_pred = log_pred,
bilog_pred = bilog_pred,
MSFE = MSFE,
MAFE = MAFE,
PIT = PIT))
}
#' @export
recursive_seperate <- function(y, t_start = 100, t_pred = 12, K, p, dist = "MST", SV = T, outname = NULL, samples = 30000, burnin = 10000, thin = 1){
t_max = nrow(y)
if (is.null(outname)) outname = paste("Recursive_", dist, "_", SV, "_", t_start+1, ".RData", sep = "")
time_current <- t_start # index of the longer y
t_current <- time_current - p # length of the shorter y
y_current <- matrix(y[c(1:time_current), ], ncol = K)
prior <- get_prior(tail(y_current, time_current - p), p = p, dist=dist, SV = SV)
inits <- get_init(prior, samples = samples, burnin = burnin, thin = thin)
if (SV) {
Chain <- BVAR.SV(y = tail(y_current, time_current - p), K = K, p = p, dist = dist,
y0 = head(y_current, p), prior = prior, inits = inits)
} else {
Chain <- BVAR.novol(y = tail(y_current, time_current - p), K = K, p = p, dist = dist,
y0 = head(y_current, p), prior = prior, inits = inits)
}
y_obs_future <- matrix(y[c((time_current+1):(time_current+t_pred)), ],ncol = K)
forecast_err <- forecast_density(Chain = Chain, y_obs_future = y_obs_future, t_current = t_current)
out_recursive <- list(time_id = time_current+1,
forecast_err = forecast_err,
dist = dist,
SV = SV,
y_current = y_current,
y_obs_future = y_obs_future)
save(out_recursive, file = outname)
return( out_recursive)
}
# recursive_forecast <- function(y, t_start = 100, t_pred = 12, K, p, dist = "MST", SV = T, core = NULL){
# t_max = nrow(y)
# max_rolling <- t_max - (t_start + t_pred) + 1
#
# log_pred <- array(NA, dim = c(t_pred, K, max_rolling))
# bilog_pred <- array(NA, dim = c(t_pred, K*(K-1)*0.5, max_rolling))
# MSFE <- array(NA, dim = c(t_pred, K, max_rolling))
# MAFE <- array(NA, dim = c(t_pred, K, max_rolling))
# PIT <- array(NA, dim = c(t_pred, K, max_rolling))
#
# for (time_id in c(1:max_rolling)){
# if (time_id %% 10 == 0 ) cat("At rolling ", time_id, " \t")
# time_current <- time_id + t_start - 1
# y_current <- y[c(1:time_current), ]
# prior <- get_prior(tail(y_current, time_current - p), p = p, dist=dist, SV = SV)
# inits <- get_init(prior, samples = 15000, burnin = 5000, thin = 1)
# if (SV) {
# Chain <- BVAR.SV(y = tail(y_current, time_current - p), K = K, p = p, dist = dist,
# y0 = head(y_current, p), prior = prior, inits = inits)
# } else {
# Chain <- BVAR.novol(y = tail(y_current, time_current - p), K = K, p = p, dist = dist,
# y0 = head(y_current, p), prior = prior, inits = inits)
#
# }
# y_obs_future <- y[c((time_current+1):(time_current+t_pred)), ]
# forecast_err <- forecast_density(Chain = Chain, y_obs_future = y_obs_future)
# log_pred[,,time_id] <- forecast_err$log_pred
# bilog_pred[,,time_id] <- forecast_err$bilog_pred
# MSFE[,,time_id] <- forecast_err$MSFE
# MAFE[,,time_id] <- forecast_err$MAFE
# PIT[,,time_id] <- forecast_err$emp_CDF
# }
# mlog_pred <- apply(log_pred, MARGIN = c(1,2), FUN = mean)
# mbLP <- apply(bilog_pred, MARGIN = c(1,2), FUN = mean)
# mMSFE <- apply(MSFE, MARGIN = c(1,2), FUN = mean)
# mMAFE <- apply(MAFE, MARGIN = c(1,2), FUN = mean)
#
# return( list( mLP = mlog_pred,
# mbLP = mbLP,
# mMSFE = mMSFE,
# mMAFE = mMAFE,
# log_pred = log_pred,
# bilog_pred = bilog_pred,
# MSFE = MSFE,
# MAFE = MAFE,
# PIT = PIT))
# }
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.