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.
Thank you for your time!
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.
 A: I assume the following from your problem, please confirm:


*

*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).

