How to Handle Many Times Series Simultaneously? I have a data set including the demand of several products (1200 products) for 25 periods and I need to predict the demand of each product for the next period. At first, I wanted to use ARIMA and train a model for each product, but because of the number of products and tunning of (p,d,q) parameters, it is so time-consuming and it is not practical. Is it recommended to use a regression where previous demands are independent variables (Autoregressive)?
Can I know if there is any method to train a single model for the demand prediction of all 1200 products? I would be thankful if you can suggest any library in Python because I am using Python.
 A: The problem with the mass-fitting packages that have been suggested is they uniformly fail to deal with latent deterministic structure such as pulses, level/step shifts, seasonal pulses and time trends or efficiently deal with user-suggested causals as per https://autobox.com/pdfs/SARMAX.pdf
Additionally the compute time can be a serious complication. AUTOBOX ( which I helped to develop) has a very sophisticated model building phase which archives models and a very quick forecasting option that reuses previously developed model reducing the forecasting time to a small fraction of the rigorous model development time while adjusting the new forecast for recent data observed after the model had been developed and stored. This was implemented for Annheuser-Busch's 600,000 store forecast project for some 50+ items taking into account Price and Weather .
Models can be updated in a rolling fashion , replacing prior models as needed.
No need for parametric restrictions OR omitting the simultaneous effect of causal variables as in VAR and VARIMA while solely relying on only the past of all series a la ARIMA . 
There is no need to have just 1 model with 1 set of parameters as models can and should be tailored/optimized to the individual series.
Unfortunately there is no Python solution yet but hope springs eternal.
A: 1200 products is the main driver of the dimensionality of your problem. Now you have only 25 periods. This is very little data, insufficient to do any kind of blanket correlation analysis. In other words you don't have data to have a simultaneous forecast of all products without reducing the dimensionality. This pretty much eliminates all VARMA and other nice theoretical models. It's impossible to deal with the coefficients of these models, there's too many of them to estimate.
Consider a simple correlation analysis. You'd need (1200x1200 + 1200)/2 cells in the covariance/correlation matrix. You have only 25 data points. The matrix will be rank defficient to enormous degree. What are you going to do? Broadly you have two simple approaches: separate forecasts and factor model. 
The first approach is obvious: you run each product independently. The variation is to group them by some feature, e.g. sector such as "mens closing".
The second approach is to represent the product demand as $d_i=\sum_jF_{j}\beta_{ji}+e_i$, where $F_j$ is a factor. What are the factors? These could be exogenous factors such as GDP growth rate. Or they could be exogenous factors , e.g. those you obtained with PCA analysis. 
If it's an exogenous factor, then you'd need to obtain betas by regressing the series on these factors individually. 
For PCA, you could do a robust PCA and get first few factors with their weights which are you betas. 
Next, you analyze the factors, and build a forecasting model to produce $\hat F_j$ and plug them back to your model to obtain forecast of product demand. You could run a time series model for each factor, even a vector model such as VARMA for several factors. Now, that the dimensionality of the problem was reduced, ou may have enough data to build time series forecasting.
A: As Ben mentioned, the text book methods for multiple time series are VAR and VARIMA models. In practice though, I have not seen them used that often in the context of demand forecasting. 
Much more common, including what my team currently uses, is hierarchical forecasting (see here as well). Hierarchical forecasting is used whenever we have groups of similar time series: Sales history for groups of similar or related products, tourist data for cities grouped by geographical region, etc...
The idea is to have a hierarchical listing of your different products and then do forecasting both at the base level (i.e. for each individual time series) and at aggregate levels defined by your product hierarchy (See attached graphic). You then reconcile the forecasts at the different levels (using Top Down, Botton Up, Optimal Reconciliation, etc...) depending on the business objectives and the desired forecasting targets. Note that you won't be fitting one large multivariate model in this case, but multiple models at different nodes in your hierarchy, which are then reconciled using your chosen reconciliation method. 

The advantage of this approach is that by grouping similar time series together, you can take advantage of the correlations and similarities between them to find patterns (such a seasonal variations) that might be difficult to spot with a single time series. Since you will be generating a large number of forecasts that is impossible to tune manually, you will need to automate your time series forecasting procedure, but that is not too difficult - see here for details. 
A more advanced, but similar in spirit, approach is used by Amazon and Uber, where one large RNN/LSTM Neural Network is trained on all of the time series at one. It is similar in spirit to hierarchical forecasting because it also tries to learn patterns from similarities and correlations between related time series. It is different from hierarchical forecasting because it tries to learn the relationships between the time series itself, as opposed to have this relationship predetermined and fixed prior to doing the forecasting. In this case, you no longer have to deal with automated forecast generating, since you are tuning only one model, but since the model is a very complex one, the tuning procedure is no longer a simple AIC/BIC minimization task, and you need to look at more advanced hyper-parameter tuning procedures, such as Bayesian Optimization.  
See this response (and comments) for additional details. 
For Python packages, PyAF is available but nor very popular. Most people use the HTS package in R, for which there is a lot more community support. For LSTM based approaches, there is Amazon's DeepAR and MQRNN models which are part of a service you have to pay for. Several people have also implemented LSTM for demand forecasting using Keras, you can look those up. 
A: Generally when you have multiple time-series you would use some kind of vector-based model to model them all simultaneously.  The natural extension of the ARIMA model for this purpose is the VARIMA (Vector ARIMA) model.  The fact that you have $1200$ time-series means that you will need to specify some heavy parametric restrictions on the cross-correlation terms in the model, since you will not be able to deal with free parameters for every pair of time-series variables.
I would suggest starting with some simple vector-based model (e.g., VAR, VMA, VARMA) with low degree, and some simple parameter restrictions for cross-correlation.  See if you can find a reasonable model that incorporates cross-correlation to at least one degree of lag, and then go from there.  This exercise will require reading up on vector-based time-series models.  The MTS package and the bigtime pacakage in R has some capabilities for dealing with multivariate time-series, so it would also be worth familiarising yourself with these packages.
A: I am not sure if you are interested in cloud-based solutions, but Amazon makes an algorithm they call "DeepAR" available through AWS SageMaker, as seen here.
This algorithm is specifically intended to be able to learn from multiple input time series in order to create forecasts, including static and dynamic features; as seen in this excerpt from the above linked page:

The training input for the DeepAR algorithm is one or, preferably, more target time series that have been generated by the same process or similar processes. Based on this input dataset, the algorithm trains a model that learns an approximation of this process/processes and uses it to predict how the target time series evolves. Each target time series can be optionally associated with a vector of static (time-independent) categorical features provided by the cat field and a vector of dynamic (time-dependent) time series provided by the dynamic_feat field.

Unfortunately, as far as I can tell, they do not make this algorithm available for offline/self-hosted use.
