DEPTS: DEEP EXPANSION LEARNING FOR PERIODIC TIME SERIES FORECASTING
Year: 2021

Abstract:
Periodic time series (PTS) forecasting plays a crucial role in a variety of industries to foster critical tasks, such as early warning, pre-planning, resource scheduling, etc. However, the complicated dependencies of PTS signals on its inherent periodicity as well as the sophisticated composition of various periods hinder the performance of PTS forecasting. In this paper, we introduce a deep expansion learning framework, DEPTS, for PTS forecasting. DEPTS starts with a decoupled formulation by introducing the periodic state as a hidden variable, which stimulates us to make two dedicated modules to tackle the aforementioned two challenges. First, we develop an expansion module on top of residual learning to perform a layer-by-layer expansion of those complicated dependencies. Second, we introduce a periodicity module with a parameterized periodic function that holds sufficient capacity to capture diversified periods. Moreover, our two customized modules also have certain interpretable capabilities, such as attributing the forecasts to either local momenta or global periodicity and characterizing certain core periodic properties, e.g., amplitudes and frequencies. Extensive experiments on both synthetic data and real-world data demonstrate the effectiveness of DEPTS on handling PTS. In most cases, DEPTS achieves significant improvements over the best baseline. Specifically, the error reduction can even reach up to 20% for a few cases. Finally, all codes for this paper will be publicly available.

heading: 1 INTRODUCTION
text: Time series (TS) with apparent periodic (seasonal) oscillations, referred to as periodic time series (PTS) in this paper, is pervasive in a wide range of critical industries, such as seasonal electricity spot prices in power industry (Koopman et al., 2007), periodic traffic flows in transportation (Lippi et al., 2013), periodic carbon dioxide exchanges and water flows in sustainability domain (Seymour, 2001; Tesfaye et al., 2006). Apparently, PTS forecasting plays a crucial role in these industries since it can foster their business development by facilitating a variety of capabilities, including early warning, pre-planning, and resource scheduling (Kahn, 2003; Jain, 2017).
Given the pervasiveness and importance of PTS, two obstacles, however, largely hinder the performance of existing forecasting models. First, future TS signals yield complicated dependencies on both adjacent historical observations and inherent periodicity. Nevertheless, many existing studies did not consider this distinctive periodic property (Salinas et al., 2020; Toubeau et al., 2018; Wang et al., 2019; Oreshkin et al., 2020). The performance of these methods has been greatly restrained due to its ignorance of periodicity modeling. Some other efforts, though explicitly introducing periodicity modeling, only followed some arbitrary yet simple assumptions, such as additive or multiplicative seasonality, to capture certain plain periodic effects (Holt, 1957; 2004; Vecchia, 1985b; Taylor & Letham, 2018). These methods failed to model complicated periodic dependencies beyond much simplified assumptions. The second challenge lies in that the inherent periodicity of a typical real-world TS is usually composed of various periods with different amplitudes and frequencies. For example, Figure 1 exemplifies the sophisticated composition of diversified periods via an eightyears hourly TS of electricity load in a region of California. However, existing methods (Taylor & Letham, 2018; Smyl, 2020) required the pre-specification of periodic frequencies before estimating other parameters from data, which attempted to evade this obstacle by transferring the burden of periodicity coefficient initialization to practitioners.
To better tackle the aforementioned two challenges, we develop a deep expansion learning framework, DEPTS, for PTS forecasting. The core idea of DEPTS is to build a deep neural network that conducts the progressive expansions of the complicated dependencies of PTS signals on periodicity to facilitate forecasting. We start from a novel decoupled formulation for PTS forecasting by introducing the periodic state as a hidden variable. This new formulation stimulates us to make more customized and dedicated designs to handle the two specific challenges mentioned above.
For the first challenge, we develop an expansion module on top of residual learning (He et al., 2016; Oreshkin et al., 2020) to conduct layer-by-layer expansions between observed TS signals and hidden periodic states. With such a design, we can build a deep architecture with both high capacities and efficient parameter optimization to model those complicated dependencies of TS signals on periodicity. For the second challenge, we build a periodicity module to estimate the periodic states from observational data. We represent the hidden periodic state with respect to time as a parameterized periodic function with sufficient expressiveness. In this work, for simplicity, we instantiate this function as a series of cosine functions. To release the burden of manually setting periodic coefficients for different data, we develop a data-driven parameter initialization strategy on top of Discrete Cosine Transform (Ahmed et al., 1974). After that, we combine the periodicity module with the expansion module to perform end-to-end learning.
To the best of our knowledge, DEPTS is a very early attempt to build a customized deep learning (DL) architecture for PTS that explicitly takes account of the periodic property. Moreover, with two delicately designed modules, DEPTS also owns certain interpretable capabilities. First, the expansions of forecasts can distinguish the contributions from either adjacent TS signals or inherent periodicity, which intuitively illustrate how the future TS signals may vary based on local momenta and global periodicity. Second, coefficients of the periodicity module have their own practical meanings, such as amplitudes and frequencies, which provide certain interpretable effects inherently.
We conduct experiments on both synthetic data and real-world data, which all demonstrate the superiority of DEPTS on handling PTS. On average, DEPTS reduces the error of the best baseline by about 10%. In a few cases, the error reduction can even reach up to 20%. Besides, we also include extensive ablation tests to verify our critical designs and visualize specific model components to interpret model behaviors.

heading: 2 RELATED WORK
text: TS forecasting is a longstanding research topic that has been extensively studied for decades. After a comprehensive review of the literature, we find three types of paradigms in developing TS models. At an early stage, researchers developed simple yet effective statistical modeling approaches, including exponentially weighted moving averages (Holt, 1957; 2004; Winters, 1960), auto-regressive moving averages (ARMA) (Whittle, 1951; 1963), the unified state-space modeling approach as well as other various extensions (Hyndman & Khandakar, 2008). However, these statistical approaches only considered the linear dependencies of future TS signals on past observations. To handle highorder dependencies, researchers attempted to adopt a hybrid design that combines statistical modeling with more advanced high-capacity models (Montero-Manso et al., 2020; Smyl, 2020). At the same time, with the great successes of DL in computer vision (He et al., 2016) and natural language processing (Vaswani et al., 2017), various DL models have also been developed for TS forecasting (Rangapuram et al., 2018; Toubeau et al., 2018; Salinas et al., 2020; Zia & Razzaq, 2020). Among them, the most representative one is N-BEATS (Oreshkin et al., 2020), which is a pure DL
architecture that has achieved state-of-the-art performance across a wide range of benchmarks. The connections between DEPTS and N-BEATS have been discussed in Section 4.2.
As for PTS forecasting, many traditional statistical approaches explicitly considered the periodic property, such as periodic ARMA (PARMA) (Vecchia, 1985a;b) and its variants (Tesfaye et al., 2006; Anderson et al., 2007; Dudek et al., 2016). However, as discussed in Sections 1 and 3, these methods only followed some arbitrary yet simple assumptions, such as additive or multiplicative seasonality, and thus cannot well handle complicated periodic dependencies in many real-world scenarios. Besides, other recent studies either followed the similar assumptions for periodicity or required the pre-specification of periodic coefficients (Taylor & Letham, 2018; Smyl, 2020). To the best of our knowledge, we are the first work that develops a customized DL architecture to model complicated periodic dependencies and to capture diversified periodic compositions simultaneously.

heading: 3 PROBLEM FORMULATIONS
text: We consider the point forecasting problem of regularly sampled uni-variate TS. Let xt denote the time series value at time-step t, and the classical auto-regressive formulation is to project the historical observations xt−L:t = [xt−L, . . . , xt−1] into its subsequent future values xt:t+H = [xt, . . . , xt+H−1]:
xt:t+H = FΘ(xt−L:t) + t:t+H , (1)
where H is the length of the forecast horizon, L is the length of the lookback window, FΘ : RL → RH is a mapping function parameterized by Θ, and t:t+H = [ t, . . . , t+H−1] denotes a vector of independent and identically distributed Gaussian noises. Essentially, the fundamental assumption behind this formulation is the Markov property xt:t+H ⊥ x0:t−L|xt−L:t, which assumes that the future values xt:t+H are independent of all farther historical values x0:t−L given the adjacent shortterm observations xt−L:t. Note that most existing DL models (Salinas et al., 2020; Toubeau et al., 2018; Wang et al., 2019; Oreshkin et al., 2020) directly follow this formulation to solve TS. Even traditional statistical TS models (Holt, 1957; 2004; Winters, 1960) are indeed consistent with that if omitting those long-tail exponentially decayed dependencies introduced by moving averages.
To precisely formulate PTS, on the other hand, this assumption needs to be slightly modified such that the dependency of xt:t+H on xt−L:t is further conditioned on the inherent periodicity, which can be anchored by associated time-steps. Accordingly, we alter the equation (1) into
xt:t+H = F ′ Θ(xt−L:t, t) + t:t+H , (2)
where other than xt−L:t, F ′
Θ : RL × R → RH takes an extra argument t, which denotes the forecasting time-step. Existing methods for PTS adopt a few different instantiations of F ′Θ. For example, Holt (1957; 2004) developed several exponentially weighted moving average processes with additive or multiplicative seasonality. Vecchia (1985a;b) adopted the multiplicative seasonality by treating the coefficients of the auto-regressive moving average process as time dependent. Smyl (2020) also adopted the multiplicative seasonality and built a hybrid method by coupling that with recurrent neural networks (Hochreiter & Schmidhuber, 1997), while Taylor & Letham (2018) chose the additive seasonality by adding the periodic forecast with other parts as the final forecast.

heading: 4 DEPTS
text: In this section, we elaborate on our new framework, DEPTS. First, we start with a decoupled formulation of (2) in Section 4.1. Then, we illustrate the proposed neural architecture for this formulation in Sections 4.2 and 4.3. Last, we discuss the interpretable capabilities in Section 4.4.

heading: 4.1 THE DECOUPLED FORMULATION
text: To explicitly tackle the two-sided challenges of PTS forecasting, i.e., complicated periodic dependencies and diversified periodic compositions, we introduce a decoupled formulation (3) that refines (2) by introducing a hidden variable zt to represent the periodic state at time-step t:
xt:t+H = fθ(xt−L:t, zt−L:t+H) + t:t+H , zt = gφ(t), (3)
where we treat zt ∈ R1 as a scalar value to be consistent with the uni-variate TS xt ∈ R1, we use fθ : RL × RL+H → RH to model complicated dependencies of the future signals xt:t+H on the local observations xt−L:t and the corresponding periodic states zt−L:t+H within the lookback and forecast horizons, and gφ : R1 → R1 is to produce a periodic state zt for a specific time-step t. The right part of Figure 2 depicts the overall data flows of this formulation, in which the expansion module fθ and the periodicity module gφ are responsible for handling the two aforementioned PTSspecific challenges, respectively.

heading: 4.2 THE EXPANSION MODULE
text: To effectively model complicated periodic dependencies, the main challenge lies in the trade-off between model capacity and generalization. To avoid the over-fitting issue, many existing PTS approaches relied on the assumptions of additive or multiplicative seasonality (Holt, 1957; Vecchia, 1985b; Anderson et al., 2007; Taylor & Letham, 2018; Smyl, 2020), which however can hardly express periodicity beyond such simplified assumptions. Lately, residual learning has shown its great potentials in building expressive and generalizable DL architectures for a variety of crucial applications, such as computer vision (He et al., 2016) and language understanding (Vaswani et al., 2017). Specifically, N-BEATS (Oreshkin et al., 2020) conducted a pioneer demonstration of introducing residual learning to TS forecasting. Inspired by these successful examples and with full consideration of PTS-specific challenges, we develop a novel expansion module fθ on top of residual learning to characterize the complicated dependencies of xt:t+H on xt−L:t and zt−L:t+H .
The proposed architecture for fθ, as shown in the middle part of Figure 2, consists of N layers in total. As further elaboration in the left part of Figure 2, each layer `, share an identical residual structure consisting of three residual branches, which correspond to the recurrence relations of z(`)t−L:t+H , x (`) t−L:t, and x̂ (`) t:t+H , respectively. Here x (`) t−L:t and z (`) t−L:t+H denote the residual terms of xt−L:t and zt−L:t+H after `-layers expansions, and x̂ (`) t:t+H denotes the cumulative forecasts after ` layers. In layer `, three residual branches are specified by two parameterized blocks, a local block f lθl(`) and a periodic block fpθp(`), where θl(`) and θp(`) are their respective parameters.
First, we present the updating equation for z(`−1)t−L:t+H , which aims to produce the forecasts from periodic states and exclude the periodic effects that have been used. To be more concrete, fpθp(`) takes in z(`−1)t−L:t+H and emits the `-th expansion term of periodic states, denoted as v (`) t−L:t+H ∈ RL+H . v (`) t−L:t+H has two components, a backcast component v (`) t−L:t and a forecast one v (`) t:t+H . We leverage v (`) t−L:t to exclude the periodic effects fromx (`−1) t−L:t and adopt v (`) t:t+H as the portion of forecasts purely from the `-th periodic block. Besides, when moving to the next layer, we exclude v(`)t−L:t+H from z (`−1) t−L:t+H as z (`) t−L:t+H = z (`−1) t−L:t+H − v (`) t−L:t+H to encourage the subsequent periodic blocks to focus on the unresolved residue z(`)t−L:t+H .
Then, since v(`)t−L:t is related to the periodic components that have been used to produce a part of forecasts, we construct the input to f lθl(`) as (x (`−1) t−L:t − v (`) t−L:t). Here the purpose is to encourage f lθl(`) to focus on the unresolved patterns within x (`−1) t−L:t. f l θl(`) emits u(`)t−L:t and u (`) t:t+H , which correspond to the local backcast and forecast expansion terms of the `-th layer, respectively. After that, we update x(`)t−L:t by further subtracting u (`) t−L:t from (x (`−1) t−L:t − v (`) t−L:t) as x (`) t−L:t = x (`−1) t−L:t − v (`) t−L:t − u (`) t−L:t. Here the insight is also to exclude all analyzed patterns of this layer to let the following layers focus on unresolved information. Besides, we update x̂(`)t:t+H by adding both u(`)t:t+H and v (`) t:t+H as x̂ (`) t:t+H = x̂ (`−1) t:t+H + u (`) t:t+H + v (`) t:t+H . The motivation of such expansion is to decompose the forecasts from the `-th layer into two parts, u(`)t:t+H and v (`) t:t+H , which correspond to the part from local observations excluding redundant periodic information and the other part purely from periodic states, respectively.
Note that before the first layer, we have x(0)t−L:t = xt−L:t, z (0) t−L:t+H = zt−L:t+H , and x̂ (0) t:t+H = 0. Besides, we collect the cumulative forecasts x̂(N)t:t+H of theN -th layer as the overall forecasts x̂t:t+H . Therefore, after stacking N layers of z(`)t−L:t+H , x (`) t−L:t, and x̂ (`) t:t+H , we have the following triply residual expansions that encapsulate the left and middle parts of Figure 2:
zt−L:t+H = z (0) t−L:t+H = N∑ `=1 v (`) t−L:t+H + z (N) t−L:t+H ,
xt−L:t = x (0) t−L:t = N∑ `=1 (u (`) t−L:t + v (`) t−L:t) + x (N) t−L:t,
x̂t:t+H = x̂ (N) t:t+H = N∑ `=1 (u (`) t:t+H + v (`) t:t+H),
(4)
where z(N)t−L:t+H and x (N) t−L:t+H are deemed to be the residues irrelevant to forecasting.
Connections and differences to N-BEATS. Our design of fθ shares the similar insight with NBEATS (Oreshkin et al., 2020), which is stimulating a deep neural network to learn expansions of raw TS signals progressively, whereas N-BEATS only considered the generic TS by modeling the dependencies of xt:t+H on xt−L:t. In contrast, our design is to capture the complicated dependencies of xt:t+H on xt−L:t and zt−L:t+H for PTS. Moreover, to achieve periodicity modeling, N-BEATS produces coefficients solely based on the input signals within a lookback window for a group of predefined seasonal basis vectors with fixed frequencies and phases. However, our work can capture diversified periods in practice and model the inherent global periodicity.
Inner architectures of local and periodic blocks. The local block f lθl(`) aims to produce a part of forecasts based on the local observations excluding redundant periodic information as (xt−L:t −∑`−1 i=1 v (i) t−L:t). Thus, we reuse the generic block developed by Oreshkin et al. (2020), which consists of a series of fully connected layers. As for the periodic block fpθp(`), which handles the relatively stable periodic states, we can adopt a simple two-layer perception. Due to the space limit, we include more details of inner block architectures in Appendix A.

heading: 4.3 THE PERIODICITY MODULE
text: To represent the sophisticated periodicity composed of various periodic patterns, we estimate zt via a parameterized periodic function gφ(t) that holds sufficient capacities to incorporate diversified periods. In this work, for simplicity, we instantiate this function as a series of cosine functions as gφ(t) = A0 + ∑K k=1Ak cos(2πFkt+Pk), where K is a hyper-parameter denoting the total number of periods, A0 is a scalar parameter for the base scale, Ak, Fk, and Pk are the scalar parameters for the amplitude, the frequency, and the phase of the k-th cosine function, respectively, and φ represents the set of all parameters. Coupling gφ with fθ illustrated in Section 4, we can effectively model the
periodicity-aware auto-regressive forecasting process in the equation (2). However, it is extremely challenging to directly conduct the joint optimization of φ and θ from random initialization. The reason is that in such a highly non-convex condition, the coefficients in φ are easily trapped into numerous local optima, which do not necessarily characterize our desired periodicity.
Parameter Initialization. To overcome the optimization obstacle mentioned above, we formalize a two-stage optimization problem based on raw PTS signals to find good initialization for φ. First, we construct a surrogate function, gMφ (t) = A0 + ∑K k=1Mk · Akcos(2πFkt + Pk), to enable the selection of a subset of periods via M = {Mk, k ∈ {1, · · · ,K}}, where each Mk ∈ {0, 1} is a mask variable to enable or disable certain periods. Note that gφ(t) is equivalent to gMφ (t) when every Mk is equal to one. Then, we construct the following two-stage optimization problem:
M∗ = arg min ‖M‖1<=J LDval(gMφ∗(t)), φ∗ = arg min φ LDtrain(gφ(t)), (5)
where LDtrain and LDval denote the discrepancy losses on training and validation, respectively; the inner stage is to obtain φ∗ that minimizes the discrepancy between zt and xt on the training data Dtrain; the outer stage is a binary integer programming on the validation data Dval to find M∗ that can select certain periods with good generalization, and the hyper-parameter J controls the maximal number of periods being selected. With the help of such two-stage optimization, we are able to estimate generalizable periodic coefficients from observational data as a good starting point for φ to be jointly optimized with θ. Nevertheless, it is still costly to perform exact optimization of equations (5) in practice. Thus, we develop a fast approximation algorithm to obtain an acceptable solution with affordable costs. Our approximation algorithm contains the following two steps:
• conducting Discrete Cosine Transform (Ahmed et al., 1974) of PTS signals on Dtrain to select top-K cosine bases with the largest amplitude as an approximated solution of φ∗;
• iterating over the selected K cosine bases from the largest amplitude to the smallest one and greedily select J periods that generalize well on the validation set.
Due to the space limit, we include more details of this approximation algorithm in Appendix B. After obtaining approximated solutions φ̃∗ and M̃∗, we fix M = M̃∗ to exclude those unstable periodic coefficients and initialize φ with φ̃∗ to avoid being trapped into bad local optima. Then, we follow the formulation (3) to perform the joint learning of φ and θ in an end-to-end manner.

heading: 4.4 INTERPRETABILITY
text: Owing to the specific designs of fθ and gφ, our architecture is born with a degree of interpretability. First, for fθ, as shown in equations (4), we decompose x̂t:t+H into two types of components, u
(`) t:t+H and v (`) t:t+H . Note that v (`) t:t+H is conditioned on zt−L:t+H and independent of xt−L:t. Thus,∑N
`=1 v (`) t:t+H represents the portion of forecasts purely from periodic states. Meanwhile, u (`) t:t+H depends on both xt−L:t and zt−L:t+H , and it is transformed by feeding the subtraction of v (`) t−L:t
from x(`−1)t−L:t into the `-th local block. Thus, we can regard ∑N `=1 u (`) t−L:t as the forecasts from the local historical observations excluding the periodic effects, referred to as the local momenta in this paper. In this way, we can differentiate the contribution to the final forecasts into both the global periodicity and the local momenta. Second, gφ, the periodicity estimation module in our architecture, also has interpretable effects. Specifically, the coefficients in gφ(t) have practical meanings, such as amplitudes, frequencies, and phases. We can interpret these coefficients as the inherent properties of the series and connect them to practical scenarios. Furthermore, by grouping various periods together, gφ provides us with the essential periodicity of TS filtering out various local momenta.

heading: 5 EXPERIMENTS
text: Our empirical studies aim to answer three questions. 1) Why is it important to model the complicated dependencies of PTS signals on its inherent periodicity? 2) How much benefit can DEPTS gain for PTS forecasting compared with existing state-of-the-art models? 3) What kind of interpretability can DEPTS offer based on our two customized modules, fθ and gφ? To answer the first two questions,
we conduct extensive experiments on both synthetic data and real-world data, which are illustrated in Section 5.1 and 5.2, respectively. Then, Section 4.4 answers the third question by comparing and interpreting model behaviors for specific cases.
Baselines. We adopt the state-of-the-art DL architecture, N-BEATS (Oreshkin et al., 2020), as our primary baseline since it has been shown to outperform a wide range of DL models, including MatFact (Yu et al., 2016), Deep State (Rangapuram et al., 2018), Deep Factors (Wang et al., 2019), and DeepAR (Salinas et al., 2020), and many competitive hybrid methods (Montero-Manso et al., 2020; Smyl, 2020). Besides, we also include PARMA as a reference to be aware of the positions of conventional statistical models. For this baseline, we leverage the AutoARIMA implementation provided by Löning et al. (2019) to search for the best configurations automatically.
Evaluation Metrics. To compare different models, we utilize the following two metrics, normalized deviation, abbreviated as nd, and normalized root-mean-square error, denoted as nrmse, which are conventionally adopted by Yu et al. (2016); Rangapuram et al. (2018); Salinas et al. (2020); Oreshkin et al. (2020) on PTS-related benchmarks.
nd =
1 |Ω| ∑ (i,t)∈Ω |xit − x̂it|
1 |Ω| ∑ (i,t)∈Ω |xit| , nrmse =
√ 1 |Ω| ∑ (i,t)∈Ω(x i t − x̂it)2
1 |Ω| ∑ (i,t)∈Ω |xit|
, (6)
where i is the index of TS in a dataset, t is the time index, and Ω denotes the whole evaluation space.

heading: 5.1 EVALUATION ON SYNTHETIC DATA
text: To intuitively illustrate the importance of periodicity modeling, we generate synthetic data with various periodic dependencies and multiple types of periods. Specifically, we generate a simulated TS signal xt by composing an auto-regressive signal lt, corresponding to the local momentum, and a compounded periodic signal pt, denoting the global periodicity, via a function f c as xt = f c(lt, pt), which characterizes the dependency of xt on lt and pt. First, we produce lt via an auto-regressive process, lt = ∑L i=1 αilt−i + l t, in which αi is a coefficient for the i-lag dependency, and the error term lt ∼ N (0, σl) follows a zero-mean Gaussian distribution with standard deviation σl. Then, we produce pt by sampling from another Gaussian distribution N (zt, σp), in which zt is characterized by a periodic function (instantiated as gφ(t) in Section 4.3), and σp is a standard deviation to adjust the degree of dispersion for periodic samples. Next, we take three types of f c(lt, pt), (lt + pt), (lt + pt)
2, and (lt + pt)3, to characterize the linear, quadratic, and cubic dependencies of xt on lt and pt, respectively. Last, after data generation, all models only have access to the final mixed signal xt for training and evaluation.
Due to the space limit, we include the main results in Figure 3 and leave finer grained parameter specifications and more experimental details to Appendix C. For each setup (linear, quadratic, cubic) in Figure 3, we have searched for the best lookback length (L) for N-BEATS and the best number of periods (J) for DEPTS on the validation set and re-run the model training with five different random seeds to produce robust results on the test set. We can observe that for all cases, even with an exhaustive search of proper lookback lengths for N-BEATS, there exists a considerable performance gap between it and DEPTS, which verifies the utility of explicit periodicity modeling. Moreover, as the periodic dependency becomes more complex (from linear to cubic), the average error reduction of DEPTS over N-BEATS keeps increasing (from 7% to 11%), which further demonstrates the importance of modeling high-order periodic effects.

heading: 5.2 EVALUATION ON REAL-WORLD DATA
text: Other than simulation experiments, we further demonstrate the effectiveness of DEPTS on realworld data. We adopt three existing PTS-related datasets, ELECTRICITY1, TRAFFIC2, and M4 (HOURLY)3, which contain various long-term (quarterly, yearly), mid-term (monthly, weekly), and short-term (daily, hourly) periodic effects corresponding to regular economic and social activities. These datasets serve as common benchmarks for many recent studies (Yu et al., 2016; Rangapuram et al., 2018; Salinas et al., 2020; Oreshkin et al., 2020). For ELECTRICITY and TRAFFIC, we follow two different test splits defined by Salinas et al. (2020) and Yu et al. (2016), and the evaluation horizon covers the first week starting from the split date. As for M4 (HOURLY), we adopt the official test set. Besides, we note that the time horizons covered by these three benchmarks are still too short, which results in very limited data being left for periodicity learning if we alter the time split too early. This drawback of lacking enough long PTS limits the power of periodicity modeling and thus may hinder the research development in this field. To further verify the importance of periodicity modeling in real-world scenarios, we construct two new benchmarks with sufficiently long PTS from public data sources. The first one, denoted as CAISO, contains eight-years hourly actual electricity load series in different zones of California.4 The second one, referred to as NP, includes eight-years hourly energy production volume series in multiple European countries.5 Accordingly, we define four test splits that correspond to all four seasons of the last year for robust evaluation.
For all benchmarks, we search for the best hyper-parameters of DEPTS on the validation set. Similar to N-BEATS (Oreshkin et al., 2020), we also produce ensemble forecasts of multiple models trained with different lookback lengths and random initialization seeds. Tables 1 and 2 show the overall performance comparisons. On average, the error reductions (nd) of DEPTS over N-BEATS on ELECTRICITY, TRAFFIC, M4 (HOURLY), CAISO, and NP are 12.5%, 3.5%, 8.7%, 13.3%, and 9.9%, respectively. Interestingly, we observe some prominent improvements in a few specific cases,
1https://archive.ics.uci.edu/ml/datasets/ElectricityLoadDiagrams20112014 2https://archive.ics.uci.edu/ml/datasets/PEMS-SF 3https://github.com/Mcompetitions/M4-methods/tree/master/Dataset/Train 4http://www.energyonline.com/Data 5https://www.nordpoolgroup.com/Market-data1/Power-system-data
such as 18.7% in ELECTRICITY (2014-09-01), 23.1% in CAISO (2020-10-01), and 15.2% in NP (2020-10-01). At the same time, we also observe some tiny improvements, such as 2.6% in TRAFFIC (2008-06-15) and 3.3% in CAISO (2020-07-01). These observations imply that the predictive abilities and the complexities of periodic effects may vary over time, which corresponds to the changes in performance gaps between DEPTS and N-BEATS. Nevertheless, most of the time, DEPTS still brings stable and significant performance gains for PTS forecasting, which clearly demonstrate the importance of periodicity modeling in practice.
Due to the space limit, we leave more details about datasets and hyper-parameters used in real-world experiments to Appendix D. Moreover, to achieve effective periodicity modeling, we have made several critical designs, such as the triply residual expansions in Section 4.2 and the composition of diversified periods in Section 4.3. We also conduct extensive ablation tests to verify these critical designs, which are included in Appendix E.

heading: 5.3 INTERPRETABILITY
text: In Figure 4, we illustrate the interpretable effects of DEPTS via two cases, the upper one from ELECTRICITY and the bottom one from TRAFFIC. First, from subplots in the left part, we observe that DEPTS obtains much more accurate forecasts than N-BEATS and PARMA. Then, in the middle and right parts, we can visualize the inner states of DEPTS to interpret how it makes such forecasts. As Section 4.4 states, DEPTS can differentiate the contributions to the final forecasts x̂t:t+H into the local momenta ∑N `=1 u (`) t:t+H and the global periodicity ∑N `=1 v (`) t:t+H . Interestingly, we can see that DEPTS has learned two different decomposition strategies: 1) for the upper case, most of the contributions to the final forecasts come from the global periodicity part, which implies that this case follows strong periodic patterns; 2) for the bottom case, the periodicity part just characterizes a major oscillation frequency, while the model relies more on the local momenta to refine the final forecasts. Besides, the right part of Figure 4 depicts the hidden periodic state zt estimated by our periodicity module gφ(t). We can see that gφ(t) indeed captures some inherent periodicity. Besides, the actual PTS signals also present diverse variations at different time, which further demonstrate the importance of leveraging fθ to model the dependencies of xt:t+H on both xt−L:t and zt−L:t+H . We include more case studies and interpretability analysis in Appendix F.

heading: 6 CONCLUSION
text: In this paper, we develop a novel DL framework, DEPTS, for PTS forecasting. Our core contributions are to model complicated periodic dependencies and to capture sophisticated compositions of diversified periods simultaneously in TS forecasting. Extensive experiments on both synthetic data and real-world data demonstrate the effectiveness of DEPTS on handling PTS. Moreover, periodicity modeling is actually an old and crucial topic for traditional TS modeling but is rarely studied in the context of DL. Thus we hope that the new DL framework together the two new benchmarks with evident periodicity and sufficiently long observations can facilitate more future research on PTS.

heading: A BLOCK ARCHITECTURES
text: As illustrated in Section 4.2, the local block f lθl(`) produces forecasts based on local observed PTS signals excluding redundant periodic effects. This goal aligns with that of N-BEATS to extract informative representations from generic TS signals. Therefore, we reuse the generic block design of N-BEATS to instantiate f lθl(`). Here we include a brief description of the local block for completeness. Please refer to Section 3.1 in (Oreshkin et al., 2020) for more details.
Local Block. The left part of Figure 5 shows the detailed architecture within a local block, where we use x̃(`)t−L:t = x (`−1) t−L:t − v (`) t−L:t to denote the portion of the local observations x (`−1) t−L:t excluding the periodic effects v(`)t−L:t for the `-th layer. After taking in x̃ (`) t−L:t, we pass it through four fullyconnected layers and then obtain the backcast coefficients c(`)b and the forecast coefficients c (`) f via two linear projections:
u (`),1 t−L:t = FC`,1(x̃ (`) t−L:t), u (`),2 t−L:t = FC`,2(u (`),1 t−L:t), u (`),3 t−L:t = FC`,3(u (`),2 t−L:t),
u (`),4 t−L:t = FC`,4(u (`),3 t−L:t), c (`) b = LINEAR b `(u (`),4 t−L:t), c (`) f = LINEAR f ` (u (`),4 t−L:t),
where FC is a standard fully-connected layer with ReLU activation (Nair & Hinton, 2010), and LINEAR denotes a linear projection function. Then, we pass these coefficients to the basis layers, hb(·) and hf (·), to obtain the backcast term u(`)t−L:t and the forecast term u (`) t:t+H , respectively. The generic choice for hb(·) can simply be another linear projection function, which is also adopted by us since it produces more competitive and stable performance on PTS-related benchmarks than other interpretable basis layers, as shown by (Oreshkin et al., 2020) in Appendix C.4.
Periodic Block. The periodic block fpθp(`) aims to extract predictive information from associated periodic states, which are relatively simple and stable compared with rapidly shifting PTS signals. Therefore, we can adopt a simple architecture while still maintain desired effects. In this work, we use one-layer standard fully-connected layer to encode z(`−1)t−L:t and leverage another two linear projection functions to obtain the backcast term v(`)t−L:t and the forecast term v (`) t:t+H as the periodic effects of the `-th layer.
v (`),1 t−L:t+H = FC`(z (`−1) t−L:t+H), v (`) t−L:t = LINEAR `(v (`),1 t−L:t), v (`) t:t+H = LINEAR `(v (`),1 t:t+H),
where FC and LINEAR share the same meanings mentioned above. Moreover, when training for multiple series simultaneously, we use a series-specific scalar parameter αi (i is the series index) to take account of differences in the strengths of periodicity by updating v(`)t−L:t+H as αi · v (`) t−L:t+H .

heading: B PARAMETER INITIALIZATION FOR THE PERIODICITY MODULE
text: As illustrated in Section 4.3, we leverage a fast approximation approach to obtain an acceptable solution of the two-stage optimization problem (5) with affordable costs in practice. Algorithm 1 summarizes the overall procedure for this fast approximation.
Algorithm 1: Parameter initialization for the periodicity module. Input: Dtrain = x0:Tv , Dval = xTv:T , K, and J Conduct DCT over x0:Tv . Sort the top-K cosine bases by amplitudes in descending order to obtain φ̃∗ = {Ã∗0} ∪ {Ã∗k, F̃ ∗k , P̃ ∗k }Kk=1. Initialize M̃∗ = 0. for j in [1, · · · ,K] do
if ‖M̃∗‖1 < J then Update M̃∗j by arg minMj∈{0,1} LDval(g Mj φ̃∗ (t)) else return φ̃∗ and M̃∗
end end Output: φ̃∗ and M̃∗
First, we divide the whole PTS signals x0:T into the training part Dtrain = x0:Tv and the validation part Dval = xTv:T , where Tv is the split time-step. Then, the inner optimization stage is to identify the optimal parameter set φ∗ that can best fit the training data:
φ∗ = arg min φ LDtrain(gφ(t)), gφ(t) = A0 + K∑ k=1 Ak cos(2πFkt+ Pk), (7)
where the hyper-parameter K controls the capacity of gφ(t) and the discrepancy training loss LDtrain can be instantiated as the mean square error ∑Tv−1 t=0 ‖gφ(t)− xt‖22. Directly optimizing (7) via gradient descent from random initialization is inefficient and time-consuming since it involves numerous gradient updates and is easily trapped into bad local optima. Fortunately, our instantiation of gφ(t) as a group of cosine functions shares the similar format with Discrete Cosine Transform (DCT) (Ahmed et al., 1974). Accordingly, we conduct DCT over x0:Tv and select top-K cosine bases with the largest amplitudes, which characterize the major periodic oscillations of this series, as the approximated solution φ̃∗ of (7).
Next, we enter the outer optimization stage to select certain periods with good generalization:
M∗ = arg min ‖M‖1<=J LDval(gMφ∗(t)), gMφ∗(t) = A∗0 + K∑ k=1 Mk ·A∗k cos(2πF ∗k t+ P ∗k ), (8)
where the hyper-parameter J further constrains the expressiveness of gMφ (t) for good generalization. Conducting exact optimization of this binary integer programming is also costly since it involves an exponentially growing parameter space. Similarly, to capture the major periodic oscillations as much as possible, we develop a greedy strategy that iterates the selected K cosine bases from the largest amplitude to the smallest and greedily assigns 1 or 0 to Mk depending on whether the k-th period further reduces the discrepancy loss on the validation data. Specifically, assuming K periods are already sorted by their amplitudes descendingly and are indexed by k, we construct another
surrogate function gMjφ∗ (t) for the j-th greedy step:
g Mj φ∗ (t) = Mj ·A ∗ j cos(2πF ∗ j t+ P ∗ j ) + [ A∗0 + j−1∑ k=1 M̃∗k ·A∗k cos(2πF ∗k t+ P ∗k ) ] , (9)
where {M̃∗k} j−1 k=1 is determined in previous steps, Mj is an integer parameter to be set in the current step. Thus, for the j-th step, we are actually updating M̃∗j by
M̃∗j = arg min Mj∈{0,1} LDval(g Mj φ∗ (t)). (10)
Besides, to tolerate the approximation errors introduced by φ̃∗, which may result in shifted periodic oscillations, we use Dynamic Time Warping to measure the discrepancy of gMjφ∗ (t) and xt on Dval. We continue this greedily updating steps until selecting J periods in total or completing the traverses of all K selected periods. Finally, we obtain an approximated solution M̃∗ of (8).
Complexity Analyses. We also provide the complexity analyses of Algorithm 1, which runs very fast in practice and takes up negligible time compared with training neural networks. Let us denote the length of training series as Lt and the length of validation series as Lv . First, the complexity of conducting DCT over training series is O(Ltlog(Lt)). Then, the complexity of selecting top-K frequencies with the largest amplitudes is O(Ltlog(K)), which can be ignored since K << Lt. Next, we need to select at most J frequencies greedily based on the generalization errors on the validation set. Since we measure the generalization errors via dynamic time warping, the total worst complexity for this selection procedure is O(KL2v). In total, the worse complexity of our approximation algorithm for a series is O(Ltlog(Lt) + KL2v). In practice, Lv , the length of the validation series, is relatively small, and K, the maximum number of frequencies, can be regarded as a constant. So, the squared complexity term O(KL2v) is not a big trouble.

heading: C MORE DETAILS ON SYNTHETIC EXPERIMENTS
text: As Section 5.1 states, we produce a TS lt via an auto-regressive process, lt = ∑L i=1 αilt−i + l t, in which αi is a coefficient for the i-lag dependence, and the error term lt ∼ N (0, σl) follows a zero-mean Gaussian distribution with standard deviation σl. Specifically, we set L as 3 and σl as 1. We leverage uniform samples from [−1, 1] to initialize {αi}3i=1 and also uniformly sample three values from [0, 5 for the initial points, l−3, l−2, and l−1. Then, we produce pt by sampling from another Gaussian distribution N (zt, σp), in which zt is characterized by a periodic function (instantiated as gφ(t) in Section 4.3), and σp is a standard deviation to adjust the degree of dispersion for periodic samples. Specifically, we also set σp as 1 and produce zt via a composition of three cosine bases, 8cos(2π(t + 2)/50), 4cos(2π(t + 3)/10), 2cos(2πt/4), and a base level, 30. These three cosine bases represent long-term, mid-term, short-term periodic effects, respectively, which are very similar to the circumstance in practice. Next, we take three types of f c(lt, pt), (lt + pt), (lt + pt)
2, and (lt + pt)3, to characterize the linear, quadratic, and cubic dependencies of xt on lt and pt, respectively. We repeat the above procedure for 5000 time steps and divide them into 4000, 100, and 900 for training, validation, and evaluation, respectively. Figure 6 shows the first 1000 time steps of these synthetic series. Note that after data generation, all models only have access to the final mixed signal xt for training and evaluation.
Moreover, as illustrated in Section 5.1, we search for the best loobkack length (L) for N-BEATS and the best number of periods (J) for DEPTS. The lookback length for DEPTS is fixed as 3, which is also determined by hyper-parameter tuning on the validation set. Figure 7 shows detailed comparisons of N-BEATS and DEPTS for different configurations of L and J . We can see that N-BEATS always needs a relatively long lookback window, such as 48 or 96 time steps, to capture those periodic patterns effectively. Besides, further increasing the lookback length will introduce more irrelevant noises, which overwhelm effective predictive signals and thus result in more worse performance. In contrast, with effective periodicity modeling, DEPTS can achieve better performance by using a short lookback window, which is also consistent with the auto-regressive process that governs the local momenta.

heading: D MORE DETAILS ON REAL-WORLD EXPERIMENTS
text: D.1 DATASETS
Table 3 includes main statistics of the five datasets used by our experiments. We can see that the existing datasets (ELECTRICITY, TRAFFIC, and M4 (HOURLY)) utilized by recent studies usually have a large number of series but with relatively short lengths. Therefore, it is hard to identify or evaluate yearly or quarterly periods on these benchmarks. In contrast, CAISO and NP contain tens of series with the lengths of several years, which can better illustrate the inherent periodicity of these series and serve as complementary benchmarks for PTS modeling.
For N-BEATS, we use its default hyper-parameters6 for ELECTRICITY, TRAFFIC, and M4 (HOURLY), and we report its hyper-parameters searched on CAISO and NP in Table 4. Besides, N-BEATS used multiple loss functions, such as sMAPE or MASE, for model training, and we also follow these setups. Tables 5 and 6 include the hyper-parameters of DEPTS for all five datasets. Note
6https://github.com/ElementAI/N-BEATS
that, all these hyper-parameters are searched on a validation set, which is defined as the last week before the test split. Moreover, for a typical dataset with multiple series, we build an independent periodicity module gφ for each series and perform respective parameter initialization procedures as described in Appendix B. Then, for all datasets (splits), we train 30 models (6 lookback lengths × 5 random seeds) for both N-BEATS and DEPTS and then produce ensemble forecasts for fair and robust evaluation.

heading: E ABLATION TESTS
text: As Figure 8 shows, we adopt three ablated variants of DEPTS to demonstrate our critical designs in the expansion module (Section 4.2):
• DEPTS-1: removing the residual connection of (x(`−1)t−L:t−v (`) t−L:t) so that the outputs of the
local block u(`)t−L:t are only conditioned on the raw PTS signals xt−L:t, which correspond to the mixed observations of local momenta and global periodicity.
• DEPTS-2: removing the residual connection of (x̂(`−1)t:t+H + v (`) t:t+H) so that the contribu-
tions to the forecasts only come from the local block, which takes in the signals excluding periodic effects progressively.
• DEPTS-3: removing the residual connection of (z(`−1)t−L:t+H)− v (`) t−L:t+H so that the inputs
to the periodic block of each layer are the same hidden variables zt−L:t+H .
We also construct another four baselines to demonstrate the importance of our customized periodicity learning:
• NoPeriod: removing the periodic blocks by directly feeding (xt − zt) to N-BEATS. • RandInit: randomly initializing periodic coefficients (φ) and directly applying the end-to-
end learning.
• FixPeriod: fixing the periodic coefficients (φ) after the initialization stage and only tuning θ during end-to-end optimization.
• MultiVar: treating zt as a covariate of xt and feeding (xt, zt) into an N-BEATS-style model via two channels.
Moreover, as illustrated in Section 4.3 and Appendix B, the maximal number of selected periods J is a critical hyper-parameter to balance expressiveness and generalization of gφ(t). Thus, we conduct experiments with different J to verify its sensitivity on different datasets. Tables 7 and 8 include experimental results of these model variants on ELECTRICITY, TRAFFIC, CAISO, and NP with different J . Since we only identify one reliable period via Algorithm 1 on M4 (HOURLY), we report its results separately in Table 9.
First, as Tables 7 and 8 show, J is a crucial hyper-parameter that has huge impacts on forecasting performance. The reason is that if J is too small, the periodicity module gφ cannot produce effective representations of the inherent periodicity to boost the predictive ability. While if it is too large, gφ has a high risk of over-fitting to the irrelevant noises contained by the training data, which also results in poor predictive performance. Moreover, the interactions between local momenta and global periodicity may vary over time. Therefore, it is critical to search for a proper J for each PTS and
each split point to pursue better performance. Fortunately, we demonstrate that the hyper-parameter tuning of J on the validation set can ensure its good generalization abilities on the subsequent test horizons.
Then, let us focus on Table 7 to compare DEPTS with DEPTS-1, DEPTS-2, and DEPTS-3. First, we can see that DEPTS-1 usually produces the worst performance in most cases, which demon-
strates that excluding periodic effects from raw PTS signals can stably and significantly boost the performance of PTS forecasting. Second, for most cases, DEPTS outperforms DEPTS-2, and the performance gaps can be remarkable, such as 0.139 vs. 0.148 on ELECTRICITY and 0.020 vs. 0.021 on CAISO. These results verify the importance of including the portion of forecasts solely from the periodicity module. Third, DEPTS-3 can produce competitive results compared with DEPTS in many cases. Nevertheless, after selecting the best J for each dataset, DEPTS still slightly outperforms DEPTS-3 in most cases. Besides, from Table 9, we also observe that DEPTS performs much better than DEPTS-3 on M4 (HOURLY). Thus we retain the residual connection to reduce the periodic effects leveraged by previous layers.
Next, let us focus on Table 8 to compare DEPTS with other four baselines, NoPeriod, RandInit, MultiVar, and FixPeriod. First, we observe that NoPeriod usually produces the worst preformance. The reason is that (xt–zt) denotes the raw time-series signal subtracting the periodic effect, so it is challenging for the model to forecast the future signals, xt:t+H , solely based on the periodicityagnostic inputs, (xt−L:t–zt−L:t). Second, RandInit also produces much more worse results than DEPTS, which demonstrate the importance of initializing periodic coefficients (Section 4.3). Third, DEPTS performs much better than FixPeriod in most cases, which demonstrate the effectiveness of fine-tuning periodic coefficients after the initialization stage. Last, we observe that sometimes MultiVar can produce comparable and even slightly better results than DEPTS. However, after selecting the best J on each dataset for these two models, we find that DEPTS still outperforms MultiVar consistently and significantly, which also demonstrates the superiority of our expansion learning. Moreover, as Table 9 shows, DEPTS outperforms all these baselines by a large margin on M4 (HOURLY), which containing very short PTS with 854 observations on average. Given limited data, all our critical designs, such as properly initializing periodic coefficients, fine-tuning periodic coefficients, and conducting expansion learning to decouple the dependencies of xt on zt, play much crucial roles in producing accurate forecasts.

heading: F MORE CASE STUDIES AND INTERPRETABILITY ANALYSES
text: In the following, we further study the interpretable effects of DEPTS with more cases. Figures 9, 10, 11 and 12 show the additional two cases on ELECTRICITY, TRAFFIC, CAISO, and NP, respectively. Following Figure 4, we compare the forecasts of N-BEATS and DEPTS in the left side, differentiate the forecasts of DEPTS into the local part (DEPTS-L) and the periodic part (DEPTS-P) in the middle side, and plot the hidden state zt together with the PTS signals in the right side. The general observations are that with the help of explicit periodicity modeling, DEPTS achieves better performance than N-BEATS in PTS forecasting, and DEPTS has learned diverse behaviors for different cases. Besides, we also include their critical periodic coefficients (amplitude Ak, frequency Fk, and phase Pk) in Tables 10, 11, 12, and 13. We find that DEPTS can learn many meaningful periods that are consistent with practical domains.

Review summaries:
This paper focuses on an important property of time series, periodicity, and mainly studies the problem of periodic time series (PTS) forecasting. This work solves two main research challenges of PTS forecasting: (1) to learn the dependencies of observation data on periodicity; (2) to learn the sophisticated compositions of various periods. The authors propose a deep expansion framework on top of residue learning for dependency learning and a parameterized surrogate function for periodicity learning. Extensive experiments on both synthetic data and 5 real-world datasets demonstrate a significant improvement on time series forecasting tasks when considering periodicity especially.

The paper addresses the problem of time series forecasting, esp. with periodic dependencies. The authors propose a model that combines a learnt one-dimensional sum of cosines as periodic signal with residual feedforward neural network (N-BEATS). They propose to learn the model in two stages: first to estimate the cosines with a discrete cosine transform and greedily selecting the K ones with largest amplitude, second to refine both the periodic time encoding and the overall forecasting model end-to-end. In experiments on synthetic, three real-life datasets from the literature and two new real-life datasets with longer training segments, that should allow to identify periodicities more easily, they show that they outperform the underlying N-BEATS model mostly consistently.  

This paper proposed a novel DL framework for periodic time series forecasting. The contributions of the paper include modeling complicated periodic dependency and capturing compositions of diversified periods. The authors conduct extensive experiments on synthetic data and real-world datasets, showing the effectiveness and interpretability of the proposed methods