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.
trend_new = np.arange(t_last +1, t_last + forcast horizon +1)
$\endgroup$