Robust time-series regression for outlier detection

I'm navigating my way through the plethora of regression models to find some form of standardized residuals that could help "score" the observations in proportion to their "outlyingness" for the purpose of anomaly detection.

The data is a (weekly) seasonal time-series of (over-dispersed) hourly call-counts, thus advocating a need for GLMs. The predictor (baseline) is the historical median for the same day,hour and it is observed that it closely follows the current day,hour call-counts. Here are the results of OLS and ANOVA:

Call:
lm(formula = observed ~ baseline, data = temp)

Residuals:
Min      1Q  Median      3Q     Max
-7183.7  -184.9    -2.9   273.2  5514.9

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 90.558071  40.358140   2.244    0.025 *
baseline     1.009935   0.005434 185.838   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 814.4 on 1338 degrees of freedom
Multiple R-squared:  0.9627,    Adjusted R-squared:  0.9627
F-statistic: 3.454e+04 on 1 and 1338 DF,  p-value: < 2.2e-16

Analysis of Variance Table

Response: observed
Df     Sum Sq    Mean Sq F value    Pr(>F)
baseline     1 2.2907e+10 2.2907e+10   34536 < 2.2e-16 ***
Residuals 1338 8.8748e+08 6.6329e+05
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Example of overdispersed counts:

10793.0
10277.0
8686.0
6189.0
4008.0
2398.0
1553.0
755.0
489.0
375.0
480.0
987.0
2423.0
4694.0
6950.0
8220.0
8904.0
9070.0
9313.0
9778.0
10129.0
10804.0
10529.0
10022.0
9838.0


You can find a sample of the data here, these are hourly counts starting from 2016-03-31 16:00:00 Pacific Time adjusted for time-zones. If you can't view it on google drive, please download the file.

Unsurprisingly, the errors are zero mean, but non-gaussian and heteroscedastic (variance increasing with predictor and fitted values).

It is also observed that the residuals are autocorrelated and the PACF suggests an AR(1) model on the residuals (rho~0.6), which is expected in case of time-series data.

I dont think seasonality plays any role since we already adjust for it using medians of the same day, hour.

The following constraints are thus at hand:

1. Over-dispersed counts data suggest the need for Negative Binomial Regression or a suitable transformation to the data, perhaps semi-parametric regression and an estimator for the conditional variance?
2. The presence of outliers indicate using robust regression methods.
3. Auto-correlation in the residuals suggest using an AR(1) model, eg. a cochrane-orkutt procedure to adjust estimates.
4. "Scores" could be pearson, deviance, anscombe residuals or perhaps outlier statistics such as influence etc.

I'd like to know if my assessment of the available options is correct, or am I missing something? It would be very helpful if someone can share their wisdom from past experiences with such data as well.

PS: I'd prefer more "general" methods that could work with potentially different distributions of time-series (non-parametric methods) and are computationally fast as there are millions of time-series.

EDIT: Thanks to @IrishStat for his useful comments and analysis, I learnt about the requisite tools and their (un)availability in the current open-source stack. Although, i'm restricted to open-source and had to take a slight detour with whatever tools that were available, nevertheless, his analysis using AUTOBOX (from which I borrowed many ingredients) was insightful of what would be a good approach for modeling such data.

• please post an example of an hourly series that is "over-dispersed" as I have never heard that term before. Use an excel file (single column) format. demand-planning.com/2010/03/18/… might be of interest to you Oct 6, 2016 at 15:42
• @IrishStat: By overdispersion I meant that variance>mean, the observed variance, mean for the hourly time-series are respectively (17770553.883233838, 6348.05). The same numbers for the daily time-series are (856976633.1585213, 149234.8596491228). I have also added a plot for better visualization. Oct 6, 2016 at 15:58
• @IrishStat: added a data-sample showing over-dispersion Oct 6, 2016 at 16:17
• We see data like this quite frequently . Some series have error variance proportional to the expected value requiring power transformations. Some series have error variance that changes at certain points over time while some series have error variance changing as a result of day effects or hourly effects.. As a quick guess yours looks like the latter where not only does the expected value change based on the day of the week but the error variance also changes. Oct 6, 2016 at 16:45
• @IrishStat Yes, thats in accordance with my observations too. So, in your experience, how would you go about doing regression analysis with the objective of outlier "scoring" on this data? Oct 7, 2016 at 16:27

I took your 90 days of data (24 hourly readings per day) and analyzed it using AUTOBOX a piece of software that I have helped develop using a 28 day forecast horizon. The documentation for the approach can be found in the User Guide available from the AFS website. I will try and give you you a general overview here. The data is analyzed in a parent-to-child approach where a model is initially developed for the daily totals and here incorporating memory,daily effects and anomalies, level shifts , local time trends etc.. This "parent model" leads to forecasts and confidence limits based upon possible different daily error variances which also differed by hour ( and the potential for anomalies to arise in the future using simulation procedures. The next step is to identify 24 causal models using the parent as a possible predictor using memory as needed , level shifts as needed while identifying and remedying possible anomalies/level shifts/local time trends. As an example of this let me show you the graphical output for hour 0,6,9,15,and 18 . As an example consider the model for hour 4 . GROUPT, the daily total series is significant along with three days of the week (2,3,and 4). Additionally 5 unusual values (outliers) were found at periods 89,55,19,78 and 2. There is a level shift (permanent change upwards in the mean value) at period 20 and a partial daily effect (day 3) starting at period 19 . Now we have 24 sets of forecasts for the children for the next 28 days and a set of forecasts for the parent for the next 28 days. We reconcile these two to obtain the final forecasts for 24 hours for the 28 day forecast. The reconciliation can be done in a parent-to-child or a child-to-parent manner. Following are snapshots of both procedures. . Clearly 90 days of data insufficient to capture , weekly effects, holiday effects , specific days of the month effects , long-weekend effects , week-in-month effects, monthly effects etc.. but if you had a longer series you can get the picture as to what might be possible. The detailed output can be made available to you or any interested party by contacting me as it is to voluminous to post. With this analysis one might want to creatively program the approach with free R tools but there are a lot of pitfalls awaiting such enterprise. Hope this helps your research into what I think is a very important statistical problem regarding detecting exceptional events and accounting for their effect in the forecast horizon.

• This is amazing1 Oct 14, 2016 at 12:14
• seemingly much more advanced methods than I had thought of ! Such data is ubiquitous nowadays, for instance count of tweets etc. and research on detecting unusual events indeed holds a lot of promise. For data size, I only had data 90 days available, however there are several such timeseries (counts of categorical sub-attributes) which follow similar patterns and could therefore be tackled similarly. I would have to dig into the documentation in order to get more detail for my own set of analyses, but this serves as a good starting point, thanks for your help. Oct 14, 2016 at 12:29