MA Component: Add Bayesian Structure Time Series Model into rapbf

Design

Overview

Add simple bayesian structure time series model into pick best forecast component.

Detailed behavior

The mathematical models behind the approach are:

• Observation Equation

\[Time Series = Trend + Seasonality + Regression + Noise\] \[Y_t = U_t + S_t + \beta^T X_t+ \epsilon_t\]

• State/Transition/Process Equation

\[U_t=U_{t-1}+\delta_{t-1}+\omega_t\] \[\delta_t=\delta_{t-1}+v_t\] \[S_t=\sum_{s=1}^{s-1} S_{t-s}+\gamma_t\] Where \(\epsilon_t\), \(\omega_t\), \(v_t\) and \(\gamma_t\) are independent components of Gaussian random noise. \(U_t\) is the current level of the trend, the current slop of the trend is \(\delta_t\). The seasonal component \(S_t\) can be thought of as a set of \(s\) of seasons. \(X_t\) is the external regressors and \(\beta^T\) is the regression coefficients.

If the univariate time series data set is used for forecasting, the Bayesian structural model can be simplified as follow, since the series doesn’t have any regressors.

• Observation Equation

\[Time Series = Trend + Seasonality + Noise\] \[Y_t = U_t + S_t + \epsilon_t\]

• State/Transition/Process Equation

\[U_t=U_{t-1}+\delta_{t-1}+\omega_t\] \[\delta_t=\delta_{t-1}+v_t\] \[S_t=\sum_{s=1}^{s-1} S_{t-s}+\gamma_t\]

Considerations

• The methodology in the prototype is simple Bayesian time series model for the univariate time series forecasting.
  • The pros of the prototyped model is that Bayesian time series models are more transparent than ARIMA model. They also facilitate better handling of uncertainty, a key feature when planning for the future. Bayesian time series model fits perfectly with sequential learning and decision making and it directly leads to exact small sample results.
  • The cons of the prototyped model is that the partial seasonality will be threw into noise/error term. The accurary of the forecasting results might be slighly worse than ARIMA model for some cases.
• What are not in the prototype ?
  • The model with regressor (i.e. external variable) is not considered.
  • The input data is partial year data is not considered.
  • Since there is no project using Bayesian time series model at RA, the model is not covered in any Atlas Toolkit so far.

Function Overview

Function

rabsts(data, season.num, season.dur, pred.horizon)

Input

•       data - Input data.

•       season.num - The number of season to be modeled. 
        i.e., for a time series with quarterly seasonality, season.num=4.

•       season.dur - The number of time periods in each season. 
        i.e., when season.num=4, 
        if it is monthly data, then season.dur=3;                
        if it is weekly data, then season.dur=13.

•       pred.horizon - The number of periods you wish to predict. 

  Note: seasonal component (season.num and season.dur) can be determined by external analysis.

Output

A data frame containing the following as columns:

•       result$mean - The posterior mean of the prediction; 

•       result$median - The posterior median of the prediction.

Required R Package

•       bsts

Prototype Code in R

Function Definition

## BEGIN: function ##
rabsts <- function(data, season.num, season.dur, pred.horizon) {
  ##get trend or/and seasonal state
  if(season.num==0){
    ss <- bsts::AddLocalLinearTrend(list(), data)}
  else{
    ss <- bsts::AddLocalLinearTrend(list(), data)
    ss <- bsts::AddSeasonal(ss, data, nseasons = season.num, season.duration = season.dur)}
  
  ##build Bayesial model
  bsts.model <- bsts::bsts(data, state.specification = ss, niter = 666, ping=0, seed=1000)
  
  ##predict
  result<-bsts::predict.bsts(bsts.model, horizon = pred.horizon, burn = SuggestBurn(0.1, bsts.model), quantiles = c(.025, .975))
  pred<-data.frame(result$mean,result$median)
  return(pred)
}
## END: function ##

Example of calling the function with sample data

First, we need to load the required R packages.

suppressMessages(library(bsts))

Now, we need some input data. The unit test data sets for rapbf package is hired for the prototype. It looks like:

Segment the data by seriesId:

#Input data
dat <- subset(inputData, seriesId=="1")
dat <- droplevels(dat)
#Holdout data
holdout.dat<-subset(holdoutData, seriesId=="1")
holdout.dat <- droplevels(holdout.dat)

Apply the Bayesial time series model for forecast:

• Seasonal series:

x<-rabsts(dat$value,season.num=12,season.dur=1,pred.horizon=3) #assume that there are 12 seasons and each season is one month.
data.frame(x,holdout.dat$value)

• Non-seasonal series:

y<-rabsts(dat$value,season.num=0,pred.horizon=3) #If no season, set season.num=0.
data.frame(y,holdout.dat$value)