Anomaly detection using vector autoregression

Question

I want to detect anomalies in multivariate time series using statistical approaches. In particular. I want to use a vector autoregression model like VAR, VARMA or VARIMA, to predict a time stamp $x_t$ given $x_{t-1},...,x_{t-lag}$. Here $x_t$ is multivariate timestamp having dimension $d$.
Then, I would use the deviation between $x_t$ and $\hat{x}_t$ ($\hat{x}_t$ is the prediction of the VAR model) to assign an anomaly score to $x_t$.
One way would be to use the euclidean distance between $x_t$ and $\hat{x}_t$. However this is critical, because it ignores the relationship between the $d$ variables. Some papers suggest to use the Mahalanobis distance: here.
I wanted to ask if anyone has experience what is an appropriate way to compute the deviation between $x_t$ and $\hat{x}_t$. I would appreciate any suggestion. It would also be very helpful if you could give a hint to the python function of it, if it does exist.

Answer 1

A principled approach is to make probabilistic predictions. This gives a probability distribution over possible values at each timepoint, rather than a simple point prediction. The methods you described can provide this. The model should capture the ordinary behavior of the system, so it should ideally be fit to data containing no anomalies. Then, a point can be identified as a possible anomaly if it's unlikely according to the model.

Given the predictive distribution at a particular timepoint, one can identify the highest density region (HDR) containing $100 (1-\alpha)\%$ of the probability mass. You can think of this as a multivariate analog of a prediction interval. But, it doesn't account for uncertainty in estimating the model parameters. Also, it may consist of multiple disjoint subregions if the predictive distribution is multimodal. Points generated according to the model should fall within the HDR with probability $1-\alpha$. If a point falls outside the HDR, call it a possible anomaly. The expected false positive rate is $\alpha$, since this is the probability that a point generated according to the model will fall outside the region.

In general, the HDR of a continuous probability distribution can be defined as the set of all possible points where the density is greater than some threshold. And, the threshold is defined such that integrating the distribution over this set yields $1-\alpha$. For details and calculation methods, see Hyndman (1996). From this definition, it's apparent that that points outside the HDR (which we're calling possible anomalies) have low predictive density.

If the predictive distribution is Gaussian (e.g. as in typical autoregressive models), the HDR will be an ellipsoid defined as the set of all $x$ where $D(x, \hat{x}_t)^2 \le \chi^2_{inv}(1-\alpha)$. Here, $D(x, \hat{x}_t)^2$ is the squared Mahalanobis distance between $x$ and the mean of the predictive distribution at time $t$. $\chi^2_{inv}$ is the inverse chi-squared CDF with degrees of freedom equal to the dimensionality of $x$. So, an observed data point $x_t$ can be identified as a possible anomaly simply by checking the Mahalanobis distance $D(x_t, \hat{x}_t)$.

Note that Euclidean distance is equivalent to Mahalanobis distance in the special case where the variance is equal in all directions. Otherwise, it's a poor choice because it fails to account for the shape of the predictive distribution. For the same reason, the Mahalanobis distance approach won't work well if the predictive distribution is non-Gaussian.

References

Hyndman, R. J. (1996). Computing and graphing highest density regions. The American Statistician, 50(2), 120-126.