Rolling forecasts

[AEdlerfi]

Hi,

I'd like to save/access the output of forecast draws at different time points. I think the model output only has the fc.ydraws etc for the last time period. Is there any way to save these output as the model is run?

Essentially, I want to do a 5 year rolling forecast starting from the start of estimation to the end.

Thanks

Adam

[pollytatouin]

Hi @AEdlerfi ,

I am myself new to the package, but I've been looking into it for the past couple of hours (hence why I'm here) and perhaps I can share with you some useful thoughts.

Is my understanding correct that by rolling forecasts you don't mean that you want to re-fit the model for every window, but rather you want to use the betas for a particular point in time and compute forecasts based on these?

If so, are you working on replicating Lubik and Matthes (2015)? This is what I'm exploring and my understanding is that they indeed take forecasts 5 years ahead at each point in time, using the relevant betas, but only estimate the model once.

I think this can be achieved using the Beta.postmean object within a loop.

Happy to discuss further

[AEdlerfi]

Hi @pollytatouin,

Yes that's correct, I'm applying the LM framework. What you describe is what I have done, but I don't think there's a method to derive the uncertainty around the mean estimate. I originally thought I could take the distribution of betas at each point and use them to construct a upper and lower bound of the forecast, but this resulted in odd results. So I've simply used a frequentist confidence interval assuming normality, which is not 100% correct given the model.

My thinking was to access the forecast distribution at each interval, but it's not saved as the model is estimated. I might see if I can alter the code and run it past @FK83 to see if he approves. That is, if what I suggest actually solves my problem.

Have you thought on this?

Also, I'd be happy to share my work so far with you.

[pollytatouin]

Yes that sounds like what we should do. I've investigated the problem a bit further and here are my thoughts:

Taking the vignette's forecasting example (that is, with usmacro.update), the estimation function provides us with 10 forecasts, of each 3 variables, for the 5000 retained iterations. In our case, we would need an object that's 10 forecasts, of each 3 variables, for 5000 iterations, at each 206 points in time. This goes for the fc.mdraws and fc.vdraws objects (not sure about fc.ydraws).

Equipped with those, we can then tweak the predictive.density function to select a particular "point in time" dimension:

predictive.density <-
function (fit, v = 1, h = 1, cdf = FALSE, point-in-time = 1) 
{
    nv <- length(fit$fc.ydraws[, 1, 1, 1]) # these are just to check lengths, so we can set "1"
    nh <- length(fit$fc.ydraws[1, , 1, 1])
    if (v > nv | h > nh | v != round(v) | h != round(h)) 
        stop("Please choose appropriate variable and horizon indices")
    v.ind <- vels(nv)[v]
    m <- fit$fc.mdraws[v, h, , point-in-time]
    s <- sqrt(fit$fc.vdraws[v.ind, h, , point-in-time])
    if (cdf == FALSE) {
        return(function(q) sapply(q, function(o) mean(dnorm(o, 
            mean = m, sd = s))))
    }
    else {
        return(function(q) sapply(q, function(o) mean(pnorm(o, 
            mean = m, sd = s))))
    }
}

Now, to get these objects in the first place, we need to tweak the estimation function. Trying to identify which part of it is relevant and assuming the parameter pdrift = 1, I think the work would happen here:

            if (pdrift == TRUE) {
                tempfc <- getfcsts(lastcol(Btdraw), lastcol(Atdraw), 
                  lastcol(Sigtdraw), Qdraw, Sdraw, Wdraw, t(y), 
                  nf, p)
                fc.m[, , aux.ct] <- tempfc$mean
                fc.v[, , aux.ct] <- tempfc$variance
                fc.y[, , aux.ct] <- tempfc$draw
            }

where instead of supplying the lastcol of the objects, we would loop across all columns of them. I'm not certain how to treat the objects that are taken as a whole (Qdraw, Sdraw...). Maybe we leave them as they are?

Please let me know your thoughts on this. Polly

[AEdlerfi]

Hi @pollytatouin ,

I think that all makes sense. There may be an issue with the amount of data being saved. Perhaps a more elegant solution would be to tweak the above to generate the required forecast and store that.

I'm pretty busy but should have some time over Christmas to dig into this further.

Adam

[FK83]

Hi @AEdlerfi and @pollytatouin,

thanks for your comments. The function get_all_forecasts in the following code chunk computes forecasts for all time periods, based on pdrift = FALSE(that is, the parameters stay constant over the forecast horizon). Please have a look and let me know if the results make sense to you.

Best regards, Fabian

rm(list = ls())
set.seed(20191216)
library(bvarsv)
data(usmacro.update)

# Estimate model
fit <- bvar.sv.tvp(usmacro.update, nrep = 10000, pdrift = FALSE)

# Helper function (select lower triangular matrix elements)
sel_lower_tri <- function(x){
  x[lower.tri(x)]
}

# Y is the data matrix
# fit is an object obtained from applying bvar.sv.tvp to Y
# nf is the maximal forecast horizon 
get_all_forecasts <- function(Y, fit, nf = 10){
  # Get parameter draws
  A_draws <- fit$A.draws
  Beta_draws <- fit$Beta.draws
  # Nr of MC draws
  mc_draws <- dim(Beta_draws)[3]
  # Lag length
  p <- fit$p 
  # Nr of variables in VAR
  M <- ncol(Y)
  # Reconstruct length of training period
  tau <- nrow(Y) - dim(Beta_draws)[2] - p
  # Effective sample size
  # note that data in Y and draws in parameter arrays are (tau+p) periods apart
  # (no parameter draws for (tau+p) pre-sample periods)
  t <- nrow(Y) - tau - p
  # Get non-redundant elements of fit$logs2.draws 
  logs2_draws <- fit$logs2.draws[,1:t,]
  # Output array for forecasts
  fc_m <- fc_y <- array(NA, c(M, t, nf, mc_draws))
  fc_v <- array(0,c(M*(M+1)*0.5,t, nf, mc_draws))
  # loop over time
  for (jj in 1:t){
    # Conditioning data
    fcdat <- Y[(tau+jj+1):(tau+jj+p), , drop = FALSE]
    # loop over MC draws
    for (ii in 1:mc_draws){
      # Parameters for current time period and MC draw
      parmat <- bvarsv:::beta.reshape(fit$Beta.draws[,jj,ii], M, p)       
      Atdraw_tmp <- sel_lower_tri(A_draws[,((jj-1)*M+1):(jj*M),ii])
      Sigtdraw_tmp <- logs2_draws[,jj,ii]
      auxb <- solve(bvarsv:::sigmahelper1(matrix(Atdraw_tmp, ncol = 1), M)) * 
        diag(exp(0.5*Sigtdraw_tmp))
      # Forecast variance
      tmpvar <- auxb %*% t(auxb) 
      # Loop over forecast horizons
      for (hhh in 1:nf){
        auxl <- bvarsv:::varfcst(parmat, tmpvar, fcdat, hhh)
        fc_m[,jj,hhh,ii] <- auxl$mean
        fc_v[,jj,hhh,ii] <- bvarsv:::vechC(auxl$variance)
        fc_y[,jj,hhh,ii] <- bvarsv:::mvndrawC(auxl$mean, auxl$variance)
      }      
    }
  }
  list(fc_m = fc_m, fc_v = fc_v, fc_y = fc_y)
}

# Apply function
fc_all <- get_all_forecasts(usmacro.update, fit)

# Last time period
t_sel <- dim(fit$Beta.postmean)[3]
# Check consistency with fit for last time period
all.equal(fc_all$fc_m[, t_sel,,], fit$fc.mdraws)
all.equal(fc_all$fc_v[, t_sel,,], fit$fc.vdraws)
library(magrittr)
fc_all$fc_y[1,t_sel,,] %>% t %>% 
  (function(x) apply(x, 2, function(z) quantile(z, c(.05, .95))))
fit$fc.ydraws[1,,] %>% t %>% 
  (function(x) apply(x, 2, function(z) quantile(z, c(.05, .95))))

[pollytatouin]

Hi @FK83 Fabian.

Thank you very much for coming back to us and providing this code snippet. I did a few tests this morning and it seems to be working. I also think that setting pdrift = FALSE is indeed more appropriate for this exercise.

One question: on the very last steps, where you compute the confidence regions, I'm getting different results between the two objects, despite setting a seed. The differences are not huge so I don't mind it, I just wanted to check if this is the expected outcome?

[FK83]

Hi @pollytatouin,

small differences in the very last steps are expected since the respective draws in fit and in fc_all are based on different calls of the function bvarsv:::mvndrawC.

Best regards, Fabian

[AEdlerfi]

Thank you so much @FK83 !!

[pollytatouin]

Hi @FK83 ,

Thank you for the details and again for taking the time to provide guidance! I now have a question about the way the forecasts are computed for nf > 1.

Looking at the snippet you provided as well as the bvar.sv.tvp function (pdrift = F), it seems to me that the forecasts are not computed recursively. In other words, the forecast at h=2 does not depend on the forecast at h=1. I'm saying that because the loop over forecasts horizons does not seem to be using outputs from the loop to calculate the next round.

If that's correct, then how would a forecast at, say, h=20 be computed? I think that for the LM framework to be valid, the process has to be recursive, such that the forecasts converge to a given point.

What prompted this questions is that I get forecasts at h=20 that are much more volatile than in LM.

Thanks you very much.

[FK83]

Hi @pollytatouin,

the forecasts could be computed recursively, but using the function bvarsv:::varfcst (as done in the code) is equivalent. The speed of convergence of the forecasts (in the forecast horizon h) depends on the persistence of the model, and the TVP-VAR can generate very persistent parameter draws.

Best, Fabian