# How to deal with large number of times series?

I have a large number of times series (think of the order of 10k), one data point per day (think of the order of 100 days), so I have about 1M data points in total. (One of the things that I will try is a vector autoregressive model, but for now I wanted to start with something simpler. Simple regression where the predictors aren't the series itself, but something else).

One predictor that I have found works quite well is the day-of-week, since the time series turn out to be highly seasonal along a week. This is a category with 7 values.

One way to organize this data is to say that each series is a different value in the category named "series" (a cateogry with 10k different values). So so far I have the "series" category and the "day-of-week" category.

My first question is generally what is a good multiple regression model for this?

I'm using python's statsmodel (but you can contribute ideas in whichever language you prefer). In statsmodel syntax, a model that I've tried is "y ~ C(series) + C(dayofweek)". This model has 7+10k+intercept variables. Doesn't work very well because according to this model the only way in which the series differ is a by a constant shift in the Y axis, which isn't great.

The next-to-simplest model is "y ~ C(series) * C(dayofweek)". This is a model with 7 * 10k + intercept number of variables. The problem is that already statsmodel can't fit this. If I pick a subset of 500 series out of the 10k, statsmodel manages to fit this and the result isn't bad at all. But for the 10k series for some reason it just can't handle it.

What should I do? Intuitively I think that something like https://stats.stackexchange.com/a/230678/188779 should do the trick. The time series have other features associated, and it I could start by grouping similar series together, say in clusters of size 20, that's 10k series divided by 20 equals 500 which is something that I already know I can fit. And after that's done I could try other features still to distinguish individual series within their cluster of size 20.

My second question is how do I implement the above idea? I believe this is called a hierarchical model and the answer linked above what they say is that hierarchical models are a good way to deal with categories with very high cardinality, which is something that makes sense. But how exactly do I implement this? What is the formula for a hierarchical model? I'm looking for either a classical statistics formula, or a statsmodel formula, either would be a step in the right direction.

EDIT (more information): I'm adding some extra information about the series because the appropriateness of a response will depend on the characteristics of the data. All series are very similar. They are all strongly seasonal (the weekends being very different from the weekdays). By doing a seasonal decomposition you see that the trend either extremely weak for some series, or non-existent for others. Seasonal is a decomposition is a good idea since the residuals look healthy. By eyeballing it seems that the biggest difference between series is the trend level and the size of the seasons (ie difference between max and min of seasons).

EDIT (response to https://stats.stackexchange.com/a/319234/188779) Thank you! So after reading your response and considering the above edit, the simple predictors that you suggested that I look for could be "trend level" and "season size"? What do you think? From knowing the data I think this would do quite well. Responding to the second part of your answer: my data is hierarchical (think country > county > city clusters) but from knowing the data I don't think this will be very predictive. But thanks for mentioning ANOVA, certainly a lot to play with there.

EDIT (again response to https://stats.stackexchange.com/a/319234/188779) I've just tried your suggestion (using trend as a predictor) and the results are excellent. To be clear, this model has only 8 parameters (one intercept, plus 7-1=6 parameters for the seasonal categorical encoding, plus trend) and fits 10k series over 100 days with R^2=0.92, and all parameters with t-test p-value=0, overall F-test p-value also =0. So I have a follow up question: this is all in-sample; now how do I predict values in the future? If you give me a date in the future I don't know the trend that date... Should I now be looking into ARIMA models of the trend? Vector auto-regressive models of the trend? Thank you.

• To the last part: for forecasting with deterministic variables like trend or season dummies you can just use the extension of the existing series, e.g. for integer trend trend_new = np.arange(t_last +1, t_last + forcast horizon +1) – Josef Dec 17 '17 at 21:42
• @user333700 Hey, thanks for the suggestion. What you described is about the simplest that can be done. What would be your next suggestion? – oneloop Dec 18 '17 at 13:28

• You want to predict the next day (t+1) for each of the 10K series.
• In your model the variable C(series) is a vector of 10K values where there are 10K-1 zeros and only one 1 for the series you are trying to predict.

In that case, I don't see the advantage of creating a huge vector with ~N*10K values for the coefficients (and associated training matrix), as they will be highly independent. You can try to fit them separately as the biggest correlation comes from the day-of-the-week feature.

Question 1 Assuming the above, clustering could work to reduce the number of predictors. Instead of having 10K predictors, if the series are statistically similar you can come up with a clever clustering scheme that would allow you to reduce the number of predictors to M << 10K. You should run a performance analysis to come up with the best value of M (where M == 10K is the best you can achieve). The features you use for clustering are also very important, you can start with first order statistical estimators (sample mean, sample variance, number of counts) and increase to more complex estimators (e.g. sample skewness).

Question 2 Clustering is not an implementation for a hierarchical model. It will treat all the series the same way. From your description one cannot say that the data is structured in hierarchical layers so you can use a HLM. If that is the case, you should describe better your series data and how the 10K series relate to each other.

Summing up:

• If the 10K series do not have a hierarchical structure you can use independent predictors instead of a massive one as the correlation comes only from the day of the week.
• If the 10K series do not have a hierarchical structure, to reduce the number of predictors and depending on the nature of the data itself you can use a clustering algorithm to preprocess the data and use M << 10K predictors instead. You will pay a "quantization noise" on it, but it can deliver acceptable results.
• If the 10K series do have a nested structure you can consider using an HLM. You should first break-down the dataset into the levels you want to model and then you can use statsmodel's ANOVA (http://www.statsmodels.org/devel/anova.html).
• Thank you so much for taking the time. I won't mark as correct yet to see if other people throw in ideas ;-) But give it a couple of days and I'll mark it as correct. Thank you! – oneloop Dec 17 '17 at 13:18
• Made another edit, please have a look :-) – oneloop Dec 17 '17 at 19:06
• Hi @oneloop, regarging your question on prediction in the 3rd EDIT: when you fit your model in the predictors (your 8 coefficients, including intercept) you should be able to predict any value (that applies to all the times series). If you want to test out-of-sample values I can suggest a classical train/test split (say 85%, 15%), don't show the test data to the fit and try to see how it predicts those values. – drublackberry Dec 18 '17 at 9:39
• Hi, thanks again. Yes! For a few days in the future this model works great out of sample. I told you I only have 100 days but in reality I have years of data. Over 100 days a trend plus a weekly season is great, but over longer periods trends move. Should I linearly extrapolate trends? Should I try an ARIMA model just on the trend component? I'm asking because even though I've seen ARIMA models applied to the data, I don't think I've seen them applied just to the trend part. Is this statistically bad practice for some reason? Or is it ok as long as the results are ok? – oneloop Dec 18 '17 at 13:19
• From what you say maybe the model of trend+weekday can work over short periods (100 days) but cannot fit well over longer periods (years). You might want to add more intelligence (e.g. year seasonality). Wrt ARIMA over trend: you can fit anything with ARIMA as long as it follows the definition of an ARIMA process, try to see if the trend residuals follow it (by looking into its autocovariance) and see what it comes back. But I would go for option#1 and add more intelligence to the model. – drublackberry Dec 18 '17 at 14:24