Background

• Ultra-long time series are increasingly accumulated in many cases.

  • hourly electricity demands
  • daily maximum temperatures
  • streaming data generated in real-time

• Forecasting these time series is challenging.

  • time-consuming training process
  • hardware requirements
  • unrealistic assumption that the DGP remains invariant over a long time interval

• Some attempts are made in the vast literature.

  • discard the earliest observations
  • allow the model itself to evolve over time
  • apply a model-free prediction
  • develop methods using the Spark’s MLlib library

Electricity load data

• The Global Energy Forecasting Competition 2017 (GEFCom2017)

• Hourly electricity load data of zones spanning New England

  • $8$ bottom zones & $2$ aggregated zones

  • without considering the hierarchical configuration

• Ranging from March 1, 2003 to April 30, 2017 ( $124,171$ time points)

  • training periods
    March 1, 2003 - December 31, 2016

  • test periods
    January 1, 2017 - April 30, 2017 ( $h=2,879$ )

Electricity load data

Distributed forecasting

• Zhu, Li & Wang (2019) tackle regression problems on distributed systems by developing a distributed least squares approximation (DLSA) method.

• Local estimators are computed by worker nodes in a parallel manner.

• Local estimators are delivered to the master node to approximate global estimators by taking a weighted average.

Distributed forecasting

Parameter estimation problem

• For an ultra-long time series $\{y_t;t=1,2,\dots ,T\}$. Define $S=\{1,2,\cdots ,T\}$ to be the timestamp sequence.

• The parameter estimation problem can be formulated as $$f(\theta ,\Sigma |y_t,t\in S).$$

• Suppose the whole time series is split into $K$ subseries with contiguous time intervals, that is $S={\cup }_{k=1}^K{S}_k$.

• The parameter estimation problem is transformed into $K$ sub-problems and one combination problem as follows: $$f(\theta ,\Sigma |y_t,t\in S)=g(f_1({\theta }_1,{\Sigma }_1|y_t,t\in S_1),\dots ,f_K({\theta }_K,{\Sigma }_K|y_t,t\in S_K)).$$

Distributed forecasting

Distributed forecasting

Why focus on ARIMA models

Advantages

• The most widely used statistical time series approachs.

• ARIMA models can handle non-stationary series via differencing and seasonal patterns

• ARIMA models frequently serve as benchmark methods for model combination because of their excellent performance in forecasting.

• ARIMA models can be converted into AR representations (linear form).

• Hyndman & Khandakar (2008) developed the automatic time series forecasting with ARIMA models to easily implement the order selection process.

Limitations

• The likelihood function of such time series model could hardly scale up.

Automatic ARIMA modeling

AR representation

• A seasonal ARIMA model is generally defined as $$\begin{array}{c}(1-\sum _{i=1}^p{\varphi }_i{B}^i)(1-\sum _{i=1}^P{\Phi }_i{B}^{im})(1-B{)}^d{(1-B^m)}^D(y_t-{\mu }_0-{\mu }_{1}t)\\ =(1+\sum _{i=1}^q{\theta }_i{B}^i)(1+\sum _{i=1}^Q{\Theta }_i{B}^{im}){\epsilon }_t.\end{array}$$

• Let the term $y_t-{\mu }_0-{\mu }_{1}t$ be denoted as $x_t$.

• By utilizing the polynomial multiplication, the seasonal ARIMA model is converted into a non-seasonal ARMA( $u,v$ ) model (possibly non-stationary) $$\begin{array}{c}(1-\sum _{i=1}^u{\varphi }_i^{\prime }{B}^i)x_t=(1+\sum _{i=1}^v{\theta }_i^{\prime }{B}^i){\epsilon }_t.\end{array}$$

AR representation

• The AR representation of the ARMA( $u,v$ ) can be obtained by long division of AR and MA polynomials.

• Given two polynomials ${\varphi }^{\prime }\left(B\right)=(1-{\sum }_{i=1}^u{\varphi }_i^{\prime }{B}^i)$ and ${\theta }^{\prime }\left(B\right)=(1+{\sum }_{i=1}^v{\theta }_i^{\prime }{B}^i)$, we have $$\pi \left(B\right)x_t=\frac{{\varphi }^{\prime }\left(B\right)}{{\theta }^{\prime }\left(B\right)}{x}_t={\epsilon }_t,$$ where $\pi \left(B\right)=(1-{\sum }_{i=1}^{\infty }{\pi }_i{B}^i)$.

• The linear representation of the original seasonal ARIMA model can be given by $$y_t={\beta }_0+{\beta }_{1}t+\sum _{i=1}^{\infty }{\pi }_i{y}_{t-i}+{\epsilon }_t,$$ where $${\beta }_0={\mu }_0(1-\sum _{i=1}^{\infty }{\pi }_i)+{\mu }_1\sum _{i=1}^{\infty }i{\pi }_i\text{and}{\beta }_1={\mu }_1(1-\sum _{i=1}^{\infty }{\pi }_i).$$

Estimators combination

• Some excellent statistical properties of the global estimator obtained by DLSA has been proved (Zhu, Li & Wang, 2019).

• We extend the DLSA method to solve time series modeling problem.

• Define $L(\theta ;y_t)$ to be a second order differentiable loss function, we have $$\begin{array}{cc}L\left(\theta \right)& =\frac{1}{T}\sum _{k=1}^K\sum _{t\in S_k}L(\theta ;y_t)\\ & =\frac{1}{T}\sum _{k=1}^K\sum _{t\in S_k}\{L(\theta ;y_t)-L({\hat{\theta }}_k;y_t)\}+c_1\\ & \approx \frac{1}{T}\sum _{k=1}^K\sum _{t\in S_k}{(\theta -{\hat{\theta }}_k)}^{\top }\ddot{L}({\hat{\theta }}_k;y_t)(\theta -{\hat{\theta }}_k)+c_2\\ & \approx \sum _{k=1}^K{(\theta -{\hat{\theta }}_k)}^{\top }\left(\frac{T_k}{T}{\hat{\Sigma }}_k^{-1}\right)(\theta -{\hat{\theta }}_k)+c_2.\end{array}$$

Estimators combination

• The global estimator calculated by minimizing the weighted least squares takes the following form $$\begin{array}{cc}\tilde{\theta }& ={\left(\sum _{k=1}^K\frac{T_k}{T}{\hat{\Sigma }}_k^{-1}\right)}^{-1}\left(\sum _{k=1}^K\frac{T_k}{T}{\hat{\Sigma }}_k^{-1}{\hat{\theta }}_k\right),\\ \tilde{\Sigma }& ={\left(\sum _{k=1}^K\frac{T_k}{T}{\hat{\Sigma }}_k^{-1}\right)}^{-1}.\end{array}$$

• ${\hat{\Sigma }}_k$ is not known and has to be estimated.

• We approximate a GLS estimator by an OLS estimator (e.g., Hyndman et al., 2011) while still obtaining consistency.

• We consider approximating ${\hat{\Sigma }}_k$ by ${\hat{\sigma }}_k^{2}I$ for each subseries.

Point forecasts

The $h$-step-ahead forecast can be calculated as $${\hat{y}}_{T+h|T}={\tilde{\beta }}_0+{\tilde{\beta }}_1(T+h)+\{\begin{array}{cc}{\sum }_{i=1}^p{\tilde{\pi }}_i{y}_{T+1-i},& h=1\\ {\sum }_{i=1}^{h-1}{\tilde{\pi }}_i{\hat{y}}_{T+h-i|T}+{\sum }_{i=h}^p{\tilde{\pi }}_i{y}_{T+h-i},& 1<h<p.\\ {\sum }_{i=1}^p{\tilde{\pi }}_i{\hat{y}}_{T+h-i|T},& h\ge p\end{array}$$

Prediction intervals

• In the global model, the standard error of $h$-step ahead forecast is formally expressed as $${\tilde{\sigma }}_h^2=\{\begin{array}{cc}{\tilde{\sigma }}^2,& h=1\\ {\tilde{\sigma }}^2(1+{\sum }_{i=1}^{h-1}{\tilde{\theta }}_i^2),& h>1\end{array},$$ where ${\tilde{\sigma }}^2=\mathrm{tr}\left(\tilde{\Sigma }\right)/p$.

• The central $(1-\alpha )\times 100\%$ prediction interval of $h$-step ahead forecast can be defined as $${\hat{y}}_{T+h|T}\pm {\Phi }^{-1}(1-\alpha /2){\tilde{\sigma }}_h.$$

Experimental setup

• Number of subseries: $150$

  • the length of each subseries about $800$
  • the hourly time series in M4 ranges from $700$ to $900$
  • the time consumed by automatic ARIMA modeling process is within 5 minutes

• AR order: $2000$

• The experiments are carried out on a Spark-on-YARN cluster

  • one master node and two worker nodes
  • Each node contains 32 logical cores, 64 GB RAM and two 80GB SSD local hard drives

Distributed forecasting results

Distributed forecasting results

Distributed forecasting results

Distributed forecasting results

Distributed forecasting results

Sensitivity analysis

Number of split subseries

Sensitivity analysis

Maximum values of model orders

Summary

• A distributed time series forecasting framework using the industry-standard MapReduce framework.

• The local estimators trained on each subseries are combined using weighted least squares to minimize a global loss function.

• Our framework

  • works better than competing methods for long-term forecasting.

  • achieves improved computational efficiency in optimizing the model parameters.

  • allows that the DGP of each subseries could vary.

  • can be viewed as a model combination approach.

Thanks!