key: cord-0820653-xpypu047
authors: Cuevas, Ángel; Ledo, Ramiro; Quilis, Enrique M.
title: Seasonal adjustment of the Spanish sales daily data
date: 2021-11-22
journal: SERIEs (Berl)
DOI: 10.1007/s13209-021-00251-7
sha: 460a204346184fe267ceb5cc1f06865c5514404a
doc_id: 820653
cord_uid: xpypu047

We present a procedure to perform seasonal adjustment over daily sales data. The model adjusts daily information from the Immediate Supply of Information System for Value Added Tax declaration forms compiled by the Spanish Tax Agency. The procedure performs signal extraction and forecasting at the daily frequency, by means of an unobserved components model. The daily information allows a permanently updated monitoring of the short-term economic conditions of the Spanish economy.

Economic data provided by tax sources are gaining popularity among economic analysts and forecasters due to its timely availability, reliability, coverage and direct economic meaning. In this way, tax-based data have a clear function in the design of a model aimed at nowcasting and short-term forecasting.

Updated and reliable forecasts play a critical role for budgetary planning and for the anticipation of risky situations due to adverse shocks. In particular, the daily information provided by the new Immediate Supply of Information system (SII, Suministro Inmediato de Información), based on the Value Added Tax (VAT) forms and developed by the Tax Agency, opens new perspectives for the integration, on a real-time basis, of tax-based reliable quantitative information with macro-data observed at lower frequencies.

This paper presents a system for seasonal adjustment of very high frequency data derived from tax sources. Specifically, the model treats the daily sales time series B Enrique M. Quilis enrique.martinq@correo.aeat.es 1 Research and Statistics Division, Tax Agency, Spain provided by the SII. Seasonal adjustment (SA) of daily data is necessary, in the same way as happens to monthly and quarterly data, in order to provide a meaningful signal of its underlying evolution. However, this task is notably more difficult than in the case of monthly or quarterly data due to the complexity of its seasonal component, formed by several subcomponents, some of them linked to fractional periodicities and its noisy nature (Ladiray et al. 2018) .

The paper is organized as follows. The second section presents the main characteristic of the data. The third section develops the econometric methodology, which is implemented in two stages. In the first one, we apply a preliminary treatment of the deterministic effects and, in the second one, we use an univariate structural model, based on an unobserved component representation, to decompose and forecast the daily sales data. The empirical results are presented in the fourth section. The fifth section is focused in predictive evaluation. The paper ends presenting the main conclusions and future developments.

Daily sales series of monthly VAT taxpayers comes from the Immediate Supply of Information System (SII, Sistema Inmediato de Información), introduced in January 2017 and officially implemented by the Spanish Tax Agency since July 2017 (Tax Agency 2017).

This system allows the exchange of tax information between the Spanish Tax Agency and taxpayers required by the SII practically in real time, by supplying the detail of the invoicing records within four days, through the electronic platform of the Spanish Tax Agency.

In this way, both tax management and tax compliance are improved (e.g., by the taxpayers comparison of the information in their books with the information provided by their customers and suppliers).

The group compulsorily included in the SII is made up of all those taxpayers whose obligation to declare VAT is monthly:

• Large Companies (turnover greater than 6,010,214.04e in the previous year).

• Those companies that pay taxes through the special regime for groups of companies. 1 • Registered, voluntarily, in the Monthly VAT Return Registry.

In addition, this system is applied to those taxpayers who voluntarily adopt it. Thus, the SII comprises about 63,000 taxpayers representing around 70% of the country's total business turnover, with a great diversity of coverage by activities (see "Annex 1" for more details).

This statistic is particularly innovative at the international level both for its frequency (daily instead of monthly) and for its compilation procedure (exhaustive information based on administrative records instead of partial sampling surveys). Figure 1 shows the daily time series, from July 1st, 2017 to December 31th, 2020. The sample comprises from the full-fledged implementation of the SII until the end of 2020. At first glance, it is possible to see the large volatility of the series, due to the different invoicing patterns, that we will comment on later, such as the large effect of the end of the month, which makes very difficult to distinguish a relevant signal of underlying evolution.

As a general rule, firms must send their invoices four days after they have been issued. However, the information received by the Tax Agency is not always the final one, as it can be completed and corrected in the following days. Experience indicates that it is necessary to wait around two weeks for the stabilization of the levels that allows us to consider the data as definitive.

Moreover, it is appropriate to show how the daily information from the SII helps to anticipate the behavior of the domestic sales that are included in the monthly VAT data. Most of the companies that belonging to the SII must present the corresponding monthly periodic self-assessment forms. These companies can submit these self-assessment forms until the 25th of the month following the month to which they refer (except for the exceptions of December, which can be submitted until January 27, and January, which can be submitted until February 23).

It should be noted that there are minor conceptual differences between the daily data and the data from self-assessment forms (mainly due to operations with reversed taxpayer, and their corresponding modifications). In order to determine whether these differences are very significant, Fig. 2 shows the evolution of the year-on-year growth rates of both sources for total domestic sales in 2019 and 2020.

We can appreciate the great similarity between both data sources, directionally and in magnitude (specifically, the correlation coefficient between both series is 0.99). So, the usefulness of the daily series is clearly demonstrated in order to anticipate its monthly counterpart.

The series are still at an experimental stage, requiring the approval by the Statistical Council in order to be included in the National Statistical Plan, in the same way as happened with the remaining statistics compiled by the Tax Agency (e.g., Sales, 

In this section, we present the econometric methodology used in the paper. Modeling economic daily time series poses several and difficult challenges due to the coexistence of multiple seasonal components linked to various frequencies, the complex structure of the calendar, 2 the strength of its irregular component and its sensitivity to exogenous factors (e.g., outliers) that distort its usual behavior (Ladiray et al. 2018) .

We use the structural model of unobserved components proposed by De Livera et al. (2011) to perform modeling and seasonal adjustment of the Spanish daily sales data. The model, called TBATS (acronym for Trigonometric seasonality, Box-Cox transformation, ARMA innovations, Trend and Seasonality), complemented with a suitable dynamic-regression preprocessing, provides a flexible although parsimonious way to handle the complex nature of daily time series. The econometric methodology has two steps:

• Preprocessing (linearization). In this step, we apply an intervention analysis by means of exogenous deterministic variables designed to control for the presence of outliers and specific calendar effects that, due to their moving nature, do not fit well into the structural representation considered by TBATS ). This intervention analysis is a preprocessing step of the observed series that renders it suitable for TBATS. • Structural decomposition using TBATS. The preprocessed time series is decomposed into trend-cycle, seasonality and irregularity. As we will expose below, the seasonal component has a complex nature due the coexistence of multiple seasonal patterns, some of them with fractional periodicities.

Let us now explain both steps with some detail.

The rank-mean analysis of the time series, as well as its wide range of variation, suggests the convenience of using a preliminary logarithmic transformation. In the same vein, the estimation of the λ parameter of the Box-Cox transformation performed by TBATS is very close to zero, 3 providing additional support for the log-transformation. The analysis of the residuals both from the structural TBATS model and from a reduced-form model (Cuevas et al., 2019) confirms the need to control for some effects linked to: (i) bank holidays, (ii) the inexact periodicity of the monthly seasonal component and (iii) its interaction with the weekly seasonal component.

These effects are represented by deterministic variables, and their impact on the observed time series is estimated using a regression model that includes a linear trend and a multiple seasonal component affine to the one used in the trigonometric seasonal representation used by TBATS. The regression model is given as follows:

Being:

• z t : (log-transformed) observed variable. It is interesting to note that the regression model (1) can be considered as a oneequation approximation to the complete structural TBATS model that will be presented below, especially due its similar treatment of the (multiple) seasonality. This similitude enhances the complementarity of steps 1 and 2.

The number of harmonics, k i , associated with each seasonal component (weekly, monthly and annual) can be determined by means of a preliminary estimate of the TBATS model applied to the original time series.

The TBATS approach is based on the representation of the unobserved components (trend, seasonality, irregularity) by means of explicit dynamic models (Harvey 1989) . 5 Following the structural approach, the model incorporates a parsimonious but rather general representation of the trend. It also includes an explicit model for the irregular component that acts as a sort of "safety valve," accommodating elements that, for whatever reason, did not find a proper fit within the basic systematic components (trend and seasonality). In this way, the plain representation of these two components does not compromise neither the fit of the model to the sample nor its forecasting performance.

The TBATS model assumes that the (possibly Box-Cox transformed) observed series (z t ) results from the aggregation of three unobserved components: trend (p t ), seasonality (s t ) and a stationary innovation (u t ). This innovation plays an additional role as the stochastic input for the other two components. In this way, both the trend and the seasonality depend on a single shock that, properly scaled and filtered, generates them. In general, we assume that u t evolves according to a stationary and invertible autoregressive and moving average (ARMA) model:

The ultimate shock e t is a Gaussian white noise:

An interesting feature of TBATS is that it can handle complex seasonal patterns, comprising both multiple periodicities (weekly, monthly and yearly) and fractional periodicities (e.g., 30.4375 days for monthly seasonality or 365.25 for annual seasonality).

This complex seasonal pattern is one of the most important differences between a daily time series and its monthly/quarterly counterpart. In the latter case, the seasonal component is unique and of integer periodicity (12 months and 4 quarters, respectively). Of course, this additional complexity requires an additional layer of specific modeling.

Assuming that there are I seasonal components of different periodicity, total seasonality is the sum of all of them:

Each seasonal subcomponent is linked to a basic frequency and with k of its harmonics, according to the following equation:

Being:

• w i,j is the frequency of the j-th harmonic linked to the i-th seasonal subcomponent.

• m i is the periodicity, in time units, of the seasonal subcomponent (e.g., 7 days for the weekly seasonality).

In this way, the seasonality associated with each basic frequency is obtained by adding the signals associated with that basic frequency and its k harmonics:

These individual terms are determined according to a bivariate vector autoregressive (VAR) process that includes S and an auxiliary factor Q. 6

Equation (7) is a deterministic Fourier expansion centered on the frequency (5), stochastically perturbed by the common innovation of the system, u t . This innovation is scaled using the parameters γ 1 and γ 2 . This representation stands out for its parsimony, since only four parameters are involved: two scale parameters and two initial conditions, regardless of the time scale of the seasonality, which can be very large: greater than 28 and 360 periods in the monthly and annual case, respectively.

In this way, for each seasonal subcomponent (e.g., weekly), k first-order VAR representations are defined, as many as harmonics are needed to represent it.

The magnitude of the scale parameters determines the proximity to a deterministic behavior of the seasonal subcomponent. At the limit, if both are zero, the component is completely deterministic.

Finally, the trend p t is a random walk, I(1), with a first-order autoregressive, AR(1), drift:

Being:

• g t : drift.

• Φ: damping parameter that controls for the impact of the drift on the trend. In general, 0 ≤ Φ ≤ 1.

• α: scale parameter that modules the impact of the innovation on the trend. As a rule, α ≥ 0.

The next equation defines the drift:

Being:

• b: location parameter that represents the steady state of the drift, provided that Φ < 1. • β: scale parameter that modules the impact of the common innovation on the drift.

In general, β ≥ 0.

Equations (8) and (9) provide a parsimonious yet flexible representation for the trend. In this way, depending on Φ, we can get an I(2) or an I(1) trend. If Φ = 1, we obtain an IMA(2,1) trend. In addition, the scale parameters determine the closeness to a deterministic behavior. As a special case, if 0 < Φ < 1 and β = 0, we get a random walk with a constant drift. If, in addition, α = 0, the trend becomes completely deterministic.

The TBATS procedure sets the model in state space form, computes its likelihood and maximizes it using the model parameters as instruments. It also determines the most appropriate Box-Cox transformation and, once applied, the proper number of harmonics for the seasonal subcomponents, starting with j = 1. In all the cases, the different combinations are ranked according to the Akaike information criterion (AIC) and the one that minimizes AIC is chosen.

Finally, TBATS performs a search for the most adequate ARMA(p,q) model for the innovation, starting with a white noise (p = q = 0). If the innovation fails to be considered as a white noise, a search along p and q is implemented, selecting the combination that minimizes the AIC.

We turn now to the empirical results derived from the application of the methodology presented in the previous section, following its two-step scheme.

The first step consists of estimating Eq. (1). Let us now describe the exogenous variables considered in x t . The bank holidays variable is based on the official working calendar, including national and regional holidays. In this way, we get 20 daily time series (1 national, 17 regional and 2 for the autonomous cities). We have built a single regression variable by combining the 20 time series according to its weight on the distribution of interior sales as reported by the 56 offices of the Tax Agency (Cuevas et al. 2019) . 7 The role of this variable is improved if we restrict it to be binary, setting 2/3 as the threshold. 7 An alternative estimate, using weights from the Spanish Regional Accounts, yields similar results. The number of harmonics associated with each seasonal component (weekly, monthly and annual) is 3, 9 and 5, respectively. These values are derived from the preliminary estimate of a TBATS model applied to the original time series

The effects linked to the deterministic part of the monthly seasonality are collected using three binary variables that separately consider whether the day is the beginning of the month, the 15th day, or the end of the month, adopting the value 1 in this case and 0 in the rest. These three effects interact with the weekly seasonality, for which three additional binary variables are considered that adopt the value of 1 if, in addition to being the beginning of the month, the 15th day or the end of the month, the day is also a weekend (Saturday or Sunday). Table 1 shows the results of the estimation of the deterministic effects, by means of the regression model (1).

It is worth noting the strength of the holiday effect, as well as the impact of the end of the month effect, in itself and in its interaction with the weekly seasonality. The beginning of the month effect and day-15th effect, although significant, are of a lesser magnitude.

The daily time series, corrected from the deterministic effects quantified in the previous table, is decomposed by means of the TBATS model. The next table presents the estimation of its parameters.

The estimated parameters indicate that all the seasonal components are quite stable, especially the one related to the weekly seasonality. This stability does not prevent the estimated components from being relatively complex, given the high number of harmonics required to represent them.

The estimated trend is a random walk whose specific innovation is much less volatile than the corresponding to the common innovation u t .

From the previous table, we can ascertain that the common innovations u t that affect the series are mainly reflected in its trend and its annual seasonality. Monthly seasonality and, above all, weekly seasonality are relatively immune to the shocks. This ranking is an alternative way to quantify proximity to a deterministic scheme and is shown in Fig. 3 .

The TBATS model, whose parameters are shown in Table 2 , allows the estimation of the unobserved components that underlie the observed series by means of the Kalman filter. These components are presented in Fig. 4 . The information provided by this decomposition makes it possible to infer, in the first place, the weekly seasonal profile, the most stable component of all. As can be seen in the following graph, the weekly pattern of daily sales shows a strong contrast between the weekend, especially Sunday, and the rest of the week (Fig. 5) .

In the same vein, we can check the mean profile of the annual seasonality. As can be seen in the next figure, its pattern is relatively complex, comprising several peaks and troughs. Among the former, the days that belong to December, June and July stand out. On the other hand, the trough days are registered in August, January and February (Fig. 6) .

To fully analyze the monthly seasonal profile, it is necessary to combine the stochastic component estimated by TBATS, as shown in Fig. 4 , with the estimation of the deterministic effects of monthly periodicity, (beginning of the month, day 15th and end of the month effects), presented in Table 1 . These effects, properly centered to average zero in the long term, are shown in the following graph (Fig. 7) .

The main features of the monthly seasonality are, first and above all, the preeminence of its deterministic component, the important role of the end of the month effect and, finally, the lesser relevance of both the beginning of the month and the 15th day effects.

Stochastic decomposition of the daily SII sales. Note level refers to trend, season1 refers to weekly seasonality, season2 refers to monthly seasonality and season3 refers to annual seasonality One additional feature of the TBATS model is that can be used to forecast the daily sales in order to have permanently filled the current month. Of course, the forecasts can be extended as much as required if, for whatever reasons, more distant forecasts are required. In Sect. 5, we present additional information about the predictive performance of the model.

As mentioned before, the firms included in the Immediate Information Supply system (SII) are classified according to the activity declared by their own companies. Specifically, the activities are classified according to the breakdown of the Economic Activity Tax (IAE, Impuesto de Actividad Económica), which has a correspondence with the standard NACE-2009 classification (see "Annex 1" for details on its coverage).

In this way, the methodological approach previously described can be applied to all the branches of activity considered. For simplicity, we are going to show the most significant results for total sales and the main branches of activity. Additionally, we are going to focus in the recent recessive episode due to the COVID-19 health crisis, because the extraordinary and unusual nature of this shock deserves special attention, showing how it has affected the different economic activities. Figures 8 and 9 show, respectively, the final series corrected from deterministic effects and seasonal and calendar components, and their corresponding year-on-year rate growth. We have used a 28 days moving average filter to render the levels comparable with those of a monthly time series.

Transport Hostelry The adjustment carried out allows us to extract a clear signal, useful for the monitoring of sales in each sector on a real-time basis. Likewise, the decline in activity caused by the impact of the COVID-19 shock has been clearly uneven by activity branch. The greatest impact caused by the lockdown and related policies have had the largest impact on the hostelry and transport sectors, with falls that remain at 60% in the first case.

Daily forecasting is an important use of the TBATS model, in addition to seasonal adjustment. In this way, for example, we can complete the current month having an updated forecast that can be plugged in nowcasting models. Daily forecasts can also be used to assess and monitor the short-term situation. The exercise proposed in this section seeks to evaluate the forecasting accuracy of this model against other alternatives that generate daily projections.

The training set to be considered will be the SII series up to 2019/12/31, while 2020 would be used as the test set. It is necessary to emphasize that information is only available since mid-2017, so in order to have one-year forecast horizon and at least 20% of the sample for the test set (Hyndman and Athanasopoulos 2021), we cannot evaluate previous years. Therefore, TBATS forecasting accuracy is checked during an atypical period due to the effect of the COVID-19 shock. Figure 10 shows 2020 the first quarter out-of-sample forecast from 2019/12/31 versus real data. As can be seen, TBATS captures well the seasonality and returns a In order to further refine the evaluation of the predictive performance, the forecasting accuracy evaluation is done through two exercises based on rolling forecast schemes with different time frames. Firstly, a rolling forecast exercise day by day until the month is completed, on the other hand, a daily forecast scheme month by month until the end of the quarter, which will be extended secondly. As alternative models, versus the TBATS model as benchmark, the following methods are proposed:

1. RW: random walk. 2. STLM: STL (Seasonal and Trend decomposition using LOESS 8 ) method. 3. ARIMA: ARIMA model obtained from AIC minimization. 9 4. MRL: linear time series with trend and multiple seasonality.

These four methods cover various projection approaches: simple (RW), nonparametric (STLM), parametric in reduced form (ARIMA), and parametric based on an essentially deterministic regression model (MRL). Moreover, all the models have been constructed with the series in logarithms and corrected from the COVID-19 shock with dummy variables.

Finally, as a measure of predictive precision, the comparison between models is based on the root of the mean square error (RMSE). Besides, the second exercise completes the analysis using the Diebold and Mariano test for the null hypothesis that SERIEs the two methods have the same forecast accuracy using the modified version proposed by Harvey et al. (1997) . Table 3 and Fig. 11 show the RMSE of the daily forecast until the end of January 2020. TBATS clearly improves the results of the RW and the ARIMA model. The MRL improves it in specific moments, while the STLM presents a lower RMSE only Table 4 and plotted in Fig. 12 . As in the daily case, the three models with the best accuracy are TBATS together with the MRL and the STLM. However, MRL improves TBATS in the first quarter of the year, coinciding with the beginning of the pandemic. In the second quarter, TBATS beats the rest of the models in June. Finally, in the remaining months from the third quarter, TBATS registers a lower RMSE than the proposed alternatives.

Complementing the analysis, the results of the Diebold and Mariano test month by month through the p values are represented in Table 5 . TBATS is the benchmark model against the proposed alternatives. It can be seen how the RW and the ARIMA models show different accuracy from the benchmark in all periods of the test sample. In the MRL case, there is a major rejection of the null hypothesis in the second semester once the model assimilates the impact of the pandemic. In the case of the STLM, the analysis of the results month by month does not provide a clear picture.

Summarizing, TBATS model provides a lower RMSE than the proposed alternatives in most of the prediction intervals considered. Besides, the Diebold and Mariano test rejects in most of the cases the null hypothesis of same forecast accuracy. The deviations from the forecast on real data occur from the COVID-19 outbreak. However, the reliability of the method is comparable to other models in this atypical period. Harvey et al.1997, version) for the null hypothesis of the equal forecast accuracy of two forecast methods. The number represents the p value *Rejection of the null hypothesis at the 10% level of significance **Rejection of the null hypothesis at the 5% level of significance ***Rejection of the null hypothesis at the 1% level of significance

The use of structural, unobserved components time series model has been successful both to extract non-seasonal information and to forecast the daily data. From a purely statistical view, the employ of tax data stands out as an inexpensive and reliable way to use the information routinely provided by the economic agents when complying with their tax duties, replacing the use of costly sample surveys or qualitative information derived from data sources whose representativeness cannot be properly assessed. Of course, part of this reliability comes from the positive reinforcement due to penalties and fines if compliance is not perfect. Few data sources are backed by this reinforcement schedule.

For future work, we plan to compare TBATS with both STLM (Ollech 2018 ) and with Prophet (Taylor and Letham 2018) in the realm of seasonal adjustment, in order to ascertain the performance of model-based methods (TBATS) versus nonparametric filtering methods (STLM) and Bayesian regression methods (Prophet).

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long SERIEs as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.

The branches of activity that can be considered, with their corresponding coverage in 2019, comprise the following sections of the NACE: 

STL: A seasonal-trend decomposition procedure based on Loess

Incorporando información fiscal de frecuencia diaria en la previsión macroeconómica a corto plazo

Forecasting time series with complex seasonal patterns using exponential smoothing

Automatic modelling of daily series of economic activity

Forecasting, structural time series models and the Kalman filter

Testing the equality of prediction mean squared errors

Modeling time series with calendar variation

Modeling considerations in the seasonal adjustment of economic time series

Seasonal adjustment of daily and weekly data

Time series analysis and forecasting. Scientific Computing Associates, River Forest Mariano R, Murasawa Y (2003) A new coincident index of business cycles based on monthly and quarterly series

Seasonal adjustment of daily time series

Forecasting at scale

Suministro Inmediato de Información (SII) del Impuesto sobre el Valor Añadido (IVA)

Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations