Joint work with

Yanfei Kang	Rob J Hyndman	Feng Li
Beihang University	Monash University	Central University of Finance and Economics

Motivation

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

  • hourly electricity demands
  • daily maximum temperatures

• 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

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

Why focus on ARIMA models

ARIMA models

• the most widely used statistical time series forecasting approaches.

• frequently serve as benchmark methods for model combination.

• handle non-stationary series via differencing and seasonal patterns.

• the automatic ARIMA modeling was developed to easily implement the order selection process (Hyndman & Khandakar, 2008).

• be converted into AR representations (linear form).

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}$$

• 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 (Distributed Least Squares Approximation) has been proved (Zhu, Li & Wang 2021, JCGS).

• 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.

Output forecasts

• The $h$-step-ahead point 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}$$

• 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.$$

  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$.

Experimental setup

• Number of subseries: $150$

  • the length of each subseries is 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

• The rationality of setting the AR order to $2000$.

• DARIMA always outperforms the benchmark method regardless of point forecasts or prediction intervals.

Distributed forecasting results

• If long-term observations are considered, DARIMA is preferable, especially for interval forecasting.

• The achieved performance improvements of DARIMA become more pronounced as the forecast horizon increases.

Distributed forecasting results

Distributed forecasting results

• Our approach has captured the decreasing yearly seasonal trend.

• Both DARIMA and ARIMA have captured the hourly seasonality, while DARIMA results in forecasts closer to the true future values than ARIMA.

Distributed forecasting results

Conclusions

• 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.

• Spark implementation: @xqnwang/darima

• Website: https://xqnwang.rbind.io

• Twitter: @Xia0qianWang