I am working with a large amount of time series. These time series are basically network measurements coming every 10 minutes, and some of them are periodic (i.e. the bandwidth), while some other aren't (i.e. the amount of routing traffic).

I would like a simple algorithm for doing an online "outlier detection". Basically, I want to keep in memory (or on disk) the whole historical data for each time series, and I want to detect any outlier in a live scenario (each time a new sample is captured). What is the best way to achieve these results?

I'm currently using a moving average in order to remove some noise, but then what next? Simple things like standard deviation, mad, ... against the whole data set doesn't work well (I can't assume the time series are stationary), and I would like something more "accurate", ideally a black box like:

double outlier_detection(double* vector, double value);

where vector is the array of double containing the historical data, and the return value is the anomaly score for the new sample "value" .

  • 1
    Just for clarity, here's the original question on SO: stackoverflow.com/questions/3390458/… – Matt Parker Aug 2 '10 at 21:42
  • 1
    I think we should encourage posters to post links as part of the question if they have posted the same question at another SE site. – user28 Aug 2 '10 at 21:47
  • yes, you're completely right. Next time I'll mention that the message is crossposted. – gianluca Aug 2 '10 at 21:53
  • I also suggest you check out the other Related links on the right hand side of the page. This is a popular question and it has come up in a variety of questions previously. If they aren't satisfactory it is best to update your question about the specifics of your situation. – Andy W Sep 3 '12 at 15:10
  • Good catch, @Andy! Let's merge this question with the other one. – whuber Sep 3 '12 at 16:08

14 Answers 14

Here is a simple R function that will find time series outliers (and optionally show them in a plot). It will handle seasonal and non-seasonal time series. The basic idea is to find robust estimates of the trend and seasonal components and subtract them. Then find outliers in the residuals. The test for residual outliers is the same as for the standard boxplot -- points greater than 1.5IQR above or below the upper and lower quartiles are assumed outliers. The number of IQRs above/below these thresholds is returned as an outlier "score". So the score can be any positive number, and will be zero for non-outliers.

I realise you are not implementing this in R, but I often find an R function a good place to start. Then the task is to translate this into whatever language is required.

tsoutliers <- function(x,plot=FALSE)
    x <- as.ts(x)
        resid <- stl(x,s.window="periodic",robust=TRUE)$time.series[,3]
        tt <- 1:length(x)
        resid <- residuals(loess(x ~ tt))
    resid.q <- quantile(resid,prob=c(0.25,0.75))
    iqr <- diff(resid.q)
    limits <- resid.q + 1.5*iqr*c(-1,1)
    score <- abs(pmin((resid-limits[1])/iqr,0) + pmax((resid - limits[2])/iqr,0))
        x2 <- ts(rep(NA,length(x)))
        x2[score>0] <- x[score>0]
        tsp(x2) <- tsp(x)
  • +1 from me, excellent. So > 1.5 X inter-quartile range is the consensus definition of an outlier for time-dependent series? That would be nice to have a scale-independent reference. – doug Aug 3 '10 at 3:06
  • The outlier test is on the residuals, so hopefully the time-dependence is small. I don't know about a consensus, but boxplots are often used for outlier detection and seem to work reasonably well. There are better methods if someone wanted to make the function a little fancier. – Rob Hyndman Aug 3 '10 at 3:45
  • Really thank you for your help, I really appreciate. I'm quite busy at work now, but I'm going to test an approach like yours as soon as possible, and I will come back with my final considerations about this issue. One only thought: in your function, from what I see, I have to manually specify the frequency of the time series (when constructing it), and the seasonality component is considered only when the frequency is greater than 1. Is there a robust way to deal with this automatically? – gianluca Aug 3 '10 at 15:59
  • Yes, I have assumed the frequency is known and specified. There are methods to estimate the frequency automatically, but that would complicate the function considerably. If you need to estimate the frequency, try asking a separate question about it -- and I'll probably provide an answer! But it needs more space than I have available in a comment. – Rob Hyndman Aug 3 '10 at 23:40
  • 2
    @Marcin, I recommend taking a stab at it yourself. Maybe paste your solution at gist.github.com and post a SO question when you're done, to have others check your work? – Ken Williams Jan 24 '13 at 16:24

A good solution will have several ingredients, including:

  • Use a resistant, moving window smooth to remove nonstationarity.

  • Re-express the original data so that the residuals with respect to the smooth are approximately symmetrically distributed. Given the nature of your data, it's likely that their square roots or logarithms would give symmetric residuals.

  • Apply control chart methods, or at least control chart thinking, to the residuals.

As far as that last one goes, control chart thinking shows that "conventional" thresholds like 2 SD or 1.5 times the IQR beyond the quartiles work poorly because they trigger too many false out-of-control signals. People usually use 3 SD in control chart work, whence 2.5 (or even 3) times the IQR beyond the quartiles would be a good starting point.

I have more or less outlined the nature of Rob Hyndman's solution while adding to it two major points: the potential need to re-express the data and the wisdom of being more conservative in signaling an outlier. I'm not sure that Loess is good for an online detector, though, because it doesn't work well at the endpoints. You might instead use something as simple as a moving median filter (as in Tukey's resistant smoothing). If outliers don't come in bursts, you can use a narrow window (5 data points, perhaps, which will break down only with a burst of 3 or more outliers within a group of 5).

Once you have performed the analysis to determine a good re-expression of the data, it's unlikely you'll need to change the re-expression. Therefore, your online detector really only needs to reference the most recent values (the latest window) because it won't use the earlier data at all. If you have really long time series you could go further to analyze autocorrelation and seasonality (such as recurring daily or weekly fluctuations) to improve the procedure.

  • 2
    This is an extraordinary answer for practical analysis. Never would have thought needed to try 3 IQR beyond the quartiles. – John Robertson Oct 9 '12 at 18:45
  • 2
    @John, 1.5 IQR is Tukey's original recommendation for the longest whiskers on a boxplot and 3 IQR is his recommendation for marking points as "far outliers" (a riff on a popular 60's phrase). This is built into many boxplot algorithms. The recommendation is analyzed theoretically in Hoaglin, Mosteller, & Tukey, Understanding Robust and Exploratory Data Analysis. – whuber Oct 9 '12 at 21:38
  • This confirms time series data I have been trying to analyse. Window average and also a window standard deviations. ((x - avg) / sd) > 3 seem to be the points I want to flag as outliers. Well at least warn as outliers, I flag anything higher than 10 sd as extreme error outliers. The problem I run into is what is an ideal window length? I'm playing with anything between 4-8 data points. – Josh Peak Jun 29 '16 at 8:00
  • 1
    @Neo Your best bet may be to experiment with a subset of your data and confirm your conclusions with tests on the remainder. You could conduct a more formal cross-validation, too (but special care is needed with time series data due to the interdependence of all the values). – whuber Jun 29 '16 at 12:10

(This answer responded to a duplicate (now closed) question at Detecting outstanding events, which presented some data in graphical form.)

Outlier detection depends on the nature of the data and on what you are willing to assume about them. General-purpose methods rely on robust statistics. The spirit of this approach is to characterize the bulk of the data in a way that is not influenced by any outliers and then point to any individual values that do not fit within that characterization.

Because this is a time series, it adds the complication of needing to (re)detect outliers on an ongoing basis. If this is to be done as the series unfolds, then we are allowed only to use older data for the detection, not future data! Moreover, as protection against the many repeated tests, we would want to use a method that has a very low false positive rate.

These considerations suggest running a simple, robust moving window outlier test over the data. There are many possibilities, but one simple, easily understood and easily implemented one is based on a running MAD: median absolute deviation from the median. This is a strongly robust measure of variation within the data, akin to a standard deviation. An outlying peak would be several MADs or more greater than the median.

There is still some tuning to be done: how much of a deviation from the bulk of the data should be considered outlying and how far back in time should one look? Let's leave these as parameters for experimentation. Here's an R implementation applied to data $x = (1,2,\ldots,n)$ (with $n=1150$ to emulate the data) with corresponding values $y$:

# Parameters to tune to the circumstances:
window <- 30
threshold <- 5

# An upper threshold ("ut") calculation based on the MAD:
library(zoo) # rollapply()
ut <- function(x) {m = median(x); median(x) + threshold * median(abs(x - m))}
z <- rollapply(zoo(y), window, ut, align="right")
z <- c(rep(z[1], window-1), z) # Use z[1] throughout the initial period
outliers <- y > z

# Graph the data, show the ut() cutoffs, and mark the outliers:
plot(x, y, type="l", lwd=2, col="#E00000", ylim=c(0, 20000))
lines(x, z, col="Gray")
points(x[outliers], y[outliers], pch=19)

Applied to a dataset like the red curve illustrated in the question, it produces this result:


The data are shown in red, the 30-day window of median+5*MAD thresholds in gray, and the outliers--which are simply those data values above the gray curve--in black.

(The threshold can only be computed beginning at the end of the initial window. For all data within this initial window, the first threshold is used: that's why the gray curve is flat between x=0 and x=30.)

The effects of changing the parameters are (a) increasing the value of window will tend to smooth out the gray curve and (b) increasing threshold will raise the gray curve. Knowing this, one can take an initial segment of the data and quickly identify values of the parameters that best segregate the outlying peaks from the rest of the data. Apply these parameter values to checking the rest of the data. If a plot shows the method is worsening over time, that means the nature of the data are changing and the parameters might need re-tuning.

Notice how little this method assumes about the data: they do not have to be normally distributed; they do not need to exhibit any periodicity; they don't even have to be non-negative. All it assumes is that the data behave in reasonably similar ways over time and that the outlying peaks are visibly higher than the rest of the data.

If anyone would like to experiment (or compare some other solution to the one offered here), here is the code I used to produce data like those shown in the question.

n.length <- 1150
cycle.a <- 11
cycle.b <- 365/12
amp.a <- 800
amp.b <- 8000

x <- 1:n.length
baseline <- (1/2) * amp.a * (1 + sin(x * 2*pi / cycle.a)) * rgamma(n.length, 40, scale=1/40)
peaks <- rbinom(n.length, 1,  exp(2*(-1 + sin(((1 + x/2)^(1/5) / (1 + n.length/2)^(1/5))*x * 2*pi / cycle.b))*cycle.b))
y <- peaks * rgamma(n.length, 20, scale=amp.b/20) + baseline
  • This is a really interesting solution and I appreciate that I can implement it without using R (just using plain JavaScript in a web application). Thanks! – hgoebl Oct 10 '15 at 7:45

If you're worried about assumptions with any particular approach, one approach is to train a number of learners on different signals, then use ensemble methods and aggregate over the "votes" from your learners to make the outlier classification.

BTW, this might be worth reading or skimming since it references a few approaches to the problem.

I am guessing sophisticated time series model will not work for you because of the time it takes to detect outliers using this methodology. Therefore, here is a workaround:

  1. First establish a baseline 'normal' traffic patterns for a year based on manual analysis of historical data which accounts for time of the day, weekday vs weekend, month of the year etc.

  2. Use this baseline along with some simple mechanism (e.g., moving average suggested by Carlos) to detect outliers.

You may also want to review the statistical process control literature for some ideas.

  • 1
    Yeah, this is exactly what I am doing: until now I manually split the signal into periods, so that for each of them I can define a confidence interval within which the signal is supposed to be stationary, and therefore I can use standard methods such as standard deviation, ... The real problem is that I can not decide the expected pattern for all the signals I have to analyze, and that's why I'm looking for something more intelligent. – gianluca Aug 2 '10 at 21:37
  • Here is a one idea: Step 1: Implement and estimate a generic time series model on a one time basis based on historical data. This can be done offline. Step 2: Use the resulting model to detect outliers. Step 3: At some frequency (perhaps every month?), re-calibrate the time series model (this can be done offline) so that your step 2 detection of outliers does not go too much out of step with current traffic patterns. Would that work for your context? – user28 Aug 2 '10 at 22:24
  • Yes, this might work. I was thinking about a similar approach (recomputing the baseline every week, which can be CPU intensive if you have hundreds of univariate time series to analyze). BTW the real difficult question is "what is the best blackbox-style algorithm for modeling a completely generic signal, considering noise, trend estimation and seasonality?". AFAIK, every approach in literature requires a really hard "parameter tuning" phase, and the only one automatic method I found is an ARIMA model by Hyndman (robjhyndman.com/software/forecast). Am I missing something? – gianluca Aug 2 '10 at 22:38
  • Please keep in mind I'm not too lazy for investigating these parameters, the point is that these values need to be set according to the expected pattern of the signal, and in my scenario I can't make any assumption. – gianluca Aug 2 '10 at 22:40
  • ARIMA models are classic time series models that can be used to fit time series data. I would encourage you to explore the application of ARIMA models. You could wait for Rob to be online and perhaps he will chime in with some ideas. – user28 Aug 2 '10 at 22:44

Seasonally adjust the data such that a normal day looks closer to flat. You could take today's 5:00pm sample and subtract or divide out the average of the previous 30 days at 5:00pm. Then look past N standard deviations (measured using pre-adjusted data) for outliers. This could be done separately for weekly and daily "seasons."

  • Again, this works pretty well if the signal is supposed to have a seasonality like that, but if I use a completely different time series (i.e. the average TCP round trip time over time), this method will not work (since it would be better to handle that one with a simple global mean and standard deviation using a sliding window containing historical data). – gianluca Aug 2 '10 at 22:02
  • 1
    Unless you are willing to implement a general time series model (which brings in its cons in terms of latency etc) I am pessimistic that you will find a general implementation which at the same time is simple enough to work for all sorts of time series. – user28 Aug 2 '10 at 22:06
  • Another comment: I know a good answer might be "so you might estimate the periodicity of the signal, and decide the algorithm to use according to it", but I didn't find a real good solution to this other problem (I played a bit with spectral analysis using DFT and time analysis using the autocorrelation function, but my time series contain a lot of noise and such methods give some crazy results mosts of the time) – gianluca Aug 2 '10 at 22:06
  • A comment to your last comment: that's why I'm looking for a more generic approach, but I need a kind of "black box" because I can't make any assumption about the analyzed signal, and therefore I can't create the "best parameter set for the learning algorithm". – gianluca Aug 2 '10 at 22:09
  • @gianluca As you have intimated the underlying ARIMA structure can mask the anomaly. Incorrect formulation pf possible cause variables such as hour of the day, day-of-the-week, holiday effects etc can also mask the anomaly. The answer is fairly clear you need to have a good eqaution to effectively detect anomalies. To quote Bacon , "For whoever knows the ways of Nature will more easily notice her deviations and, on the other hand, whoever knows her deviations will more accurately describe her ways." – IrishStat Sep 6 '12 at 17:03

An alternative to the approach outlined by Rob Hyndman would be to use Holt-Winters Forecasting . The confidence bands derived from Holt-Winters can be used to detect outliers. Here is a paper that describes how to use Holt-Winters for "Aberrant Behavior Detection in Time Series for Network Monitoring". An implementation for RRDTool can be found here.

Spectral analysis detects periodicity in stationary time series. The frequency domain approach based on spectral density estimation is an approach I would recommend as your first step.

If for certain periods irregularity means a much higher peak than is typical for that period then the series with such irregularities would not be stationary and spectral anlsysis would not be appropriate. But assuming you have identified the period that has the irregularities you should be able to determine approximately what the normal peak height would be and then can set a threshold at some level above that average to designate the irregular cases.

  • 2
    Could you explain how this solution would detect "local irregularities"? Presenting a worked example would be extremely helpful. (To be honest, I'm suggesting you do this because in carrying out such an exercise I believe you will discover that your suggestion is not effective for outlier detection. But I could be wrong...) – whuber Sep 3 '12 at 15:03
  • 1
    @whuber The spectral analysis will only identify where all the peaks are. The next step would be to fit a yime series model using sine and cosine terms with the frequencies determined from the spectral analysis and the amplitudes estimated from the data. If irregularities mean peaks with very high amplitudes then I think a threshold on the amplitude would be appropriate. If local irregularities means that for a period the amplitude sometimes is significantly larger than other then the series is not stationary and spectral analysis would not be approriate. – Michael Chernick Sep 3 '12 at 15:35
  • 1
    I don't follow the conclusion about lack of stationarity. For instance, the sum of a regular sinusoidal waveform and a marked Poisson point process would be stationary, but it would not exhibit any of the periodicity you seek. You would nevertheless find some strong peaks in the periodogram, but they would tell you nothing relevant to the irregular data peaks introduced by the Poisson process component. – whuber Sep 3 '12 at 16:07
  • 1
    A stationary time series has a constant mean. If the peak for a periodic component can change over time it can cuase the mean to change over time and hence the seires would be nonstationary. – Michael Chernick Sep 3 '12 at 16:16

Since it is a time series data, a simple exponential filter http://en.wikipedia.org/wiki/Exponential_smoothing will smoothen the data. It is a very good filter since you don't need to accumulate old data points. Compare every newly smoothed data value with its unsmoothed value. Once the deviation exceeds a certain predefined threshold (depending on what you believe an outlier in your data is), then your outlier can be easily detected.

In C I will do the following for a real-time 16 bit sample (I believe this is found somewhere here < Explanation - https://dsp.stackexchange.com/questions/378/what-is-the-best-first-order-iir-approximation-to-a-moving-average-filter>)

#define BITS2 2     //< This is roughly = log2( 1 / alpha ), depending on how smooth you want your data to be

short Simple_Exp_Filter(int new_sample) 
{static int filtered_sample = 0;
long local_sample = sample << 16; /*We assume it is a 16 bit sample */
filtered_sample += (local_sample - filtered_sample) >> BITS2;   
return (short) ((filtered_sample+0x8000) >> 16); //< Round by adding .5 and truncating.   

int main()
newly_arrived = function_receive_new_sample();
filtered_sample = Simple_Exp_Filter(newly_arrived);
if (abs(newly_arrived - filtered_sample)/newly_arrived > THRESHOLD)
 return 0;   

You could use the standard deviation of the last N measurements (you have to pick a suitable N). A good anomaly score would be how many standard deviations a measurement is from the moving average.

  • Thank you for your response, but what if the signal exhibits a high seasonality (i.e. a lot of network measurements are characterized by a daily and weekly pattern at the same time, for example night vs day or weekend vs working days)? An approach based on standard deviation will not work in that case. – gianluca Aug 2 '10 at 20:57
  • For example, if I get a new sample every 10 minutes, and I'm doing an outlier detection of the network bandwidth usage of a company, basically at 6pm this measure will fall down (this is an expected an totaly normal pattern), and a standard deviation computed over a sliding window will fail (because it will trigger an alert for sure). At the same time, if the measure falls down at 4pm (deviating from the usual baseline), this is a real outlier. – gianluca Aug 2 '10 at 20:58

what I do is group the measurements by hour and day of week and compare standard deviations of that. Still doesn't correct for things like holidays and summer/winter seasonality but its correct most of the time.

The downside is that you really need to collect a year or so of data to have enough so that stddev starts making sense.

  • Thank you, that's exactly what I was trying to avoid (having a lot of samples as baseline), because I would like a really reactive approach (e.g. online detection, maybe "dirty", after 1-2 weeks of baseline) – gianluca Aug 10 '10 at 14:49

I suggest the scheme below, which should be implementable in a day or so:


  • Collect as many samples as you can hold in memory
  • Remove obvious outliers using the standard deviation for each attribute
  • Calculate and store the correlation matrix and also the mean of each attribute
  • Calculate and store the Mahalanobis distances of all your samples

Calculating "outlierness":

For the single sample of which you want to know its "outlierness":

  • Retrieve the means, covariance matrix and Mahalanobis distances from training
  • Calculate the Mahalanobis distance "d" for your sample
  • Return the percentile in which "d" falls (using the Mahalanobis distances from training)

That will be your outlier score: 100% is an extreme outlier.

PS. In calculating the Mahalanobis distance, use the correlation matrix, not the covariance matrix. This is more robust if the sample measurements vary in unit and number.

For the case where one has to compute the outliers quickly, one could use the idea of Rob Hyndman and Mahito Sugiyama ( https://github.com/BorgwardtLab/sampling-outlier-detection , library(spoutlier), function qsp ) to compute the outliers as follows:

rapidtsoutliers <- function(x,plot=FALSE,seed=123)
    x <- as.numeric(x)
    tt <- 1:length(x)
    qspscore <- qsp(x)
    limit <- quantile(qspscore,prob=c(0.95))
    score <- pmax((qspscore - limit),0)
        x2 <- ts(rep(NA,length(x)))
        x2[score>0] <- x[score>0]
        tsp(x2) <- tsp(x)

anomaly detection requires the construction of an equation which describes the expectation. Intervention Detection is available in both a non-causal and causal setting . If one has a predictor series like price then things can get a little complicated. Other responses here don't seem to take into account assignable cause attributable to user specified predictor series like price and thus might be flawed. Quantity sold may well depend on price , perhaps previous prices and perhaps quantity sold in the past. The basis for the anomaly detection ( pulses,seasonal pulses, level shifts and local time trends ) is found in https://pdfs.semanticscholar.org/09c4/ba8dd3cc88289caf18d71e8985bdd11ad21c.pdf

  • The link is not working, could you please fix it. Thanks – Pankaj Joshi Aug 7 at 7:02
  • done .................. – IrishStat Aug 7 at 10:51

Your Answer


By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.