On projection methods for functional time series forecasting

Two nonparametric methods are presented for forecasting functional time series (FTS). The FTS we observe is a curve at a discrete-time point. We address both one-step-ahead forecasting and dynamic updating. Dynamic updating is a forward prediction of the unobserved segment of the most recent curve. Among the two proposed methods, the first one is a straightforward adaptation to FTS of the $k$-nearest neighbors methods for univariate time series forecasting. The second one is based on a selection of curves, termed \emph{the curve envelope}, that aims to be representative in shape and magnitude of the most recent functional observation, either a whole curve or the observed part of a partially observed curve. In a similar fashion to $k$-nearest neighbors and other projection methods successfully used for time series forecasting, we ``project'' the $k$-nearest neighbors and the curves in the envelope for forecasting. In doing so, we keep track of the next period evolution of the curves. The methods are applied to simulated data, daily electricity demand, and NOx emissions and provide competitive results with and often superior to several benchmark predictions. The approach offers a model-free alternative to statistical methods based on FTS modeling to study the cyclic or seasonal behavior of many FTS.


Introduction
In the last few decades, technological advancements have simplified and decreased the cost of data collecting and storing processes. This new paradigm has made data scientific to confront not only big data sets but also complex data structures. Functional data analysis (FDA) arises naturally in this context to exploit the information recorded over a continuum such as time or space. The literature has experienced a prominent development of techniques and tools to take advantage of this kind of data. While the monographs of [23,38,60] provide complete picture of the methodological and theoretical contributions, recent advances in the field can be found in survey papers [1,18,27,79].
This article addresses the problem of forecasting functional time series (FTS). These are time series {y k , k ∈ N} where each y k is a random function t → y k (t), t ∈ [a, b]. We refer the reader to Hörmann and Kokoszka [28] for theory and examples of FTS. One-step-ahead forecasting is one of the most important problems in FTS analysis, and it has been addressed by diverse prominent literature [see, e.g., 2,3,7,9,11,30,31,34,59]. Also, forward prediction of a curve which is already being observed, termed dynamic updating by Shang and Hyndman [66], a topic of current interest [64,67]. As many classical methods for time series forecasting, for example, those based on the ARIMA family models, the benchmark methods for FTS forecasting are based on fitting some statistical model of some approximated data generation. The most common is to fit ARIMA or VAR models to the Functional Principal Components Scores (FPC) of the functions, see FAR model [9,34]. Klepsch and Klüppelberg [37] considered Functional Moving Average (FMA) and FARMA. [42] and [41] consider functional autoregressive fractionally integrated moving average (ARFIMA). Here we discuss an alternative model-free approach that addresses both one-step-ahead forecasting and dynamic updating in a unified way.
The literature of artificial intelligence and machine learning for time series forecasting has gained considerable prominence over the last decade [see, e.g., 51]. Besides the lack of an underlying model, many of these methods exhibit interesting features compared with classical statistical methods, including their ability to capture nonlinearities [53]. Among the most powerful machine learning technics used for time series forecasting, we must mention artificial neural networks [82] and k-nearest neighbors (KNN). The simplicity of implementing the algorithm and its high accuracy have popularized KNN. It has demonstrated to be a strong contender in forecasting competitions [52]. The method can be schematized as follows: 1. Consider T observations, x 1 , . . . , x T , from a time series and a prediction horizon h. Our goal is to forecast p = (x T +1 , . . . , x T +h ), i.e. the time series h periods ahead. 2. Consider all the possible observed (E + h)-tuples (x i−(E−1) , . . . , x i , x i+1 , . . . , x i+h ). This means E ≤ i ≤ T − h. 3. Denote x i = (x i−(E−1) , . . . , x i ) and its projection P(x i ) = (x i+1 , . . . , x i+h ). Note that p is the projection of x T . 4. Find the k-nearest neighbors to x T from {x i , E ≤ i ≤ T − h}. 5. Predict p by a convex combination of the projections of the k-nearest neighbors.
In addition to the distance function used to find the nearest neighbors and the convex combination used for predicting, KNN involves setting the sometimes called embedding dimension E and the number of neighbors denoted by k. The simplex projection [71] and the S-map projection [70] are methods related to KNN that provide skillful forecasts, free of k.
Moreover, these model-free methods have outperformed any forecast based on models [58]. Surprisingly, functional versions have not been developed for FTS. In contrast, in the context of FDA, the KNN method has been explored and theoretically supported towards regression problems [10,12,40,43] or functional classification problems [13,29]. Biau et al. [10] introduces asymptotic results when the response variable has a finite dimension. Kara et al. [36] provides asymptotic results that, exceptionally, allows for an automatic choice data-driven of the number of neighbors for several operators, among them functional regression with a scalar response. Lian [43] proves the consistency of the KNN regression estimate when both dependent and independent variables are functions. Recently, [22] use the nearest neighbours method for bandwidth selection on a problem of scalar-on-function local linear regression. Remarkably, [83] addresses the problem of forecasting final prices of auctions via functional KNN. Roughly, their forecasting is based on the k-nearest price curves to the item being predicted and their corresponding final prices. We consider dynamic updating in line with these authors. For example, if we predict afternoon electricity demand, we use afternoon data from days with morning data near to the morning of the day being predicted. For one-step-ahead forecasting, we mimic the projection procedure inherent to KNN. For example, for tomorrow's forecasting electricity demand, we consider next-day demands of daily curves near to today. Moreover, we take ideas of [71], and we do not only consider the nearest daily curves but those that surround the most recent data from above and below, when this is possible, to select a neighbor that is representative of its shape and its magnitude. This is the key to produce coherent band predictions, in which the observed part of the day being predicted is inside the band delimited by the curves that we are using for predicting. Whereas the curve to predict may be entirely outside of the band delimited by its k-nearest curves. In this vein, we combine two concepts, nearness and centrality, commonly termed depth in the functional context. Functional depth [25,56] has received a great deal of attention since it was introduced by Fraiman and Muniz [24]. It has been used for several applications; including classification [16,17,19,29,46], outlier detection [8,21,35,54,55], populations comparison [47,48,74], and clustering [68,75]. Functional versions of boxplots and other graphical tools based on different depths have also been proposed for visualizing curves to discover features from a sample that might not be apparent using other methods [20,32,62,72]. However, up to our knowledge, the concept of depth has not been used for forecasting.
Here we introduce a depth-based projection method for forecasting that constitutes a contribution to the literature of nonparametric FDA, an area popularized by [23] and that it is actively growing, see the review by Ling and Vieu [44]. Also, we consider an adaptation to FTS of the KNN briefly schematized above and we follow ideas from [36] to make it automatic. We also consider several benchmark methods in the field for comparisons. For that, we use daily functions of electricity demand and nitrogen monoxide (NOx) emissions. Both are FTS that have already been considered by the literature. Also, both exhibit a strong calendar-effect [78], making them appropriate to illustrate the excellent performance of the method discussed here. As we discuss below, the projection methods have been intended to forecast FTS with seasonal or cyclical patterns. Through a simulation study, we investigate the impact on forecast accuracy from the departure of stationarity. We remark that, although we do not assume any stationarity of first or second order, we do not consider applications to FTS with a significant trend. Specifically, we suppose the curves of the FTS fluctuate in a finite band.
This paper is organized as follows. Section 2 describes the projection methods considered here. Section 3 provides results of comparing the proposed methods with benchmark methods. Section 4 concludes with a general discussion.

Projection methods for forecasting
Let {y i , i ∈ N} be a FTS [28]. Here, we assume that each y i is a univariate random function and, without loss of generality, we suppose that its domain is [0, 1].
• Dynamic updating; when we have observed y 1 , ..., y n−1 on [0, 1] and y n (t) on [0, q], for a given 0 q < 1, and we want to predict {y n (t) : t ∈ (q, 1]}. In both scenarios, we refer to the most recently observed functional datum, either {y n (t) : t ∈ [0, 1]} or {y n (t) : t ∈ [0, q]}, as the focal curve. The latter will be denoted by f and its observation domain by D f , regardless of whether it is [0, 1] or [0, q]. For each i, 1 ≤ i ≤ n − 1, we consider the restriction of y i to D f , this is the curve with graph equal to {(t, y i (t)) : t ∈ D f }. We denote this restriction by y i| f and we refer to Y n = {y 1| f , . . . y n−1| f } as the past-focal-curves sample. Note that in one-step-ahead forecasting y i| f ≡ y i .
Denote by p the curve or curve segment to predict, this is {y n+1 (t) : t ∈ [0, 1]} or {y n (t) : t ∈ (q, 1]}. The domain of p, either [0, 1] or (q, 1], is denoted by D p . Then, for each y i| f ∈ Y n , we consider the projection in dynamic updating.
Note that P f = p and that, although not explored it in this article, h-step-ahead projection could be achieved by Many FTS regularly repeat dependency patterns between past focal curves and their projections. The change in shape and magnitude of the curves over time of some FTS can be visualized from the rainbow plots of Hyndman and Shang [32] and other related plots [65].
In dynamic updating, we observe that y i| f and Py i| f are a partition of y i in two curve segments. The dependence patterns between past focal curves and their projections can be observed from groups of curves with similar colors in the rainbow plots already mentioned. Similar curves may be distant in time due to seasons and cycles, making specific patterns be repeated over time. It is not necessarily an easy task to visualize dependency patterns between focal curves and their projections; it depends on the complexity of the FTS and the curves' shapes, and the number of curves involved. Nonetheless, the cyclic nature of many FTS suggests that there may be statistical regularity between the morphology of focal curves and the morphology of their projections. The core idea of the methods discussed here is based on the existence of this kind of regularity. Thus, if f is surrounded by past focal curves that capture both its shape and magnitude, we could predict p from these past focal curves' projections. Of course, when the FTS has a strong global trend, the future values are always out of the historical data magnitude. Therefore, the projection methods discussed here are not suitable without detrending in such cases.

The focal-curve envelope
The straightforward approach for surrounding f with past focal curves is by using KNN. However, as we already commented, f does not have to be necessarily centred in the band delimited by its k-nearest curves. Even f may be completely outside of this band. To produce bands where f is central, which is natural for band estimation, we use the concept of functional depth. In particular, we use the modified band depth of López-Pintado and Romo [47] with bands formed by two curves. This is a widely used depth measure that allows a million curves to be ranked in only tenths of seconds [73] making our algorithm efficient even for large data sets. Strictly speaking, for any J ⊂ Y n with |J| ≥ 2 and y ∈ J ∪ { f }, we consider In line with Sun and Genton [72], who consider central regions based on the modified band depth for visualizing the spread of the deepest curves, we consider the Focal Central Region (FCR) delimited by J k . This is the region So that R(J k ) can capture the shape and magnitude of f , R(J k ) must be roughly centered at f and its mean width should be small. Having these objectives in mind, we look for a set J of past focal curves, as large as possible in order of not restricting the range of k, such that: 2) f is enveloped by J as much as possible. Here, we say f is more enveloped by J than by J if and only if 3) f is surrounded by near curves of J, as many as possible.
Algorithm 1 provides a set of past focal curves with the three features above that we call the focal-curve envelope and denote by J from now on. In a nutshell, the algorithm iteratively selects as many past curves as possible. From the nearest to the farthest to f , we obtain an envelope of f with an increasing depth.
Let y be the nearest curve to f from Y and N = {y } for y ∈ Y \ {y }, from the nearest curve to the farthest from f , do if f is more enveloped by N ∪ {y} than by N then As an illustration of how Algorithm 1 works, see top panels of Figure 1, where 100 curves from a simulated FTS are considered. The left panels show the twelve curves selected on the first iteration of Algorithm 1 (top panel) and the 12-nearest neighbors (bottom panel). For measuring nearness, hereafter, we use L 2 distance. Besides this is the most widely used distance in the KNN context, we base our choice by looking for a small mean squared error on the prediction method discussed in Subsection 2.3. We emphasize the intervals where the nearest curves do not envelope the focal curve. Right panels show the band delimited by the focal-curve envelope obtained after five iterations and formed by forty-six curves and the corresponding band from the 46-nearest curves. As a general rule, Algorithm 1 selects near curves that are enveloping the focal until it is the deepest one. In contrast, nearest neighbors neither necessarily envelope the focal curve nor make it the deepest one.

Prediction
We make predictions by projecting the k-nearest curves to the focal and projecting the curves that belong to the focal-curve envelope. We compute weighted means of these projections for punctual prediction, assigning more weights to the curves closer to the focal. Following the literature on KNN [36,52], and simplex projections or Smap [58], we consider two types of weights: On one hand, weights inversely proportional to the square distance to the focal, which we identify in the sequel by 1/d weights. On the other hand, exponential weights. Specifically, let f 1 , . . . , f k be the k-nearest curves to f , d i = f i − f 2 , and more general d y = y − f 2 , let J be the focal-curve envelope of f and θ > 0, what we do is to consider the following estimators of p = P f : • functional KNN projections (fKNN) with 1/d and exponential weightŝ • Envelope projections (EP) with 1/d and exponential weightŝ For fixed θ, k in (1) will be chosen by minimizing the prediction Mean Square Error (MSE) on a rolling forecasting origin procedure, as it is customary in the literature of KNN for time series. We denote by k * the value that minimizes the MSE. Similarly, we can select θ in (2). Let θ * be the resulting value of minimizing MSE by rolling origin. We discard selection of optimal pairs (θ, k) forp k,θ in (1) by a high computational cost involved in this two-dimensional optimization problem. To compare the methods, on one hand we computep * k andp. On the other hand, we fix θ = 1 and computep k * ,θ=1 andp θ=1 . In addition, we computep θ * .
For band prediction, we consider the region delimited by the projections of the k deepest curves of J. This is the band This band prediction corresponds to a central functional region such as those defined by [72]. As with the fKNN projections, k is chosen by minimizing a rolling forecasting origin procedure's cost function. For that, we consider the Winkler score. This is a measure that penalizes both bandwidth and non-coverage for a given α ∈ (0, 1), penalizing even more heavily when observations are remarkably outside the band the confidence levels are higher [26,80]. Thus, for a given (1 − α) mean coverage, i.e., the proportion of time that p is in the band, k is chosen by minimizing the score. Mimicking this procedure, we also consider the region delimited by the projection of the k nearest curves for band prediction based on fKNN. Figure 2 shows an example of point and band prediction based on 1000 curves from a simulated FTS corresponding to the model discussed in Section 3.1 with µ = 0. Both one-step-ahead forecasting and dynamic updating provided by EP withp θ are illustrated.

Results
For one-period-ahead forecasting, we compare fKNN and EP estimators with benchmark predictions for stationary FTS. In particular, we implement a variant of the method proposed by Hyndman and Ullah [34] and Hyndman and Shang [31], who suggest functional prediction based on modeling the functional principal component scores as independent time series. Following Aue et al. [9], we will consider the first d functional scores as a d-variate time series since the score vectors could have non-diagonal autocorrelations. The number of scores d is chosen as the smallest integer such that the first d principal components explain at least 80% of the variance of the data. Then, we fit a d-dimensional autoregressive model [50] to the scores' time series. This method is termed FPCF. For validating FPCF as a benchmark method, we repeated the simulation study conducted by Aue et al. [9] and obtained competitive results with the best of their methods (see Appendix 5.1).
For dynamic updating, we compare fKNN, and EP with the method of [39] for predicting partially observed functional data. We also tested the well-known method of Yao et al. [81] but Kraus [39] achieved better results in our experiments. When the FTS is constructed by slicing a periodic univariate time series, we also consider the Block-Moving (BM) approach by Shang and Hyndman [66]. The BM approach is designed to re-arrange the time series so that all the FTS are ultimately observed and, then, to apply FPCF is possible. The literature devoted to dynamic updating has successfully proposed methodologies based on ordinary least squares (OLS) regression, ridge regression (RR), penalized least squares (PLS) regression, or functional linear regression [64,66,67]. However, we do not compare with them because we deal with large data sets where these methods might be computationally inefficient.
In other circumstances, they certainly are the benchmark of the field. Finally, to understand the scale of how superior a method is with respect to the other, we also report about the two most basic estimators. They are the historical mean and the naive prediction. This is just the most recent observed datum corresponding to the prediction range. In our notation, the naive predictor of Py n is Py n−1 .
The implementation of the Envelope Projection method is available at https://github.com/aefdz/nnFTS. The benchmark methods for Functional Time Series are available in the ftsa R package by Hyndman and Shang [33]. Finally, the method for estimating partially observed data was implemented by [39] with the code available from https://www.davidkraus.net/?page_id=349.

Stationary FTS with random shocks
To simulate FTS with the features discussed at the end of Section 2.1, we consider slices of a scalar process with complex periodic rhythm that artificially corrupt with shocks. This process is denoted Y µ , being µ ∈ [0, 1) the expected proportion of affected periods. The process is the sum of three independent components corresponding to: If µ = 0, there are no shocks and the corresponding FTS is stationary. With increasing µ, we explore different non-stationary regimes. Figure 3 shows two random trajectories of Y 0.2 and their respective FTS.
We performed 300 realizations of Y µ , 100 for each level of contamination considered. They were µ = 0.0, 0.2 and 0.4. For each realization, we considered 1000 periods and, from the corresponding 1000 slices, the exercise of one-step-ahead prediction and dynamic updating on the last 100 slices. For dynamic updating, we chose q = 0.5, corresponding to half slice. The resulting MSEs are shown in Table 1. In summary, we found the proposed methods were superior to the benchmark method. They were even better at higher contamination levels, and fKNN is often superior to EP.

Spanish electricity demand
Functional methods for electricity demand forecasting have been widely used [4-6, 14, 57, 59, 63, 76, 77]. These methods often consider particular features of the data, such as weekly seasonality, the effect of holidays, and outliers' presence. Thus, researchers process weekdays separately from Saturdays and Sundays and apply methods for replacing outliers. In addition to this, Raña et al. [59] use additional data from functional exogenous covariates in their Functional Additive Models (FAM) for next-day forecasting. Certainly, they provide the benchmark predictions in functional Spanish electricity demand forecasting.
We consider the challenge of predicting the daily curves of a whole year (2018) without considering any specific peculiarity of the electricity demand by using raw data. For dynamic updating, we considered the practical case of half-day ahead prediction. Our study is based on electricity demand at every 10-minute interval, from January 1, 2014, to December 31, 2018. The data set is available at http://www.ree.es/es/. Following Raña et al. [59], we compute the mean absolute percentage error (MAPE) in addition to the MSE. Viewing the weekly seasonality of electricity demand and the naive prediction corresponding to yesterday for predicting today, we also consider the naive prediction corresponding to the same day one week ago, identified below as seasonal naive (Naive in the tables). The results are shown in Table 2. From this table, we remark that FPCF is slightly inferior to the naive method and significantly inferior to seasonal naive. However, the projection methods were superior to the benchmark competitors for both one-step-ahead forecasting and dynamic updating, and EP slightly superior to fKNN.
Although they are different studies, MAPE of Raña et al. [59] (in the last four rows of the last column of their Table  1) might reference the expected magnitude of the forecasting performance. The naive method in both studies produces the same MAPE (6.383 vs. 6.39). However, we remark that MAPE computed by Raña et al. [59] corresponds to 2012 prediction instead of 2018, using demand each hour, instead of every 10 minutes, and therefore conclusions might be taken with caution. From this extrapolation, FPCF in Table 2 seems to be by far inferior to FAM. This points out the convenience of removing outliers and taking into account the data's weekly seasonality to improve predictions as suggested by Raña et al. [59]. Finally, MAPE of fKNN and EP turned out to be by far smaller than those of FAM, despite being based on data without preprocessed, without splitting the data into Weekdays, Saturdays, and Sundays, and without covariates.
Notably, Raña et al. [59] also reports coverage and mean Winkler scores of prediction bands, providing an approximate order of magnitude of the results to compare them with the prediction bands based on projection methods introduced at the end of Subsection 2.3. Table 3 displays the scores resulting from projection methods for α = 0.05, 0.10 and 0.20. The first noteworthy result is that the mean coverages obtained are often above (1 − α). Second, the mean Winkler scores are often smaller than that those previously reported. Also, we want to remark that the projection methods take into account the weekly and yearly periodicity of the data without any prior information. For illustrating this, we plot the frequency distribution of the number of days between predicted days and days used in the EP prediction in Figure 4. The exponential weights properly averaged and scaled are also plotted in the same figure. This figure reveals the most relevant past information to forecast one-period-ahead and the electricity demand's seasonal structure. In particular, we observe that the most recent day contributes the most to the point forecast, followed by the same weekday but one week before (highest bars at the bottom panel). This weekly contribution vanishes until we approach the values of one year ago (highest bars at the top panel). For example, if we aim to forecast the demand for a given Wednesday, then the envelope would consider data from the last Tuesday, past Wednesdays, and same calendar days of other years in the same season, in order of importance.

NOx emissions in Madrid city center
Air quality is nowadays in the spotlight due to its impact on human health [45]. Many cities gather almost continuous data of different air quality measures, such as nitrogen monoxide (NOx). The Cityhall of Madrid has recorded data since 2001 at different city points and has made them publicly available at https://datos.madrid. es/portal/site/egob/. In this section, we analyze the data from a station located in Plaza de España, Madrid downtown, because it has been functioning uninterrupted since the very beginning. By discarding days with fails on record, we end up with 6187 daily functions from 2001 to 2017 observed each hour.
Following [60], we evaluate these functions every 15 minutes by applying a positive smoothing with K = 11 Fourier basics functions (periodic data) and a λ obtained from minimizing the generalized cross-validation coefficient. As we did with the Spanish electricity demand, we forecast the last recorded year (2017), predicting half a day for dynamic updating exercise. The results are also shown in Table 4.
Results for BM are not reported for the high computational cost involved. Neither MAPE is reported for avoiding divisions by zero. Once again, projection methods were superior to benchmark competitors in real case studies, as shown in Table 4.

Conclusion
The functional projection methods offer an empirical insight into forecasting functional time series of cyclic nature that repeat curve patterns and dependency patterns between contiguous curves. Both one-step-ahead prediction and dynamic updating are addressed in a unified way, without attending to any statistical model that could explain data generation's mechanism.  The projection methods were superior to benchmark methods from the functional time series literature for simulated and real data. The Spanish electricity demand curves illustrated how projection methods take into account structural features of the underlying functional time series, such as weekly and yearly periodicity, crucial in the literature of electricity demand forecasting. Particular characteristics of the electricity demand related to the calendar effect and outliers' presence were anticipated by the method for making accurate predictions. Although these peculiarities are known and often used to model and forecast electricity demand, they could be unknown for other data sets. In these cases, the projection methods may be particularly attractive.
According to the data set, the two projection methods were competitive with each other, being on average one better than the other. By construction, they are similar if the k-nearest curves surround the most recent functional datum (termed here focal curve), making it deep. However, they differ if this curve is outlying with respect to their nearest curves. The method based on curve envelope projections performs better in the simulation and empirical data analyses in such cases. The superiority is because both punctual and band predictions are computed from curves surrounding the focal from below and above. However, for focal curves inlying with respect to their nearest curves, the method based on k-nearest neighbors projections can be even superior, making both methods similar in average, as we observed with the data sets considered in the previous section. This result reveals the common statistical tradeoff between efficiency and robustness that merits a separate study between both methods. Nevertheless, our main goal is to propose a parameter-free forecasting method for functional time series. We remark that the envelope method with θ = 1 presents a competitive and entirely data-driven method. Based on our implementations and executions of the algorithms, the projection methods discussed here lead to substantive improvements in the computational time with respect to the benchmark. Additionally, the method could easily accommodate explanatory variables by-products of a suitable distance between multivariate functions and a multivariate functional depth, already available in the literature [15,29,35,49]. Finally, we note that the methodology is intuitive and data-driven, making the valuable approach to a broad audience, even without any prior knowledge of functional time series.  Table 5: MSE based on functional average, FPCF, EPF and EPU 0.25 , for the two vector standard deviations, (σ1) and (σ2), chosen by Aue et al. [9]. The first row of each setting (κ 1 , κ 2 ) corresponds to small n and the second row to large n.
• MSE based on DB-U 0.25 was by far the smallest.
We emphasize that FAR processes are not necessarily the kind of FTS that we have discussed at the end of Subsection 2.1, and consequently, neither EPF nor EPU should work well, although their performance is reasonably good. For the pair (0.8, 0.0), EPF was competitive. For these values, the future realization is generated only with the focal curve and innovations. Unlike (0.0, 0.8), the focal curve is not involved, making EPF a little competitive. For (0.4, 0.4), the difference of MSE based on FPCF and EPF was roughly 0.1, below to 7% of relative difference. For (0.2, 0.0), the FTS is mainly noise, and the average is the best prediction.

Stationary FTS with random shocks
The noise X is generated with a stationary Gaussian process with zero mean and squared exponential autocovariance function Cov(X(t), X(s)) = exp(− |t − s| 2 2l 2 X ).
This is the default autocovariance function in Gaussian processes simulation [61]. The parameter l X determines the length of the 'wiggles' in the trajectories.
The random function f is a trajectory of a stationary Gaussian process with zero mean and periodic autocovariance function Cov( f (t), f (s)) = exp(− 2 sin 2 (π|t − s|) The parameter l f determines lengthscale in the same way as in the squared exponential autocovariance function. We emphasize that the trajectories of f are periodic functions of period 1. The random shocks process S µ is based on two independent processes. First, an occurrence process {x i , i ≥ 0} being a two states Markov chain, with x 0 = 0 and transition probabilities Second, a random periodic function for modeling irregular shock patterns. Specifically, we use Z(t) = σ 2 g (g 2 (t)−g 2 (0)), g being an independent copy of f . The parameter σ g determines the average distance of g away from zero. The squares are taking for generating a non-zero mean pattern, and g 2 (0) is subtracted for beginning to shock from zero. The random shock at time t is then for t < u Z(t − u)χ ρ (t − u), for t ≥ u u being the minimum time that can elapse before the first shock and ρ = µ/(2 − µ). In our simulations, u is uniformly distributed on [0, 1]. Since the stationary distribution of {x i } is P(x i = 1) = ρ/(1 + ρ) and each shock affects two contiguous periods with probability equals to 1, the proportion of periods affected by shocks, when t → ∞, is µ.
The parameters considered for examples and simulations are l f = 1, l X = 0.2, and σ 2 g = 10.