# Estimating same model over multiple time series

I have a novice background in time series (some ARIMA estimation/forecasting) and am facing a problem I don't fully understand. Any help would be greatly appreciated.

I am analyzing multiple time series, all over the same time interval and all of the same frequency, all describing a similar type of data. Each series is just one variable, there are no other corresponding predictors that I'm looking at.

I have been asked to estimate a single model that describes ALL the series - for example, imagine I could find one ARIMA (p,d,q) with the same orders, coefficients, etc. that could fit all the series. My supervisor does not want me to separately estimate each series, nor does he want me to do some sort of VAR model with dependencies between the series.

My question is: what would I even call such a model, and how might I go about estimating / forecasting it? If it's easier for you to use code examples, I speak both SAS and R.

• I'm trying to do the same thing. Apparently, there is something called a 'Multivariate AutoRegressive' model out there. I've found reference to it, but not how to use it. Based on the linked paper, I presume it has been implemented in R. journal.r-project.org/archive/2012-1/…
– Mox
Dec 19 '14 at 6:34
• The standard approach is vector autoregression and there is an R package called var. Dec 19 '14 at 7:43
• Is a vector autoregession different than autoregression on panel data? Or is it a matter of different fields, different names? The plm package was suggested for panel data. cran.r-project.org/web/packages/plm/vignettes/plm.pdf clidyn.ethz.ch/papers/arfit.pdf
– Mox
Dec 19 '14 at 20:43

You could do a grid search: start with ARIMA(1,0,0) and try all the possibilities up to ARIMA(5,2,5) or something. Fit the model to each series, and estimate a scale-independent error measurement like MAPE or MASE (MASE would probably be better). Choose the ARIMA model with the lowest average MASE across all your models.

You could improve this procedure by cross-validating your error measurement for each series, and also by comparing your results to a naive forecast.

It might be a good idea to ask why you're looking for a single model to describe all of the series. Unless they're generated by the same process, this doesn't seem like a good idea.

• Thank you - I will try this. I agree that this doesn't seem like the best idea. The argument I got was that each series does not have enough observations (~28) for a good estimation and that it would be more robust to estimate over all the series. I'm not sure I agree with that argument. Feb 17 '12 at 20:42

One way to do that is to construct a long time series with all of your data, and with sequences of missing values between the series to separate them. For example, in R, if you have three series (x, y and z) each of length 100 and frequency 12, you can join them as follows

combined <- ts(c(x,rep(NA,56),y,rep(NA,56),z,rep(NA,56)),frequency=12)


Notice that the number of missing values is chosen to ensure the seasonal period is retained. I've padded out the final year with 8 missing values and then added four missing years (48 values) before the next series. That should be enough to ensure any serial correlations wash out between series.

Then you can use auto.arima() to find the best model:

library(forecast)
fit <- auto.arima(combined)


Finally, you can apply the combined model to each series separately in order to obtain forecasts:

fit.x <- Arima(x,model=fit)
fit.y <- Arima(y,model=fit)
fit.z <- Arima(z,model=fit)

• +1, neat trick. Judging by the OP's comment to another answer, I was going to suggest some sort of panel data model, but this is way better. Feb 20 '12 at 11:47
• This is a huge help, much appreciated. Can you explain more re: the number of missing values is chosen to ensure the seasonal period is retained? Sorry I didn't fully follow - thanks. Feb 23 '12 at 23:26
• If there is seasonality in the data (as there often is with monthly observations), you want the long series to still have the Januaries a multiple of 12 apart, the Februaries a multiple of 12 apart, and so on. Then, when the model is being chosen, the seasonality can be modelled appropriately. Feb 24 '12 at 0:32
• To add to this trick - you could add external regressors (xreg) indicating category membership. This would account for separate means for the different series, while still keeping other coefficients in common. Dec 18 '14 at 12:49

Estimating single model for multiple time series is the realm of panel data econometrics. However in your case with no explanatory variable @Rob Hyndman answer is probably the best fit. However if it turns out that the means of time series are different (test it, since in this case @Rob Hyndman's method should fail!), but ARMA structure is the same, then you will have to use Arellano-Bond type estimator. The model in that case would be:

$$y_{it}=\alpha_i+\rho_1 y_{i,t-1}+...+\rho_p y_{i,t-p}+\varepsilon_{it}$$

where $$i$$ indicates different time series and $$\varepsilon_{it}$$ can have the same covariance structure across all $$i$$.

• Really appreciate your solution and the others as well. You mention that: However if it turns out that the means of time series are different (test it, since in this case @Rob Hyndman's method should fail!) Can you explain more about why this is? Thanks. Feb 23 '12 at 23:25
• @sparc_spread, suppose it's just two series. One is centered at about 0 with variance 1 and the other is centered at 1000 with variance 1. Then If both series are fit using the same coefficients, that means we're constraining alpha_1 = alpha_2, so the predictions for both series would be around 500, horribly off. Basically, treating all series as belonging to the same model may require some recentering/normalization prior to fitting the joint model. Mar 26 '17 at 11:18

An alternative to Rob Hyndman's approach, to make a single data series, is to merge the data. This might be appropriate if your multiple time series represent noisy readings from a set of machines recording the same event. (If each time series is on a different scale you need to normalize the data first.)

NOTE: you still only end up with 28 readings, just less noise, so this may not be appropriate for your situation.

t1=xts(jitter(sin(1:28/10),amount=0.2),as.Date("2012-01-01")+1:28)
t2=xts(jitter(sin(1:28/10),amount=0.2),as.Date("2012-01-01")+1:28)
t3=(t1+t2)/2


• This will only work if all the signals are similar in nature (e.g. periodic) are all in phase - in your example if two of the sine waves were 180 degrees out of phase they would completely cancel!
– tdc
Feb 24 '12 at 9:43
• Yes, averaging your data is only appropriate when each date series is supposed to represent the same data, and (you are happy to assume that) they only differ in that each has different noise. Feb 26 '12 at 1:34

I would look at hidden Markov models and dynamic Bayesian networks. They model time series data. Also they are trained using multiple time series instances e.g. multiple blood pressure time series from various individuals . You should find packages in Python and R to build those. You might have to define structure for these models.