Preprocessing: stationarity, transforms and denoising

1. Why raw data is not enough

The deck opens by stating the problem without softening it: raw data could be unsuitable for forecasting algorithms. Two reasons are given, and they are different in kind.

• Some algorithms require data with specific characteristics. This is a hard constraint: an AR model is only feasible for stationary processes, full stop.
• Some algorithms perform better when data is presented in some specific format. This is a soft one: an MLP with sigmoid activations will train, but badly, on unscaled inputs.

So data preprocessing has the objective of making raw data better compatible with the selected algorithm. Note the direction of the arrow: preprocessing is chosen for a model, not in the abstract. And then the warning that should be printed on the inside of every student's eyelids: there are many different preprocessing options, and no deterministic rule on what to do.

What preprocessing is aiming at is stated in one line under a picture of white noise: that is the easiest time series to forecast, and data preprocessing tries to revert to this case. It sounds paradoxical — flat noise is unforecastable — but it is exactly right. Every structure you can strip out and model explicitly (trend, seasonality, changing variance) is structure the fitting algorithm no longer has to discover for itself. What should be left over is the part nobody can predict.

2. The three fundamental moments

A moment is a specific quantitative measure of the shape of a function. The statistical description of a time series rests on three of them.

These three quantities are the vocabulary in which stationarity is about to be defined, and the reason the definition has exactly three clauses.

3. Stationarity, in three conditions

A time series is stationary if:

1. its mean is constant over time (stationarity in means);
2. its variance is constant over time (stationarity in variance);
3. the relation between values separated by k time periods depends only on k, not on time (stationarity in covariance).

Or, compressed into the sentence worth memorising: a series is stationary if its distribution is invariant with respect to translations in time. Slide the window anywhere along the series and it looks statistically the same.

The slides give two contrasting examples. Daily total female births in California, 1959 is a stationary process: no trends, stationary in mean and variance, no seasonality or cyclicity. Accidental deaths in the U.S.A. (1973-1978) is stationary in mean and variance but with seasonality — a reminder that the three conditions fail independently, and that a series can be well behaved in two of them and useless in the third.

4. The difference transform

We try to make a series stationary in mean by differencing: subtracting from the current value the value of the previous time moment.

first order:     ∇yₜ = yₜ - yₜ₋₁
second order:    ∇²yₜ = ∇yₜ - ∇yₜ₋₁ = yₜ - 2yₜ₋₁ + yₜ₋₂
seasonal (L):    ∇ₜyₜ = yₜ - yₜ₋₄

If the series has a linear trend, the resulting differenced series is stationary in means. More generally, the difference transform subtracts from each entry the entry back-shifted by a constant time lag: a simple way for removing a systematic structure from the series. The process can be repeated — difference the differenced series — to remove second-order trends. A seasonal structure is removed by subtracting the observation from the prior season: e.g. 12 time steps before, for monthly data with a yearly seasonal structure.

On the fil rouge

Applying a first difference to the TV sales series gives, for t = 2 … 20:

The slides annotate this table with two words: removed the linear trend! The differenced series oscillates around a roughly constant level instead of climbing. What it has not removed is the seasonal alternation, still visible as the regular sign pattern, nor the growth of the swing amplitude.

The invertibility is the point. difference() throws away the first interval observations and produces increments; invert_difference() replays them from the original head. Every forecast produced on differenced data must go back through this function before anyone sees it.

5. The log transform

Where differencing fixes the mean, the log transform attacks the variance — and, combined with diff, both.

Why the logarithm turns growth into addition

The slides derive the key approximation carefully. The function y = ln(x) is tangent at (1, 0) to the line y = x - 1, so ln(1+r) ≈ r when r is small with respect to 1 (say r ≤ 0.2). Now suppose x has a small percent increase, for example 5%: it changes from x to (1+r)x with r = 0.05. Then

ln((1+r)x) = ln(x) + ln(1+r) ≈ ln(x) + r

If x increases by 5%, the log of x increases from ln(x) to ln(x) + 0.05. Since a percent variation of a series is (y_t - y_t-1) / y_t-1, the conclusion is that a percent variation is almost equal to ln(y_t) - ln(y_t-1) — the log-difference is the growth rate.

The slides check it numerically on the first airline passenger records:

The last two columns agree to three decimals, and disagree slightly more when r is largest — exactly as the approximation predicts.

Three properties

The two-step ladder on the fil rouge

Taking logs first and then differencing gives the pair of columns the slides annotate as removed the trend on variance! and removed the trend on means and variance!

6. Power transforms and Box-Cox

The log is one member of a family. A power transform is a family of functions that create a monotonic transformation of data using power functions, used to stabilise variance and make the data more normal-distribution-like. On a time series, this can remove a change in variance over time. Popular examples: the log transform, the square root transform, and the generalised Box-Cox transform, which transforms positive data according to

          ⌈ (xᵥ - 1) / λ     if λ ≠ 0
BC(x) =   ⌊ ln(x)             if λ = 0

and where the λ value is chosen to maximise the log-likelihood function, making the data as normally distributed as possible. So Box-Cox with λ = 0 is the log transform; λ = 0.5 is essentially a square root; λ = 1 leaves the shape alone. Choosing λ automatically is what makes it attractive: you do not have to guess how strong the variance stabilisation should be.

from math import log, exp
from scipy.stats import boxcox
import pandas as pd

def invert_boxcox(value, lam):
    if lam == 0:                       # log case
        return exp(value)
    return exp(log(lam * value + 1) / lam)   # general case

data = pd.read_csv('FilRouge.csv', usecols=[1]).values.flatten()
transformed, lmbda = boxcox(data)
print(f"data={transformed}, lmbda={lmbda}")
inverted = [invert_boxcox(x, lmbda) for x in transformed]

Note the shape of that snippet: transform, use, invert. It is the same three-beat rhythm as the difference transform, and it is the rhythm the exam project has to implement for whatever preprocessing it chooses.

7. Standardization and normalization

Two rescalings that are constantly confused, so the slides define them separately.

xₛ = (xᵢ - x̄) / s

from sklearn.preprocessing import StandardScaler
data = pd.read_csv('FilRouge.csv', usecols=[1]).values   # column vector!
transformer = StandardScaler()
transformer.fit(data)
transformed = transformer.transform(data)     # mean 0, std 1
inverted    = transformer.inverse_transform(transformed)

xₙ = (xᵢ - xᵐᵢₙ) / (xᵐₐₓ - xᵐᵢₙ)

from sklearn.preprocessing import MinMaxScaler
data = pd.read_csv('FilRouge.csv', usecols=[1]).values
transformer = MinMaxScaler()
transformer.fit(data.T)
transformed = transformer.transform(data)     # min 0, max 1
inverted    = transformer.inverse_transform(transformed)

8. Missing values

The slides refuse to be encouraging here, and they are right to: there is NO good way to deal with missing data, and imputation does not necessarily give better results. Imputation algorithms divide into simple but fast approaches, like mean imputation, and more advanced algorithms that need more computation time, like Kalman smoothing.

ds  = pd.Series(y)
ds1 = ds.ffill()              # LOCF
ds2 = ds.bfill()              # NOCB
ds3 = ds.fillna(ds.mean())    # mean imputation
ds4 = ds.interpolate()        # linear interpolation

9. Outliers, dimensionality reduction, augmentation

From here the deck marks its own content as advanced topics: sketches of subjects that deserve much more grounding, offered as pointers to issues that could need to be addressed in order to enhance forecasting effectiveness. These are also the "advanced extensions" that lift an exam project into the top band.

Outliers, or anomaly detection

An outlier is an observation that significantly differs from other observations of the same feature. Time series with outliers have three characteristics: there is a systematic pattern (deterministic) and some variation (stochastic); only a few data points are outliers; and outliers are significantly different from the rest of the data.

An anomaly detection algorithm must either label each time point as anomaly / not anomaly, or forecast a signal for some point and test whether the actual value deviates from the forecast enough to be deemed an anomaly. Two decomposition-based approaches are given:

• STL decomposition: split the signal into seasonal, trend and residue; analyse the deviation of the residue and introduce a threshold on it, and you have an anomaly detector. Similarly with ARIMA or exponential smoothing.
• Classification and regression trees: use supervised learning to teach trees to classify anomaly and non-anomaly points.

from statsmodels.tsa.seasonal import seasonal_decompose

df     = pd.read_csv("traffico16.csv")
col    = "ago2"
series = df[col].fillna(value=df[col].mean())
result = seasonal_decompose(series, model='additive', period=7)
Resid  = result.resid
std    = Resid.std()

plt.plot(Resid, "o", label="datapoints")
plt.hlines(0, 0, len(Resid))
plt.hlines( 1.5*std, 0, len(Resid), color="red", label="std limits")
plt.hlines(-1.5*std, 0, len(Resid), color="red")
plt.legend(); plt.show()

Dimensionality reduction

In multivariate analysis, datasets can be defined over too many dimensions — attributes can number in the hundreds or thousands. Five reasons are given for reducing them: less storage space; less computation and training time; some algorithms do not perform well in large dimensions; it takes care of multicollinearity by removing redundant features; and it helps in visualising data.

It can be done in two different ways: (1) keeping only the most relevant variables from the original dataset — properly called feature selection; (2) finding a smaller set of new variables, each a combination of the input variables, containing basically the same information — sometimes called dimensionality reduction proper.

Data augmentation and bootstrapping

Data augmentation means generating artificial data in order to reduce the variance of the output (to reduce the number of errors). Widely used in image analysis, but useful also when dealing with too short time series — which, given the fil rouge has twenty points, is a live concern. Common techniques for time series classification: extrapolation (fill relevant fields with values based on heuristics — a trend-seasonality model can compute predicted values for instants not covered by the series); tagging (common records tagged to a group); aggregation (values estimated using averages and means); probability (values populated based on the probability of events).

Bootstrapping is useful both for augmentation and for uncertainty estimation: it generates many pseudo-samples from your observed series to empirically estimate the distribution of a statistic without assuming a parametric model. In forecasting, this translates into confidence intervals derived from data alone.

Why it matters, in the slides' own list: no distributional assumptions; works with non-linear models; forecast interval estimation; model selection validation; and it handles short series well.

10. Noise reduction and spectral analysis

Denoising can be heavily borrowed from signal processing techniques to improve the signal-to-noise ratio (SNR). The general recipe is analysis in the frequency domain: apply a Fourier transform (or a wavelet transform) to the series, and then an appropriate filter.

Fourier, DFT, FFT

The premise is the classic one: any periodic function (or function defined on an interval) can be reconstructed by adding an appropriate, infinite number of sinusoids with appropriate amplitudes and frequencies. Spectral analysis is a technique to discover periodicities underlying time series, and to that end we first transform data from the time domain to the frequency domain.

The Fourier transform of a function of time is a complex-valued function of frequency, whose magnitude represents the amount of that frequency present in the original function, and the transformation is invertible. For data series, the Discrete Fourier Transform (DFT) converts a finite sequence of equally spaced samples into a same-length sequence of equally spaced samples of a complex-valued function of frequency. The Fast Fourier Transform (FFT) is an algorithm that computes the DFT efficiently, in O(n log n).

Two facts about the output that you must not forget when coding it:

• Symmetry. For a sequence of real values, the transform satisfies FT(k) = FT(n-k) — symmetry around n/2 — so only half of the transform data is needed.
• Nyquist-Shannon. A signal containing frequency components up to frequency F must be sampled at a rate of at least 2F to be fully reconstructed. Conversely, sampling at frequency F can reconstruct the signal only up to the F/2 component — its Nyquist frequency.

And the sentence that explains why any of this is a preprocessing technique: noise has high frequency; cutting out high frequencies prevents from representing noise. A second, less obvious use is noted too: spectral decomposition can be used as externals — as exogenous variables for SARIMAX — because sinusoids extend into the future; usually a low number of frequencies is used.

Doing it, step by step

The pipeline needs a stationary signal first: detrend the series and keep only the seasonal component, either via seasonal_decompose(...).seasonal or via scipy.signal.detrend(). Only half of the transformed values are used, because of FFT symmetry. And remember that periods (wavelengths) are the inverse of frequencies.

detrended_series = signal.detrend(orig_values, type="linear")
fft = np.fft.fft(detrended_series)
n   = len(fft)
one_sided_fft = fft[:int(n/2)+1]        # conjugate symmetry: half is enough
Fs      = 1                              # 1 observation per unit of time
nyquist = Fs/2                           # only frequencies up to Fs/2 matter
power   = np.abs(one_sided_fft)**2 / n   # magnitude squared over length

freq_sample = np.array(range(n)) / n     # normalized frequency grid
x_freq      = freq_sample[:int(n/2)+1]
frequency_max = freq_sample[np.where(power == np.amax(power))]
print("Max period {}".format(1/frequency_max))

# ---------- low pass filter
f_cut = 1/4      # keep only periods longer than the denominator
mask  = np.logical_or(freq_sample <= f_cut, freq_sample >= 1 - f_cut)

fft_lowpass = fft * mask                 # reconstruction needs the FULL fft array
reconstructed_detrended = np.real(np.fft.ifft(fft_lowpass))

t     = np.arange(len(orig_values))       # put the trend back
coeff = np.polyfit(t, orig_values, 1)
trend = coeff[0] * t + coeff[1]
reconstructed_series = reconstructed_detrended + trend

The mask line deserves a second look: it keeps frequencies below the cut-off and above 1 - f_cut, because the second half of the FFT array is the mirror image of the first. Masking only the left half would destroy the symmetry and give a complex reconstruction. A convenience alternative is np.fft.rfft / np.fft.irfft, which keep only half the frequencies for you, or scipy.signal.periodogram, which returns the spectrum directly.

11. Autoencoders and Kalman filters

Autoencoders: a neural approach to noise filtering

An autoencoder is a simple neural network trained with the same dataset both as input and as output of the network, where the network has fewer hidden neurons than the input/output ones. The bottleneck is the whole trick: the network cannot copy its input through, so it must learn a compressed description. Since noise has much higher dimensions than regular data, this process reduces the noise. The slides note that this is similar to PCA and is a general dimensionality reduction technique.

class Autoencoder(nn.Module):
    def __init__(self, input_dim, latent_dim):
        super(Autoencoder, self).__init__()
        self.encoder = nn.Sequential(
            nn.Linear(input_dim, 64), nn.ReLU(),
            nn.Linear(64, latent_dim), nn.ReLU()
        )
        self.decoder = nn.Sequential(
            nn.Linear(latent_dim, 64), nn.ReLU(),
            nn.Linear(64, input_dim), nn.Sigmoid()   # output in [0, 1]
        )
    def forward(self, x):
        latent = self.encoder(x)
        return self.decoder(latent)

Three further uses are listed beyond denoising. Anomaly detection: autoencoders reconstruct "normal" patterns well but struggle with anomalies, so the reconstruction errors correspond to outliers, which can be removed or flagged before the forecasting model. Handling missing data: they can be trained to impute gaps, giving a cleaner, complete dataset. Dimensionality reduction and feature extraction: rather than feeding raw, high-dimensional series into a forecasting model, use the encoder output as features — these learned features tend to capture temporal structure (trends, seasonality, correlations) more effectively than hand-crafted ones.

Kalman filters

The Kalman filter is an algorithm that provides estimates of some unknown variables given their measurements observed over time. It is called a filter because it filters out noise. The precise statement of its optimality is worth quoting: if all noise is Gaussian, the Kalman filter minimises the mean square error of the estimated parameters; if the noise is NOT Gaussian, the Kalman filter is the best linear estimator, but non-linear estimators may be better.

Its practical appeal is architectural: Kalman filters have a relatively simple form and require small computational power — they do not need to keep any history other than the previous state — and they are very fast, which makes them well suited to real-time problems and embedded systems.

The model assumes the true state at time t evolves from the state at t-1:

xₜ = Fₜ xₜ₋₁ + Bₜ uₜ + wₜ        state transition
zₜ = Hₜ xₜ + vₜ                  observation

The slides add the practical footnote: matrices need to be identified for the specific application. That identification is the real work.

Predict and update

The algorithm consists of two stages, cycled for each new measurement — and the slides note that the terms prediction and update are often called propagation and correction in different literature.

1. Prediction phase. The filter predicts the next value in the series from the current value and its parameters, and predicts the uncertainty (error, variance) in that prediction according to the process noise present in the system.
2. Update phase. The filter reads the actual measurement, computes the difference between predicted and measured value, and decides which value to keep — both being unreliable — by calculating the Kalman Gain. A new value and a new uncertainty follow, and are fed back into the prediction step.

The gain itself has a memorable form:

K = Error in Prediction / (Error in Prediction + Error in Measurement)

Its value ranges from 0 to 1: near 0, the predicted value is close to the actual value; near 1, the measured value is close to the actual value. The cycle continues until the difference between predicted and measured value converges to zero.

from pykalman import KalmanFilter

measurements = y.to_numpy()
kf = KalmanFilter(transition_matrices=[1],
                  observation_matrices=[1],
                  initial_state_mean=measurements[0],
                  initial_state_covariance=1,
                  observation_covariance=10,
                  transition_covariance=10)
state_means, state_covariances = kf.filter(measurements)

The slides report the outcome on their example series as randomness reduced by 50%, and add the sentence that ties this whole chapter back to the rest of the course: once you have the filtered series, you can make forecasts by training the model you have chosen using the "clean" data.

Check your understanding