Discrete time Infinite Impulse Response low-pass filters are widely used in many fields such as engineering, physics and economics. Once applied to a given time series, they have the ability to pass low frequencies and attenuate high frequencies. As a result, the data are expected to be less noisy. A properly filtered signal, is generally more informative with positive repercussions involving qualitative aspects – e.g. visual inspection and interpretation – as well as quantitative ones, such as its digital processing and mathematical modelling. In order to effectively disentangle signal and noise, the filter smoothing constant, which controls the degree of smoothness in First Order Discrete Time Infinite Impulse Response Filters, has to be carefully selected. The proposed method conditions the estimation of the smoothing parameter to a modified version of the information criterion of the type Hannan - Quinn which in turns is built using the Estimated Log Likelihood Function of a model of the class SARIMA (Seasonal Auto Regressive Moving Average). Theoretical evidences as well as an empirical study conducted on a particularly noisy time series will be presented.

Keywords: big data, discrete time infinite impulse response filters, denoising, hannan quinn information criterion, sarima models, time series

Among the many denoising methods and techniques successfully employed for univariate time series – e.g. based on regression,¹ Kalman filter,^2,3 decomposition,⁴ wavelet^5,6 and non-linear method⁷– those based on algorithms of the type Infinite Impulse Response (IIR) exponential filters have been massively used, given their satisfactory performances (see, for example,⁸ and, more recently⁹). Such methods are useful for their ability to maximize the amount of relevant information that can be extracted from “real life” time series. In fact, regardless the scientific field time dependent data are collected for (e.g. engineering, economics, physics, environmental), they can never be error–free. In spite of all of the efforts and precautions one might take in order to provide clean data – e.g. robust data acquisition methods, reliable routine checks, sophisticated procedures for error correction, fail safe data storage and data communication lines – reality is way too complex for such procedures to be completely reliable.

Noise, in fact, is simply an ubiquitous entity able to affect virtually all the stages a given analysis of a time series can be broken into, showing uncountable expressions that can be only partially controlled, never fully removed nor exactly pinpointed. Many are the fields where noise reduction methods are employed: data mining, satellite data, radio communications, radar, sonar and automatic speech recognition, are just a few of them. Often, the treatment of noisy data is critical, as in the case of bio–medical signals, tracking systems for defense purposes or economics, where the trend or the seasonal patterns of vital economic variables can be obscured or distorted by noise components. However, the type of time series likely to be easily affected by noise – as well as by sudden or unpredictable high frequency variations – can be found in the economic and social fields. A telling example, is the great deal of data generated by web–based services (e.g. Google or Twitter), which are commonly and publicly available, in many instances free of charge. However, statistical estimation procedures based on this type of data can embody many sources of uncertainty and instability. They are related, for example, to technical (computer alphabetization) or psychological (e.g. emotional driven spikes and other form of data perturbations) reasons. Even when data are collected and validated by Official Organisms (e.g. national and super-national statistical offices or central banks), strong psychological components can play a significant role in determining erratic and/or noisy behaviors. This is the case, for instance, of the Economic Sentiment Indicators – provided by many National Institute of Statistics – which are purposely designed to capture the amount of optimism (pessimism) towards the future behavior of a set of economic variables. Being these data able to reflect, at least to some extent, the future decision making strategies of a population of reference, they can show dangerous instabilities and irregularities. In other words, when an opinion is in the process to form or change then many analysts are mostly interested in its future developments; however, this is usually the points in time where the data show complex structures, which are usually hard to capture.

The need to conduct more precise model–based investigations under noisy data structures has been the motivating factor of the proposed method. As it will be seen, it delivers interesting performances and is reasonably fast, the latter characteristic making the proposed method a viable option with high dimensional data sets, e.g. of the type big data. It is also risk–free, i.e. when there are no filtering–related benefit – according to a suitable objective function – the method is designed to simply leave the data unfiltered (the $\lambda$ –smoothing parameter is set to $1$ ). The paper proceeds as follows: in Section 2 the problem is defined along with the statistical model assumed for the underlying data generating mechanism, whereas in Section 3 the employed low–pass filter is illustrated, with a special focus on its calibration issues. The method is explained in Section 4 and the step-by-step implementation of the algorithm is given in paragraph 4.3. Finally, an empirical experiment, conducted on a time series provided by the search engine Google and recovered in the repository called Google Trends (for Italy the address is https://trends.google.it/trends/?geo=IT), is presented in section 5.

Throughout the paper, the time series (also defined as signal) of interest is intended to be a real–valued, uniformly sampled, sequence of data points of length $T$ , denoted as:

$x_t\mathrm{<mo>:</mo><mo>=</mo>}\left\{{\mathrm{<mo\; stretchy="false">(</mo>}{x}_t\mathrm{<mo\; stretchy="false">)</mo>}}_{t\in {\mathbb{Z}}^{+}}^T\right\}\mathrm{.}$     (1)

Each datum $x_t$ , being observed with an error, say ${\epsilon }_t$ , is called $x_t^o$  in contrast with the “pure’, noise–free, unobservable one $x_t^u$ . Therefore, the formal set up is of the type

$x_t^o\mathrm{<mo>=</mo>}{x}_t^u+{\epsilon }_t\mathrm{.}$

It is worth stressing that the theoretical framework we are dealing with, does not envision ${\epsilon }_t$  as being the standard $0$ –mean constant variance ${\sigma }^2$  idiosincratic component. In fact, ${\epsilon }_t$  can also embody all the unpredictable components arising when “real–life” signals come into play and which, in general, are able to generate unwanted phenomena, like sudden outburst of energy (outliers) and noise. As a consequence, at least two detrimental effects are embedded in the signal and can affect its statistical properties: accuracy reduction and generation of heavy noisy components. The former relates to the addiction (subtraction) of significant quantities to (from) $x_u$ , so that we have $x^o\mathrm{<mo>=</mo>}k{x}^u\mathrm{<mo>;</mo>}k\ne 1$ , while the latter is a subtler, but nonetheless dangerous, disturbing phenomenon. It arises, for example, when in a survey some statistical units are not longer available and must be replaced by “similar” ones. When such a corrupted information set is applied as an input to a mathematical model for any purposes, e.g. simulation or prediction, the related outcomes are likely to be affected.

The data generating process

The observed signal is supposed to be a realization of a statistical model of the class $SARIMA$  (Seasonal Auto Regressive Integrated Moving Average), which is a generalization of the $SARIMA$  (Auto Regressive Integrated Moving Average) class proposed by George et.,¹⁰ have been introduced to model complex dynamics i.e. of the type stochastic seasonal. It can be expressed as:^10,11

${\varphi }_p\mathrm{<mo\; stretchy="false">(</mo>}B\mathrm{<mo\; stretchy="false">)</mo>}{\Phi }_P\mathrm{<mo\; stretchy="false">(</mo>}{B}^S\mathrm{<mo\; stretchy="false">)</mo>}{\nabla }_S^D{\nabla }^d{X}_t\mathrm{<mo>=</mo>}{\theta }_q\mathrm{<mo\; stretchy="false">(</mo>}B\mathrm{<mo\; stretchy="false">)</mo>}{\Theta }_Q\mathrm{<mo\; stretchy="false">(</mo>}{B}^S\mathrm{<mo\; stretchy="false">)</mo>}{\epsilon }_t$

where $B$  denotes the backward shift operator, $d$  and $D$  denote the non – seasonal and seasonal difference operator respectively and:

${\nabla }^d\mathrm{<mo>=</mo>1}-B^d$ ;

${\varphi }_p\mathrm{<mo\; stretchy="false">(</mo>}B\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>1}-{\varphi }_{1}B-{\varphi }_2{B}^2-\mathrm{....}-{\varphi }_p{B}^p$ ;

${\theta }_q\mathrm{<mo\; stretchy="false">(</mo>}B\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>1}-{\theta }_{1}B-{\theta }_2{B}^2-\mathrm{....}-{\theta }_q{B}^p$ ;

${\Phi }_P\mathrm{<mo\; stretchy="false">(</mo>}{B}^S\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>1}-{\Phi }_1{B}^S-{\Phi }_2{B}^{2S}-\mathrm{....}-{\Phi }_P{B}^{PS}$ ;

${\Theta }_Q\mathrm{<mo\; stretchy="false">(</mo>}{B}^S\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>1}-{\Theta }_1{B}^S-{\Theta }_2{B}^{2S}-\mathrm{....}-{\Theta }_q{B}^{QPS}$ ;

with $\varphi$ , $\theta$ , $\Phi$ , $\Theta$ , respectively, the non seasonal autoregressive and moving average parameters and the seasonal autoregressive and moving average parameters. Finally ${\epsilon }_t$  is a $0$ -mean and finite variance white noise. The model can be estimate when the stationary and invertibility conditions are meet for both the autoregressive and moving average polynomials respectively, that is ${\varphi }_P\mathrm{<mo\; stretchy="false">(</mo>}B\mathrm{<mo\; stretchy="false">)</mo>}{\Phi }_p\mathrm{<mo\; stretchy="false">(</mo>}B\mathrm{<mo\; stretchy="false">)</mo>}$ =0 and ${\Theta }_Q\mathrm{<mo\; stretchy="false">(</mo>}B\mathrm{<mo\; stretchy="false">)</mo>}{\theta }_q\mathrm{<mo\; stretchy="false">(</mo>}B\mathrm{<mo\; stretchy="false">)</mo>}$ =0 have their root lying outside the unit circle. Generally, this model is abbreviated as $SARIMA\mathrm{<mo\; stretchy="false">(</mo>}p\mathrm{,}d\mathrm{,}q\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}P\mathrm{,}D\mathrm{,}Q{\mathrm{<mo\; stretchy="false">)</mo>}}^S$ . When the process is stationary and no seasonal patterns are detected, the model collapses to a pure $ARMA\mathrm{<mo\; stretchy="false">(</mo>}p\mathrm{,}q\mathrm{<mo\; stretchy="false">)</mo>}$  model. The explanation of the ARIMA theory is outside the scope of this paper, therefore the interested reader is referred to Makridakis⁸ and Stock.¹²

Under Gaussianity, the SARIMA parameter vector

$\text{G}\equiv \mathrm{<mo\; stretchy="false">(</mo>}{\varphi }_p\mathrm{,}{\theta }_q\mathrm{,}{\Phi }_P\mathrm{,}{\Theta }_Q\mathrm{,}{\sigma }^2\mathrm{<mo\; stretchy="false">)</mo>}$         (2)

can be estimated via Conditional $MLE$  (Maximum Likelihood Estimation) method, which:

1. In the version employed to run the empirical experiment (Section 5), uses the joint probability density function for ${\epsilon }_1\mathrm{,}\dots \mathrm{,}{\epsilon }_T$ , i.e.
2. $f\mathrm{<mo\; stretchy="false">(</mo>}{\epsilon }_1\mathrm{,}\dots \mathrm{,}{\epsilon }_T\mathrm{<mo>|</mo>}\text{G}\mathrm{<mo>=</mo><mo\; stretchy="false">(</mo>2}\pi {\sigma }^2{\mathrm{<mo\; stretchy="false">)</mo>}}^{-T\mathrm{<mo>/</mo>2}}\exp \mathrm{<mo>\{</mo>}-\frac{1}{2{\sigma }^2}{\displaystyle \sum _{t\mathrm{<mo>=</mo>1}}^T}{\epsilon }^2\mathrm{<mo>\}</mo><mo>;</mo>}$

and

3. Maximizes the function given by

$\log f\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo>|</mo>}{\epsilon }_0\mathrm{<mo>=</mo>0<mo>;</mo>}\text{G}\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}-\frac{T}{2}\log 2\pi {\sigma }^2-\frac{1}{2{\sigma }^2}{\displaystyle \sum _{t\mathrm{<mo>=</mo>1}}^T}{\epsilon }_t^2\mathrm{,}$

with $\mathrm{<mo\; stretchy="false">[</mo>}{\epsilon }_t\mathrm{<mo>;</mo>}t\mathrm{<mo>=</mo>1,}\dots \mathrm{,}T\mathrm{<mo\; stretchy="false">]</mo>}$  being recursively estimated using

${\epsilon }_t\mathrm{<mo>=</mo>}{x}_t-\text{G}{\epsilon }_{t-1}\mathrm{.}$        (3)

Equation (3) holds by equating the first $\mathrm{<mo\; stretchy="false">(</mo>}P+Q\mathrm{<mo\; stretchy="false">)</mo>}$  observations to 0, i.e. $y_0\mathrm{,}{y}_{-1}\mathrm{,}\dots \mathrm{,}{y}_{-\mathrm{<mo\; stretchy="false">(</mo>}P+Q\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo>=</mo>0}$  and setting the first $Q-1$  innovations to 0: ${\epsilon }_0\mathrm{,}{\epsilon }_{-1}\dots \mathrm{,}{\epsilon }_{-\mathrm{<mo\; stretchy="false">(</mo>}Q-\mathrm{1<mo\; stretchy="false">)</mo>}}\mathrm{<mo>=</mo>0}$ .

By the Central Limit Theorem, by “properly” averaging sub-sequences of a noisy signal, a better approximation of the “true” ones is obtained as a result.¹³ Consistently, one strategy for noise control relies on a filter of the type discrete time Infinite Impulse Reponse.¹⁴ This class of filters is a powerful tool for handling noise in many fields: in¹⁵ its performances has been compared with many other methods on environmental data (other than simulated ones) whereas its application for standard control-chart procedures is documented in Layth et al.,¹⁶ and Stuart Hunter.¹⁷ Its usefulness for interactive systems and information extraction for complex turbulent flows has been discussed respectively in Casiez¹⁸ and Adrien Cahuzac,¹⁹ whereas the effectiveness of smoothing–driven approaches in economic time series forecasting is documented in the excellent book from Hyndman et al.,²⁰ and in Codrut.²¹

IIR filters envision the level of the stochastic process (1) at a given time $t$ , say ${\xi }_t$  to be dependent to the present and past observations according to

${\xi }_t\mathrm{<mo>=</mo>}\lambda x_t+\lambda \mathrm{<mo\; stretchy="false">(</mo>1}-\lambda \mathrm{<mo\; stretchy="false">)</mo>}{x}_{t-1}+\lambda {\mathrm{<mo\; stretchy="false">(</mo>1}-\lambda \mathrm{<mo\; stretchy="false">)</mo>}}^2{x}_{t-2}+\dots \mathrm{<mo>=</mo>}\lambda {\displaystyle \sum _{j\mathrm{<mo>=</mo>0}}^{\infty }}{\mathrm{<mo\; stretchy="false">(</mo>1}-\lambda \mathrm{<mo\; stretchy="false">)</mo>}}^j{x}_{t-j}\mathrm{.}$

Using the lag operator $L$ , such that $L^k{x}_t\mathrm{<mo>=</mo>}{x}_{t-k}$ , (4) can be re-expressed as

${\xi }_t\mathrm{<mo>=</mo>}\lambda \mathrm{<mo\; stretchy="false">[</mo>1}+\mathrm{<mo\; stretchy="false">(</mo>1}-\lambda \mathrm{<mo\; stretchy="false">)</mo>}B+{\mathrm{<mo\; stretchy="false">(</mo>1}-\lambda \mathrm{<mo\; stretchy="false">)</mo>}}^2{B}^2+\dots \mathrm{<mo\; stretchy="false">)</mo><mo>=</mo><mo\; stretchy="false">[</mo>1}-\mathrm{<mo\; stretchy="false">(</mo>1}-\lambda \mathrm{<mo\; stretchy="false">)</mo>}B{\mathrm{<mo\; stretchy="false">]</mo>}}^{-1}$         (4)

and rewritten as

${\xi }_t\mathrm{<mo>=</mo><mo\; stretchy="false">[</mo>1}-\mathrm{<mo\; stretchy="false">(</mo>1}-\lambda \mathrm{<mo\; stretchy="false">)</mo>}B\mathrm{<mo\; stretchy="false">]</mo>}{\xi }_t\mathrm{<mo>=</mo>}\lambda x_t\mathrm{.}$     (5)

From (5) it is evident that ${\xi }_t$  is a weighted average of $x_t$  (current datum) and the previous level ${\xi }_{t-1}$ , whose weight sequence is determined by the smoothing parameter $\lambda$ . Finally, (4) can be expressed, by skipping the lag operator, in a recurrence form as

${\xi }_t\mathrm{<mo>=</mo>}\lambda x_t+\mathrm{<mo\; stretchy="false">(</mo>1}-\lambda \mathrm{<mo\; stretchy="false">)</mo>}{\xi }_{t-1}\mathrm{<mo>;</mo>}{\xi }_0\mathrm{<mo>=</mo>}{x}_0\mathrm{,}$         (6)

or in a error correction form, i.e. determined by the equation $e_t\mathrm{<mo>=</mo><mo\; stretchy="false">(</mo>}{x}_t-{\lambda }_{t-1}\mathrm{<mo\; stretchy="false">)</mo>}$ , as follows:

${\xi }_t\mathrm{<mo>=</mo>}{\xi }_{t-1}+\lambda e_t\mathrm{.}$

The form (6) is the one used in the sequel. Here, the first term represents the contribution added by each and every datum updating our time series whereas the second one accounts for the inertia from the previous observations. Therefore, the amount of noise left in the signal as well as the system lag (the responsiveness of the system output to changes in the input) are both proportional to $\lambda$ .

It is worth emphasizing that the order of a low–pass filter reflects the amount of strength it uses to attenuate a predetermined frequency band. First order filters of the type (4) reduce the signal amplitude by half for each frequency doubling while higher order filters amplify this effect proportionally with their order. For the purpose pursued here, the order 1 appeared to be appropriate since the filter is:
[(a)]

1. “Only” required to denoise (e.g. no forecasting purposes are pursued);
2. Easier to tune;
3. More interpretable.

However, properly tuning the constant parameter is not a trivial task nor something that – at least in general – can be left to subjective judgment. In fact, a too high $\lambda$  can leave significant amount of idiosincratic components in the signal while too low values determine the speed at which the older data are dampened to be too slow (over smoothing).

The smoothing constant and the objective function

As already pointed out, the smoothing parameter (also called smoothing factor or smoothing constant), controls the amount of past information contributing to the formation of the present signal level. Formally, a point in time distant $\mathrm{<mo\; stretchy="false">(</mo>}t+\mathrm{1<mo\; stretchy="false">)</mo>}-k\mathrm{<mo>=</mo>}m$  lags influences the level of the signal at time $t+1$  by an amount given by ${\xi }_{t+1}\mathrm{<mo>=</mo>}\lambda {\mathrm{<mo\; stretchy="false">(</mo>1}-\lambda \mathrm{<mo\; stretchy="false">)</mo>}}^{t+1-m}$ . Even though the proposed method is designed to provide an estimate of $\lambda$ , it might be interesting to investigate the frequency property of the filter, in order to obtain meaningful information from the selected $\lambda$ . The related cut–off frequency can be derived^15,22 by $z$ –transorming (6), i.e.

$\widehat{\xi }\mathrm{<mo\; stretchy="false">(</mo>}z\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\frac{\lambda }{1-\mathrm{<mo\; stretchy="false">(</mo>1}-\lambda \mathrm{<mo\; stretchy="false">)</mo>}{z}^{-1}}\widehat{x}\mathrm{<mo\; stretchy="false">(</mo>}z\mathrm{<mo\; stretchy="false">)</mo>,}$

which can be regarded as the digital version of the analog filter (6). By defining ${\Delta }_t$  as the sampling period and equating $z$  to the frequency $\omega$ , i.e. $z\mathrm{<mo>=</mo>}{e}^{\mathrm{<mo\; stretchy="false">(</mo>}i\omega {\Delta }_t\mathrm{<mo\; stretchy="false">)</mo>}}$ ; $i\mathrm{<mo>=</mo>}\sqrt{\mathrm{<mo\; stretchy="false">(</mo>}}-\mathrm{1<mo\; stretchy="false">)</mo>}$ , the power spectra of ${\xi }_t$  and $x_t$ , respectively called $\mathrm{<mo>|</mo>}\widehat{\xi }\mathrm{<mo\; stretchy="false">(</mo>}\omega {\mathrm{<mo\; stretchy="false">)</mo><mo>|</mo>}}^2$  and $\mathrm{<mo>|</mo>}\widehat{x}\mathrm{<mo\; stretchy="false">(</mo>}\omega {\mathrm{<mo\; stretchy="false">)</mo><mo>|</mo>}}^2$  are related to each other through

$\mathrm{<mo>|</mo>}\widehat{\xi }\mathrm{<mo\; stretchy="false">(</mo>}\omega {\mathrm{<mo\; stretchy="false">)</mo><mo>|</mo>}}^2\mathrm{<mo>=</mo>}\frac{{\lambda }^2}{\mathrm{<mo>|</mo>}{e}^{\mathrm{<mo\; stretchy="false">(</mo>}i\omega {\Delta }_t\mathrm{<mo\; stretchy="false">)</mo>}}-\mathrm{<mo\; stretchy="false">(</mo>1}-\lambda {\mathrm{<mo\; stretchy="false">)</mo><mo>|</mo>}}^2}\mathrm{<mo>|</mo>}\widehat{x}\mathrm{<mo\; stretchy="false">(</mo>}\omega {\mathrm{<mo\; stretchy="false">)</mo><mo>|</mo>}}^2\mathrm{.}$    (7)

By equating the first term in (7) to $\frac{1}{2}$ , the cut–off frequency at which the amplitude is reduced by a factor of 2, say $\tilde{F}$ , is found, i.e.

$\lambda \mathrm{<mo>=</mo>}\frac{1}{3}\sqrt{\mathrm{6<mo\; stretchy="false">(</mo>1}-\cos \mathrm{<mo\; stretchy="false">(</mo>}\tilde{w}{\Delta }_t\mathrm{<mo\; stretchy="false">)</mo>}+\mathrm{<mo\; stretchy="false">(</mo>1}-\cos {\mathrm{<mo\; stretchy="false">(</mo>}\tilde{w}{\Delta }_t\mathrm{<mo\; stretchy="false">)</mo>}}^2-\mathrm{<mo\; stretchy="false">(</mo>1}-\cos \mathrm{<mo\; stretchy="false">(</mo>}\tilde{w}{\Delta }_t\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}}\approx \frac{2\pi \tilde{F}{\Delta }_t}{\sqrt{3}}\approx 3.62579\tilde{F}{\Delta }_t\mathrm{.}$

and therefore

$\lambda \mathrm{<mo>=</mo>}\frac{1}{3}\sqrt{\mathrm{6<mo\; stretchy="false">(</mo>1}-\cos \mathrm{<mo\; stretchy="false">(</mo>}\tilde{w}{\Delta }_t\mathrm{<mo\; stretchy="false">)</mo>}+\mathrm{<mo\; stretchy="false">(</mo>1}-\cos {\mathrm{<mo\; stretchy="false">(</mo>}\tilde{w}{\Delta }_t\mathrm{<mo\; stretchy="false">)</mo>}}^2-\mathrm{<mo\; stretchy="false">(</mo>1}-\cos \mathrm{<mo\; stretchy="false">(</mo>}\tilde{w}{\Delta }_t\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}}\approx \frac{2\pi \tilde{F}{\Delta }_t}{\sqrt{3}}\approx 3.62579\tilde{F}{\Delta }_t\mathrm{.}$

The above approximation is possible being in practice $\tilde{w}{\Delta }_t\mathrm{<mo><</mo><mo><</mo>1}$ , for the temporal integration to be possible.

The Hannan–Quinn information criterion

As it will be outlined in the sequel, the smoothing constant is the minimiser of a penalized Log-Likelihood function, estimated on a model of the class SARIMA applied to both the original and filtered time series ${\widehat{\xi }}_t$ . The objective function chosen belongs to the class of the information criterion, i.e. the Hannan-Quinncriterion.²³ This order selector has been constructed from the law of iterated algorithm and shows a penalty function growing at a very slow rate, as the samples size increases. It is defined as follows:

$HQc\mathrm{<mo>=</mo>}\log L\mathrm{<mo\; stretchy="false">(</mo>}\text{<mover accent="true"> G <mo>^</mo> </mover>}\mathrm{<mo>|</mo>}{X}_t\mathrm{<mo\; stretchy="false">)</mo>}+2k\log \mathrm{<mo\; stretchy="false">(</mo>}\log \mathrm{<mo\; stretchy="false">(</mo>}T\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>,}$

being $\text{<mover accent="true"> G <mo>^</mo> </mover>}\equiv \mathrm{<mo\; stretchy="false">(</mo>}\text{f}\mathrm{,}\text{F}\mathrm{,}\text{q}\mathrm{,}\text{Q}\mathrm{,}\sigma \mathrm{<mo\; stretchy="false">)</mo>}$  the vector of the estimated SARIMA parameters of (??) whereas $L$  is the Likelihood function.

By estimating the log–likelihood function $\log L\mathrm{<mo\; stretchy="false">(</mo>}\text{<mover accent="true"> G <mo>^</mo> </mover>}\mathrm{<mo>|</mo>}{X}_t\mathrm{<mo\; stretchy="false">)</mo>}$  with $\mathrm{<mo\; stretchy="false">(</mo>}{\widehat{\sigma }}^2\mathrm{<mo>|</mo>}\text{<mover accent="true"> G <mo>^</mo> </mover>}\mathrm{<mo\; stretchy="false">)</mo>}$  and writing the penalty term as $\rho \mathrm{<mo>=</mo>2}$ k $\log \mathrm{<mo\; stretchy="false">(</mo>}\log \mathrm{<mo\; stretchy="false">(</mo>}T\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}$ , the Hannan Quinn information criterion can be written as

$HQc\mathrm{<mo>=</mo>}f\mathrm{<mo\; stretchy="false">(</mo>}{\widehat{\sigma }}_o^2\mathrm{,}\rho \mathrm{<mo\; stretchy="false">(</mo>}\text{G}\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo><mo>|</mo>}{X}_t^o\mathrm{<mo\; stretchy="false">)</mo>,}$

where also the subscript $\mathrm{<mo\; stretchy="false">(</mo>}{}_o\mathrm{<mo\; stretchy="false">)</mo>}$ , stands for “observed”.

The standard identification procedure for the selection of the best SARIMA model order under this information criterion, follows the general rule as the other Information Criteria and is denoted as Minimum Information Criterion Estimation ( $MICE$ ). In essence, it is based on the minimization of the information criterion itself. Here, the model $M_0$  minimizing the HQ criterion, is the winning one, that is:

$M_0^o\mathrm{<mo>:</mo><mo\; stretchy="false">(</mo>}\text{<mover accent="true"> G <mo>^</mo> </mover>}\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\arg \underset{\text{<mover accent="true"> G <mo>^</mo> </mover>}}{{\displaystyle \min }}f\mathrm{<mo\; stretchy="false">(</mo>}{\widehat{\sigma }}^2\mathrm{,}\text{G}\mathrm{<mo>|</mo>}{X}_t^o\mathrm{<mo\; stretchy="false">)</mo>,}$            (8)

with $\Gamma$  as defined in (2).

The proposed method can be regarded as an optimal averaging procedure under a penalized log–likelihood function. “Properly” averaging a signal as a noise reduction method is consistent to the principle that, if repeated measures are conducted many times on a given signal, part will tend to accumulate but the noise will be irregular and tend to cancel itself. This result is connected to the fact that the mean standard deviation of N measurements is smaller by a factor of N than the standard deviation of a single measurement. Practically, if one compute the average of “many” samples of a noisy signal the noise–induced random fluctuations tend to die out and the signal is therefore stronger.

In more details, when one records the signal $x_t^o$  (1) $t$  times (the sampling rate) over a predefined time interval, say $\tau$ , being $\tau \le T$ , he can compute $l$  averages each of length $\frac{l}{\tau }$ . Two critical points arise at this time, i.e. how many data points should be averaged and in which manner. The latter question has been already answered (Section 3) whereas the former – which is strictly related with the smoothing parameter $\lambda$  – controls the number of observations of $x_t^o$  to be averaged (and to what extent), so that the noise which randomly fluctuate with equal probability above or below the “true” level $x_t^u$  tends to cancel out while the signal builds up.

The $\lambda$  parameter is chosen as the minimizer of (8) once a “great” number of candidates ${\lambda }_j\mathrm{<mo>;</mo>}j\mathrm{<mo>=</mo>1,}\dots \mathrm{,}\Lambda$  is tried. Therefore, (8) is modified as

$M_0^{\xi }\mathrm{<mo>=</mo><mo\; stretchy="false">(</mo>}\text{<mover accent="true"> G <mo>^</mo> </mover>}\mathrm{,}\widehat{\lambda }\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\arg \underset{\text{<mover accent="true"> G <mo>^</mo> </mover>}\mathrm{,}\widehat{\lambda }}{{\displaystyle \min }}f\mathrm{<mo\; stretchy="false">(</mo>}{\widehat{\sigma }}^2\mathrm{<mo\; stretchy="false">(</mo>}{\epsilon }^o\mathrm{<mo\; stretchy="false">)</mo>,}\text{G}\mathrm{<mo>|</mo>}{\xi }_t\mathrm{<mo\; stretchy="false">(</mo>}\lambda \mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>,}$           (9)

being ${\xi }_t\mathrm{<mo>=</mo>}f\mathrm{<mo\; stretchy="false">(</mo>}{x}^o\mathrm{<mo\; stretchy="false">(</mo>}\lambda \mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}$ .

What presented in (9) is a slightly modified version of the Hannan Quinn criterion, called here $mHQ$  (which will be justified below), whose minimum value delivers the optimal smoothing parameter ${\lambda }^{\ast }$ .

From the last equation, it is clear that both the optimal model and the optimal $\lambda$  are simultaneously estimated. However, while the $ML$  parameters’ estimation is conducted on the filtered versions ${\xi }_t$  of the observed signal, the variance parameter is estimated using the observed data $x^o$ . Such an operation is possible assuming $\mathbb{E}\mathrm{<mo\; stretchy="false">[</mo>}{\epsilon }_t\mathrm{<mo\; stretchy="false">]</mo><mo>=</mo>0}$  and serves the purpose to “force” the filtered output to be “not too far” from the observed signal $x_t^o$ . In practice, the Information Criterion selects an “optimal” model – built on a ${\lambda }^{\ast }$ –filtered version of the data (which therefore are used to build the ML estimations of the $SARIMA$  parameters) whose goodness of fit is computed with respect the observed (noisy) signal $x_t^o$ . In other words, mHQ criterion has a stronger penalty than the standard one. Such an additional weight (the penalty associated with the residual variance computed on $x^o$ ) is proportional to $\mathrm{<mo>|</mo>}{\sigma }_o^2-{\sigma }^2\mathrm{<mo\; stretchy="false">(</mo>}{\xi }_t\mathrm{<mo\; stretchy="false">)</mo><mo>|</mo>}$  and serves the purpose to prevent the H-Q criterion from choosing $\lambda$  values that can lead to over smoothing. In fact, the greater the value of $\lambda$  the greater the distance between $x_t^o$  and ${\xi }_t$  and therefore the greater the absolute values of the difference $\mathrm{<mo>|</mo>}{\sigma }_o^2-{\sigma }^2\mathrm{<mo\; stretchy="false">(</mo>}{\xi }_t\mathrm{<mo\; stretchy="false">)</mo><mo>|</mo>}$ . When no benefit can be obtained by applying the filter, we have $M^{\xi }\mathrm{<mo\; stretchy="false">(</mo>}HQ\mathrm{<mo\; stretchy="false">)</mo>}\ge M_0$  therefore $x^o$  and $\xi$  coincide and $\mathrm{<mo>|</mo>}{\sigma }_o^2-{\sigma }^2\mathrm{<mo\; stretchy="false">(</mo>}{\xi }_t{\mathrm{<mo\; stretchy="false">)</mo><mo>|</mo><mo>=</mo>0}}_t$ , i.e. the best SARIMA model is built on $x_t^o$  and its order selection follows the MICE rule according to equation (8). On the contrary, when a benefit from filtering is detected – i.e. $mHQ\mathrm{<mo><</mo>}HQ$  – therefore the best model is the one built on the filtered time series according to a $\lambda$  value ( ${\lambda }^{\ast }$ ) at which the information criterion (9) is minimized and its order selection follows the MICE rule. Again, such a conclusion is justified as the error term is “average out” by the filter. At this point, we have two different SARIMA models fitting $x^o$  and ${\xi }_t$  and generating two error vectors. Recalling (3), and consistently with the notation so far used, it will be

${\epsilon }_t^o\mathrm{<mo>=</mo>}{x}_t^o-{\text{G}}^o{\epsilon }_{t-1}\mathrm{<mo>;</mo>}{x}^o\mathrm{<mo>=</mo>1,}\dots \mathrm{,}T\mathrm{<mo>;</mo>}{\epsilon }_0\mathrm{<mo>=</mo>0.}$               (10)

and ${\epsilon }^{\xi }$  found in a similar way, replacing the superscript identified by the letter $o$  with the Greek letter $\xi$ .

Moreover, by virtue of Granger theorem,²⁴ it is always possible for an estimated SARIMA model order to depart from the “true” one when the input is affected by white noise components. In such a case, the CSS estimates resulting from the function $\log Lik\mathrm{<mo\; stretchy="false">(</mo>}\widehat{\Gamma }\mathrm{<mo>|</mo>}{x}^o\mathrm{<mo\; stretchy="false">)</mo>}$  can be regarded as converging to the pseudo–true value of the parameters’ vector and therefore the order detected by means of the $mHQ$  – which employs the function $\log Lik\mathrm{<mo\; stretchy="false">(</mo>}\widehat{\Gamma }\mathrm{<mo>|</mo>}x{i}_t\mathrm{<mo\; stretchy="false">)</mo>}$  for the estimation of the parameters – criterion will be closer to the “true” model order.

Signal improvement assessment
At this point, it is clear that for ${\Gamma }^0\mathrm{<mo>=</mo>}\Gamma \mathrm{<mo\; stretchy="false">(</mo>}\xi \mathrm{<mo\; stretchy="false">)</mo>}$  the lower bound in the gain delivered by the procedure is attained (i.e. no gain and $\lambda \mathrm{<mo>=</mo>1}$ ). In this section, a simple way for assessing the performances of the filter (if delivered) is presented. At this point it might be worthwhile emphasizing that the method presented in this paper is not designed to improve the forecasting performance of SARIMA models. In fact, the observations coming from the future, under an unchanged data generation mechanism, will be affected, in average, by the same noise components affecting the previous observations. However, projecting the filtered observations into the future is a useful exercise for at least two reasons: i) to verify that the forecasting performances delivered by the two SARIMA models are “not too far” from each other (meaning that the ${\xi }_t$  is a “bona fide” version of the original one; ii) in many cases, improvements in the forecasting performances have been observed. This can be explained by the fact that the order selection and estimation procedures tend to work better under less noisy signals.

With that said, the signal-to-noise ratio (SNR) has been used as an indicator of the quality of the filtered signal. It is defined as the ratio of the power of the signal (the meaningful information) to the power of the noise (the unwanted portion of the signal). In its “pure” form it is $SNR\mathrm{<mo>=</mo>}\frac{P_{signal}}{P_{noise}}$ , $P$  being the average power. In the framework at hand, both the power levels are estimated by their respective variances $\frac{{\widehat{\sigma }}_{x_t}^2}{{\widehat{\sigma }}_{{\epsilon }_t}^2}$ , being the error terms obtained using (10). Practically, by comparing $SN{R}^o\mathrm{<mo>=</mo>}\frac{{\widehat{\sigma }}_{x_t^o}^2}{{\widehat{\sigma }}_{{\epsilon }_t^o}^2}$  with $SN{R}^{\xi }\mathrm{<mo>=</mo>}\frac{{\widehat{\sigma }}_{{\xi }_t}^2}{{\widehat{\sigma }}_{{\epsilon }_t^{\xi }}^2}$ , one is able to assess how much of the unwanted components has been removed under ${\lambda }^{\ast }$

The algorithm

In what follows, the method is presented in the form of a step–by–step algorithm. In particular, steps 2-6 are conducted on the training set whereas the remaining 7–8 use the test set.

1. The time series under investigation is split in two sets: training and test;
2. A grid of tentative smoothing constant ${\lambda }_j\mathrm{<mo>;</mo>}j\mathrm{<mo>=</mo>1,}\dots \Lambda \in \mathfrak{L}$  is built;
3. A SARIMA model is estimated on the original time series $x^o$  using the H–Q criterion;
4. $\Lambda$  SARIMA models are estimated on $\Lambda$  filtered time series ${\xi }_j^{\mathrm{<mo\; stretchy="false">(</mo>}\lambda \mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo>=</mo>}{x}^o\mathrm{<mo\; stretchy="false">(</mo>}{\lambda }_j\mathrm{<mo\; stretchy="false">)</mo><mo>;</mo>}j\mathrm{<mo>=</mo>1,}\dots \mathrm{,}\Lambda$ ;
5. If $mHQ\ge HQ$  stop and use the unfiltered time series otherwise proceed;
6. The ${\lambda }^{\ast }$  smoothing constant associated with the winner model (according to the modified H–Q criterion mHQ) is then chosen;
7. Horizon $h\mathrm{<mo>;</mo>}h\mathrm{<mo>=</mo>1,}\dots \mathrm{,}H$  forecast mean square errors are computed for both the filtered and unfiltered series;
8. Both the statistics are compared for consistency.

Clearly, if one is not interested in forecast comparison – e.g. to speed up the computations and have faster outcomes – then steps 1, 7 and 8 can be skipped.

The experiment presented in this paragraph is aimed at showing the capabilities of the proposed method. A time series – collected by the popular search engine Google and recovered in an “ad hoc” repository called Google Trends – will be employed. It refers to the relative number of Italian Internet users which inputted the keyword “caviar” (caviale, in Italian) into the search engine Google, between January 2004 and June 2018. The data have been download on June 6th 2018 (h 03:29 PM GMT). Time series length is 174, split in training set – ending in December 2014 (132 data) – and test set (42 data), ending in June 2018. The last datum has been purposely left in the data set, for it is an even more unstable one (in fact, it is supported only by a small portion of the month of reference). A visual inspection of the series, depicted in Figure 1, shows many irregularities, a pronounced, non stationary seasonality as well as a non stationary trend. Eteroschedasticity is also an issue. Such a behavior might be considered typical for a luxury good especially when observed, as in the present case, for a long period of time.

Figure 1 Relative number of the keyword “caviar” searched by the Italian Internet users. Training set (continuous line) and test set (dashed line).

The whiteness of the residuals has been tested using a test of the type “Portmanteau’, i.e. the Ljung-Box test,²⁵ which is designed to account for the sum of the autocorrelations (denoted by the Greek letter $\rho$ ) until a predetermined, arbitrary lag $K$ , i.e.
Publication

$Q\mathrm{<mo>=</mo>}T\mathrm{<mo\; stretchy="false">(</mo>}T+\mathrm{2<mo\; stretchy="false">)</mo>}{\displaystyle \sum _{k\mathrm{<mo>=</mo>1}}^K}\frac{{\widehat{\rho }}_k{}^2}{T-K}\stackrel{H_0}{\sim }{\chi }_{\mathrm{<mo\; stretchy="false">(</mo>}K\mathrm{<mo\; stretchy="false">)</mo>}}^2\mathrm{<mo>;</mo>}{Q}^{rej}\mathrm{<mo>></mo>}{\chi }_{\mathrm{<mo\; stretchy="false">(</mo>1}-\alpha \mathrm{,}K\mathrm{<mo\; stretchy="false">)</mo>}}^2\mathrm{.}$

Here $Q^{rej}$  is, for significance level $\alpha$ , the critical region for rejection of the hypothesis of randomness under the ${\chi }^2$  distribution.

The estimated SARIMA order in the case of $x^o$  is $\mathrm{<mo\; stretchy="false">(</mo>1,1,0<mo\; stretchy="false">)</mo>}\times \mathrm{<mo\; stretchy="false">(</mo>2,0,0<mo\; stretchy="false">)</mo>}$  with an HQ criterion of 4.97. In the case of ${\xi }_t$ , the estimated SARIMA order is $\mathrm{<mo\; stretchy="false">(</mo>0,0,1<mo\; stretchy="false">)</mo>}\times \mathrm{<mo\; stretchy="false">(</mo>0,1,1,<mo\; stretchy="false">)</mo>}$ , for a modified HQ criterion equal to 4.05. Therefore, the filter is applied. The effects of its application are presented in Figure 2, where a more regular signal can be noticed. Such an impression is confirmed by looking at the SNRs, which is equal to 6.88 and 7.06 respectively for the unfiltered and filtered time series. The optimal smoothing constant ${\lambda }^{\ast }\mathrm{<mo>=</mo>.805}$ , has been obtained at iteration number 82 (see Figure 3), using a grid of 121 $\lambda$  values (equally spaced from 0.400 to 1; incremental step= 0.005). The residuals (computed in both the cases on the original time series) pass the ‘Portmanteau’ test (computed for the lags 1–24) only in the case of ${\xi }_t$  (p value $\mathrm{<mo>=</mo>0.382}$ ), whereas for the raw signal we have p value $\approx 0$ . Such a conclusion is consistent with the patterns shown by the empirical autocorrelation functions (see Figure 4) for $x^0$  (graph a) and ${\xi }_t$  (graph b). In the same graph, the residual density distributions for $x^0$  (graph c) and ${\xi }_t$  (graph d) confirm a better behavior when the proposed method is applied. Finally, in Figure 5, the prediction performances obtained using ${\xi }_t$  appears to be slightly better until lag 8.

Figure 2 Unfiltered (continuous line) and filtered (dashed line) time series.

Figure 3 Values of the modified Hanna–Quinn information criterion for each iteration (attempted $\lambda$  value). The vertical line shows the iteration number at which the minimum mHQ value is found.

Figure 4 Autocorrelation functions for $x_t^o$  (graph a), ${\xi }_t$  (graph b) and residual density for $x_t^o$  (graph c) and ${\xi }_t$  (graph d).

Figure 5 Mean Absolute Error of the predictions obtained from the two SARIMA models applied on the raw time series (continuous line) and the filtered time series (dashed line). Time horizon 1–12.

One effective way to reduce noise components affecting our time series is to use a First Order Discrete Time Infinite Impulse Response Filter. As it is well known, its filtering performances are critically dependent on the smoothing constant, which therefore has to be carefully fine–tuned. The method proposed in this paper provides the practitioner with a useful – and reasonably fast – tool to filter out such unwanted components. However, that the objective function used – based on the Hannan Quinn information criterion – has given good results, does not imply that is the best possible one. For example, with small data set it could be worth considering different information criteria, specifically designed for such a case, e.g. the second-order AIC.²⁶ Finally, different stochastic models can be always considered, according to the underlying data generating process. For instance, whenever the time series shows memory characteristics of the type long range, the proposed procedure can be associated to a suitable model (e.g. of the class Autorgressive Fractional Integrated Moving Average).

None.

Author declares that there is no conflict of interest.

1. Richard T O'Connell, Anne B Koehler. Forecasting, time series, and regression: An applied approach, volume 4. USA: South–Western Pub; 2005.
2. Ruey S Tsay. Analysis of Financial time series. Volume 543. USA: John Wiley & Sons; 2005.
3. Harold Wayne Sorenson. Kalman ltering: theory and application. IEEE. 1985.
4. Phillip G Gould, Anne B Koehler, J Keith Ord, et al. Forecasting time series with multiple seasonal patterns. European Journal of Operational Research. 2008;191(1):207–222.
5. David L Donoho. De–noising by soft–thresholding. IEEE transactions on information theory. 1995;41(3):613–627.
6. Albert Cohen, Ingrid Daubechies, Bjorn Jawerth, et al. Multiresolution analysis, wavelets and fast algorithms on an interval.  Comptes Rendus de l Académie des Sciences. 1993;316(5):417–421.
7. Holger Kantz, Thomas Schreiber. Nonlinear time series analysis, volume 7. USA: Cambridge university press; 2004.
8. S Makridakis. Forecasting: Methods and applications. 1998.
9. Torsten kohler, Dirk Lorenz. A comparison of denoising methods for one dimensional time series. Bremen: Germany, University of Bremen; 2005.
10. George EP Box, Gwilym M Jenkins, Gregory C Reinsel, et al. Time series analysis: forecasting and control. USA: John Wiley & Sons; 2015.
11. Peter J Brockwell, Richard A Davis. Time series: theory and methods. Springer Science & Business Media. 2013.
12. JH Stock, MW Watson. Introduction to econometrics. Boston: Addison–wesley; 2003.
13. MJ Evanich, AO Newberry, LD Partridge. Some limitations on the removal of periodic noise by averaging. Journal of Applied Physiology. 33(4):536–541;1972.
14. CC Holt. Forecasting trends and seasonals by exponentially weighted averages. carnegie institute of technology. Technical report, Pittsburgh ONR memorandum. 1957.
15. Adrien Cahuzac, Jerome Boudet, Pierre Borgnat, et al. Smoothing algorithms for mean–ow extraction in large–eddy simulation of complex turbulent ows. Physics of uids. 2010;22(12):125104.
16. Layth C Alwan, Harry V Roberts. Time series modelling for statistical process control. Statistical Process Control in Automated Manufacturing. 1988;15:45–65.
17. J Stuart Hunter. The exponentially weighted moving average. Journal of quality technology. 1986;18(4):203–210.
18. Gery Casiez, Nicolas Roussel, Daniel Vogel. 1 euro Filter: a simple speed–based low–pass Filter for noisy input in interactive systems. In: Proceedings of the SIGCHI Conference on Human Factors in Computing Systems. ACM. 2012;2527–2530.
19. Adrien Cahuzac, Jerome Boudet, Pierre Borgnat, et al. Smoothing algorithms for mean–ow extraction in large–eddy simulation of complex turbulent ows. Physics of uids. 2010;22(12):125104.
20. Rob Hyndman, Anne B Koehler, J Keith Ord, et al. Forecasting with exponential smoothing: the state space approach. Springer Science & Business Media. 2008.
21. Fat Codrut, A Maria, Dezsi Eva, et al. Exchange–rates forecasting: Exponential smoothing techniques and arima models. Annals of Faculty of Economics. 2011;1(1):499–508.
22. Athanasios Papoulis. Signal analysis. volume 191. New York: McGraw–Hill; 1977.
23. Edward J Hannan, Barry G Quinn. The determination of the order of an autoregression. Journal of the Royal Statistical Society. Series B (Methodological). 1979; 190–195.
24. Clive WJ Granger. Essays in econometrics: collected papers of Clive WJ Granger. Volume 32. USA: Cambridge University Press; 2001.
25. Greta M Ljung, George EP Box. On a measure of lack of t in time series models. Biometrika. 1978;65(2):297–303.
26. Kenneth P Burnham, David R Anderson. Model selection and multimodel infer–ence: a practical information–theoretic approach. Springer Science & Business Media. 2003.