A Composite Quantile Fourier Neural Network for Multi-Step Probabilistic Forecasting of Nonstationary Univariate Time Series
Kostas Hatalis Shalinee Kishore
stat.ML
































kmh511@lehigh.edu


Department of Electrical and Computer Engineering, 
Lehigh University, Bethlehem, PA 18105, USA

[label1]The research reported in this manuscript was supported by the National Science Foundation; grant number 1442858
[cor1]Corresponding author




	Point forecasting of univariate time series is a challenging problem with extensive work having been conducted. However, nonparametric probabilistic forecasting of time series, such as in the form of quantiles or prediction intervals is an even more challenging problem. In an effort to expand the possible forecasting paradigms we devise and explore an extrapolation-based approach that has not been applied before for probabilistic forecasting. We present a novel quantile Fourier neural network is for nonparametric probabilistic forecasting of univariate time series. Multi-step predictions are provided in the form of composite quantiles using time as the only input to the model. This effectively is a form of extrapolation based nonlinear quantile regression applied for forecasting. Experiments are conducted on eight real world datasets that demonstrate a variety of periodic and aperiodic patterns. Nine naive and advanced methods are used as benchmarks including quantile regression neural network, support vector quantile regression, SARIMA, and exponential smoothing. The obtained empirical results validate the effectiveness of the proposed method in providing high quality and accurate probabilistic predictions. 




probabilistic forecasting artificial neural networks time series multi-step prediction quantile regression









§ INTRODUCTION


Univariate time series based deterministic or point forecasting is a well-studied field that has numerous usages. Examples of applications include finance <cit.>, topic behavior <cit.>, traffic flow <cit.>, and renewable power <cit.>. There are several approaches to forecasting with different classes of methods. Approaches include having a sliding window of past data to predict future data, recurrent models, and extrapolation based regression such as signal approximation <cit.>. In all the approaches, methods are divided into two classes as linear or nonlinear. In the first class, methods include linear regression, AutoRegression (AR), AutoRegressive–Moving-Average (ARMA) models, and exponential smoothing. The second class of methods is nonlinear models which are predominantly machine learning based such as support vector regression, nonlinear autoregression neural networks, and recurrent neural networks. A thorough overview of time series and machine learning based deterministic forecasting can be found in <cit.>.

Despite the wide use of deterministic forecasting, it does have a significant disadvantage in that it can result in individual observation errors which can be significant. Additionally, deterministic forecasting lacks information on associated uncertainty. A solution to these problems is Probabilistic Forecasting (PF) where the goal is to produce fully probabilistic predictions that derive quantitative information on the uncertainty. A PF takes the form of a predictive probability distribution over future time horizons and aims to maximize the sharpness of predictive densities while subject to reliability. Sharpness refers to the concentration of the predictive distribution and reliability refers to the accuracy of the forecasted probability in conveying the actual probability of events.

A popular application of PF is in the fields of renewable energies and power systems. A probabilistic forecast is vital for different operations to renewable energy farms. This includes managing the optimal level of generation reserves <cit.>, optimizing production <cit.>, and bidding strategies for electricity markets <cit.>. Applications to the power grid include load analysis <cit.>, smart meters <cit.>, scheduling <cit.>, system planning <cit.>, unit commitment <cit.>, and energy trading <cit.>. A thorough overview of probabilistic wind and solar power forecasting is provided in <cit.> and <cit.>. Further popularity of PF in energy systems can be seen by the large number of submissions in the Global Energy Forecasting Competition of 2014 <cit.>.

Additional motivation for this work comes from the lack of studies conducted for nonparametric PF of univariate time series, particularly for multi-step prediction. Most nonparametric PF work utilize a multivariate regression based framework, while most univariate time series PFs are made using parametric assumptions, for instance Gaussian intervals. There is such little understanding of prediction intervals for time series forecasting, that the popular M4 Competition for time series forecasting <cit.> included a PF track for producing prediction interval forecasts specifically to study their performance in more detail. There own findings demonstrate that standard forecasting methods fail to estimate the uncertainty in time series forecasts properly, and therefore that more research is required in time series PF.

There are several main classes in the type of PF models which include if they are parametric or nonparametric, direct or indirect, and the type of inputs they use for prediction. In PF, we are first trying to predict one of two types of density functions, either parametric or nonparametric. When the future density function is assumed to take a specific distribution, such as the Normal distribution, then this is known as parametric probabilistic forecasting. For processes where no assumption is made about the shape of the distribution, a nonparametric probabilistic forecast can be made. Nonparametric predictions can be made in the form of quantiles, prediction intervals, or full density functions. For example, nonlinear and non-stationary data, such as wind speeds or stocks, may not correspond to fixed or known distributions. When in need of forecasting such data it can be more beneficial to apply a nonparametric probabilistic forecast to estimate the distribution instead then assume it is shaped.

The second classification of PF models is whether they are combined with point forecasts or not. For instance in <cit.> a deterministic and probabilistic forecast for wind power is combined. This approach is known as an indirect PF method. First, a point forecast is made such as with support vector regression, and then prediction intervals for point forecast values are obtained with a PF method such as quantile regression. On the other hand, when a PF method estimates future quantiles or prediction intervals without using as input point forecasts, this is known as direct forecasting. The last distinction to be made with PF models is if past lagged data are used as inputs to the forecast model or if future exogenous variables are used too. For instance in renewable power PF, if Numerical Weather Predictions (NWP) exist for each forecasting horizon that we are interested in, then those exogenous NWPs are used as inputs to provide a PF in that prediction horizon. When NWPs are not given then lagged past time values of renewable power can be used for prediction. 

We introduce a new approach to developing a nonparametric direct PF where the input to the model is neither exogenous variables nor past data but instead treats the series as a signal. This approach is motivated by Fourier extrapolation which is the process by which a Fourier transform is applied to a data set to decompose it into a sum of sinusoidal components thus interpreting it as a signal. In time series analysis this is related to Harmonic regression. In accounting for periodic and non-periodic aspects of a signal such as a trend, Fourier Neural Networks (FNN) have been proposed. 

FNNs are feedforward neural networks with sinusoidal activation functions that model the Fourier transform. While other machine learning methods such as convolutional or recurrent neural networks have shown success in time series forecasting <cit.>, these architectures are not well suited for modeling Fourier extrapolation. Thus, most recently, a new FNN called Neural Decomposition (ND) was proposed in <cit.> that can decompose a signal into a sum of its constituent parts, model trend, and reconstruct a signal beyond the training samples. ND can provide a prediction by having time as its only input similarly to an inverse Fourier transform. We propose a PF model motivated by the ND model. 

Several works have explored Fourier extrapolation based deterministic forecasting with sinusoidal neural networks <cit.>, but none have yet explored it for probabilistic forecasting. We are the first to introduce an FNN for forecasting composite quantiles that we dub the quantile Fourier neural network. The main contributions of our approach can be summarized as follows:


	
  * We show how to estimate composite quantiles using a Fourier neural network using a smoothed pinball loss function.
	
  * We demonstrate how this extrapolation based quantile forecasting is able to model periodic and non-periodic components of nonstationary time series.
	
  * We demonstrate how to conduct multi-step probabilistic forecasts in the form of quantiles or prediction intervals.
	
  * We design experiments to validate our approach for direct probabilistic forecasting and provide insight how this method can generalize modeling uncertainty on real-world datasets.


The contents of the paper are: in Section <ref> we provide the mathematical background on probabilistic forecasting, quantile regression, and evaluation methods. In Section <ref> we review existing architectures of FNNs, go over our model, its architecture, training, and weighting initialization scheme. Results and discussion of validating our method are presented in Section <ref>. We conclude the paper and review future research directions in Section <ref>.



§ PROBABILISTIC FORECASTING
 

This section highlights the underlying mathematics in probabilistic forecasting, overviews linear quantile regression which forms the foundation of our proposed method, and summarizes the main evaluation metric for density forecasts. Given a random variable Y_t at time t, its probability density function is defined as f_t and its the cumulative distribution function as F_t. If F_t is strictly increasing, the quantile q_t^(τ) at time t of the random variable Y_t with nominal proportion τ is uniquely defined on the value x such that P(Y_t < x) = τ or equivalently as the inverse of the distribution function q_t^(τ) = F_t^-1(τ). A quantile forecast q̂_t+z^(τ) is an estimate of the true quantile q_t+z^(τ) for the lead time t+z, given a predictor values. 

Prediction intervals are another type of probabilistic forecast and give a range of possible values within which an observed value is expected to lie with a certain probability β∈ [0,1]. A prediction interval Î^(β)_t+z produced at time t for future horizon t + z is defined by its lower and upper bounds, which are the quantile forecasts Î^(β)_t+z = [ q̂_t+z^(τ_l) ,q̂_t+z^(τ_u)] = [ l_t^(β),u_t^(β)] whose nominal proportions τ_l and τ_u are such that τ_u - τ_l = 1-β. A nonparametric probabilistic forecast f̂_t+z <cit.> can be made of the density function by gathering a set of M quantiles forecasts such that f̂_t+z = {q̂_t+z^(τ_m) ,m=1,...,M|0≤τ_1 < ... < τ_M ≤ 1 } with chosen nominal proportions spread on the unit interval.

Quantile regression is a popular approach for nonparametric probabilistic forecasting. Koenker and Bassett <cit.> introduce it for estimating conditional quantiles and is closely related to models for the conditional median <cit.>. Minimizing the mean absolute function leads to an estimate of the conditional median of a prediction. By applying asymmetric weights to errors through a tilted transformation of the absolute value function, we can compute the conditional quantiles of a predictive distribution. The selected transformation function is the pinball loss function as defined by

    ρ_τ(u) = {[    τ u  u ≥ 0; (τ-1)u  u < 0 ].,


where 0 < τ < 1 is the tilting parameter. The pinball loss function penalizes low-probability quantiles more for overestimation than for underestimation and vice versa in the case of high-probability quantiles. Given a vector of predictors X_t where t = 1,...,N, a vector of weights W and intercept b coefficient in a linear regression fashion, the conditional τ_i quantile, for i ∈ [0,1], is given by q̂_t^(τ_i ) = W^⊤ X_t +b. To determine estimates for the weights and intercept for composite quantile regression we solve the following minimization problem

    min_W,b1/NM∑_i=1^M∑_t=1^Nρ_τ_i (y_t-q̂_t^(τ_i)),


where y_t is the observed value of the predictand and M is the number of quantiles we are estimating. The formulation above in Eq. (<ref>) can be minimized by a linear program. Here are many variations of QR. In <cit.> local QR is applied to estimate different quantiles, while in <cit.> a spline-based QR is used to estimate quantiles of wind power. In <cit.> quantile loss gradient boosted machines are used to estimate many quantiles, and in <cit.> multiple quantile regression is used to predict a full distribution with optimization achieved by using the alternating direction method of multipliers. Quantile regression forests <cit.> are another approach in forecasting which are an extension of regression forests based on classification and regression trees. Due to their flexibility in modeling elaborate nonlinear data sets, artificial neural networks are another dominant class of machine learning algorithms that can be used to enhance QR. Taylor <cit.> implemented and demonstrated a quantile regression neural network (QRNN) method, combining the advantages of both QR and a neural network. In <cit.> an autoregressive version of QRNN is used for applications to evaluating value at risk, and <cit.> implements the QRNN model in R as a statistical package. We extend the QRNN model to form our quantile Fourier neural network.



 §.§ Evaluation Methods


In probabilistic forecasting it is essential to evaluate the quantile estimates and derived predictive intervals. To evaluate quantile predictions, one can use the pinball function directly for evaluation called the quantile score (QS). We choose QS as our main quantile measure for the following reasons. When averaged across many quantiles it can evaluate full predictive densities, it is found to be a proper scoring rule <cit.>, and it is related to the continuous rank probability score. QS calculated overall N test observations and M quantiles is defined as


    QS = ∑_t=1^N∑_m=1^Mρ_τ_m (y_t - q̂_t^(τ_m))


where y_t is an observation used to forecast evaluation. To evaluate full predictive densities, this score is averaged across all target quantiles for all look ahead time steps using equal weights. A lower QS indicates a better forecast. In some applications it may be needed to have prediction intervals (PIs) and as such, we look at two PI evaluation measures: reliability and sharpness. Reliability, also known as calibration, is a measure which states that over an evaluation set the observed and nominal probabilities should be as close as possible, and the empirical coverage should ideally equal the preassigned probability. Sharpness is a measure of the width of prediction intervals, defined as the difference between the upper u_t^β_i and lower l_t^β_i interval values. 

For interval reliability we use the average coverage error (ACE) metric <cit.> and for measuring interval sharpness we use the sharpness score (SS). For measuring reliability, PIs show where future observations are expected to lie, with an assigned probability termed as the PI nominal confidence (PINC) 100(1-β_i)%. Here i = 1...M/2 indicates a specific coverage level. The coverage probability of estimated PIs is expected to eventually reach a nominal level of confidence over the test data. A measure of reliability which shows target coverage of the PIs is the PI coverage probability (PICP), which is defined by

    PICP_i=1/N∑_t=1^N1{y_t ∈ I_t^β_i(x_t)}.


For reliable PIs, the examined PICP should be close to its corresponding PINC. A related and easier to visualize assessment index is the average coverage error (ACE), which is defined by

    ACE = ∑_ i = 1^M/2 |PICP_i-100(1-β_i)|.


This assumes calculation across all test data and coverage levels. To ensure PIs have high reliability, the ACE should be as close to zero as possible. A high reliability can be easily achieved by increasing or decreasing the distance between lower and upper interval bounds. Thus, the width of a PI can also influence its quality. For measuring the effective width of PIs we use the sharpness score proposed by <cit.> which measures how wide PIs are by focusing on the mean size of the intervals only. We define q̂^u_t-q̂^l_t as the size of the central interval forecast with nominal coverage rate (1 - β). For lead times t = 1...N_test, a measure of sharpness for PIs is then given by the mean size of the intervals

    SS=1/N_test∑_t=1^N_test(q̂^u_t-q̂^l_t).


A lower sharpness score is considered more ideal, but too small and the PIs would not cover enough of the observed data. Thus sharpness is typically a measure to be considered along with reliability and a skill score.








§ PROPOSED METHODOLOGY
 

Fourier analysis examines the approximation of functions through their decomposition as a sum or product of trigonometric functions, while Fourier synthesis focuses on the reconstruction of a signal from its decomposed oscillatory components. These well-known processes have a large utility in time series analysis. By decomposing a time series into its frequencies one could then interpolate missing time values by reconstructing the original signal. Further applications include modeling seasonality and even prediction of a time series through extrapolation of an approximated signal. In the application of Fourier analysis for time series analysis, an important method is the discrete Fourier transform (DFT), which converts a series into its frequency domain representation, and the inverse discrete Fourier transform (iDFT) which maps the frequency representation back to the time domain. The transforms can be expressed as either a summation of complex exponentials or sines and cosines by Euler’s formula. In this section we explore existing works on Fourier networks that directly use iDFT in their operation or mimic it, then we describe our proposed FNN methodology for quantile forecasting.



 §.§ Fourier Neural Networks


Neural networks with sine as an activation function are difficult to train in theory and when initialized randomly yield poor results <cit.>. Thus, few works have attempted to explore Fourier analysis with sinusoidal neural networks. We highlight most of the works here. One of the first FNNs was introduced by Adrian Silvescu <cit.> who developed Fourier-like neurons for learning boolean functions. The FNN model used the units of the network to approximate a DFT in its output. Similar in spirit to an FNN was a Fourier transform neural network introduced in <cit.> that uses the Fourier transform of the data as input to an artificial neural network. FNNs have since been used for stock prediction <cit.>, aircraft engine fault diagnostics <cit.>, harmonic analysis <cit.>, and extensions include a single input multiple outputs based FNNs that can turn nonlinear optimization problems into linear ones <cit.>, FNNs for output feedback learning control schemes <cit.>, and deep FNNs for lane departure prediction <cit.>.

There are two recent works that study FNNs for time series prediction that use the Fourier transform of the data as weights. The first is an FNN presented by Gashler and Ashmore in <cit.>. Their technique uses the fast Fourier transform (FFT) to approximate the DFT and then uses the obtained values to initialize the weights of the neural network. Their model uses a combination of sinusoid, linear, and softplus activation units for modeling periodic and non-periodic components of a time series. However, their trained models were slightly out of phase with their validation data. The second study on FNNs for time series prediction is presented by Godfrey and Gashler <cit.> who proposed a similar model to <cit.> called neural decomposition (ND), except that they do not use the Fourier transform to directly initialize any weights.  

The ND model is inspired by the inverse discreet Fourier transform where given time t as input it attempts to model the signal x(t). However, there are some distinctions between ND and iDFT. First ND allows sinusoid frequencies to be trained and second ND can also model non-periodic components in a signal such as trend. With the ability to train the frequencies, ND learns the actual period of a signal whereas iDFT assumes that the underlying function always has a period equal to the size of the samples it represents. ND is a feedforward neural network with a single hidden layer with N nodes and has one input and one output node. Hidden nodes are composed of sinusoid units for capturing the periodic component in an underlying signal and other activation functions, such as linear or sigmoid units, for capturing the non-period component. Parameters of ND are initialized in such a way so as to mimic the iDFT. ND is then trained with stochastic gradient descent with backpropagation. ND was applied to time series deterministic forecasting and showed very promising results across different data sets, often beating state-of-the-art methods such as LSTM, SVR, and SARIMA.



 §.§ Quantile Fourier Neural Network


Inspired by the ND model we propose a new forecasting method which we call the quantile Fourier neural network (QFNN).  Unlike ND and other FNNs our QFNN  model is trained to extrapolate quantiles of an underlying time series, rather than point estimates. We apply sinusoid activation functions to allow QFNN to fit periodic data, and coupled with an augmentation function it is able to probabilistically forecast a time series that is made up of non-periodic components too. The model is defined as follows. Let each a_jk represent an amplitude, each ω_k represents a frequency, each ϕ_k represents a phase shift, and  b^[2]_τ and b^[1]_k represent bias terms for the quantile signal representation. Let f(t) be an augmentation function that represents the non-periodic components of the signal. QFNN then can be defined by

    q_t^τ = f(t) + b^[2]_τ + ∑_k=1^N(  a_τ, k·cos( ω_k t + ϕ_k ) + b^[1]_k)


where given time t as the input, it attempts to predict the τ-level quantile. QFNN is loosely modeling a time series as a partial Fourier cosine series

    x(t) = A_0 + ∑_n=1^N A_n cos(n ω_0 t + ϕ_n )


where ω_0 = 2 π/T, T is the period of the signal x(t), and A_0, A_n, and ϕ_n are real numbers. The main difference between Eq. <ref> and Eq. <ref> is that QFNN does not fix the period of the signal to a predetermined size T, it allows for bias terms, it has an augmentation function to represent non-periodic components of the signal, and it learns the frequencies ω_k versus keeping them at a fixed size. The bias terms in the output layer of the network are important because it allows to shift the level of each quantile appropriately.

The hidden layer of QFNN is composed of N units with a sinusoid activation function and an arbitrary number of units with other activation functions to calculate f(t). The output layer is composed of M number of linear units that represent quantiles. The parameters a_jk, being the weights between the hidden and output layers allows us to model different amplitudes for composite quantiles while simultaneously learning the frequency and phases for all quantiles in the hidden layer. Utilizing conventional neural network notation W^[1] is a matrix of the f(t) unit parameters and the frequency parameters in Eq. 3. The b^[1] vector represents the phases of the sinusoidal components, W^[2] is a parameter matrix of the amplitudes, and we also add additional bias terms to the output nodes for each quantile with the b^[2] vector.

To estimate quantiles we need to solve the minimization problem described in Eq. <ref>. However, the pinball function ρ  in Eq. <ref>. is not differentiable at the origin, x = 0. The non-differentiability of ρ makes it difficult to apply gradient-based optimization methods in fitting the quantile regression model. Gradient-based methods are preferable for training neural networks since they are time efficient, easy to implement and can yield a local optimum. Therefore, we need a smooth approximation of the pinball function that allows for the direct application of gradient-based optimization. A smooth approximation of the pinball function in Eq. (<ref>)  is proposed by Zheng in <cit.> as

    S_τ,α(u) = τ u + αlog( 1 + exp(-u/α) ),


where  α>0 is a smoothing parameter and τ∈ [0,1] is the quantile level we're attempting to estimate. In Fig. <ref> we see the pinball function with τ = 0.5 as the red line and the a smooth approximation as the blue line with α = 0.2. Zheng proves <cit.> that in the limit as α→ 0^+ then S_τ,α(u) = ρ_τ(u). With this smooth approximation we can then define the cost minimization problem for QFNN as

    E = 1/NM∑_t=1^N∑_m = 1^M[ τ_m (y_t - q̂_t^(τ_m)) + αlog( 1 + exp( - y_t-q̂_t^(τ_m)/α) ) ].

where M number of τ's we are trying to estimate in the output layer. The input to hidden neurons is calculated, in vectorization notation, by Z_t^[1] =  W^[1] t + b^[1], the output of the hidden layer then uses the logistic activation function H_t = cos(  Z_t^[1]), f(  Z_t^[1]). The input to output neurons is then calculated by Z_t^[2] =  W^[2] H_t + b^[2], and the output layer uses the identity activation function Q̂_t =  Z_t^[2] where Q̂_t is a vector output of the estimated composite quantiles. An architectural view of the QFNN is shown in Fig. <ref>.



 §.§ Implementation Details


The proposed QFNN model is trained using gradient descent with backpropagation. The training process allows the model to learn better frequencies and phase shifts so that the sinusoid units more accurately represent the seasonality of an underlying time series. Since frequencies and phase shifts can change, the model learns a more reliable periodicity of the underlying series rather than assuming the period is of a predetermined size. Training also tunes the weights of the augmentation function which estimates the non-periodic component of the time series. 

In addressing the problem of overfitting, the ND model proposes the use of L1 regularization on the output weights to promote sparsity and shrink the less essential sinusoidal components. However, ND is designed for estimating a single output - the expected value of the time series. In our QFNN model we are estimating multiple quantile outputs where each cosine component may play a different role in the estimation of each quantile. Thus we do not want to shrink cosine components, but are instead interested in learning different possible models during the training phase to see how each cosine component affects individual estimates. To achieve such regularization we use dropout <cit.>. 

Dropout is a recently proposed regularization technique that randomly drops units, along with their connections, from the neural network during training. This can significantly reduce overfitting and gives major improvements over other regularization methods. Our main proposal for the usage of dropout in QFNN vs L1 regularization is to efficiently approximate in training all possible sub-models of a given architecture and then take their average. No regularization was applied for the augmentation function f(t).

Next, for parameter initialization, there is a considerable distinction in how QFNN is initialized compared to other FNNs. Instead of initializing parameters to mimic the iDFT we randomly set all bias terms near 0, the output weights W^[2] which represent the amplitudes are initialized randomly close to 1, and the input weights W^[1] which represent the frequencies are set to multiples of π k where k is a specific hidden node. The input weight parameters of the augmentation function f(t) are initialized randomly to close to 1 and its bias terms randomly set near 0. 

Before training starts, the input data is preprocessed in the same fashion as in <cit.> to improve learning. First, the time associated with each training sample is normalized between 0 (inclusive) and 1 (exclusive) on the time axis. By doing this normalization testing data points will have a time value greater than or equal to 1. With this normalization the 1/N term in the frequencies is taken into account by transforming t into t/N. Next, if the max value in the training set is greater then 10, then the training set is scaled between 0 and 10. This normalization avoids any look ahead bias by using on the max value in the training set only. Both these preprocessing steps expedite learning and help prevent the model from falling into local optimums.

We present the full proposed QFNN methodology architecture in Fig. <ref>. Summarizing our methodology, we first partition a given time series into training and testing sets. Preprocessing of the training set is then conducted which includes applying a logarithmic filter if multiplicative seasonality is present. A multiplicative time series can be identified by visually inspecting a plot of the series or by analyzing it's autocorrelation plot. This is also considered a key step in the Box–Jenkins decomposition method <cit.>. This preprocessing step is needed since QFNN is an additive model. Next if the training data has points above 10, the data is normalized between 0 and 10. Next, the training and testing times are normalized so that training times are between 0 and 1. Parameters of the QFNN are then initialized as described in the previous paragraphs. Training of the model is conducted using batch gradient descent. After the max number of training epochs is reached the QFNN model is ready to be used for testing. Forecasts can be provided for a multi-step period of indefinite time steps. After a test set prediction is made, preprocessing steps are reversed if any were conducted. Preprocessing steps may include scaling a time series back to its original scale or removing the logarithmic filter from the outputs by exponentiating the predictions. We implement QFNN in Python 3.6 using the Keras/tensorflow framework <cit.>.







§ VALIDATION
 

In this section, we evaluate the effectiveness of the QFNN method for the estimation of quantiles and prediction intervals on 8 diverse datasets, described in Section <ref>, and compare results with 9 benchmark methods described in Section <ref>. For each case (dataset), we conduct three kind of studies to see how the proposed algorithm handles. We analyze the ability to estimate single median percentile values, multiple percentiles (which are used combined to form multiple prediction intervals), and extremely low/high percentiles. For experimentation in each case study, we use 50% of the data for training and the other 50% for testing. All models only saw observations in the training sets. Test data were never presented to the models and were used just to calculate evaluation metrics. The use of 50% of the time series for testing was done to achieve the goal of long-term multi-step PF. 

For the analysis of the estimation of median percentiles, of nominal level τ = 0.5, we use QS for evaluation. For the analysis of estimating multiple percentiles, we look at 100 quantiles whose nominal values are equally spaced between 0 and 1 (0.0099, 0.01988, 0.0297, ..., 0.9702, 0.9801, 0.9900). These quantiles can be combined to form 50 prediction intervals with nominal coverage rates of 0.95% to 98.01%. We also analyze extreme valued percentiles and prediction intervals whose nominal coverage is very low or high. We estimate ten low and high valued quantiles (with levels of 0.005, 0.010, 0.015, 0.02, 0.025, 0.975, 0.98, 0.985, 0.99, 0.995) which combined together form five prediction intervals with nominal coverage of 95%, 96%, 97%, 98%, and 99%. For analyzing multiple and extreme percentile estimates we use QS, and for analyzing prediction interval performance we use ACE and SS.

In each of the experiments, we use one linear activation function for the augmentation function of QFNN to capture the trend. We found that using more than one augmentation function or using nonlinear activation functions such as tanh, sigmoid, and rectified linear units did not provide any significant improvements in trend estimation. The dropout regularization rate can have a large impact on the performance of QFNN in terms of overfitting. Therefore, we apply grid search using the training data for validation. We pick an optimal value of dropout, which we search in a range between 5% to 60%, based on the lowest QS. Grid search is applied for each case study. An example grid search is seen in Fig. <ref> for choosing the dropout rate for estimating the median percentile for each dataset. The QS metric is standardized to allow for comparison between datasets.

For all experiments, we use a maximum training iteration of 10,000 for QFNN and a smoothing rate of 0.01. The learning rate of 0.5 was used in all experiments, except when estimating the single median value for the Air Passengers dataset which failed to converge during training. In that case, a lower learning rate of 0.1 was used. Through initial analysis of training results, we found that altering the training iterations, smoothing rate, or learning rate had little impact, so no form of hyperparameter tuning was used for them. 




 §.§ Case Studies
 

We carry out experiments on eight nonstationary univariate time series datasets, seven being real-world case studies and one synthetic case. These datasets were explicitly picked because they display a diverse set of periodic and aperiodic patterns such as trend, additive and multiplicative seasonality, multiple seasonality, cycles, and irregular patterns. Table <ref> the characteristics of all the datasets. The first case study is the classical Air Passengers time series <cit.> which is composed of 144 samples of the number of passengers flying each month from January 1949 to December 1960. It has a positive linear trend and multiplicative seasonality.

The second case study is the yearly mean and monthly smoothed total sunspot numbers from 1700 to 2017 <cit.>. It consists of 318 samples with a time granularity of one year. This time series includes an unstable (non-constant) seasonal patterns over time. Case study 3 is the load demand from ISO New England <cit.>. Its time series is composed of 744 samples for January 2017. Target values represent real-time demand in MW for wholesale market settlement from revenue quality metering. This case study displays seasonal and cyclical patterns. Internet traffic data in bits from a private ISP with centers in 11 European cities is used for the 4rth case study <cit.>, which exhibits multiple seasonality. We use the data that corresponds to a transatlantic link and was collected from on June 18 to July 16, 2005. 

The highly random movements of the stock market are almost impossible to predict but some stocks may exhibit unseen cycles or trends over more extended periods of time <cit.>. To examine such possible patterns we use the closing stock prices of Apple Inc. over five years from 2012 to the beginning of April 2018 <cit.>. The next two case studies are normalized solar and wind power for September 2012 and January 2012. These two datasets come from the Global Energy Forecasting Competition of 2014 <cit.>. Wind power is highly chaotic and is very difficult to forecast from univariate time series.

The last case study looks at ocean wave elevation, the main motivation for using such data is the irregular sinusoidal nature of waves. Due to the difficulty of finding high resolution deep ocean wave elevation measurements we construct a synthetic dataset. For simulating of ocean waves we focus on vertical sensors for predicting irregular wave formations. Under generally well accepted assumptions <cit.>, the wave elevation for sensor locations (x,y) on the ocean surface for all times t the exact time waveform which would be observed at a particular point in the ocean can be described by

    H (x,y,t)= ∑_i=1^L A_i cos(ω_i^2/g(xcos(β_i)+y sin(β_i))- ω_i t + ϕ_i ),


which has the parameters A for the amplitude, ω for the frequency measured in radians per second (rads/s), β for the wave angular direction in radians measured relative to the x-axis, and ϕ for the phase in radians. We chose to estimate waves at the origin with L = 2, and for each parameter we arbitrarily chose the values A = [1,1.5], ω = [0.5π, 0.1 π], and ϕ = [1.2, 1.4]. To each observation we also include additive white Gaussian noise which we assume come from the sensors.



 §.§ Benchmark Methods
 

To thoroughly examine the forecasting accuracy of our QFNN method we compare it to nine simple and state-of-the-art probabilistic forecasting methods. These include two naive approaches, the uniform and persistence methods. Three-time series models which are the autoregressive integrated moving average model, the seasonal autoregressive integrated moving average model, and exponential smoothing with trend and seasonality. Lastly, we use four advanced PF methods: linear quantile regression, polynomial quantile regression, composite support vector quantile regression, and a composite quantile regression neural network. All methods are implemented in Python 3.6.

The uniform method (UM), commonly used in wind power PF studies<cit.>, is the simplest of all the methods. UM assumes that any observation in the time series has equal probability to occur at any time step. The support of the UM is defined by the parameters a and b which are the minimum and maximum values of the training set for each case study. Quantiles are then defined by F^-1(τ) = (1-τ)a + τ b for τ∈ [0,1]. For deterministic forecasting, the persistence forecast method is a very popular benchmark and is known to be hard to outperform for single point or short look-ahead forecasts. We use the persistence method (PM) <cit.> for PF as a benchmark where the forecast error is assumed to be random and normally distributed, it's mean and variance are computed by the latest observations. For our experiments, we use the last S observations from the training set to calculate the moments of the PM distribution where S corresponds to the size of the seasonality derived from the autocorrelation function (ACF). To ensure that UM and PM can estimate appropriate multi-step forecasts we extend both of them by adding an estimated linear trend component from the training data. 

The next three benchmark methods are well established time series models. We use the autoregressive integrated moving average (ARIMA) model, seasonal ARIMA (SARIMA) model, and exponential smoothing with trend and seasonality (ETS) model, also known as the Holt-Winters seasonal method. We choose ARIMA because it can eliminate non-stationarity through an initial differencing step to better fit time series for prediction, and we select SARIMA and ETS to capture periodic patterns better. Parameters of ARIMA and SARIMA are selected using grid search with the Akaike information criterion and the application of the parsimony principle to prevent over-fitting. The seasonal parameter S for SARIMA and ETS is chosen using the ACF. Quantiles are estimated for ARIMA, SARIMA, and ETS assuming the normality assumption <cit.>. ARIMA and SARIMA are implemented in Python using the sarimax function from the statsmodels package <cit.>, and ETS is implemented in Python from a holtwinters package <cit.>. 

Assuming a normal distribution for quantile prediction with ARIMA, SARIMA, and ETS is a parametric PF approach and can be somewhat restrictive and may not appropriately estimate the forecast distribution. Therefore, we use four advanced nonparametric PF methods linear quantile regression (QR), polynomial QR (PQR), composite support vector quantile regression (SVQR) <cit.>, and composite quantile regression neural network (QRNN) <cit.>. All four of these advanced benchmarks use time as their input for direct comparison to QFNN. QR, PQR, and QRNN are implemented in Keras/tensorflow <cit.>. We implement composite SVQR using the liquidSVM package with quantile regression loss and Gaussian kernel <cit.>. Our last benchmark, QRNN, is similar in structure to our QFNN model where we use one hidden layer, rectified linear units in the hidden layer, and uses a smooth approximation of the pinball loss function. The main differences in QRNN are that L2 regularization is used, no augmentation function is used, and all parameters are randomized during initialization. For QRNN we use a maximum training iteration of 10,000, a learning rate of 0.1, and a smoothing rate of 0.01, and a L2 regularization rate of 0.001 for all experiments.



 §.§ Results and Discussion




  §.§.§ Median Percentile Analysis


We first analyze the performance of QFNN and benchmark methods for predicting the median percentile for each case study. We utilize QS for evaluation between benchmarks (on the test component of each dataset), and standardized QS to see result comparisons across case studies. From Fig. <ref> we see the standardized QS for every method. QFNN tied closely to ETS with a nonstandard QS of 22.25 for QFNN and 20.82 for ETS. QFNN resulted in the lowest QS for every dataset except for solar and wind power where PM and ETS had better performance. For figures <ref> to <ref> we plot the estimated median percentiles where red dots represent the underlying case study time series. The colored curves represent median forecasts of QFNN and 3 benchmark methods that were able to capture the most information of periodic and aperiodic patterns visually. We found that only the benchmark methods SARIMA and ETS for case studies 1 to 3 and 7 to 8 in Table <ref> were able to fit meaningful periodic patterns and thus we only display these three benchmarks. We also include QRNN, a quantile regression based benchmark for comparison to QFNN with the same loss function. Most benchmark methods yielded poor fits for case studies 4, 6, and 7; therefore we do not show these plots. We discuss a visual evaluation of individual case studies in Section <ref>.



  §.§.§ Multiple and Extreme Percentile Estimation


To establish QFNN's ability to estimate composite quantiles for multi-step probabilistic forecasting we conduct an experiment to predict 100 percentiles with nominal coverage rates equally space between 0 and 1. These quantile can be combined together to form 50 PIs with equally spaced PINC. We use QS to evaluate the 100 percentiles, and ACE and SS to evaluate the 50 PIs. Figures <ref> to <ref> showcase QFNN estimation of 50 prediction intervals across test and training data from all 8 case studies. The red line in these figures represents the time series observations. Standardized QS for QFNN and benchmark methods across the 8 datasets for estimating 100 quantiles are presented in Fig. <ref>. We showcase standardized scores so as to make comparisons between datasets. 

QFNN yields the lowest scores for each dataset except for solar and wind power where ETS and PM had the lowest scores. ACE scores are presented in Fig. <ref> in the form of a reliability plot. QFNN averaged a 23% ACE rate across the 8 datasets, and most methods yielded a high ACE for the stock price data indicating poor coverage. The lowest ACE across all methods was the wave elevation time series. ARIMA had the best coverage for the wave elevation time series. Standardized SS is presented in Fig. <ref> as a sharpness plot. QFNN yields the lowest sharpness for every dataset except for the Air Passengers time series where the persistence method had the narrowest mean intervals due to its variance estimation from the training data. The narrow QFNN intervals explain the poor coverage highlighted by the ACE scores.

In many cases we are interested in extreme percentile values such as in anomaly detection or risk management in finance. Having well estimated extreme valued quantiles can have a big impact in these fields. Therefore, as over-viewed in the beginning of this section, we analyze ten low and high valued quantiles for QFNN, which combined form 5 high PINC values. As with multiple quantile evaluation above, we utilize QS to measure the performance of extreme percentiles, and ACE and SS to evaluate the PIs. In Fig. <ref> QS plot is presented for QFNN across the 8 datasets for estimating extreme valued percentiles. QS for most time series is lower for extreme nominal rates such as 0.5% or 99.5%, and higher for percentiles with coverage of 2.5% and 97.5%. For PIs, we see in Fig. <ref> the reliability plot for QFNN. The air passengers, sunspots, and wave elevation time series showed the best coverage of observations by the PINC rates, while wind power, solar power, and stock price time series resulted in poorer fits with high ACE scores up to 64% for wind power. Lastly, a sharpness plot for QFNN is presented in Fig. <ref>. As expected, as the PINC increases for PIs, the their width increases.



  §.§.§ Individual Case Study Analysis
 

In our first individual case study analysis, we inspect the air passengers dataset, where the first six years of data (72 samples) are used for training QFNN and each benchmark method. The next six years were used for prediction. An ACF evaluation finds that the time series has a season of S = 12, and grid search found ARIMA(2,1,3) and SARIMA(1,0,0)(1,0,1)[12] to be the best hyperparameters for ARIMA and SARIMA respectively. We see in Fig. <ref> QRNN estimates the trend but not the seasonality. ETS and SARIMA estimate both trend and seasonality well, but the median forecasts fall below the test data. QFNN learns the shape of the data to capture the median appropriately. Fig. <ref> shows prediction intervals fitting the test data well. In the last 3 seanal periods, QFNN however does estimate percentiles above the observed data. For the remaining experiments, ARIMA and SARIMA hyperparameters are displayed in Table <ref>.

Our second experiment demonstrates the power of QFNN in modeling non-constant seasonal patterns. The sunspots case study is used which has seasonal patterns of varying amplitudes. The years from 1700 to 1858 were used for training, and the years 1859 to 2017 were used for testing. We see in Fig. <ref> that QRNN fails to capture any meaningful pattern in its prediction. SARIMA captures a seasonal pattern that is out of phase with the sunspot test series, and ETS shots off in the test set with a positive trend. QFNN captures a seasonal pattern that is a bit more in phase with the number of sunspots over the years and is also able to learn multiple patterns of the sunspot time series. Fig. <ref> shows prediction intervals fitting the test data surprising well, QFNN captures higher peaks around 1943 and lower peaks around 1903. 

The third experiment uses the real-time load demand case study. We use the first 372 hours in the time series for training. In the median plot of Fig. <ref> SARIMA captures the daily seasonality but fails to capture any cycles in the test set, and ETS goes up in the test set with a positive trend. QFNN learns both the seasonal and cyclical pattern of the load demand. The capture of the cyclical pattern in the load data by QFNN is better presented in Fig. <ref> where we see a tight fit of the daily seasonal and weekly cyclical pattern. The only deviation being around January 21, 2017, which shows a lower observed demand in load possibly due to a warmer weekend and less power needed for heating. The fourth experiment uses the solar power case study and uses the first 380 hours of training. The training set includes samples from both sunny and non-sunny days where solar power is lower than average. QFNN learns a constant daily quantile pattern for all its prediction intervals in Fig. <ref>. 

The fifth experiment using the closing stock prices of the Apple corporation is considered a fascinating case study due to the highly random nature of stock movements. For training, the first 2.5 years of closing prices are used, and testing is composed of closing prices up until the start of April 2018. In Fig. <ref> we see a long-term positive trend and possibly a cycle in the stock price of Apple across the five years. In the plot, SVQR fits the linear trend of the stock series but nothing else. ETS learns a non-existing seasonal pattern, and SARIMA doesn't seem to capture any meaningful pattern. With QFNN we see that it learns the cyclic and positive trend of the stock price which follows the test data better than any other method. This is also demonstrated in Fig. <ref> that up until the end of 2017 QFNN follows the trend and cycle, but then in 2018, the price of Apple jumps higher than the prediction. 

The prediction intervals forecasted by QFNN for the simulated wave elevation case study is shown in Fig. <ref> and the median predictions in Fig. <ref>. It is not a surprise that since ocean waves can be modeled by a sum of sinusoids that QFNN can estimate well the amplitudes, frequencies, and phases of the irregular periodic patterns. For the median percentile predictions, QRNN captures no meaningful pattern, and ETS and SARIMA both capture a periodic pattern whose phase and frequency do not visually match the observed wave elevations.

The remaining experiments are the internet and wind case studies. Median plots of these experiments are not shown since the benchmark methods performed poorly in capturing the multiple or irregular seasonal patterns in these time series. We present the prediction intervals by QFNN in Fig. <ref> where the first 343 hours are used for training. We see that QFNN can learn the multiple seasonal patterns of internet traffic data. The last experiment conducted is on the wind power dataset. In Fig. <ref> we see no identifiable periodic or aperiodic patterns in the wind time series training set. This explains why QFNN has a hard time forecasting the test set. We do see a few peaks predicted by QFNN such as on January 27 and the 30th, but overall the PF is very poor. The wind experiment demonstrates that nonstationary time series without periodic patterns or trend can not be predicted by QFNN. 



§ CONCLUSION
 

Probabilistic predictions can provide a much better analysis of uncertainty then point forecasts. In this paper, a novel approach for probabilistic forecasting is presented called the quantile Fourier neural network. The proposed approach uses a smooth approximation to the pinball ball loss function for estimating composite quantiles. Furthermore, the proposed model provides forecasts using extrapolation based regression instead of autoregression. Extrapolation based regression has not been studied before for probabilistic forecasting. Empirical results on real world univariate time series showcase that our model is able to appropriately capture periodic and aperiodic components to provide high-quality probabilistic predictions. 

From the thorough skill, reliability, and sharpness analysis and visual inspection of median percentile, 100 percentile, and extreme value percentile predictions we can form the following conclusions on which circumstances QFNN performs well and why. QFNN is able to appropriately capture the phase, frequency, and amplitude of periodic components of seasonal time series. This is evident in its quantile estimates for seasonal time series such as wave elevation and load demand. QFNN is also able to capture multiple seasonality such as in the fit for the time series of sunspots and internet traffic. Lastly, QFNN performs well when linear trend is present such as in the air passengers time series, and is even able to loosely capture long term cycles such as in the stock price time series. However, the shortcomings of QFNN are that it is unable to capture patterns in more chaotic data such as wind power and short term stock prices. Another shortcoming is that QFNN is designed for signal extrapolation, and thus it will repeat the same patterns it learned during training on testing data. This can make it ideal for time series imputation. But it is unable to predict unforeseen patterns or patterns dependent on other variables beyond time. For instance, QFNN is not able to predict non-sunny days for solar power, it performs bad on wind power, and it is not able to predict the daily values of stock prices. 

Given the novelty of our approach, more research needs to be conducted to assess its application to more domains and under different scenarios. Further studies could also look at other cost functions such as using the interval score that measures overall skill of PIs and sharpness <cit.>. This could provide appropriate prediction intervals that may have even higher reliability and while maintaining narrow intervals. Furthermore, the influence of exogenous variables as additional inputs to the quantile Fourier neural network could be explored.



§ REFERENCES


elsarticle-num













	Datasets used in the experiments.
	
	
		Case Study     Target     Samples     Time Granularity     Reference 

		1     Air Passengers     144     Month     <cit.> 

		2     Sunspots     318     Year     <cit.> 

		3     Real-Time Load Demand     744     Hour     <cit.> 

		4     Internet Traffic Data (in bits)     686     Hour     <cit.> 

		5     Apple Closing Stock Price     1581     Day     <cit.> 

		6     Solar Power     760     Hour     <cit.> 

		7     Wind Power     744     Hour     <cit.> 

		8     Ocean Wave Elevation     400     Second     (simulated) 

		











































































































