Anomaly Generation and Detection
Large spikes in data series have typically been associated with the term outliers. These lead to contamination of normal process variation under most types of dynamics and is therefore the type of anomaly which is most general across different types of data generating mechanisms.
While such spikes usually decay almost instantaneously, others may decay at a slower pace (temporary shifts) and some may even be permanent (level shifts). One generating mechanism for these anomalies is through a recursive AR(1) filter. The figure below shows an example of how such anomalies can be generated with different decay rates. The half-life of the decay of such shocks can be easily derived from the theory on ARMA processes and is equal to \(-\text{ln}\left(2\right)/\text{ln}\left(\alpha\right)\), where \(\alpha\) is the AR(1) coefficient. Thus, the closer this gets to 1, the closer the shock becomes permanent whereas a coefficient close to 0 is equivalent to an instantaneous decay.
Code
x <- rep(0, 100)
x[50] <- 1
oldpar <- par(mfrow = c(1,1))
par(mar = c(2,2,1,1), cex.main = 0.8, cex.axis = 0.8, cex.lab = 0.8)
plot(filter(x, filter = 0, method = "recursive", init = 0), ylab = "", xlab = "")
lines(filter(x, filter = 0.5, method = "recursive", init = 0), col = 2, lty = 2)
lines(filter(x, filter = 0.8, method = "recursive", init = 0), col = "orange", lty = 2)
lines(filter(x, filter = 1, method = "recursive", init = 0), col = 3, lty = 3, lwd = 2)
grid()
legend("topleft", c("Additive Outlier (alpha = 0.0)",
"Temporary Change (alpha = 0.5)",
"Temporary Change (alpha = 0.8)",
"Permanent Shift (alpha = 1.0)"), cex = 0.8,
lty = c(1,2,2,3), bty = "n", col = c(1,2,"orange",3))
Identifying different types of anomalies is usually undertaken as an integrated part of the assumed dynamics of the process. An established and flexible framework for forecasting economic time series is the state space model, or one of its many variants. For example, the basic structural time series model of (Harvey 1990) is a state of the art approach to econometric forecasting, decomposing a series into a number of independent unobserved structural components:
\[ \begin{aligned} y_t &= \mu_t + \gamma_t + \varepsilon_t,\quad &\varepsilon_t\sim N(0,\sigma)\\ \mu_t &=\mu_{t-1} + \beta_{t-1} + \eta_t,\quad &\eta_t\sim N(0,\sigma_\eta)\\ \beta_t &= \beta_{t_1} + \zeta_t,\quad &\zeta_t\sim N(0,\sigma_\zeta)\\ \gamma_t &= - \sum^{s-1}_{j=1}\gamma_{t-j} + \omega_t,\quad &\omega_t\sim N(0,\sigma_\omega)\\ \end{aligned} \]
which uses random walk dynamics for the level, slope and seasonal components, each with independent noise components. Variations to this model include the innovations state space model which assumes fully correlated components and is observationally equivalent to this model but does not require the use of the Kalman filter since it starts from a steady state equilibrium value for the Kalman Gain (G) and hence optimizing the log likelihood is somewhat simpler and faster. In this structural time series representation, additive outliers can be identified by spikes in the observation equation; level shifts by outliers in the level equation; and a large change in the slope as an outlier in the slope equation (and similar logic applies to the seasonal equation). Therefore, testing for such outliers in each of the disturbance smoother errors of the models can form the basis for the identification of anomalies in the series (see (De Jong and Penzer 1998)).
Direct equivalence between a number of different variations of this
model and the more well known ARIMA model have been established. In the
ARIMA framework, a key method for the identification of anomalies is
based on the work of (Chen and Liu 1993),
and is the one we use in our approach. This approach uses multi-pass
identification and elimination steps based on automatic model selection
to identify the best anomaly candidates (and types of candidates) based
on information criteria. This is a computationally demanding approach,
particularly in the presence of multiple components and long time series
data sets. As such, we have implemented a slightly modified and faster
approach, as an additional option, which first identifies and removes
large spikes using a robust STL decomposition and MAD criterion, then
de-seasonalizes and detrends this cleaned data again using a robust STL
decomposition, before finally passing the irregular component to the
method of (Chen and Liu 1993), as
implemented in the tsoutliers
package. Finally, the outliers from the first stage are combined with
the anomalies identified in the last step (and duplicates removed),
together with all estimated coefficients on these anomalies which can
then be used to clean the data. The tsaux package
provides 2 functions: auto_regressors automatically returns
the identified anomalies, their types and coefficients, whilst
auto_clean automatically performs the data cleaning step.
The functions allow for the use of the modified approach as well as the
direct approach of passing the data to the tso
function.
One shortcoming of the current approach is that for temporary change identification, the AR coefficient defaults to 0.7 and there is no search performed on other parameters. Adapting the method to search over different values is left for future research.
We start with an illustration of data simulation and contamination. Starting with a simulated monthly series of level + slope + seasonal, we then proceed to contaminate it with outliers, temporary changes and level shifts.
Code
set.seed(154)
r <- rnorm(300, mean = 0, sd = 6)
d <- seq(as.Date("2000-01-01"), by = "months", length.out = 300)
mod <- initialize_simulator(r, index = d, sampling = "months", model = "issm")
mod <- mod |> add_polynomial(order = 2, alpha = 0.22, beta = 0.01, b0 = 1, l0 = 140)
mod <- mod |> add_seasonal(frequency = 12, gamma = 0.7, init_scale = 2)
y <- xts(mod$simulated, mod$index)
oldpar <- par(mfrow = c(1,1))
par(mfrow = c(2, 2), mar = c(2,2,1,1), cex.main = 0.8, cex.axis = 0.8, cex.lab = 0.8)
plot(as.zoo(y), main = "Original", ylab = "", xlab = "", col = "blue")
grid()
y_contaminated <- y |> additive_outlier(time = 25, parameter = 1)
plot(as.zoo(y_contaminated), main = "+Additive Outlier[t=25]", ylab = "", xlab = "", col = "green")
grid()
y_contaminated <- y_contaminated |> temporary_change(time = 120, parameter = 0.75, alpha = 0.7)
plot(as.zoo(y_contaminated), main = "+Temporary Change[t=120]", ylab = "", xlab = "", col = "purple")
grid()
y_contaminated <- y_contaminated |> level_shift(time = 200, parameter = -0.25)
plot(as.zoo(y_contaminated), main = "+Level Shift[t=200]", ylab = "", xlab = "", col = "red")
grid()
We now pass the contaminated series to the anomaly identification
function auto_regressors. The tso function
directly identifies the correct number, timing and types of anomalies.1
Code
xreg <- auto_regressors(y_contaminated, frequency = 12,
return_table = TRUE, method = "full")
#> Warning in locate.outliers.oloop(y = y, fit = fit, types = types, cval = cval,
#> : stopped when 'maxit.oloop = 4' was reached
xreg[,time := which(index(y_contaminated) %in% xreg$date)]
print(xreg)
#> Index: <type>
#> type date filter coef time
#> <char> <Date> <num> <num> <int>
#> 1: AO 2002-01-01 0.0 171.4351 25
#> 2: TC 2009-12-01 0.7 184.6145 120
#> 3: LS 2016-08-01 1.0 -103.4125 200While we can call directly the auto_clean function,
we’ll instead show how to use the returned table to decontaminate the
data.
Code
x <- cbind(additive_outlier(y_contaminated, time = xreg$time[1], add = FALSE),
temporary_change(y_contaminated, time = xreg$time[2], alpha = xreg$filter[2],
add = FALSE),
level_shift(y_contaminated, time = xreg$time[3], add = FALSE))
x <- x %*% xreg$coef
y_clean <- y_contaminated - x
oldpar <- par(mfrow = c(1,1))
par(mfrow = c(2,1), mar = c(2,2,1,1), cex.main = 0.8, cex.axis = 0.8, cex.lab = 0.8)
plot(as.zoo(y_contaminated), main = "Anomaly Decomtanination", xlab = "", ylab = "", ylim = c(90, 700))
lines(as.zoo(y_clean), col = 2, lty = 2, lwd = 2)
lines(as.zoo(y), col = 3, lty = 3, lwd = 2.5)
grid()
legend("topleft", c("Contaminated", "Clean", "Original"), col = 1:3, bty = "n", lty = 1:3, cex = 0.8)
plot(zoo(x, index(y)), main = "Anomalies", xlab = "", ylab = "")
As can be observed, the routine does a very good job of identifying
the 3 types of anomalies and their timing. Additionally, the function
has an argument h which allows to project the regressors
into the future. This is particularly useful for location shifts which
propagate indefinitely as a vector of ones, and temporary shifts which
have a defined decay patterns (towards zero).
