Seasonal Adjustment of Infra-Monthly Time Series with JDemetra+

2.1. Data Peculiarities

Infra-monthly economic time series often display dynamics that cannot be seen in data recorded at traditional sub-annual intervals, such as months and quarters. Therefore, we briefly summarize their key distinctive features; detailed discussions can be found in de Livera et al. (2011), Ollech (2023), Proietti and Pedregal (2023), and Webel (2022), amongst others.

The seasonal dynamics of infra-monthly time series are typically more complex than those of monthly and quarterly data. The main reason is the potential coexistence of multiple seasonal patterns that can even have fractional periodicities (the origin of this phenomenon is explained in Remark 1 below). To exemplify this increased complexity, Figure 1 shows realized hourly electricity consumption in Germany over the course of April 2023. Both infra-daily and infra-weekly repetitive patterns are clearly visible in this series: electricity consumption is usually lower during nighttime, sharply increases in the morning, peaks around lunchtime, displays a minor dip in the mid-afternoon and eventually decreases until bedtime. This infra-daily pattern can be seen on any given day of the week; however, its average level and wingspan are noticeably smaller on Saturdays and especially on Sundays, creating a characteristic infra-weekly pattern. In addition, a U -shaped infra-yearly pattern is also present in the data, see Figure 11 below and the detailed analysis of the series in Subsection 8.2. Seasonal movements of daily data are often equally complex; however, those of weekly data are further complicated by a year-to-year phase shift in the position of the weeks within the year. For example, the first full Monday through Sunday sequence of the year can be anywhere between January 1 to 7 (as in 2024) and January 7 to 13 (as in 2019); see, for example, Pierce et al. (1984) for more details on this matter.

OPEN IN VIEWER

Figure 1.

Realized hourly electricity consumption in Germany in units of gigawatt hours (GWh) from Saturday, April 1, 2023 (00:00 AM) through Sunday, April 30, 2023 (11:00 PM). Gray area marks the Easter period (April 7–10).

Compared with monthly and quarterly data, the calendar-related dynamics of infra-monthly time series are often more complex as well since data granularity enables more accurate measurements of effects related to both fixed and moving holidays, including potential interactions with the seasonal patterns. During Easter 2023, for example, it is seen that hourly electricity consumption exhibits quite peculiar swings on Good Friday and Easter Monday whereas the series displays relatively normal infra-daily and infra-weekly dynamics on Easter Saturday and Easter Sunday (Figure 1). Data granularity also facilitates measurements of anticipatory and catch-up effects as well as of effects of holiday-related events, such as bridging days, which are usually unquantifiable for monthly and quarterly economic time series and hence not recommended for consideration in official seasonal adjustments of those data. In some countries, secular and religious activities are governed by different calendars, so entangled dual-calendar effects may also show up in some national economic data, such as private consumption.

Besides the aforementioned complex forms of recurrent variation, infra-monthly time series are also susceptible to new, or additional, types of irregularities. Some of them are reminiscent of issues that often need to be dealt with in small-sample survey data, with prime examples being missing, zero, or otherwise meager observations. Irregular spacing is another likely facet of infra-monthly time series. It is typically observed in (infra-)daily monetary data recorded only on banking days. As regards outliers, additive point outliers and structural level shifts that are well-known from the modeling of monthly and quarterly economic data are likely to be present in infra-monthly data as well. However, data granularity also entails a higher probability for the occurrence of new types of transitional outliers with ramp-like, triangular, wavy, or otherwise smooth signatures.

The above discussion shows that complex time series models will be needed to adequately capture all the specifics of infra-monthly data. Model calibration will thus require detailed knowledge of the data generating process; in particular, selecting appropriate calendar regression variables can be very time-consuming, especially when modeling multiple series at once. Therefore, sacrificing some model accuracy in favor of model parsimony (and computing time) will sometimes be inevitable.

2.2. Basic Unobserved Components Framework

Let { y t } denote an observed infra-monthly time series and assume that it can be decomposed additively, maybe after taking logs, into a trend-cyclical component { t t } , a seasonal component { s t } , a calendar component { c t } , and an irregular component { i t } , resulting in the standard unobserved components (UC) decomposition

y t = t t + s t + c t + i t .

(1)

The seasonal variation in Equation (1) is decomposed further according to

s t = ∑ τ ∈ S s t ( τ ) ,

(2)

where S = { τ 1 , τ 2 , … } with τ i ∈ [ 2 , ∞ ] for all i is a finite set of seasonal periodicities and { s t ( τ ) } denotes the seasonal pattern associated with periodicity τ . At this point, there is no need to formulate an explicit model for any UC in Equations (1) and (2). Further specifications will follow below when needed for pretreatment (Section 3) and model-based signal extraction (Sections 4 and 7), bearing in mind that the non-parametric X-11 and STL methods do not require any UC specification at all. Either way, it should be noted that the seasonal patterns in Equation (2) are often estimated sequentially, usually from the linearized observations and in ascending order of their associated periodicities in order to mitigate risks related to confounding.

F ( S ) = ⋃ τ ∈ S { ω k ( τ ) | k ∈ { 1 , … , ⌊ τ / 2 ⌋ } } ,

(3)

where ω k ( τ ) = 2 π k / τ ∈ ( 0 , π ] is the k -th seasonal frequency of { s t ( τ ) } and ⌊ x ⌋ ∈ Z is the largest integer not exceeding x . For example, S = { 24 , 168 , 8765 . 82 } for the hourly electricity consumption series (see Figures 1 and 11), so | F ( S ) | = 12 + 84 + 4 , 382 = 4 , 478 in this case. The denseness of (3) for daily data is nicely visualized in Proietti and Pedregal (2023, Fig. 2).

OPEN IN VIEWER

Figure 2.

Squared gain (8) for τ = 52 . 18 . Gray vertical line marks the chief week-of-the-year frequency. Dashed horizontal line marks 0 . 1 .

3.1. Linear Regression Model

Data pretreatment, or linearization, seeks to solve two problems prior to the application of any seasonal adjustment method: permanent removal of the calendar variation in Equation (1) and temporary outlier correction. Both problems are typically tackled jointly through running a time series regression of the form

y t = x t T β + η t ,

(4)

where x t is a k -dimensional vector of exogenous regression variables associated with calendar events and outliers at time t , β is a k -dimensional vector of unknown calendar and outlier effects, and { η t } is a residual process. To formulate a valid distributional assumption for the residuals in Equation (4), it is noted that they should take the trend-cyclical and seasonal variation present in the observations into account. For monthly and quarterly data, they are therefore modeled usually as a seasonal ARIMA process. For infra-monthly data, however, this is infeasible when S contains fractional periodicities or when | S | > 1 in Equation (2). The feasible alternative implemented in JD+ is a twofold generalization of the famous Airline model that has turned out to fit many monthly and quarterly time series reasonably well despite its rather simple two-parameter structure formed by combining first-order non-seasonal and seasonal differencing and moving average (MA) polynomials. More precisely, this alternative that will be referred to as the extended Airline model (EAM) allows for multiple first-order seasonal polynomials as well as fractional powers of the backshift operator. It is given by

δ 1 ( B ) Π τ ∈ S δ τ ( B ) η t = θ 1 ( B ) Π τ ∈ S θ τ ( B ) ε t ,

(5)

where B is the backshift operator, that is, B η t = η t − 1 , δ 1 ( B ) = 1 − B is the non-seasonal differencing operator, δ τ ( B ) = 1 − B τ is the lag- τ seasonal differencing operator, θ 1 ( B ) = 1 − θ 1 B is the non-seasonal MA operator, θ τ ( B ) = 1 − θ τ B τ is the lag- τ seasonal MA operator and { ε t } is zero-mean white noise (WN) with finite variance σ ε 2 > 0 . All MA operators are assumed to be invertible, that is, | θ 1 | < 1 and | θ τ | < 1 for all τ ∈ S . Fractional powers of B are defined through the first-order Taylor approximation at unity, using the fact that B k = B ⌊ k ⌋ B α k holds for any k ∈ R , where α k = k − ⌊ k ⌋ ∈ [ 0 , 1 ) is the fractional remainder of k . Now, B α k ≈ ( 1 − α k ) + α k B and, hence, B k is essentially approximated with a weighted average of B ⌊ k ⌋ and B ⌊ k ⌋ + 1 :

B k ≈ ( 1 − α k ) B ⌊ k ⌋ + α k B ⌊ k ⌋ + 1 ,

(6)

which is plugged into both the seasonal differencing and MA operators in Equation (5). For weekly data, for instance, we have S = { 52 . 18 } and seasonal differencing becomes

δ 52 . 18 ( B ) η t = ( 1 − B 52 . 18 ) η t = η t − ( 0 . 82 η t − 52 + 0 . 18 η t − 53 ) .

(7)

Note, however, that Equation (6) has no effect when k ∈ Z since α k = 0 in this case.

3.2. Model Properties

Approximation Equation (6) constitutes a fundamental principle of dealing with fractional periodicities. Therefore, it is interesting to study its effect on seasonal differencing from a theoretical standpoint. Moving to the frequency domain, the squared gain of δ τ ( B ) can be written as

| δ τ ( e − i ω ) | 2 = 1 − 2 c τ ( ω ) + [ c τ 2 ( ω ) + s τ 2 ( ω ) ] ,

(8)

where

c τ ( ω ) = ( 1 − α τ ) cos ( ⌊ τ ⌋ ω ) + α τ cos [ ( ⌊ τ ⌋ + 1 ) ω ]

and s τ ( ω ) is defined analogously based on the sine instead of cosine function, which simplify to the well-known standard cases c τ ( ω ) = cos ( τ ω ) and s τ ( ω ) = sin ( τ ω ) when τ ∈ N ∩ [ 2 , ∞ ] . Figure 2 shows the squared gain Equation (8) for the seasonal differencing operator applied in Equation (7), revealing that this operator does not completely annihilate week-of-the-year (WOY) dynamics as the squared gain displays dips at the WOY harmonics but does not reach zero at those dips (in fact, it remains even larger than 0 . 1 at higher WOY harmonics). This flaw is common to all operators δ τ ( B ) with a fractional τ and essentially a consequence of the (first-order) truncation in Equation (6), which results in δ τ ( B ) having only a single unit root at unity. To see this, note that this operator factorizes as

δ τ ( B ) = δ 1 ( B ) S τ ( B ) with S τ ( B ) = 1 + B + ⋯ + B ⌊ τ ⌋ − 1 + α τ B ⌊ τ ⌋ ,

(9)

where S τ ( B ) is the fractional aggregation operator associated with τ . The bottom line is that for fractional periodicities the roots of S τ ( B ) lie neither at the exact seasonal frequencies nor on the unit circle. However, they are still sufficiently close to the seasonal unit root case and therefore relatable to seasonal dynamics.

The principle in Equation (6) also facilitates the derivation of the pseudo-spectral density of the regression residuals { η t } in Equation (4). It follows directly from Equations (5) and (6) that

f η ( ω ) = σ ε 2 2 π × | θ 1 ( e − i ω ) | 2 Π τ ∈ S | θ τ ( e − i ω ) | 2 | δ 1 ( e − i ω ) | 2 Π τ ∈ S | δ τ ( e − i ω ) | 2 , ω ∈ ( 0 , π ] ,

(10)

where the squared gains of the seasonal differencing operators are given by Equation (8) and those of the other involved operators read

| θ 1 ( e − i ω ) | 2 = 1 − 2 θ 1 cos ( ω ) + θ 1 2 , | θ τ ( e − i ω ) | 2 = 1 − 2 θ τ c τ ( ω ) + θ τ 2 [ c τ 2 ( ω ) + s τ 2 ( ω ) ] , | δ 1 ( e − i ω ) | 2 = 2 [ 1 − cos ( ω ) ] .

Figure 3 shows the pseudo-spectral density Equation (10) for the weekly EAM and four combinations of MA parameters. It is seen that, as for the classic Airline model, the non-seasonal MA parameter governs the shape of Equation (10) as ω → 0 (higher/lower values of θ 1 correspond to steeper/flatter increases) whereas the seasonal MA parameter governs the shape of Equation (10) in the vicinity of the twenty-six WOY frequencies defined in Equation (3) (higher/lower values of θ 52 . 18 correspond to narrower/wider peaks). However, in contrast to the classic Airline model, the magnitude of the seasonal peaks decreases dramatically at the higher WOY harmonics due to the aforementioned fact that δ 52 . 18 ( B ) does not entirely annihilate WOY dynamics, recalling that its squared gain appears in the denominator of Equation (10).

OPEN IN VIEWER

Figure 3.

Pseudo-spectral density (10) with S = { 52 . 18 } , θ 1 = 0 . 75 (top row), θ 1 = 0 . 35 (bottom row), θ 52 . 18 = 0 . 85 (left column), θ 52 . 18 = 0 . 45 (right column), and σ ε 2 = 1 . Gray vertical line marks the chief week-of-the-year frequency.

3.3. Estimation and Forecasting

Estimation of pretreatment model Equations (4) and (5) is largely based upon the procedure described in Gómez and Maravall (1994). We briefly point out the basic principles; technical details can be found in Supplemental Appendix A.

When the observations are free of missing values, model Equations (4) and (5) is estimated by maximum likelihood (ML) with the differenced data. To this end, both the generalized least squares estimator of β and the ML estimator of σ ε 2 are expressed as a function of the MA parameters θ = ( θ 1 , θ τ 1 , … , θ τ | S | ) T , so the concentrated log-likelihood can be optimized over θ alone, which is done numerically with the Levenberg-Marquardt algorithm. Plugging the solution θ ^ into the aforementioned functions then provides the ML estimates of β and σ ε 2 . When the observations contain missing values, the “corrected additive outlier” approach developed in Gómez et al. (1999) is used, which requires some minor data manipulations. In essence, a tentative interpolation is calculated for each missing value, and an additive outlier is added in Equation (4) at the corresponding position. Then, model estimation proceeds as before, using appropriate determinantal corrections in the involved likelihood functions.

Model estimation also includes an optional automatic search for additive outliers, level shifts and lag- 1 switch interventions (a special case of reallocation outliers discussed in Wu et al. 1993). The search routine closely follows the robust iterative two-step procedure suggested by Gómez and Maravall (2001a). Starting from the model with only the user-defined regression variables, its first step is a sequential forward addition of the single most significant outlier, including parameter re-estimation. The second step consists of running a multiple regression with all outliers identified in the first step to check if the least significant one of them can be dropped from the model, in which case the first step restarts with the reduced model. This procedure continues until a stable model is found.

Given estimates β ^ and θ ^ , model Equations (4) and (5) can be used to produce forecasts for the observed time series. Let T be the sample size, h > 0 the forecast horizon, and Ω t the information set at time t . Under the assumption that the pretreatment model is correctly specified, the h -step-ahead forecast of the observations follows from Equation (4) as

y ^ T + h = x T + h T β ^ + η ^ T + h ,

(11)

where

η ^ T + h = E ( η T + h | Ω T )

(12)

is the MSE-optimal h -step-ahead linear forecast of the regression residuals derived from Equation (5). For example, when S = { 7 } and h = 1 , then Equation (12) in theory yields

η ^ T + 1 = η T + η T − 6 − η T − 7 − θ 1 ε T − θ 7 ε T − 6 + θ 1 θ 7 ε T − 7 .

(13)

In practice, θ ^ and the estimated innovations are plugged into the right-hand side of Equation (13).

3.4. Generalizations

A legitimate criticism of the EAM Equation (5) is the potential for excessive non-seasonal differencing introduced through the factorization Equation (9) even when the cardinality of S is moderate. To shield the pretreatment model against such over-differencing, JD+ also implements a generalized version of Equation (5) in which the full differencing polynomial is replaced with the tunable autoregressive (AR) integrated polynomial

( 1 − ϕ B ) δ 1 d ( B ) Π τ ∈ S S τ ( B ) ,

(14)

where ϕ ∈ ( 0 , 1 ) to ensure stationarity of the non-seasonal AR ( 1 ) polynomial and d ∈ { 1 , … , 1 + | S | } . The inclusion of the AR( 1 ) factor is optional; if this option is taken, then the non-seasonal MA(1) factor in Equation (5) will be dropped. As regards the user specification of d , Section 8.1 contains a brief discussion about the parameter’s empirical relevance for seasonal adjustment (see Remark 7).

The principle in Equation (6) may provoke another EAM-related criticism as it sacrifices the seasonal unit root properties of S τ ( B ) when τ is fractional (see Remark 4). McElroy and Livsey (2022, App. B) and an anonymous referee suggested using an alternative aggregation operator that retains those properties. This operator can be constructed via an a priori selection of the relevant seasonal frequencies and an atomic bottom-up multiplication of the respective minimal degree polynomials, yielding a complete factorization of each seasonal differencing operator δ τ ( B ) . Adopting notations from Equation (14), the entire differencing polynomial would take the form

δ 1 d ( B ) Π τ ∈ S S τ , h τ ( B ) with S τ , h τ ( B ) = Π j = 1 h τ [ 1 − 2 cos ω j ( τ ) B + B 2 ] ,

(15)

where h τ ≤ ⌊ τ / 2 ⌋ is the number of harmonics to be considered for modeling the seasonal pattern with periodicity τ . The MA polynomial in Equation (5) could be approximated in similar fashion through the operator

θ 1 ( B ) Π τ ∈ S θ τ , h τ ( B ) with θ τ , h τ ( B ) = ( 1 − λ τ B ) Π j = 1 h τ [ 1 − 2 λ τ cos ω j ( τ ) B + λ τ 2 B 2 ] ,

(16)

where λ τ = θ τ 1 / τ (assuming θ τ ≥ 0 ). This operator implies for each seasonal pattern a homogenous treatment of its variance components around the associated seasonal frequencies. As remarked by an Associate Editor, a more flexible parametrization in the spirit of the frequency-specific models discussed in Aston et al. (2007) could lead to further improvements.

5.1. Trend-Cycle Filters

5.1.1. Preliminary Trend-Cycle Estimation

Preliminary trend-cycle extraction in the classic X-11 method is carried out with centered symmetric 2 × τ moving averages for τ ∈ { 4 , 12 } . The extended X-11 approach utilizes a generalized version of those crude filters that allows for any seasonal periodicity τ . Let N odd be the set of odd (positive) integers and ⌈ k ⌉ odd denote the smallest odd (positive) integer not smaller than k , that is, ⌈ k ⌉ odd = min { l ≥ k | l ∈ N odd } . Then, the generalized 2 × τ moving average has length l τ = ⌈ τ ⌉ odd and its weights are given by

w i = { 1 / τ , i ∈ { 1 , 1 τ } with τ ∈ ℕ odd ( I { τ even } + α τ ) / ( 2 τ ) , i ∈ { 1 , 1 τ } with τ ∉ ℕ odd 1 / τ , i ∈ { 2 , … , 1 τ − 1 } ,

(17)

where I { · } is the indicator function. It is easily seen that ∑ i = 1 l τ w i = 1 holds for any choice of τ , see also Table 1 for some examples. Asymmetric variants are not used during preliminary trend-cycle extraction so that ( l τ − 1 ) / 2 observations are temporarily lost at the beginning and end of the filter output.

OPEN IN VIEWER

Table 1.

Weights (17) of Generalized 2 × τ Moving Averages for Selected Seasonal Periodicities (Rounded to Four Decimal Places).

Seasonal periodicity τ	Filter length l τ	End-points weights	Inner weights
7	7	0 . 1429	0 . 1429
30 . 4369	31	0 . 0236	0 . 0329
52 . 1775	53	0 . 0113	0 . 0192
365 . 2425	367	0 . 0003	0 . 0027

5.1.2. Refined and Final Trend-Cycle Estimation

The refined and final trend-cycle extraction filters of the classic X-11 method are essentially a set of prespecified weights for symmetric m -term Henderson filters (Henderson 1916) with m ∈ { 3 , … , 101 } ∩ N odd and their asymmetric Musgrave surrogates (Musgrave 1964; see also Doherty 2001). The extended X-11 approach implements the generalized approach to trend-cycle estimation derived in Proietti and Luati (2008), alongside three ways of calculating asymmetric variants.

Trend-cycle estimation according to Proietti and Luati (2008) is founded on applying local polynomial regressions to the input series. Let b denote the local bandwidth, or half-length, of the trend-cycle extraction filter and d the polynomial degree which the filter should be able to preserve. Also, let q ∈ { 0 , … , b } be the number of future observations available after the current time point t . The local polynomial regression model is then given by

z t = X q ζ + ε t ,

(18)

where z t = ( z t − b , … , z t + q ) T is the local span of the input series,

X q = ( 1 − b b 2 ⋯ ( − b ) d 1 − ( b − 1 ) ( b − 1 ) 2 ⋯ [ − ( b − 1 ) ] d ⋮ ⋮ ⋮ ⋱ ⋮ 1 q q 2 ⋯ q d )

is the time-constant design matrix, ζ = ( ζ 0 , … , ζ d ) T is a ( d + 1 ) -dimensional vector of unknown regression coefficients, and ε t = ( ε t − b , … , ε t + q ) T is a local vector of zero-mean Gaussian WN processes that are both mutually and serially uncorrelated. A weighted least squares (WLS) estimator of ζ is obtained through minimization of the objective function

∑ j = − b q κ j ( z t + j − ζ 0 − ζ 1 j − ⋯ − ζ d j d ) 2 ,

(19)

where { κ j } is a set of non-negative weights defined through a prespecified kernel. The WLS minimizer of (19) is found as

ζ ^ = ( X q T K q X q ) − 1 X q T K q z t

with K q = diag ( κ − b , … , κ q ) . Since the local trend-cycle estimator is given by ζ ^ 0 = e 1 , d + 1 T ζ ^ , where e k , n is the k -th unit vector of length n , we finally have ζ ^ 0 = w q T z t with

w q = K q X q ( X q Τ K q X q ) − 1 e 1 , d + 1 .

(20)

Assuming 2 b ≥ d , the weights of the symmetric kernel-based trend-cycle extraction filters result from setting q = b in (20). Figure 4 shows the weights w b and squared gains of the cubic filters derived through the Henderson, Epanechnikov, and unit-variance Gaussian kernels for bandwidths b ∈ { 4 , 6 } , which correspond to the lengths of the classic 9-term and 13-term trend-cycle extraction filters (see, e.g., Wand and Jones 1995 and Ghosh 2018 for a general introduction to kernel smoothing). For either bandwidth, the squared gains of the three filters are shaped very similarly, although the pass-band becomes slightly narrower when the Epanechnikov and Gaussian kernels are utilized in place of the Henderson kernel. Still, any of these filters qualifies for trend-cycle extraction from daily data with pronounced day-of-the-week dynamics.

OPEN IN VIEWER

Figure 4.

Coefficient curves (left panel) and squared gains (right panel) of symmetric kernel-based Henderson (solid lines), Epanechnikov (dashed lines), and unit-variance Gaussian (dotted lines) trend-cycle extraction filters (20) with q = b = 4 (black lines) and q = b = 6 (gray lines) and d = 3 . Gray vertical line (right panel) marks the chief day-of-the-week frequency.

The weights of the corresponding asymmetric variants can be obtained for any choice of q ∈ { 0 , … , b − 1 } in Equation (20). Those variants are conceptually consistent with the symmetric filters w b —and therefore called “direct” asymmetric filters—but often yield strongly localized and hence volatile trend-cycle estimates. The “cut-and-normalize” approach (Gasser and Müller 1979) is a rather simple alternative according to which the asymmetric variants are found through dropping the unneeded weights from the symmetric filter w b and dividing the remaining weights by their sum, that is

w q = w b ( q ) × [ 1 T w b ( q ) ] − 1 ,

where w b ( q ) contains only the first b + q + 1 elements of w b and 1 is a ( b + q + 1 ) -dimensional vector of ones. A more sophisticated alternative is the “minimum mean squared revision error” (MMSRE) approach also discussed in Proietti and Luati (2008). Aiming at a bias-variance trade-off, its key idea is to purposely introduce a slight bias into the trend-cycle estimators by assuming that the polynomial degree is of some lower order d ¯ < d at the boundary and by seeking asymmetric trend-cycle filters that preserve polynomials up to an even lower degree d ¯ < d ¯ subject to the reproduction constraints dictated by the symmetric filter w b . The Musgrave surrogates, for example, can be replicated within this framework through the linear-constant case ( d ¯ , d ¯ ) = ( 1 , 0 ) if the Henderson kernel is used to calculate w b according to Equation (20). The corresponding objective function, that is, the analog to the matrix representation of Equation (19), is then composed of the variance and squared bias of the trend-cycle revision error plus a vector of Lagrange multipliers linked to the reproduction constraints. The current implementation of the MMSRE approach utilizes a generalized minimization criterion derived in Grun-Rehomme et al. (2018), which essentially translates the aforementioned objective function to the frequency domain where it is decomposable into effects related to the gain and phase shift properties of the involved trend-cycle extraction filters and hence capable of balancing certain fidelity, smoothness, and timeliness criteria.

5.2. Seasonal Filters

The classic X-11 method implements a set of predefined symmetric and asymmetric 3 × k seasonal moving averages, where k ∈ { 1 , 3 , 5 , 9 , 15 } , in order to extract estimates of the seasonal component from the detrended data (alias the seasonal-irregular component). The extended X-11 approach adopts the exact same seasonal filters but allows for more combinations of different filters when extracting the initial and final seasonal factors in iterations B through D. Seasonal estimates are normalized with the generalized symmetric 2 × τ moving average with weights given in Equation (17) and its asymmetric variants defined implicitly through padding the beginning (end) of a seasonal estimate with ( l τ − 1 ) / 2 copies of its first (last) element. Principle (6) is used again in the presence of fractional periodicities, that is, weighted averages of the two corresponding adjacent integer-valued lagged values of the seasonal-irregular component constitute the input to the seasonal filter. Considering the symmetric 3 × 3 seasonal filter, for example, the estimate of the seasonal pattern is obtained from the detrended data { ( si ^ ) t } according to

s ^ t = 1 9 [ α 2 τ ( si ^ ) t − ⌊ 2 τ ⌋ − 1 + ( 1 − α 2 τ ) ( si ^ ) t − ⌊ 2 τ ⌋ ] + 2 9 [ α τ ( si ^ ) t − ⌊ τ ⌋ − 1 + ( 1 − α τ ) ( si ^ ) t − ⌊ τ ⌋ ] + 3 9 ( si ^ ) t + 2 9 [ ( 1 − α τ ) ( si ^ ) t + ⌊ τ ⌋ + α τ ( si ^ ) t + ⌊ τ ⌋ + 1 ] + 1 9 [ ( 1 − α 2 τ ) ( si ^ ) t + ⌊ 2 τ ⌋ + α 2 τ ( si ^ ) t + ⌊ 2 τ ⌋ + 1 ] .

(21)

Figure 5 shows the squared gain of the (effective) seasonal extraction filter Equation (21) for τ = 30 . 44 , that is, of the filter that extracts the day-of-the-month (DOM) pattern from detrended daily data (other examples are shown in Supplemental Appendix B). The squared gain is almost one at the chief DOM frequency but from the first harmonic onward the DOM peaks start to decline visibly in almost linear fashion toward 0 . 25 . At the same time, the average level of the squared gain at the infra-DOM frequencies is visibly different from zero between the higher DOM harmonics (with initially small infra-DOM peaks reaching almost 0 . 25 between the thirteenth and fourteenth harmonic). A similar effect has already been seen for the squared gains of seasonal differencing operators with fractional periodicities (see, e.g., Figure 2). Asymmetric seasonal filters are also obtained from applying the weights of the corresponding classic X-11 filters to the averaged seasonal-irregular component. For example, the concurrent variant of Equation (21) reads

s ^ t = 5 27 [ α 2 τ ( si ^ ) t − ⌊ 2 τ ⌋ − 1 + ( 1 − α 2 τ ) ( si ^ ) t − ⌊ 2 τ ⌋ ] + 11 27 [ α τ ( si ^ ) t − ⌊ τ ⌋ − 1 + ( 1 − α τ ) ( si ^ ) t − ⌊ τ ⌋ ] + 11 27 ( si ^ ) t ,

(22)

the squared gain of which is shown in Figure 6 for τ = 30 . 44 (left panel). Compared with the squared gain of the symmetric 3 × 3 filter, the latter is similar in shape but never drops below 0 . 02 —an effect similar to the one shown in Figure 2—and has slightly higher and broader peaks at the DOM and “in-between-DOM” frequencies. Figure 6 (right panel) also reveals that the concurrent 3 × 3 seasonal extraction filter Equation (22) introduces a phase delay of up to twenty-four days into the low-frequency variation of the filtered series, which may be neglected since the filter is applied to detrended data. Apart from that, no serious phase delays or advances are introduced beyond the chief DOM frequency.

OPEN IN VIEWER

Figure 5.

Squared gains of the symmetric 3 × 3 seasonal extraction filter Equation (21) for τ = 30 . 44 (solid line) and τ = 31 (dotted line). Gray vertical line marks the chief day-of-the-month frequency.

OPEN IN VIEWER

Figure 6.

Squared gain (left panel) and phase shift in days (right panel) of the concurrent 3 × 3 seasonal extraction filter Equation (22) for τ = 30 . 44 . Gray vertical line marks the chief day-of-the-month frequency. Dashed horizontal line (left panel) marks 0 . 02 .

8.1. Daily Births in France

Figure 7 displays the daily BIRTHS series, covering the sample period of January 1, 1968 through December 31, 2020 (two top panels). The seasonal profile of the data is assessed through a sequence of point-wise generalized Canova-Hansen (CH) test statistics (Busetti and Harvey 2003; Canova and Hansen 1995) for the null hypothesis of no seasonality against the alternative of either deterministic or non-stationary stochastic trigonometric seasonality (row- 3 panel) and through grouped boxplots of the raw BIRTHS series (bottom panels). The generalized CH statistics peak clearly at 7 days and around 365 days, suggesting the presence of day-of-the-week (DOW) and day-of-the-year (DOY) patterns. The grouped boxplots confirm this suggestion: birth counts are relatively stable from Monday through Friday but visibly lower on the weekends, especially on Sundays; they also seem to be somewhat higher during spring and summer compared to fall and winter. Based on this visual evidence, we set S = { 7 , 365 . 2425 } in Equation (2) in our subsequent analysis.

OPEN IN VIEWER

Figure 7.

Seasonal profile of the BIRTHS series: raw data over the full span (top panel) and 2015 to 2020 span (row 2), point-wise Canova-Hansen test statistics for τ ∈ { 2 , … , 367 } over full and shortened vertical axis (row 3), and boxplots of the raw data grouped by the day of the week (bottom left panel) and by the month of the year (bottom right panel).

8.1.1. Data Pretreatment and AMB Signal Extraction

We assume a multiplicative UC decomposition in Equation (1) and utilize pretreatment model Equations (4) and (5) to correct the BIRTHS series for the effects of the eleven public holidays in France, including automatic detection of additive outliers and lag- 1 switch interventions with series-length-adjusted 0 . 01 -level critical t -values obtained from the U.S. Census Bureau’s modification of the asymptotic formula derived in Ljung (1993). Table 3 reports the estimated calendar effects along with standard errors and t -statistics. Overall, each public holiday has a strong dampening effect on the birth counts, which amounts, for example, to roughly − 20 % on Christmas Day. In addition, a total of forty-four outliers has been detected automatically but the estimated effects are not reported here for the sake of brevity. The estimated EAM is given by

OPEN IN VIEWER

Table 3.

Estimated Calendar Effects for the BIRTHS Series.

Event	Date	Estimate	S.E.	t -Value
New Year’s Day	Jan 1	− 0 . 163	0 . 0057	− 28 . 5965
Easter Monday	Moving	− 0 . 179	0 . 0043	− 41 . 6279
Labor Day	May 1	− 0 . 116	0 . 0055	− 21 . 0909
Victory Day	May 8	− 0 . 129	0 . 0063	− 20 . 4762
Ascension Day	Moving	− 0 . 159	0 . 0044	− 36 . 1364
Whit Monday	Moving	− 0 . 172	0 . 0043	− 40 . 0000
Bastille Day	July 14	− 0 . 126	0 . 0056	− 22 . 5000
Assumption Day	Aug 15	− 0 . 112	0 . 0055	− 20 . 3636
All Saints’ Day	Nov 1	− 0 . 136	0 . 0055	− 24 . 7273
Armistice Day	Nov 11	− 0 . 140	0 . 0056	− 25 . 0000
Christmas Day	Dec 25	− 0 . 201	0 . 0057	− 35 . 2632

δ ( B ) w t = ( 1 + 0 . 262 ( 0 . 006 ) B ) ( 1 − 0 . 823 ( 0 . 006 ) B 7 ) ( 1 − 0 . 828 ( 0 . 005 ) B 365 . 2425 ) ε t ,

(28)

where δ ( B ) = δ 1 2 ( B ) S 7 ( B ) S 365 . 2425 ( B ) , { w t } denotes the linearized logged BIRTHS series, and standard errors are given in parentheses underneath the MA parameter estimates.

Model Equation (28) also lays the foundation for the sequential extraction of the DOW and DOY patterns with the extended AMB approach. To minimize risks associated with confounding, the DOW pattern is extracted first from the linearized logged BIRTHS series. Afterward, the DOY pattern is extracted from the DOW-adjusted linearized logged BIRTHS series. The estimated DOW pattern has a distinct time-varying amplitude that is largest during the mid-1990s whereas the estimated DOY pattern has a relatively constant amplitude but apparently undergoes some gradual structural infra-yearly changes after 1990 (see Figure 10 below). Both could be a consequence of the attempt to capture the changing volatility in the BIRTHS series during the 1980 to 2010 period. This possible explanation is supported by the fact that the volatility of the seasonally adjusted BIRTHS series remains relatively constant throughout the entire data span (Figure 8).

OPEN IN VIEWER

Figure 8.

Raw (gray line) and seasonally adjusted BIRTHS series obtained with the extended AMB approach (black line).

Remark 7. An Associate Editor suggested to elaborate on the potential excessive non-seasonal differencing that may result in the EAM version Equation (5) and on its empirical relevance for seasonal adjustment. To this end, we compare the AMB trend-cycle estimates obtained from using d ∈ { 1 , 2 , 3 } in Equation (14). Figure 9 shows that the three estimates evolve without remarkable differences over the entire sample span (left panel); some minor deviations, however, can be spotted around turning points and at the end of the sample span where the estimated first-order-integrated trend-cycle differs notably from the other two estimates (right panel). In fact, the pairwise root mean squared differences between the estimated trend-cycles amount to 5 . 282 counts ( d = 1 vs. d = 2 ), 5 . 881 counts ( d = 1 vs. d = 3 ), and 2 . 455 counts ( d = 2 vs. d = 3 ) over the entire sample span. Those similarities might be explained by the fact that each seasonal MA polynomial in EAM version Equation (5) admits a factorization similar to that in Equation (16), so the ( 1 − λ τ B ) factor could approximately cancel with a δ 1 ( B ) operator when λ τ is close to 1 . It follows from (28) that λ 7 = 0 . 823 1 / 7 ≈ 0 . 973 and λ 365 . 2425 = 0 . 828 1 / 365 . 2425 ≈ 0 . 999 when d = 2 , which could justify using d = 1 in Equation (14) instead. An analogous point can be made when running pretreatment with d = 3 in Equation (14) as the resulting estimates of the seasonal MA parameters are in this case θ ^ 7 = 0 . 955 and θ ^ 365 . 2425 = 0 . 887 . Since the choice of d affects forecasting, we also evaluate the out-of-sample forecasting performance of the three EAM variants. To this end, we drop the 366 observations recorded in 2020 from the BIRTHS series, rerun data pretreatment with d ∈ { 1 , 2 , 3 } in Equation (14) on the shortened series and generate from those models h -step-ahead forecasts via Equations (11) and (12) for h ∈ { 1 , … , 366 } . Using the dropped counts as a benchmark, the root mean squared forecast errors over the 2020 horizon are given by 147 . 069 counts ( d = 1 ) , 92 . 232 counts ( d = 2 ) , and 99 . 446 counts ( d = 3 ) . Overall, d = 2 seems to be the most appropriate choice for the BIRTHS series; nevertheless, using EAM-based pretreatment with full non-seasonal differencing ( d = 3 ) apparently does not dramatically deteriorate the final UC estimates obtained with the extended AMB approach.

OPEN IN VIEWER

Figure 9.

Raw BIRTHS series (gray line) and AMB trend-cycle estimates obtained with d = 1 (solid black line), d = 2 (dashed black line), and d = 3 (dotted black line) in (14) over the full span (left panel) and 2020 span (right panel).

8.1.2. Comparison with STS Approach

The classic Airline model for monthly and quarterly data is well-known to provide a good approximation to the reduced form of the so-called basic structural model (BSM) under certain parameter constraints (see, e.g., Maravall 1985). We are interested in finding out whether this property carries on to the EAM and specify the following additive BSM for the linearized logged BIRTHS series: the trend-cyclical component follows the local linear trend Equations (23) and (24); the DOW and DOY patterns are modeled according to the West-Harrison representation Equation (25) with trigonometric terms; finally, the irregular component in Equation (1) is assumed to be Gaussian WN with finite variance σ i 2 and to be uncorrelated with all UC innovations. When estimating this BSM with the KFS using a diffuse square-root initialization, the smoothed DOY estimates turned out to evolve rather erratically over time and even displayed some phases of explosive behavior during the second half of the sample span. We experimented with shorter sample spans and reduced sets of DOY harmonics; however, none of those modifications provided substantially more meaningful results. Therefore, we finally decided to replace the West-Harrison model Equation (25) for the DOY pattern with a (slightly) damped trigonometric cycle at the chief DOY frequency. Using notations from Equations (2) and (3), this DOY cycle is modeled in analogy to Equation (27) as

( s t ( τ ) s ~ t ( τ ) ) = ρ ( cos ω 1 ( τ ) sin ω 1 ( τ ) − sin ω 1 ( τ ) cos ω 1 ( τ ) ) ( s t − 1 ( τ ) s ~ t − 1 ( τ ) ) + ( κ t κ ~ t ) ,

where τ = 365 . 2425 , ρ = 0 . 98 , and { κ t } and { κ ~ t } are two uncorrelated Gaussian WN processes with zero means and common finite variance σ κ 2 .

The estimated UC innovation variances for this modified BSM are given by σ ^ ξ 2 = 6 . 63 × 10 − 4 , σ ^ χ 2 = 0 (which turns the local linear trend into a local level with drift), σ ^ ω 2 = 1 . 40 × 10 − 5 , σ ^ κ 2 = 3 . 19 × 10 − 2 , and σ ^ i = 0 . 317 . Accordingly, the q -ratios, that is, those estimates expressed as a fraction of their largest value, are 0 . 002 for the level innovations, zero for the slope innovations, 4 . 37 × 10 − 5 for the DOW innovations, 0 . 101 for the two DOY innovations, and 1 for the WN irregular component. The smoothed UC estimates obtained for the linearized logged BIRTHS series are finally rescaled using the exponential function. The AMB and BSM estimates of the DOW pattern are almost identical (Figure 10, left panels); their mean squared difference calculated over the entire sample span is 4 . 86 × 10 − 8 . In contrast, the AMB and BSM estimates of the DOY pattern differ noticeably. In particular, the smoothed BSM estimate displays slightly more time variation in its amplitude during the pre-1990 phase but not the aforementioned structural infra-year changes observable in the AMB estimates during the post-1990 span. The ratio between the two DOY estimates occasionally amounts up to 10 % in absolute size (Figure 10, right panels), yet their mean squared difference is 4 . 20 × 10 − 4 .

OPEN IN VIEWER

Figure 10.

AMB (top panels) and STS (middle panels) signal estimates for the BIRTHS series, and their ratios (bottom panels): day-of-the-week pattern from July through December 2020 (left panels) and day-of-the-year pattern (right panels) including annual averages (gray line).

8.2. Hourly Realized Electricity Consumption in Germany

The hourly ELCON series is considered in units of gigawatt hours (GWh) for the period of January 1, 2015 (00:00 AM) through June 30, 2023 (11:00 PM). It covers electricity supplied to the network for the general supply, excluding electricity supplied to the railroad network and to internal industrial and closed distribution networks as well as electricity consumed by the producers. Since daylight saving time (DST) is in effect in Germany, the series is irregularly spaced, albeit slightly: twenty-four observations are available for each non-DST day; twenty-three observations are available for each DST starting day, that is, the last Sunday in March when the 02:00 AM record is skipped; twenty-five observations are available for each DST ending day, that is, the last Sunday in October when the 02:00 AM value is recorded twice. To regularize the ELCON series, we interpolate the missing 02:00 AM value for each DST starting day by the average of the adjacent 01:00 and 03:00 AM records and average the two 02:00 AM values for each DST ending day.

Figure 11 shows the seasonal profile of the regularized ELCON series. Electricity consumption tends to be higher in winter than in summer, leading to a U -shaped infra-yearly pattern that is accompanied by V -shaped dips around the turn of the years (top panel). There are two noteworthy deviations from this regular behavior: the first one is the unusually rapid spring decline in 2020 initiated by the start of the COVID-19 lockdown on March 23, 2020; the second one is the unusually low level of electricity consumption during the 2022/2023 winter that was stimulated by a nationwide energy saving campaign in order to balance reduced gas imports from Russia. The ELCON series also exhibits strong infra-daily and infra-weekly dynamics (see also Figure 1). During a given Monday through Friday sequence, the median hourly consumption settles between 45 and 50  GWh before 05:00 AM, then raises up to 70  GWh until lunchtime, and eventually decreases gradually over the afternoon and less gradually after 08:00 PM (Figure 11, middle panel). The same tendency can be seen on Saturdays and Sundays, although its overall curvature is lower in level and noticeably flatter. The existence of those two patterns is confirmed by the periodogram of the differenced logged ELCON series that displays visible peaks at the corresponding fundamental frequencies and their harmonics (bottom panel).

OPEN IN VIEWER

Figure 11.

Seasonal profile of the ELCON series: raw data (top panel), boxplots grouped by the hour of the day and the day of the week (middle panel) and periodogram of the differenced logged data (bottom panel). Gray vertical lines (bottom panel) mark the chief hour-of-the-day (solid line) and hour-of-the-week (dashed line) frequencies.

8.2.1. Periodic Data Pretreatment

Preliminary empirical studies not discussed here failed to provide sufficient evidence in favor of the assumption of constant hourly calendar effects within each day in the logged regularized ELCON series. To allow for varying hourly calendar effects, we could use an appropriate set of hourly dummy regression variables in Equation (4) and set S = { 24 , 168 , 8765 . 82 } in Equation (14). However, the former set would be almost prohibitively large since nearly 400 hourly dummy variables would be required to account for calendar variation alone. Consequently, such an approach would be computationally inefficient and most likely lead to rather unstable estimates. For that reason, we adopt a periodic regression approach instead, in line with several other studies of hourly time series (e.g., Cancelo et al. 2008; Ramanathan et al. 1997). To this end, we create twenty-four daily sub-series { y t ( h ) } , where the h -th sub-series contains logged regularized electricity consumption observed during hour h ∈ { 0 , … , 23 } for each day in the sample. Each daily sub-series is then linearized separately, using in Equation (4) the automatic outlier detection routine and the common set of seventeen daily (non-DST) calendar regression variables that have been applied in an earlier study of the companion daily ELCON series (Webel 2022). More specifically, the calendar regression variables account for the effects of the following ten fixed and five moving holidays in Germany, plus two associated bridging Fridays: New Year’s Day (January 1), Epiphany (January 6), Good Friday, Easter Monday, Labor Day (May 1), Ascension, Ascension Friday, Pentecost Monday, Corpus Christi, Corpus Christi Friday, German Unification Day (October 3), 500 -th Reformation Day (October 31, 2017), All Saints’ Day (November 1), Christmas Eve (December 24), Christmas Day (December 25), Boxing Day (December 26), and New Year’s Eve (December 31). Long-term means have been removed from the moving-holiday regression variables to avoid confounding with annual seasonality. The residuals of each daily sub-series are modeled according to Equation (14) with d = 2 and S = { 7 , 365 . 2425 } , so that the entire periodic pretreatment model is given by

δ ( B ) w t ( h ) = ( 1 − θ 1 ( h ) B ) ( 1 − θ 7 ( h ) B 7 ) ( 1 − θ 365 . 2425 ( h ) B 365 . 2425 ) ε t ( h ) , h ∈ { 0 , … , 23 } ,

(29)

where δ ( B ) was defined in Equation (28) and { w t ( h ) } is the linearized h -th daily sub-series given by

w t ( h ) = y t ( h ) − x t T β ( h ) .

(30)

Note that in Equations (29) and (30) infra-daily seasonality is implicitly accounted for by the varying seasonal parameters across the twenty-four daily sub-models. Figure 12 plots the estimates of four calendar effects contained in β ( h ) and of the two seasonal MA parameters in Equation (29) against h ∈ { 0 , … , 23 } , along with point-wise standard errors (the time-varying estimates of all other effects are shown in Supplemental Appendix C; detailed results can also be retrieved from the associated R script). The estimated calendar effects display similar infra-daily curvatures that in most cases deviate visibly from the constant benchmark effects estimated for the companion daily ELCON series. The largest (absolute) estimates can usually be observed between 06:00 and 09:00 AM and—to a lesser extent—in the afternoon roughly around 03:00 PM whereas unsurprisingly the lowest (absolute) estimates tend to be recorded in the evening and during nighttime. It should be kept in mind, however, that these estimated periodic holiday effects have been obtained from only eight to nine observations for each holiday, with the obvious exception of the 500 -th Reformation Day. The estimated seasonal MA parameters remain relatively constant throughout the day, being larger than 0 . 8 for each hour of the day and hence indicating presence of strong DOW and DOY patterns in the twenty-four daily sub-series { y t ( h ) } . Their linearized versions { w t ( h ) } are then mingled across time according to

… , w t − 1 ( 0 ) , w t − 1 ( 1 ) , … , w t − 1 ( 23 ) , w t ( 0 ) , w t ( 1 ) , … , w t ( 23 ) , w t + 1 ( 0 ) , w t + 1 ( 1 ) , … , w t + 1 ( 23 ) , …

and finally transformed back to the original ELCON scale to form the precleaned hourly ELCON series.

OPEN IN VIEWER

Figure 12.

Selected estimated calendar effects and MA parameters from periodic pretreatment model Equations (29) and (30): Easter Monday (top left panel), Labor Day (top right panel), Christmas Day (middle left panel), New Year’s Eve (middle right panel), θ 7 ( h ) (bottom left panel), and θ 365 . 2425 ( h ) (bottom right panel). Shaded areas mark point-wise ± 1 S.E. intervals. Dashed horizontal lines mark parameter estimates obtained for the companion daily ELCON series (Webel 2022).