# Find the order of an ARMA model (q & p )

I fit an ARMA model in Matlab and before I calculate the predicted value with the prediction error I set the order $(p,q)$ to some random value. But how can I determine the number of AR (p) and MA lags (q)?

• 2 Comments: 1st: I edited your question, I hope I got the meaning right. 2nd: you really take random p and q? In which range? If you don't know how to fit them then what about using simple models first (AR(1), ARMA(1,1), ...)?
– Ric
Dec 10 '14 at 8:22
• A very good approach can be found here: otexts.org/fpp/8/7
– Ric
Dec 10 '14 at 14:41

The Autocorrelation Function (ACF) $\rho_k=Corr(y_t,y_{t-k})$ expresses the strength of linear dependency between the $k$-lagged realizations and hence represents an important tool for identification of the lag orders of ARMA and GARCH processes: $$\rho_k:=Corr(y_t,y_{t-k})=\frac{\gamma_k}{\gamma_0},\,\,k\in\mathbb{Z}$$ where the Autocovariance $\gamma_k$ is defined as $\gamma_k:=Cov(y_t,y_{t-k})=E((y_t-\mu)(y_{t-k}-\mu))$ with $\gamma_0=Var(y_t)$ and $\mu=E(y_t)$. The estimated Sample Autocorrelation Function (SACF) for a sample $\hat{y}=(y_1,...,y_T)$ with mean $\hat{\mu}$ is given by: $$\hat{\rho}_k:=\dfrac{\sum_{i=k+1}^T(y_i-\hat{\mu})(y_{i-k}-\hat{\mu})}{\sum_{i=1}^T (y_i-\hat{\mu})^2},\,\,k\in\mathbb{Z}$$ For weakly stationary processes, autocorrelation is symmetric as $\rho_k=\rho_{-k}$, $\gamma_k=\gamma_{-k}$. For an alternating series between positive and negative values, $\rho_k$ will alternate aswell.

For AR(p), the autocorrelations follow a weighted mean aswell, described by the {Yule Walker equations}: $$\gamma_k=a_1\gamma_{k-1}+a_2\gamma_{k-2}+...+a_p\gamma_{k-p},\,k\in\mathbb{Z}$$ $$\rho_k=a_1\rho_{k-1}+a_2\rho_{k-2}+...+a_p\rho_{k-p},\,k\in\mathbb{Z}$$ It can be shown that an AR(p) process has nonzero autocorrelation up to lag order $p$, and else decay to $0$ exponentially fast with increasing $k$, as the linear dependence of $y_t$ beyond lags $p$ is zero (only the incorporated past innovations for orders $k>p$ remain). Yule Walker equations also allow estimation of the AR(p) coefficients $a_i,i=1,...,p$ from a given sample. However, for non-stationary processes the equations tend to become unstable.

For an MA(q) process, one can show: $$\gamma_k=\left\{\begin{array}{cl} \sigma_\epsilon^2\sum_{i=k}^qb_ib_{i-k}, &k=0,1,...,q\\ 0, & k>q \end{array}\right.$$ such that autocorrelation theoretically cuts off to $0$ directly at order $q$ (without exponential decay).

For ARMA(p,q) processes, one may apply Extended Yule Walker equations to show that ACF after $q$ and $p$ decays exponentially as $k$ increases. The GARCH process orders can be identified same way by ACF of the squared residuals $\epsilon_t^2$.

Empirical sample autocorrelations however often do not cut off or decline to $0$ exponentially, as they are only estimated from real samples. The significant order with $\hat{\rho}_k=0$ is hence determined by confidence intervals. SACF is normally distributed with mean zero such that the 95% confidence interval for $\hat{\rho}_k=0$ is given by $$[\pm1.96\sqrt{Var(\hat{\rho}_k})]$$ with $Var(\hat{\rho}_k)=\frac{1}{T}\left(1+2\sum_{i=1}^T \hat{\rho}_i^2\right),k>T$.

You want to compute the BIC (Bayessian Information Criterion) or the AIC (Akaike information criterion) for different (p,q) pairs.

Here is a wikipedia article with information on how to interpret those criteria in practice.

Here is a mathworks page with detailed instructions on how to perform this task within Matlab.

Keep in mind that in practice and depending on your data you may want to choose low values for the (p,q) parameters to avoid overfitting.

Another possible solution is the EACF of Tsay and Tiao (1984) where the idea is that if the order of the AR process is known the MA can be inferred. The output is a table where the first left corner 0 is taken to be the order of the ARMA(p, q) model.

Use acf and pacf as to determine AR and MA parts. Use the position of last significant value for the two tests as the AR and MA terms respectively. or use autoarima if matlab has one with AIC or BIC coefficients. AIC returns a more general model (all possible values) while BIC results in a more constrained one (simpler).

In addition to previous answer which suggest AIC/BIC and using properties of acf/pacf plots (and their hypothesis tests) -- I would like to add that one further (more advanced) method could be to use cross validation for a time series model.

• what is cross validation!? Dec 10 '14 at 10:34

I go through the following steps when determing the ARMA(p,q) order for my data.

First, you must determine whether any transformation is needed, such as a degree of differencing or taking the log. Remember the data must be stationary to fit an ARMA model. This can be done in a number of ways, but my preference is to run an augmented-dickey fuller test on the data to see if there is a unit root, if so the data needs to be differenced. Additionally, I run an ARMA model with no p/q terms and only a constant and different orders of differencing (start with 0 and 1), and then look at the ACF plot of these models. The model with lowest standard deviation is often the optimal level of differencing. Another concern is, if you want to capture a multiplicative pattern, try taking the log (and then potentially differencing as well). In any of these attempts, if there is a wave pattern in the residuals or another nonrandom pattern, you may want to increase your order of differencing.

Secondly, after selecting your data transformations, look at ACF and PACF plot for the series. The general rule is if the ACF cuts off sharply at lag k and the PACF decays more slowly you have an MA signature; and conversely if the ACF decays more slowly and PACF cuts off you have an AR signature. Remember, if you have 2 or more orders of differencing, the general practice is to not include a constant. If you have one order of differencing, the constant is your average trend, so it is up to you whether you think that is a justifiable assumption. If you still have spikes in either the ACF or PACF at a low-order lag, you need to add another AR or MA term. If your are deciding between two models, all else being equal, the model with the lower order of differencing is generally preferred.

Another option, is to develop an algorithm which tests, on your stationarised series, the model with lowest AIC is considered the best. I don't like doing this in isolation from the the other steps, as it ignores other concerns.

Hope that helps