Seasonal Adjustment of Infra-Monthly Time Series with JDemetra+

Abstract

Both the ongoing digital transformation in many statistical agencies around the world and the COVID-19 pandemic outbreak in 2020 have increased the recognition of and demand for infra-monthly economic time series in official statistics during the past decade. Infra-monthly data often display complex forms of seasonality, such as superimposed seasonal patterns with potentially fractional periodicities, that prevent the application of traditional modeling and seasonal adjustment approaches. For that reason, JDemetra+, the official software for harmonized seasonal adjustments of monthly and quarterly data in Europe, has been augmented recently with several methods tailored to the specifics of infra-monthly data. This includes a modified TRAMO-like routine for data pretreatment and extended versions of the ARIMA model-based and X-11 seasonal adjustment approaches alongside the classic STL method. These methods can be accessed through either the graphical user interface or an R package suite, which also provides additional routines for structural time series modeling. We give a comprehensive description of all those methods and discuss the theoretical properties of their key modifications. Selected capabilities are then illustrated using daily births in France, hourly realized electricity consumption in Germany, and weekly initial claims for unemployment insurance in the United States.

1. Motivation

Many macroeconomic time series compiled and published by statistical agencies at regular sub-annual intervals display seasonal fluctuations, that is, movements that recur each year in the same season with similar intensity. The presence of such repetitive movements often obscures underlying economic trends and other non-seasonal dynamics that are usually of more interest for the majority of data users. For that reason, seasonal adjustment, that is, the process that seeks to remove seasonal fluctuations from observed time series as best as possible, has been an integral part of data production in official statistics for many decades. Fostering harmonized production processes across statistical agencies, both Eurostat, which is the statistical office of the European Union, and the European Central Bank have recommended JDemetra+ (JD+) as the time series software for conducting seasonal adjustment of official statistics within the European Statistical System and the European System of Central Banks since February 2015 (see Eurostat and European Central Bank 2015). The JD+ program is mainly developed at the National Bank of Belgium in cooperation with the Deutsche Bundesbank and the French National Institute for Statistics (Insee). It implements the popular X-11 and autoregressive integrated moving average (ARIMA) model-based (AMB) seasonal adjustment approaches for monthly and quarterly data under version branch 2.x, including their respective regARIMA and TRAMO time series regression frameworks for data pretreatment (see Grudkowska 2017 and Smyk et al. 2024 for a documentation of JD+).

A new chapter has begun with the inaugural JD+ release from the version branch 3.x in May 2023, as version 3.0 contains generalized modeling and seasonal adjustment methods for time series of any seasonal periodicity. The addition of those methods was primarily motivated by an increasing popularity of and demand for infra-monthly time series in official statistics over the past decade. This evolution has been largely fostered by the emergence of new digital data sources and automated data collection methods that provide easy and sometimes almost real-time access to such data. As a consequence, many statistical agencies have started to explore ways to replace or complement traditional measurement models with new models that are capable of blending survey data with administrative and digital data sources (cf. Jarmin 2019; Radermacher 2019). In 2020, the demand for daily and weekly data soared immediately after the outbreak of the COVID-19 pandemic as many institutional data users asked for more timely indicators to monitor current economic developments, often with the aid of so-called experimental statistics such as dashboards as well as early warning and sentiment indicators (see, e.g., Eraslan and Götz 2021; Keane and Neal 2021; Lewis et al. 2022). Although being strongly seasonal in many cases, infra-monthly economic data do not lend themselves to being modeled with classic Box-Jenkins approaches and seasonally adjusted with traditional methods due to a variety of peculiarities not seen in monthly and quarterly data. Prime examples include the coexistence of multiple seasonal patterns with potentially fractional periodicities, direct measurability of granular calendar variation and a high susceptibility to various forms of irregular variation (see the discussion in Section 2).

In general, there are three possible strategies in official statistics to resolve those issues: first, proper data regularization so that traditional modeling and seasonal adjustment approaches become applicable to the regularized infra-monthly data; second, proper modifications to traditional approaches so that they become applicable to the unaltered infra-monthly data; third, application of proper approaches not considered in official statistics so far. Following the second strategy for the most part, JD+ 3.0 implements a modified TRAMO-like pretreatment routine as well as extended versions of the AMB and X-11 seasonal adjustment approaches. Those methods are accompanied by the classic STL approach that comes without any modification for dealing with fractional periodicities. The corresponding routines have been implemented in Java and can be accessed easily through both the graphical user interface and an R package suite, which also provides additional routines for estimating structural time series models.

The purpose of this paper is to elaborate on the main modifications to the traditional modeling and seasonal adjustment approaches, to give additional theoretical details, and to discuss and illustrate selected properties of the key underlying models and filtering methods. Our target audience consists of both experienced and potentially new JD+ users who wish to find an updated documentation of the software’s current developmental stage, and—more generally—of anyone interested in recent advances in seasonal adjustment methods for infra-monthly data. The latter type of reader is also referred to the broader overviews provided in Evans et al. (2021), Ladiray et al. (2018), and Webel (2022). We start with a brief review of the main characteristics of infra-monthly time series, including the introduction of our baseline unobserved components model (Section 2). We then discuss data pretreatment through a time series regression model in which the residuals follow an Airline-type ARIMA model that allows for multiple first-order seasonal differencing and moving average polynomials as well as fractional powers of the backshift operator (Section 3). The same model is then used to introduce proper modifications to the traditional AMB approach (Section 4). The extended X-11 and classic STL approaches are presented in Sections 5 and 6, respectively, followed by a brief description of structural time series models (Section 7). Selected capabilities of these methods are then illustrated using daily births in France, hourly realized electricity consumption in Germany, and weekly initial claims for unemployment insurance in the United States (Section 8). R code snippets that can be used to rerun these examples are contained in the Supplemental Material. We conclude with some final remarks and suggestions for future developments (Section 9).

2. Infra-Monthly Time Series

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.

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=\sum _{\tau \in \mathbb{S}}{s}_t^{(\tau )},$$

(2)

where $\mathbb{S}=\{{\tau }_1,{\tau }_2,\dots \}$ with ${\tau }_i\in [2,\infty ]$ for all $i$ is a finite set of seasonal periodicities and $\{s_t^{(\tau )}\}$ denotes the seasonal pattern associated with periodicity $\tau$ . 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.

• Remark 1. The seasonal periodicities in Equation (2) are given by the long-term average numbers of observations per reference period of interest, such as a week, month, year, etc. Accordingly, some of them are fractional due to the existence of leap years and the varying lengths of months in the Gregorian calendar, which will be considered throughout this paper. More precisely, the $400$ -year cycle of the Gregorian calendar is composed of $303$ non-leap years with $365$  days and $97$ leap years with $366$  days, so the average number of days per year is $(303\times 365+97\times 366)/400=365.2425$ . Consequently, the average numbers of days per month and weeks per year are $365.2425/12=30.43688$ and $365.2425/7=52.1775$ , respectively.

• Remark 2. The time-domain complexity of the seasonal component Equation (2) translates into a dense set of seasonal frequencies in the frequency domain. To see this, note that this set is given by

$$\mathcal{F}(\mathbb{S})=\bigcup _{\tau \in \mathbb{S}}\{{\omega }_k^{(\tau )}|k\in \{1,\dots ,\lfloor \tau /2\rfloor \}\},$$

(3)

where ${\omega }_k^{(\tau )}=2\pi k/\tau \in (0,\pi ]$ is the $k$ -th seasonal frequency of $\{s_t^{(\tau )}\}$ and $\lfloor x\rfloor \in \mathbb{Z}$ is the largest integer not exceeding $x$ . For example, $\mathbb{S}=\{24,168,8765.82\}$ for the hourly electricity consumption series (see Figures 1 and 11), so $|\mathcal{F}(\mathbb{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).

• Remark 3. Regularization provides an alternative strategy for handling data with fractional seasonal periodicities. Instead of incorporating fractional periodicities directly into a model for the original data, regularization intends to construct modified data with a constant integer number of (partly artificial) observations per reference period. For example, when setting up a structural time series model for daily Dutch tax revenues, Koopman and Ooms (2003, 2006) transform the time axis and introduce artificial mid-month missing values so that each month consists of twenty-three banking days; those artificial missing values are then imputed straightforwardly during model estimation with the Kalman filter and smoother. In a similar way, Ollech (2021) uses cubic spline interpolation when modeling daily currency in circulation: each month is temporarily treated as if it had thirty-one days; leap years are temporarily shortened by skipping the records made on February, 29 (the resulting gaps are then filled with interpolations during the final calculation of the seasonally adjusted data).

Figure 2.

Squared gain (8) for $\tau =52.18$ . Gray vertical line marks the chief week-of-the-year frequency. Dashed horizontal line marks $0.1$ .

3. Data Pretreatment

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={\mathit{x}}_t^{\mathrm{T}}\mathbf{\beta }+{\eta }_t,$$

(4)

where ${\mathit{x}}_t$ is a $k$ -dimensional vector of exogenous regression variables associated with calendar events and outliers at time $t$ , $\mathbf{\beta }$ is a $k$ -dimensional vector of unknown calendar and outlier effects, and $\{{\eta }_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 $\mathbb{S}$ contains fractional periodicities or when $|\mathbb{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

$${\delta }_1(B)\underset{\tau \in \mathbb{S}}{\Pi }{\delta }_{\tau }\left(B\right){\eta }_t={\theta }_1(B)\underset{\tau \in \mathbb{S}}{\Pi }{\theta }_{\tau }\left(B\right){\epsilon }_t,$$

(5)

where $B$ is the backshift operator, that is, $B{\eta }_t={\eta }_{t-1}$ , ${\delta }_1\left(B\right)=1-B$ is the non-seasonal differencing operator, ${\delta }_{\tau }\left(B\right)=1-B^{\tau }$ is the lag- $\tau$ seasonal differencing operator, ${\theta }_1\left(B\right)=1-{\theta }_{1}B$ is the non-seasonal MA operator, ${\theta }_{\tau }\left(B\right)=1-{\theta }_{\tau }{B}^{\tau }$ is the lag- $\tau$ seasonal MA operator and $\{{\epsilon }_t\}$ is zero-mean white noise (WN) with finite variance ${\sigma }_{\epsilon }^2>0$ . All MA operators are assumed to be invertible, that is, $|{\theta }_1|<1$ and $|{\theta }_{\tau }|<1$ for all $\tau \in \mathbb{S}$ . Fractional powers of $B$ are defined through the first-order Taylor approximation at unity, using the fact that $B^k=B^{\lfloor k\rfloor }{B}^{{\alpha }_k}$ holds for any $k\in \mathbb{R}$ , where ${\alpha }_k=k-\lfloor k\rfloor \in [0,1)$ is the fractional remainder of $k$ . Now, $B^{{\alpha }_k}\approx (1-{\alpha }_k)+{\alpha }_{k}B$ and, hence, $B^k$ is essentially approximated with a weighted average of $B^{\lfloor k\rfloor }$ and $B^{\lfloor k\rfloor +1}$ :

$$B^k\approx (1-{\alpha }_k)B^{\lfloor k\rfloor }+{\alpha }_k{B}^{\lfloor k\rfloor +1},$$

(6)

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

$${\delta }_{52.18}\left(B\right){\eta }_t=(1-B^{52.18}){\eta }_t={\eta }_t-(0.82{\eta }_{t-52}+0.18{\eta }_{t-53}).$$

(7)

Note, however, that Equation (6) has no effect when $k\in \mathbb{Z}$ since ${\alpha }_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 ${\delta }_{\tau }\left(B\right)$ can be written as

$${\left|{\delta }_{\tau }\left(e^{-i\omega }\right)\right|}^2=1-2c_{\tau }\left(\omega \right)+[c_{\tau }^2\left(\omega \right)+s_{\tau }^2\left(\omega \right)],$$

(8)

where

$$c_{\tau }\left(\omega \right)=(1-{\alpha }_{\tau })\cos (\lfloor \tau \rfloor \omega )+{\alpha }_{\tau }\cos \left[(\lfloor \tau \rfloor +1)\omega \right]$$

and $s_{\tau }\left(\omega \right)$ is defined analogously based on the sine instead of cosine function, which simplify to the well-known standard cases $c_{\tau }\left(\omega \right)=\cos \left(\tau \omega \right)$ and $s_{\tau }\left(\omega \right)=\sin \left(\tau \omega \right)$ when $\tau \in \mathbb{N}\cap [2,\infty ]$ . 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 ${\delta }_{\tau }\left(B\right)$ with a fractional $\tau$ and essentially a consequence of the (first-order) truncation in Equation (6), which results in ${\delta }_{\tau }\left(B\right)$ having only a single unit root at unity. To see this, note that this operator factorizes as

$${\delta }_{\tau }\left(B\right)={\delta }_1(B)S_{\tau }\left(B\right)\mathrm{with}{S}_{\tau }\left(B\right)=1+B+\cdots +B^{\lfloor \tau \rfloor -1}+{\alpha }_{\tau }{B}^{\lfloor \tau \rfloor },$$

(9)

where $S_{\tau }\left(B\right)$ is the fractional aggregation operator associated with $\tau$ . The bottom line is that for fractional periodicities the roots of $S_{\tau }\left(B\right)$ 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.

• Remark 4. When $\tau \in \mathbb{N}\cap [2,\infty )$ and hence ${\alpha }_{\tau }=0$ , then $S_{\tau }\left(B\right)$ carries $\lfloor (\tau -1)/2\rfloor$ pairs of complex conjugate unit roots that are equally spaced around the unit circle at the seasonal frequencies given in Equation (3) and accompanied by an additional single real unit root at $-1$ if $\tau$ is even. However, when $\tau \in \overline{\mathbb{N}}\cap [2,\infty )$ and hence ${\alpha }_{\tau }>0$ , then $S_{\tau }\left(B\right)$ must contain one additional root since the polynomial order of Equation (9) has been increased by one compared to the aforementioned case. In fact, the real root at $-1$ that is present in Equation (9) when $\tau \in \overline{\mathbb{N}}\cap [2,\infty )$ is even turns into a complex root with an associated frequency close to $\pi$ when $\tau \in \overline{\mathbb{N}}\cap [2,\infty )$ and $\lfloor \tau \rfloor$ is even. Similarly, an additional negative real root is introduced when $\tau \in \overline{\mathbb{N}}\cap [2,\infty )$ and $\lfloor \tau \rfloor$ is odd compared to the case in which $\tau \in \overline{\mathbb{N}}\cap [2,\infty )$ is odd. Numerical root calculations for fractional periodicities based upon the Jenkins-Traub algorithm (Jenkins and Traub 1972) suggest that the roots of $S_{\tau }\left(B\right)$ are associated with frequencies slightly larger than the seasonal ones and have moduli slightly larger than $1$ . For example, when $\tau =52.18$ , the differences in radians between the frequencies of the complex roots of $S_{\tau }\left(B\right)$ and the respective WOY frequencies are found to increase from $5.2\times {10}^{-7}$ at ${\omega }_1^{(\tau )}$ to $9.8\times {10}^{-3}$ at ${\omega }_{26}^{(\tau )}$ while the corresponding moduli gradually increase from $1.000021$ to $1.008634$ . Conceptually, the principle in Equation (6) sacrifices the seasonal unit root properties of $S_{\tau }\left(B\right)$ in Equation (9) when $\tau$ is fractional.

The principle in Equation (6) also facilitates the derivation of the pseudo-spectral density of the regression residuals $\left\{{\eta }_t\right\}$ in Equation (4). It follows directly from Equations (5) and (6) that

$$f_{\eta }\left(\omega \right)=\frac{{\sigma }_{\epsilon }^2}{2\pi }\times \frac{{\left|{\theta }_1\left(e^{-i\omega }\right)\right|}^2\underset{\tau \in \mathbb{S}}{\Pi }{\left|{\theta }_{\tau }\left(e^{-i\omega }\right)\right|}^2}{{\left|{\delta }_1\left(e^{-i\omega }\right)\right|}^2\underset{\tau \in \mathbb{S}}{\Pi }{\left|{\delta }_{\tau }\left(e^{-i\omega }\right)\right|}^2},\omega \in (0,\pi ],$$

(10)

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

$$\begin{array}{c}{\left|{\theta }_1\left(e^{-i\omega }\right)\right|}^2=1-2{\theta }_1\cos \left(\omega \right)+{\theta }_1^2,\hfill \\ {\left|{\theta }_{\tau }\left(e^{-i\omega }\right)\right|}^2=1-2{\theta }_{\tau }{c}_{\tau }\left(\omega \right)+{\theta }_{\tau }^2[c_{\tau }^2\left(\omega \right)+s_{\tau }^2\left(\omega \right)],\hfill \\ {\left|{\delta }_1\left(e^{-i\omega }\right)\right|}^2=2[1-\cos \left(\omega \right)].\hfill \end{array}$$

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 $\omega \to 0$ (higher/lower values of ${\theta }_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 ${\theta }_{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 ${\delta }_{52.18}\left(B\right)$ does not entirely annihilate WOY dynamics, recalling that its squared gain appears in the denominator of Equation (10).

Figure 3.

Pseudo-spectral density (10) with $\mathbb{S}=\{52.18\}$ , ${\theta }_1=0.75$ (top row), ${\theta }_1=0.35$ (bottom row), ${\theta }_{52.18}=0.85$ (left column), ${\theta }_{52.18}=0.45$ (right column), and ${\sigma }_{\epsilon }^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 $\mathbf{\beta }$ and the ML estimator of ${\sigma }_{\epsilon }^2$ are expressed as a function of the MA parameters $\mathbf{\theta }=({\theta }_1,{\theta }_{{\tau }_1},\dots ,{\theta }_{{\tau }_{|\mathbb{S}|}}{)}^{\mathrm{T}}$ , so the concentrated log-likelihood can be optimized over $\mathbf{\theta }$ alone, which is done numerically with the Levenberg-Marquardt algorithm. Plugging the solution $\hat{\mathbf{\theta }}$ into the aforementioned functions then provides the ML estimates of $\mathbf{\beta }$ and ${\sigma }_{\epsilon }^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 $\hat{\beta }$ and $\hat{\theta }$ , 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 ${\Omega }_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

$${\hat{y}}_{T+h}={\mathit{x}}_{T+h}^{\mathrm{T}}\hat{\beta }+{\hat{\eta }}_{T+h},$$

(11)

where

$${\hat{\eta }}_{T+h}=\mathbb{E}({\eta }_{T+h}|{\Omega }_T)$$

(12)

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

$${\hat{\eta }}_{T+1}={\eta }_T+{\eta }_{T-6}-{\eta }_{T-7}-{\theta }_1{\epsilon }_T-{\theta }_7{\epsilon }_{T-6}+{\theta }_1{\theta }_7{\epsilon }_{T-7}.$$

(13)

In practice, $\hat{\theta }$ 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 $\mathbb{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-\varphi B){\delta }_1^d\left(B\right)\underset{\tau \in \mathbb{S}}{\Pi }{S}_{\tau }\left(B\right),$$

(14)

where $\varphi \in (0,1)$ to ensure stationarity of the non-seasonal AR $\left(1\right)$ polynomial and $d\in \{1,\dots ,1+|\mathbb{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_{\tau }\left(B\right)$ when $\tau$ 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 ${\delta }_{\tau }\left(B\right)$ . Adopting notations from Equation (14), the entire differencing polynomial would take the form

$${\delta }_1^d\left(B\right)\underset{\tau \in \mathbb{S}}{\Pi }{S}_{\tau ,h_{\tau }}\left(B\right)\mathrm{with}{S}_{\tau ,h_{\tau }}\left(B\right)=\underset{j=1}{\overset{h_{\tau }}{\Pi }}[1-2\cos {\omega }_j^{(\tau )}B+B^2],$$

(15)

where $h_{\tau }\le \lfloor \tau /2\rfloor$ is the number of harmonics to be considered for modeling the seasonal pattern with periodicity $\tau$ . The MA polynomial in Equation (5) could be approximated in similar fashion through the operator

$${\theta }_1\left(B\right)\underset{\tau \in \mathbb{S}}{\Pi }{\theta }_{\tau ,h_{\tau }}\left(B\right)\mathrm{with}{\theta }_{\tau ,h_{\tau }}\left(B\right)=(1-{\lambda }_{\tau }B)\underset{j=1}{\overset{h_{\tau }}{\Pi }}[1-2{\lambda }_{\tau }\cos {\omega }_j^{(\tau )}B+{\lambda }_{\tau }^2{B}^2],$$

(16)

where ${\lambda }_{\tau }={\theta }_{\tau }^{1/\tau }$ (assuming ${\theta }_{\tau }\ge 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.

4. Extended ARIMA Model-Based Approach

Classic AMB seasonal adjustment for monthly and quarterly time series as implemented in the SEATS program is essentially organized in three steps; see, for example, Gómez and Maravall (2001b) and Maravall (1995) for extensive overviews and, in particular, the underlying modeling assumptions. In the first step, a seasonal ARIMA model is identified for the linearized observations. The default choice is the model found for the ARIMA errors by the TRAMO pretreatment routine but sometimes minor modifications, such as changes to the polynomial orders or variance inflation, need to be made. In a second step, the requested UC models are derived from the data model by means of the canonical decomposition (Hillmer and Tiao 1982), which basically postulates the uniqueness constraint that all (additive) WN extractable from the UC models is assigned to the WN irregular component, maximizing the latter’s variance amongst all admissible decompositions. To this end, the AR and differencing operators are first factorized through polynomial division and the radian frequencies of the roots obtained this way dictate their allocation to the UC models. Afterward, the pseudo autocovariance generating function (ACGF) of the data model is decomposed using a partial fraction expansion. In a third step, drawing on ideas of Burman (1980), Wiener-Kolmogorov (WK) filters are formed from the canonically decomposed model to finally calculate the minimum mean squared error (MMSE) estimates for the signals of interest; see Wiener (1949) for a general exposition of WK theory.

The extension of the classic AMB method follows the same key steps, using either the EAM Equation (5) or its generalization Equation (14) as a starting point. Thus, the unit roots at unity will be assigned to the trend-cycle while the non-stationary $({\alpha }_{\tau }=0)$ or stationary $({\alpha }_{\tau }>0)$ roots of each fractional aggregation operator will be assigned to the corresponding seasonal pattern. Using Equation (14), the derived UC models will therefore take the general form

$$(1-\varphi B){\delta }_1^d\left(B\right)t_t={\theta }^{(t)}\left(B\right){\epsilon }_t^{(t)}$$

and

$$S_{\tau }\left(B\right)s_t^{(\tau )}={\theta }^{(\tau )}\left(B\right){\epsilon }_t^{(\tau )},\tau \in \mathbb{S},$$

where the orders of the MA polynomials ${\theta }^{(t)}\left(B\right)$ and ${\theta }^{(\tau )}\left(B\right)$ match the orders of their corresponding combined AR and differencing operators, bearing in mind that the inclusion of the $(1-\varphi B)$ factor in Equation (14) and hence its presence in the derived trend-cycle model is optional. For large $\tau$ , however, the involved polynomials are very long, which leads to computational problems when performing the canonical decomposition with the legacy algorithms. The crucial step that requires a more stable solution is the factorization of a symmetric polynomial (to be associated with the ACGF of an MA process) into its roots; the implemented modification utilizes the eigenvalues of the polynomial’s companion matrix (cf. Aurentz et al. 2015) and applies Leja ordering (Leja 1957) to retrieve more stable coefficients from the found roots.

In some cases, the decomposition of the data model can be inadmissible in the sense that not all components have a non-negative spectrum for all frequencies. Hillmer and Tiao (1982) prove that an admissible decomposition exists if and only if the sum of the spectral minima of the components obtained through a unique partial fraction decomposition of the data model’s pseudo ACGF, where uniqueness is achieved by imposing specific constraints on the polynomial orders of the UC models, is non-negative. Inadmissibility can thus be seen to be equivalent to a negative variance of the WN irregular component in this unique decomposition. To solve this issue, the classic AMB approach implements several strategies including variance inflation, which is also used here. This strategy sets the variance of the WN irregular component in the aforementioned unique partial fraction decomposition to zero, essentially eliminating this component. As a result, the sum of the remaining trend-cyclical and seasonal components yields a new and decomposable data model from which the canonical decomposition and the requisite WK filters can be computed. Note that such a collapse of the WN irregular component may also happen in other strategies for overcoming inadmissibility; a prime example is the solution based upon the specification of top-heavy UC models (Fiorentini and Planas 2001).

Another modification to the classic AMB approach concerns MMSE estimation. WK filtering is predicated on the assumption of an invertible data model; however, the long MA polynomials that arise from large seasonal periodicities are typically prone to contain near-seasonal unit root factors that cause problems with respect to the convergence of the WK filters. Therefore, the extended AMB approach utilizes the state space representation of the canonically decomposed EAM and the Kalman filter and smoother (KFS), which avoids such non-invertibility issues and, in theory, yields the same results as the WK filter. The smoothing algorithm is implicitly selected by the user through the specification as to whether or not standard deviations of the components should be provided: if so, the diffuse square-root KFS is applied; otherwise, the fast disturbance smoother (Koopman 1993) with a diffuse initialization is used. Either way, the filtering recursions produce back- and forecasts for the UC estimates, which can then be aggregated to obtain alternative model-based back- and forecasts for the observed time series.

Signal extraction with the extended AMB approach can be performed either simultaneously or sequentially. A simultaneous extraction produces MMSE estimates for each UC in a single run; in this case, the application of the more robust square-root KFS is recommended, especially when $d\ge 3$ in Equation (14). A sequential extraction targets one seasonal pattern at a time, implying $d=2$ (and the absence of the non-seasonal AR(1) factor) in Equation (14) for each run. This option can be advantageous in terms of RAM usage and computing times as the involved state space models are less complex.

5. Extended X-11 Approach

The classic X-11 seasonal adjustment method derived in Shiskin et al. (1967) and reviewed in great detail in Ladiray and Quenneville (2001) is based upon an iterative application of predefined linear filters, which produces a sequence of preliminary, refined and final estimates for the trend-cyclical, seasonal and irregular components in Equation (1). During each of the X-11 iterations (labeled B through D), the following key four-step principle is applied twice. First, the trend-cycle is extracted from the input series through the application of a $2\times \tau$ moving average or a Henderson filter, depending on the position in the entire filtering process. Second, the estimated trend-cycle is removed from the input series. Third, the detrended input series is separated into normalized seasonal and irregular components through the period-wise application of $3\times k$ seasonal moving averages. At some positions in the X-11 routine, this step also includes an automatic detection and down-weighting of extreme values in the irregular component based upon the famous $\sigma$ -limits. Fourth, the estimated seasonal component is removed from the input series. Once those iterations are completed, the final trend-cycle and irregular estimates are obtained from splitting the final estimate of the seasonally adjusted series with a Henderson filter.

The extended X-11 approach essentially adopts this iterative filtering process but incorporates some generalizations and modifications tailored to the peculiarities of infra-monthly data. The current implementation facilitates a sequential extraction of the seasonal patterns in Equation (2) without an automatic selection of trend-cycle and seasonal filters based on the $I/C$ and $I/S$ ratios and without any forecasts, even naive ones, of the seasonal patterns.

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\times \tau$ moving averages for $\tau \in \{4,12\}$ . The extended X-11 approach utilizes a generalized version of those crude filters that allows for any seasonal periodicity $\tau$ . Let ${\mathbb{N}}_{\mathrm{odd}}$ be the set of odd (positive) integers and $\lceil k{\rceil }^{\mathrm{odd}}$ denote the smallest odd (positive) integer not smaller than $k$ , that is, $\lceil k{\rceil }^{\mathrm{odd}}=min\{l\ge k|l\in {\mathbb{N}}_{\mathrm{odd}}\}$ . Then, the generalized $2\times \tau$ moving average has length $l_{\tau }=\lceil \tau {\rceil }^{\mathrm{odd}}$ and its weights are given by

$$w_i=\{\begin{array}{l}1/\tau ,i\in \left\{{1,1}_{\tau }\right\}\mathrm{with}\tau \in {\mathbb{N}}_{\mathrm{odd}}\\ \left(\mathbb{I}\left\{\tau \mathrm{even}\right\}+{\alpha }_{\tau }\right)/\left(2\tau \right),\mathrm{i}\in \left\{{1,1}_{\tau }\right\}\mathrm{with}\tau \notin {\mathbb{N}}_{\mathrm{odd}}\\ 1/\tau ,i\in \left\{2,\dots ,1_{\tau }-1\right\}\end{array},$$

(17)

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

Table 1.

Weights (17) of Generalized $2\times \tau$ Moving Averages for Selected Seasonal Periodicities (Rounded to Four Decimal Places).

Seasonal periodicity $\tau$	Filter length $l_{\tau }$	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\in \{3,\dots ,101\}\cap {\mathbb{N}}_{\mathrm{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\in \{0,\dots ,b\}$ be the number of future observations available after the current time point $t$ . The local polynomial regression model is then given by

$${\mathbf{z}}_t={\mathbf{X}}_q\mathbf{\zeta }+{\mathbf{\epsilon }}_{\mathbf{t}}\mathbf{,}$$

(18)

where ${\mathbf{z}}_t={(z_{t-b},\dots ,z_{t+q})}^{\mathrm{T}}$ is the local span of the input series,

$${\mathbf{X}}_q=\left(\begin{array}{c}1-b{b}^2\cdots {(-b)}^d\hfill \\ \multicolumn{1}{c}{1-(b-1){(b-1)}^2\cdots {[-(b-1)]}^d}\\ \multicolumn{1}{c}{⋮⋮⋮\ddots ⋮}\\ \multicolumn{1}{c}{1q{q}^2\cdots q^d}\end{array}\right)$$

is the time-constant design matrix, $\mathbf{\zeta }\mathbf{=}\mathbf{(}{\mathbf{\zeta }}_0\mathbf{,}\dots \mathbf{,}{\mathbf{\zeta }}_{\mathbf{d}}{\mathbf{)}}^{\mathrm{T}}$ is a $(d+1)$ -dimensional vector of unknown regression coefficients, and ${\mathbf{\epsilon }}_{\mathbf{t}}=({\epsilon }_{t-b},\dots ,{\epsilon }_{t+q}{)}^{\mathrm{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 $\mathbf{\zeta }$ is obtained through minimization of the objective function

$$\sum _{j=-b}^q{\kappa }_j{(z_{t+j}-{\zeta }_0-{\zeta }_{1}j-\cdots -{\zeta }_d{j}^d)}^2,$$

(19)

where $\left\{{\kappa }_j\right\}$ is a set of non-negative weights defined through a prespecified kernel. The WLS minimizer of (19) is found as

$$\hat{\zeta }={\left({\mathbf{X}}_q^{\mathrm{T}}{\mathbf{K}}_q{\mathbf{X}}_q\right)}^{-1}{\mathbf{X}}_q^{\mathrm{T}}{\mathbf{K}}_q{\mathbf{z}}_t$$

with ${\mathbf{K}}_q=\mathrm{diag}({\kappa }_{-b},\dots ,{\kappa }_q)$ . Since the local trend-cycle estimator is given by ${\hat{\zeta }}_0={\mathit{e}}_{1,d+1}^{\mathrm{T}}\hat{\zeta }$ , where ${\mathit{e}}_{k,n}$ is the $k$ -th unit vector of length $n$ , we finally have ${\hat{\zeta }}_0={\mathbf{w}}_q^{\mathrm{T}}{\mathbf{z}}_t$ with

$$w_q=K_q{X}_q{(X_q^{{\rm T}}{K}_q{X}_q)}^{-1}{e}_{1,d+1}.$$

(20)

Assuming $2b\ge d$ , the weights of the symmetric kernel-based trend-cycle extraction filters result from setting $q=b$ in (20). Figure 4 shows the weights ${\mathbf{w}}_b$ and squared gains of the cubic filters derived through the Henderson, Epanechnikov, and unit-variance Gaussian kernels for bandwidths $b\in \{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.

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\in \{0,\dots ,b-1\}$ in Equation (20). Those variants are conceptually consistent with the symmetric filters ${\mathbf{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 ${\mathbf{w}}_b$ and dividing the remaining weights by their sum, that is

$${\mathbf{w}}_q={\mathbf{w}}_b^{(q)}\times {\left[1^{\mathrm{T}}{\mathbf{w}}_b^{(q)}\right]}^{-1},$$

where ${\mathbf{w}}_b^{(q)}$ contains only the first $b+q+1$ elements of ${\mathbf{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 $\overline{d}<d$ at the boundary and by seeking asymmetric trend-cycle filters that preserve polynomials up to an even lower degree $\underset{¯}{d}<\overline{d}$ subject to the reproduction constraints dictated by the symmetric filter ${\mathbf{w}}_b$ . The Musgrave surrogates, for example, can be replicated within this framework through the linear-constant case $(\overline{d},\underset{¯}{d})=(1,0)$ if the Henderson kernel is used to calculate ${\mathbf{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\times k$ seasonal moving averages, where $k\in \{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\times \tau$ 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_{\tau }-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\times 3$ seasonal filter, for example, the estimate of the seasonal pattern is obtained from the detrended data $\{{(\widehat{\mathrm{si}})}_t\}$ according to

$$\begin{array}{c}\multicolumn{1}{c}{{\hat{s}}_t=\frac{1}{9}[{\alpha }_{2\tau }{(\widehat{\mathrm{si}})}_{t-\lfloor 2\tau \rfloor -1}+(1-{\alpha }_{2\tau }){(\widehat{\mathrm{si}})}_{t-\lfloor 2\tau \rfloor }]}\\ \multicolumn{1}{c}{+\frac{2}{9}[{\alpha }_{\tau }{(\widehat{\mathrm{si}})}_{t-\lfloor \tau \rfloor -1}+(1-{\alpha }_{\tau }){(\widehat{\mathrm{si}})}_{t-\lfloor \tau \rfloor }]}\\ \multicolumn{1}{c}{+\frac{3}{9}(\widehat{\mathrm{si}}{)}_t}\\ \multicolumn{1}{c}{+\frac{2}{9}[(1-{\alpha }_{\tau }){(\widehat{\mathrm{si}})}_{t+\lfloor \tau \rfloor }+{\alpha }_{\tau }{(\widehat{\mathrm{si}})}_{t+\lfloor \tau \rfloor +1}]}\\ \multicolumn{1}{c}{+\frac{1}{9}[(1-{\alpha }_{2\tau }){(\widehat{\mathrm{si}})}_{t+\lfloor 2\tau \rfloor }+{\alpha }_{2\tau }{(\widehat{\mathrm{si}})}_{t+\lfloor 2\tau \rfloor +1}].}\end{array}$$

(21)

Figure 5 shows the squared gain of the (effective) seasonal extraction filter Equation (21) for $\tau =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

$$\begin{array}{c}\multicolumn{1}{c}{{\hat{s}}_t=\frac{5}{27}[{\alpha }_{2\tau }{(\widehat{\mathrm{si}})}_{t-\lfloor 2\tau \rfloor -1}+(1-{\alpha }_{2\tau }){(\widehat{\mathrm{si}})}_{t-\lfloor 2\tau \rfloor }]}\\ \multicolumn{1}{c}{+\frac{11}{27}[{\alpha }_{\tau }{(\widehat{\mathrm{si}})}_{t-\lfloor \tau \rfloor -1}+(1-{\alpha }_{\tau }){(\widehat{\mathrm{si}})}_{t-\lfloor \tau \rfloor }]}\\ \multicolumn{1}{c}{+\frac{11}{27}(\widehat{\mathrm{si}}{)}_t,}\end{array}$$

(22)

the squared gain of which is shown in Figure 6 for $\tau =30.44$ (left panel). Compared with the squared gain of the symmetric $3\times 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\times 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.

• Remark 5. Seasonal filters that incorporate principle Equation (6) operate at the exact non-integer seasonal periodicities but inevitably need to compromise some extraction power at the higher seasonal harmonics. Other seasonal adjustment approaches for infra-monthly time series possess a high extraction power at all seasonal harmonics at the expense of operating at slightly inaccurate integer periodicities, which is often achieved through temporary data regularization (see Remark 3). For example, the STL-based approach for daily data implemented in the {dsa} R package (Ollech 2021) utilizes spline interpolation in order to temporarily stretch each month to thirty-one days. Setting $\tau =31$ instead of $\tau =30.44$ in Equation (21) mimics the corresponding 5-term seasonal LOESS smoother, and the aforementioned trade-off effect is visualized by the dotted squared gain in Figure 5.

Figure 5.

Squared gains of the symmetric $3\times 3$ seasonal extraction filter Equation (21) for $\tau =30.44$ (solid line) and $\tau =31$ (dotted line). Gray vertical line marks the chief day-of-the-month frequency.

Figure 6.

Squared gain (left panel) and phase shift in days (right panel) of the concurrent $3\times 3$ seasonal extraction filter Equation (22) for $\tau =30.44$ . Gray vertical line marks the chief day-of-the-month frequency. Dashed horizontal line (left panel) marks $0.02$ .

6. STL Approach

The classic STL method developed in R. B. Cleveland et al. (1990) has been conceptualized originally as an alternative to X-11 and, therefore, also rests on the iterative application of linear filters. The filtering operations are organized in the so-called inner and outer loops; broadly speaking, each run through the inner loop produces updated estimates of the trend-cyclical and seasonal components and of the seasonally adjusted data whereas each pass of the outer loop generates an estimate of the irregular component as well as robustness weights.

The inner loop is essentially modeled on X-11’s basic four-step filtering principle outlined at the beginning of Section 5 to obtain the aforementioned refined estimates. However, STL’s fundamental modification is the replacement of X-11’s predefined trend-cycle and seasonal moving averages with linear extraction filters whose weights are derived through locally weighted scatterplot smoothing (LOESS), which is a specific form of a non-parametric univariate regression. The general LOESS objective function is similar to Equation (19) with $d\in \{1,2\}$ and the kernel weights being replaced by the so-called neighborhood weights obtained through the tricube kernel. Another key modification to X-11 is the removal of low-frequency variation from the estimated seasonal component through the additional application of a convoluted MA and LOESS filter. This is done to restrain the final STL trend-cycle and seasonal smoothers from competing for the same variation in the data. In X-11 terminology, the convolution consists of a $3\times \tau \times \tau$ MA filter trailed by a LOESS smoother with a recommended bandwidth given by $\lceil \tau {\rceil }^{\mathrm{odd}}$ . As a result, the squared gain of the final seasonal LOESS smoother remains close to zero in the vicinity of $\omega =0$ , as opposed to the X-11 counterparts (cf. Figure 5). Finally, it is also worth noting that STL does not require predefined asymmetric variants of its trend-cycle and seasonal filters since smoothing near the end-points is handled internally through an appropriate modification of the neighborhood weights in the LOESS regression. The bandwidth of the latter remains fixed at $2b+1$ in the terminology of Equation (19) whereas in X-11 the filter length shrinks gradually from $2b+1$ in the symmetric case to $b+1$ in the concurrent case.

The optional outer loop yields an estimate of the irregular component, calculated as a residual by removing from the observed time series the trend-cycle and seasonal estimates obtained in the preceding inner loop. It also provides robustness weights based on the biweight kernel for the irregular component, which is conceptually akin to X-11’s extreme-value detection. During the next run through the inner loop, the combined neighborhood and robustness weights are then used in the LOESS regressions, that is, the principle of iterated WLS estimation is applied in Equation (19).

The current JD+ implementation is mostly in line with the practical recommendations given in R. B. Cleveland et al. (1990) that target “near certainty of convergence” of the STL trend-cycle and seasonal estimates. By default, two runs through the inner loop are carried out when no robustness weights are needed whereas fifteen runs through the outer loop, each time with a single run through the inner loop, are carried out otherwise. The polynomials in the LOESS objective function have orders $d=1$ for trend-cycle and $d=0$ for seasonal LOESS smoothers. If unspecified, the bandwidth of the trend-cycle LOESS smoother when extracting a seasonal pattern with periodicity $\tau$ is automatically set to $\lceil l_t(\tau ){\rceil }^{\mathrm{odd}}$ with

$$l_t\left(\tau \right)=\lfloor \frac{1.5\lfloor \tau \rfloor }{1-1.5\lfloor l_s\left(\tau \right){\rfloor }^{-1}}\rfloor ,$$

where $l_s\left(\tau \right)$ is the bandwidth of the corresponding seasonal LOESS smoother. The latter also affects the bandwidth of the LOESS smoother in the convoluted low-pass filter, which is given by $\lfloor l_s(\tau ){\rfloor }^{\mathrm{odd}}$ . The numerator of the above formula also alludes to the fact that dealing with fractional periodicities in STL amounts to always using $\lfloor \tau \rfloor$ in place of $\tau$ . Recalling the weighted-averaging principle Equation (6) employed in the extended AMB and X-11 approaches, this pragmatic solution essentially sets ${\alpha }_{\tau }=0$ regardless of $\tau$ . The size of the distortions introduced into the final UC estimates will thus increase with the true value of ${\alpha }_{\tau }$ , which is clearly unsatisfactory from a theoretical standpoint. However, ${\alpha }_{\tau }$ is relatively small for most fractional periodicities relevant to infra-monthly data, so one may wonder if those distortions should really be a matter of concern in practice.

7. Structural Time Series Models

The pretreatment and seasonal adjustment methods discussed in Sections 3 to 6 perform data linearization and signal extraction in two separate steps. Structural time series (STS) models enable simultaneous estimation of all UCs in Equations (1) and (2); however, as opposed to the other methods, they require an a priori specification of the corresponding UC models, the aggregate of which then naturally constitutes the model for the observations (see Harvey 1989 as a standard reference).

In JD+, the most general trend-cycle specification within the STS framework is the local linear model given by

$$t_t=t_{t-1}+{\nu }_{t-1}+{\xi }_t,$$

(23)

$${\nu }_t={\nu }_{t-1}+{\chi }_t,$$

(24)

where $\{{\xi }_t\}$ and $\{{\chi }_t\}$ are mutually uncorrelated Gaussian WN processes with zero means and finite variances ${\sigma }_{\xi }^2$ and ${\sigma }_{\chi }^2$ , respectively. Other popular trend-cycle models can be derived directly from this general form. For example, the local level model with drift, which is characterized by a constant slope $\nu$ in Equation (23), is obtained from setting ${\sigma }_{\chi }^2=0$ in Equation (24). Similarly, the so-called smooth trend is obtained from setting ${\sigma }_{\xi }^2=0$ in (23).

The seasonal component in (1) is specified through the West-Harrison representation (West and Harrison 1997). Assuming $|\mathbb{S}|=1$ and $\tau \in \mathbb{N}\cap [2,\infty )$ in Equation (2) for the ease of exposition (subsequent notations naturally extend to the case of multiple seasonal patterns with potentially fractional periodicities, see Proietti and Pedregal 2023 for more details), the general form reads

$$s_t={\mathbf{e}}_{1,\tau -1}^{\mathrm{T}}{\tilde{\mathbf{s}}}_t,$$

where $\left\{{\tilde{\mathbf{s}}}_t\right\}$ is the $(\tau -1)$ -dimensional reduced form of the $\tau$ seasonal effects given by

$${\tilde{\mathbf{s}}}_t={\mathbf{D}}^{-}{\mathbf{P}}^{t-1}{\mathbf{s}}_t,$$

where ${\mathbf{s}}_t^{\mathrm{T}}=(s_{1,t},\dots ,s_{\tau ,t})$ is the vector of the $\tau$ seasonal effects at time $t$ , $\mathbf{P}$ and $\mathbf{D}$ are permutation (or circular shift) and dimension reduction matrices given by

$$\mathbf{P}=\left(\begin{array}{c}{0}_{\tau -1}{\mathbf{I}}_{\tau -1}\hfill \\ \multicolumn{1}{c}{10_{\tau -1}^{\mathrm{T}}}\end{array}\right)\mathrm{and}\mathbf{D}=\left(\begin{array}{c}{\mathbf{I}}_{\tau -1}\hfill \\ \multicolumn{1}{c}{-1_{\tau -1}^{\mathrm{T}}}\end{array}\right)$$

with $0_n$ and $1_n$ being $n$ -dimensional column vectors of zeros and ones, respectively, and ${\mathbf{I}}_n$ being the $n$ -dimensional identity matrix, and ${\mathbf{D}}^{-}={\left({\mathbf{D}}^{\mathrm{T}}\mathbf{D}\right)}^{-1}{\mathbf{D}}^{\mathrm{T}}$ is the Moore-Penrose inverse of $\mathbf{D}$ . The reduced-form seasonal effects are then assumed to follow the first-order vector autoregressive model

$${\tilde{\mathbf{s}}}_t=\mathbf{T}{\tilde{\mathbf{s}}}_{t-1}+{\omega }_{\mathbf{t}},\mathbf{T}={\mathbf{D}}^{-}\mathbf{PD}=(\begin{array}{c}{0}_{\tau -2}{\mathbf{I}}_{\tau -2}\hfill \\ \multicolumn{1}{c}{-1_{\tau -1}^{\mathrm{T}}}\end{array}),$$

(25)

where $\left\{{\omega }_t\right\}$ is a vector Gaussian WN process with zero mean and non-singular covariance matrix ${\Sigma }_{\mathbf{\omega }}={\sigma }_{\omega }^2\Omega$ , where $\Omega$ is a model-specific matrix that explicitly incorporates a stochastic zero-sum restriction in the $\tau$ seasonal effects. Its exact form is given in Proietti (2000) for the crude, dummy, Harrison-Stevens and trigonometric seasonal models. For example, $\Omega ={\mathbf{D}}^{-}\mathbf{H}{\mathbf{H}}^{\mathrm{T}}{\left({\mathbf{D}}^{-}\right)}^{\mathrm{T}}$ where ${\mathbf{H}}^{\mathrm{T}}=({\mathbf{h}}_1,\dots ,{\mathbf{h}}_{\tau })$ and

$${\mathbf{h}}_k={(\cos \left[{\omega }_1^{\mathbf{(}\mathbf{\tau }\mathbf{)}}\mathbf{k}\right],\sin \left[{\omega }_1^{(\tau )}k\right],\dots ,\cos \left[{\omega }_{\lfloor \tau /2\rfloor }^{(\tau )}k\right],\sin \left[{\omega }_{\lfloor \tau /2\rfloor }^{(\tau )}k\right])}^{\mathrm{T}},k\in \{1,\dots ,\tau \},$$

(26)

for the trigonometric seasonal model, noting that the last sine term in Equation (26) is dropped when $\tau$ is even.

The entire STS model is put into a univariate linear Gaussian state space form, which also enables the inclusion of the calendar variation from Equation (1) by adding the requisite exogenous regression variables from (4) to the system matrix of the observation equation and the associated effects contained in $\mathbf{\beta }$ to the state vector. Outliers can be dealt with similarly, and an automatic detection procedure (Grassi et al. 2018) is available, which is essentially a forward-addition-backward-deletion algorithm based upon the point-wise maximum ${\tau }_t^{{✶}^2}$ -statistics derived in de Jong and Penzer (1998). The full state space model is finally estimated with the KFS with hyper-parameters being obtained through likelihood maximization, the implementation of which closely follows the algorithms described in Durbin and Koopman (2012). When automatic outlier detection is not required, the user can customize model estimation with respect to KFS initialization and likelihood-related aspects (i.e., use of concentration, marginalization, and numerical optimization routine); see, for example, Francke et al. (2010) for further details on those concepts. When automatic outlier detection is required, the Levenberg-Marquardt algorithm is always used for maximizing the profile likelihood function. The latter is favored over the diffuse likelihood mainly due to faster computing times, especially in the presence of additional prespecified regression variables. Moreover, the outliers identified under usage of the profile and diffuse likelihood functions are the same in almost all cases, and the difference in their estimated effects is typically negligible.

• Remark 6. The above West-Harrison representation of the seasonal dynamics in Equation (1) might come across as somewhat unintuitive and opaque. Following the suggestion of an anonymous referee, we note that a simpler and more intuitive specification for each seasonal pattern in Equation (2) is given by $s_t^{(\tau )}=\sum _{j=1}^{h_{\tau }}{s}_{j,t}^{(\tau )}$ , where $h_{\tau }$ is defined below Equation (15),

$$\left(\begin{array}{c}{s}_{j,t}^{(\tau )}\hfill \\ \multicolumn{1}{c}{{\tilde{s}}_{j,t}^{(\tau )}}\end{array}\right)=\left(\begin{array}{c}\cos {\omega }_j^{(\tau )}\sin {\omega }_j^{(\tau )}\hfill \\ \multicolumn{1}{c}{-\sin {\omega }_j^{(\tau )}\cos {\omega }_j^{(\tau )}}\end{array}\right)\left(\begin{array}{c}{s}_{j,t-1}^{(\tau )}\hfill \\ \multicolumn{1}{c}{{\tilde{s}}_{j,t-1}^{(\tau )}}\end{array}\right)+\left(\begin{array}{c}{\kappa }_{j,t}\hfill \\ \multicolumn{1}{c}{{\tilde{\kappa }}_{j,t}}\end{array}\right),$$

(27)

and $\left\{{\kappa }_{j,t}\right\}$ and $\left\{{\tilde{\kappa }}_{j,t}\right\}$ are two mutually and serially uncorrelated Gaussian WN processes with zero means and common finite variance ${\sigma }_{{\kappa }_j}^2$ . In addition, this structural representation is in total agreement with the ARIMA-based bottom-up approach proposed in McElroy and Livsey (2022) and mentioned in Section 3.4 since each harmonic component of the form (27) corresponds to one factor of the seasonal aggregation operator in (15).

8. Applications

We consider three real-world examples to illustrate selected capabilities of the various JD+ methods discussed in Sections 3 to 7. Daily counts of births in France (BIRTHS) are used to compare data pretreatment and model-based signal extraction with the extended AMB and STS approaches for a multiplicative UC decomposition. Hourly realized electricity consumption in Germany (ELCON) is considered to illustrate periodic data pretreatment and to compare non-parametric signal extraction with the extended X-11 and STL methods. Finally, weekly initial claims for unemployment insurance in the United States (CLAIMS) are used for exploring differences between the STS and STL estimates for an additive UC decomposition, focusing on the COVID-19 pandemic phase. Throughout this section, data linearization will be based on the EAM generalization Equation (14) with $d=2$ and without the optional non-seasonal AR( $1$ ) factor whenever $|\mathbb{S}|\ge 2$ . Each example was run on a $64$ -bit Windows OS with 32 GB RAM and different CPU and clock rate configurations; approximate computing times and details about the required R packages are summarized in Table 2. Code snippets, including data and further information about the public data sources, can be found in the Supplemental Material. It is noted that we strongly advise to operate on a system with at least 16 GB RAM when applying the STS approach to seasonal series of the given lengths and periodicities.

Table 2.

Overview of Computational Aspects for the Three Examples.

Series	Sample size	Method	R package	Reference	Computing time
BIRTHS	$19,359$	EAM	{rjd3highfreq}	Palate (2024a)	3–6 minutes
		AMB	{rjd3highfreq}	Palate (2024a)	12–15 seconds
		STS	{rjd3sts}	Palate (2024c)	1–3 seconds
ELCON	$74,472$	Periodic EAM	{rjd3highfreq}	Palate (2024a)	10–15 minutes
		X-11	{rjd3x11plus}	Palate (2024d)	1–2 minutes
		STL	{rjd3stl}	Palate (2024b)	<1 second
CLAIMS	$2,947$	STS	{rjd3sts}	Palate (2024c)	3–5 minutes
		EAM	{rjd3highfreq}	Palate (2024a)	3–5 seconds
		STL	{rjd3stl}	Palate (2024b)	<0.1 second

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 $\mathbb{S}=\{7,365.2425\}$ in Equation (2) in our subsequent analysis.

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 $\tau \in \{2,\dots ,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

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$

$$\delta \left(B\right)w_t=(1+\underset{(0.006)}{0.262}B)(1-\underset{(0.006)}{0.823}{B}^7)(1-\underset{(0.005)}{0.828}{B}^{365.2425}){\epsilon }_t,$$

(28)

where $\delta (B)={\delta }_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).

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\in \{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-{\lambda }_{\tau }B)$ factor could approximately cancel with a ${\delta }_1(B)$ operator when ${\lambda }_{\tau }$ is close to $1$ . It follows from (28) that ${\lambda }_7=0.{823}^{1/7}\approx 0.973$ and ${\lambda }_{365.2425}=0.{828}^{1/365.2425}\approx 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 ${\hat{\theta }}_7=0.955$ and ${\hat{\theta }}_{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\in \{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\in \{1,\dots ,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.

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 ${\sigma }_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

$$\left(\begin{array}{c}{s}_t^{(\tau )}\hfill \\ \multicolumn{1}{c}{{\tilde{s}}_t^{(\tau )}}\end{array}\right)=\rho \left(\begin{array}{c}\cos {\omega }_1^{(\tau )}\sin {\omega }_1^{(\tau )}\hfill \\ \multicolumn{1}{c}{-\sin {\omega }_1^{(\tau )}\cos {\omega }_1^{(\tau )}}\end{array}\right)\left(\begin{array}{c}{s}_{t-1}^{(\tau )}\hfill \\ \multicolumn{1}{c}{{\tilde{s}}_{t-1}^{(\tau )}}\end{array}\right)+\left(\begin{array}{c}{\kappa }_t\hfill \\ \multicolumn{1}{c}{{\tilde{\kappa }}_t}\end{array}\right),$$

where $\tau =365.2425$ , $\rho =0.98$ , and $\{{\kappa }_t\}$ and $\{{\tilde{\kappa }}_t\}$ are two uncorrelated Gaussian WN processes with zero means and common finite variance ${\sigma }_{\kappa }^2$ .

The estimated UC innovation variances for this modified BSM are given by ${\hat{\sigma }}_{\xi }^2=6.63\times {10}^{-4}$ , ${\hat{\sigma }}_{\chi }^2=0$ (which turns the local linear trend into a local level with drift), ${\hat{\sigma }}_{\omega }^2=1.40\times {10}^{-5}$ , ${\hat{\sigma }}_{\kappa }^2=3.19\times {10}^{-2}$ , and ${\hat{\sigma }}_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\times {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\times {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\times {10}^{-4}$ .

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

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 $\mathbb{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\in \{0,\dots ,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 $\mathbb{S}=\{7,365.2425\}$ , so that the entire periodic pretreatment model is given by

$$\delta (B)w_t^{(h)}=(1-{\theta }_1^{(h)}B)(1-{\theta }_7^{(h)}{B}^7)(1-{\theta }_{365.2425}^{(h)}{B}^{365.2425}){\epsilon }_t^{(h)},h\in \{0,\dots ,23\},$$

(29)

where $\delta (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)}-{\mathit{x}}_t^{\mathrm{T}}{\mathbf{\beta }}^{\mathbf{(}\mathbf{h}\mathbf{)}}.$$

(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 ${\mathbf{\beta }}^{\mathbf{(}\mathbf{h}\mathbf{)}}$ and of the two seasonal MA parameters in Equation (29) against $h\in \{0,\dots ,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

$$\dots ,w_{t-1}^{(0)},w_{t-1}^{(1)},\dots ,w_{t-1}^{(23)},w_t^{(0)},w_t^{(1)},\dots ,w_t^{(23)},w_{t+1}^{(0)},w_{t+1}^{(1)},\dots ,w_{t+1}^{(23)},\dots$$

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

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), ${\theta }_7^{(h)}$ (bottom left panel), and ${\theta }_{365.2425}^{(h)}$ (bottom right panel). Shaded areas mark point-wise $\pm 1$ S.E. intervals. Dashed horizontal lines mark parameter estimates obtained for the companion daily ELCON series (Webel 2022).

8.2.2. Non-Parametric Signal Extraction

Based on visual evidence (Figure 11), we now set $\mathbb{S}=\{24,168,8765.82\}$ for the precleaned ELCON series—bearing in mind that the infra-daily and infra-weekly periodicities are not prime—and extract the three patterns sequentially (again in ascending order) with the extended X-11 and STL approaches, assuming a multiplicative UC decomposition. Each X-11 run utilizes a cubic Henderson trend-cycle filter, the length of which is set to $\lceil \tau {\rceil }^{\mathrm{odd}}$ hours, cut-and-normalize asymmetric variants, $3\times 9$ seasonal filters when extracting the hour-of-the day and hour-of-the-week patterns and $3\times 5$ seasonal filters when extracting the hour-of-the-year pattern, and default $\sigma$ -limits. For STL, we specify trend-cycle and seasonal LOESS smoothers of the same lengths as the corresponding X-11 filters and disregard robustness weights. Figure 13 shows that the estimated X-11 and STL seasonal patterns are quite similar, although the combined adjustment factors occasionally differ by up to 9% in absolute size. Those large yet rare differences may also be spotted in two versions of the seasonally adjusted hourly ELCON series (Figure 14). The effects of the COVID-19 lockdown and the campaigned reduction in electricity consumption during the 2022/2023 winter are also clearly visible.

Figure 13.

X-11 (black line) and STL (gray line) signal estimates for the ELCON series: hour-of-the-day pattern over June 24 to 30, 2023 (top left panel), hour-of-the-week pattern from May through June 2023 (top right panel), hour-of-the-year pattern (bottom left panel), and ratio of the combined adjustment factors (bottom right panel).

Figure 14.

Raw (light gray line) and seasonally adjusted ELCON series obtained through the extended X-11 (black line) and STL (gray line) approaches over the full span (top panel) and 2022 to 2023 span (bottom panel).

8.3. Weekly Initial Claims for Unemployment Insurance in the United States

The weekly CLAIMS series is considered as of 1967 W01 up to 2023 W25, with economic activity being measured from Sunday through Saturday. Shorter versions of this series have also been analyzed in earlier studies on weekly data (e.g., Abeln and Jacobs 2022; Cleveland et al. 2018; Cleveland and Scott 2007; Evans et al. 2021; Proietti and Pedregal 2023).

The seasonal profile of the CLAIMS series is shown in Figure 15. The claims meander roughly between $200,000$ and $900,000$ cases, exhibiting periodic phases of higher counts during the colder months of November through February (top panels). The series also displays some sharp increases associated with tumultuous economic times, such as the global financial crisis (2008 Q3 through 2009 Q2), and unprecedented high counts immediately after the COVID-19 pandemic outbreak in March 2020. The latter peak slightly above $6.1$  million claims in 2020 W14, temporarily plateau around $1$  million claims as of 2020 W31, and eventually get back to pre-pandemic standards after a gradual decline over the course of the year 2021. The point-wise generalized CH test statistics clearly peak at thirteen and around fifty-three weeks, signaling presence of strong WOY fluctuations (bottom left panel). Grouped weekly boxplots, which have been calculated from the pre-pandemic data span, confirm this impression (bottom right panel). Their $W$ -shaped nature also reveals mild mid-year increases in the CLAIMS series that have been partly explained with annual model change-over practices in the automotive industry (W. P. Cleveland and Scott 2007). Overall, the seasonal profile of the CLAIMS series is characterized by a strong WOY pattern with dominant week-of-the-quarter dynamics, so that we could set either $\mathbb{S}=\{13,52.18\}$ or $\mathbb{S}=\{52.18\}$ in Equation (2); for the sake of sparsity, we opt for the latter in our subsequent analysis.

Figure 15.

Seasonal profile of the CLAIMS series (thousands): raw data over the full span (top left panel) and pre-pandemic span (top right panel), point-wise Canova-Hansen test statistics for $\tau \in \{2,\dots ,60\}$ (bottom left panel) and boxplots grouped by the week of the year (bottom right panel).

8.3.1. Basic Structural Model

We specify an additive BSM for the untransformed CLAIMS series, noting that this choice of UC decomposition is in line with official seasonal adjustment practices at the U.S. Bureau of Labor Statistics (BLS) in the wake of the COVID-19 pandemic outbreak (in fact, the decomposition scheme has been changed temporarily from multiplicative to additive between March 2020 and June 2021). The trend-cycle is modeled as a local linear trend according to Equations (23) and (24), the WOY pattern follows the trigonometric seasonal model Equation (26) using the full set of twenty-six harmonics, and dynamics associated with federal holidays are captured by a set of eleven weekly dummy regression variables (with the effect of New Year’s Day being assigned to W02 as in W. P. Cleveland et al. 2018). Besides an automatic search for additive outliers and level shifts with series-length-adjusted $0.01$ -level critical $t$ -values, two level shift sequences are prespecified to capture the atypical pandemic phase in 2020 (W12 through W24 and W28 through W32). Table 4 (columns 3–5) reports the estimated calendar effects for the CLAIMS series. Thanksgiving and New Year’s Day have the strongest effects whereas the effect of Independence Day is barely significant at the 5% level of significance, and that of Christmas Day is even insignificant. The estimated calendar component is shown in Figure 16, along with the smoothed estimate of the trigonometric WOY pattern (top left panel). Automatic outlier detection identified twelve additive outliers and four additional level shifts. Seven of those outliers lie in the colder months of the pre-1990 era and might be related to exceptionally cold, or warm, winters. One level shift coincides with the maximum CLAIMS counts during the global financial crisis (2009 W02), and four outliers can be associated with the post-pandemic recovery phase (top right panel). The $q$ -ratios are given by $1$ for the level innovations $\{{\xi }_t\}$ in Equation (23), zero for the slope innovations $\left\{{\chi }_t\right\}$ in Equation (24) (turning the local linear trend into a local level with drift as for the BIRTHS example), $1.74\times {10}^{-3}$ for the WOY innovations $\left\{{\omega }_t\right\}$ in (25) and $0.137$ for the WN irregular component in Equation (1). The seasonally adjusted CLAIMS series is also shown in Figure 16 (bottom panels), revealing some movements buried in the unadjusted data. A prime example is the short-lived increase to about $370,000$ seasonally adjusted claims on average during the 2005 W35 to W40 phase, that is, an approximate 10% increase compared to 2005 W34, that results from Hurricane Katrina having struck the U.S. Gulf Coast on August 29, 2005 (see also the discussion in Cleveland et al. 2018).

Table 4.

Estimated Calendar Effects for the CLAIMS Series.

Event	Date	BSM			EAM
		Estimate	S.E.	$t$ -Value	Estimate	S.E.	$t$ -Value
New Year’s Day	Jan 1	$73.030$	$4.4301$	$16.4850$	$75.000$	$4.2672$	$17.5758$
Martin Luther King Day	Moving	$-47.753$	$5.4454$	$-8.7694$	$-50.069$	$5.1051$	$-9.8077$
President’s Day	Moving	$-24.587$	$4.0781$	$-6.0289$	$-20.706$	$4.0181$	$-5.1531$
Easter	Moving	$8.608$	$2.1545$	$3.9952$	$8.457$	$2.4355$	$3.4723$
Memorial Day	Moving	$-36.010$	$4.4547$	$-8.0836$	$-35.477$	$4.2046$	$-8.4376$
Independence Day	July 4	$-10.224$	$4.4282$	$-2.3089$	$-5.674$	$4.2046$	$-1.3492$
Labor Day	Moving	$-34.185$	$4.3589$	$-7.8426$	$-31.340$	$4.1567$	$-7.5396$
Columbus Day	Moving	$-22.980$	$4.3594$	$-5.2713$	$-27.986$	$4.1583$	$-6.7300$
Veterans Day	Nov 11	$-33.414$	$4.0879$	$-8.1737$	$-35.982$	$3.9099$	$-9.2028$
Thanksgiving	Moving	$-82.731$	$4.4711$	$-18.5036$	$-86.892$	$4.2062$	$-20.6582$
Christmas Day	Dec 25	$-7.802$	$4.3650$	$-1.7875$	$-0.015$	$4.1816$	$-0.0037$

Figure 16.

Smoothed estimates for various unobserved components in the CLAIMS series (thousands): calendar component (black line) and trigonometric WOY pattern (gray line; top left panel), outlier component (top right panel), raw (gray line) and seasonally adjusted (black line) data over the full span (bottom left panel) and pre-pandemic span (bottom right panel).

• Remark 8. The main purpose of the CLAIMS example is to illustrate the STS modeling framework currently implemented in JD+. It should be noted, however, that our underlying acceptance of a time series model with constant (regression) parameters may be questioned for such a long data sequence. As pointed out by an anonymous referee, the data generating process of the CLAIMS series has undergone several changes over time. For example, every U.S. federal state has different laws on eligibility for unemployment insurance claims; these laws have changed a lot. In addition, claims needed to be filed in local offices in the distant past whereas online submissions are the norm nowadays. Hence, some holiday effects may have become less significant, or even insignificant, in recent years. Despite these valid concerns, we believe that the above BSM specification is good enough for serving our primary goal here; however, we agree that a more nuanced modeling of the CLAIMS series is advisable when targeting other purposes, such as providing guidance for policy makers. In such a case, separate treatments of shorter non-overlapping data segments, constant-parameter models with additional holiday regression variables for exceptionally early or late occurrences of some moving holidays in some years and day-of-the-week interactions of some fixed holidays, or time series models with time-varying parameters could be considered. Using temporary changes, or other types of transitional outliers, in place of the two level shift sequences could especially improve the modeling of the COVID-19 pandemic phase; in fact, the BLS has implemented such replacements during recent regular revisions of its seasonal adjustment specifications (see U.S. Bureau of Labor Statistics 2024 for more details).

8.3.2. COVID-19 Modeling and Comparison with STL Approach

Another anonymous referee pointed out that according to Abeln and Jacobs (2022) the robust STL method does not require the two level shift sequences, or other forms of outlier modeling, for taking the COVID-19 pandemic phase into account. To check this claim, we compare the smoothed BSM estimates of the seasonally adjusted CLAIMS series to their (additive) non-robust and robust STL counterparts. Utilizing bandwidths over fifty-five and eleven observations for the trend-cycle and seasonal LOESS smoothers, respectively, both STL estimates are calculated from the holiday-corrected CLAIMS series that is still contaminated with all outliers present in the raw data. Note, however, that pretreatment model Equations (4) and (5) is fed not only with the eleven user-defined holiday variables but also with the two level shift sequences and the same parameter choices for automatic outlier detection that have been specified for the BSM approach. This is done to avoid distortions to the estimated holiday effects, which are likely to arise from a willful ignorance of other known exogenous effects, or any other model misspecification. The estimated EAM-based holiday effects are in general very similar to the BSM estimates, with the more obvious insignificance of the Independence Day and Christmas Day effects being the most notable differences (see columns 6–8 of Table 4).

Mirroring Abeln and Jacobs (2022, Fig. 3), Figure 17 shows that the seasonally adjusted CLAIMS series obtained from the BSM and robust STL approaches are visually almost indistinguishable during the COVID-19 pandemic phase (top panel). The only exception are some minor deviations around the 2020/2021 turn of the year, which may be traced back to the fact that in model (4) and (5) three outliers have been automatically detected between 2020 W50 and 2021 W08 that were not identified by the BSM routine (both automatic detection routines find virtually the same outliers in the pre-pandemic span). Emphasizing the importance of applying STL with robustness weights when no outlier pretreatment is conducted, Figure 17 also highlights severe biases in the seasonally adjusted CLAIMS series obtained from the non-robust STL approach (bottom panel). Not only does the level of this seasonally adjusted series exhibit short-lived downward distortions by up to half a million cases in the March and April figures (which even become negative in 2019), it displays very long periods of visible upward biases in the other months as well. These weird and clearly undesired effects extend beyond the 2019 to 2021 span shown in Figure 17 but eventually fade out. Overall, our findings support the above claim made by Abeln and Jacobs (2022) with respect to the robust STL method, although it is not obvious to us how these authors have dealt with the calendar variation in the CLAIMS series.

Figure 17.

Raw (gray line) and seasonally adjusted CLAIMS series (thousands) obtained from the STS (solid black line) and STL (dashed black line) approaches applied to the EAM-based holiday-corrected data: robust (top panel) versus non-robust (bottom panel) STL approach.

9. Summary

We gave an elaborate description of the modeling and signal extraction methods for infra-monthly time series implemented in JDemetra+ 3.0. In particular, we highlighted the key modifications that needed to be made to official statistics’ established pretreatment and seasonal adjustment approaches in order to take the peculiarities of such granular data into account. Pretreatment for outliers and calendar variation is handled through a time series regression in which the disturbances are assumed to follow an extension of the classic Airline model, which enables the incorporation of multiple seasonal patterns and fractional seasonal periodicities. Using a first-order Taylor expansion, fractional powers of the backshift operator are approximated by a weighted average of the two adjacent integer-valued powers of the backshift operator. The same Airline-type model provides the foundation for an extended ARIMA model-based seasonal adjustment approach based on the classic concept of a canonical model decomposition. MSE-optimal estimates for the requisite unobserved components are then obtained through the Kalman filter and smoother applied to the state space representation of the decomposed extended Airline model. The above Taylor approximation is also adopted in an extended X-11 seasonal adjustment approach to retain applicability of the classic $3\times k$ seasonal filters to time series with fractional periodicities. This X-11 method also contains a rich set of kernel-based trend-cycle extraction filters derived from local polynomial regressions that encompasses X-11’s pioneering Henderson filters and Musgrave surrogates. Alternative non-parametric seasonal adjustments that utilize trend-cycle and seasonal LOESS smoothers can also be conducted, the implementation of which essentially emulates the classic STL method except for some minor modifications. Sequential estimation of prior adjustment and seasonal factors is therefore a common strategy for these extended seasonal adjustment methods. However, JDemetra+ 3.0 also implements routines for structural time series modeling that facilitate a simultaneous estimation of these factors. We illustrated selected capabilities of all those methods using around $20,000$ observations for daily birth counts in France, around $75,000$ observations for hourly realized electricity consumption in Germany, and around $3,000$ observations for weekly initial claims for unemployment insurance in the United States.

The methodological toolbox of JDemetra+ 3.0 is already quite rich but still leaves room for enhancements. First of all, some classic concepts that have proven to be useful for the seasonal adjustment of monthly and quarterly data should be made (more easily) accessible for the treatment of infra-monthly data as well. Prime examples are bias corrections for non-linear data transformations and the use of forecast extensions within non-parametric seasonal adjustment methods to enable the application of less asymmetric, or even symmetric, linear extraction filters near the end-points (which is likely to reduce real-time revisions in signal estimates). The addition of temporary changes and/or otherwise smooth transitional outliers, log/level and seasonality tests as well as tailored diagnostics (e.g., McElroy 2021) would definitely help to automate the entire pretreatment process for infra-monthly data. Several alternatives, such as models Equations (15) and (16) with either the full or a reduced set of seasonal harmonics as well as (structural) models with time-varying regression parameters, provide room for further improvements.

Apart from that, future research could also pay attention to adequate data-driven selections of trend-cycle and seasonal extraction filters for non-parametric seasonal adjustment methods. For example, the framework of Proietti and Luati (2008), which is already utilized in the extended X-11 approach, enables the identification of the optimal bandwidth for trend-cycle extraction through cross-validation given a prespecified polynomial degree and set of kernel weights. As regards model-based signal extraction, saturated seasonal models, such as the West-Harrison representation Equation (25), could be over-parametrized and computationally infeasible for large seasonal periodicities (see the BIRTHS example discussed in Section 8.1). Therefore, structural time series models with an alternative representation of seasonality, such as the more intuitive trigonometric specification (27), and/or with seasonal interactions embedded in a correlated unobserved components framework (cf. Hindrayanto et al. 2019; Koopman and Lee 2009) could be more useful in practice and may be added in a future version of JDemetra+. Periodic cubic splines (Harvey and Koopman 1993; Harvey et al. 1997) seem to be another promising alternative; however, a strategy for automatically picking both the number and positions of the constituting knots still needs to be worked out. Some inspiration might be drawn from the LOESS-based approach discussed in Proietti and Pedregal (2023).

Besides methodological improvements, it would also be interesting to compare the signal extraction methods currently implemented in JDemetra+ 3.0 with other methods. Using daily and weekly economic time series, we are currently evaluating the real-time revisions in concurrent seasonal adjustments carried out with the extended ARIMA model-based and X-11 approaches and with the STL-based approach available in the {dsa} R package (Ollech 2021) in an ongoing research project. But more empirical case studies on this subject, which could also utilize additional approaches as well as simulated time series, are certainly needed.

Supplemental Material

sj-pdf-1-jof-10.1177_0282423X241277602 – Supplemental material for Seasonal Adjustment of Infra-Monthly Time Series with JDemetra+

Supplemental material, sj-pdf-1-jof-10.1177_0282423X241277602 for Seasonal Adjustment of Infra-Monthly Time Series with JDemetra+ by Karsten Webel and Anna Smyk in Journal of Official Statistics