Lets say I have an ARMA with an AR(1,4) and an MA(1).

Update To be more specific the model to be estimated is:

$$y_t = \rho_1 y_{t-1} + \rho_4 y_{t-4} + \epsilon_t + \theta_1 \epsilon_{t-1}$$

test <- arima.sim(n=1000, list(ar=c(0.3,0,0,0.5), ma=c(0.3)))

How can this data be estimated? For now I have tried:

arima(test, c(4,0,1)

But the results are obviously not correct:

## ============================
##                 Model 1     
## ----------------------------
## ar1                 0.26 ***
##                    (0.06)   
## ar2                 0.18 ** 
##                    (0.06)   
## ar3                -0.27 ***
##                    (0.05)   
## ar4                 0.63 ***
##                    (0.04)   
## ma1                 0.93 ***
##                    (0.06)   
## intercept          -0.21    
##                    (0.32)   
## ----------------------------
## AIC              2878.68    
## BIC              2913.03    
## Log Likelihood  -1432.34    
## Q(1mon)             0.00    
## Q(6mon)             0.00    
## Q(12mon)            0.00    
## Q(24mon)            0.00    
## ============================
## *** p < 0.001, ** p < 0.01, * p < 0.05
  • 2
    $\begingroup$ This question appears to be off-topic because it is mainly about howto do somethingh in R $\endgroup$ – kjetil b halvorsen Dec 5 '14 at 17:52
  • 1
    $\begingroup$ Can you clarify if you are asking for R code, or if you are asking about the underlying issue? $\endgroup$ – gung - Reinstate Monica Dec 5 '14 at 17:58
  • $\begingroup$ I'm asking for R code. $\endgroup$ – DJJ Dec 5 '14 at 19:56
  • 2
    $\begingroup$ While the question is asking for code, a misunderstanding seems to be the cause. $y_t = \rho_1 y_{t-1} + \rho_4 y_{t-4} + \epsilon_t + \theta_2 \epsilon_{t-1}$ is an ARMA(4,1) and your code is simulating an ARMA(4,2), as in lanenok's answer. $\endgroup$ – Affine Dec 5 '14 at 20:19

As far as I see from your line arima.sim(n=1000,list(ar=c(0.3,0,0,0.5),ma=c(1,0.3)))
you are simulating ARIMA(4,0,2) process. Thus I would try to estimate it with arima(test,c(4,0,2)). By the way, what do you mean by AR(1,4) in your question? Misprint?

To answer your comment: your analytical formula is ARIMA(4,0,1) and your $\theta_2$ is conventionally denoted as $\theta_1$. In this case, I would say, you simulate with


and estimate with

arima(test,c(4,0,1), fixed=c(NA,0,0,NA,NA,NA))    
  • $\begingroup$ I hope have made my question more precise. Yes I'm simulating an ARIMA(4,0,2) but that's what all that I've found for now to simulate the above model. $\endgroup$ – DJJ Dec 5 '14 at 17:09

As already appointed by others, your process is an ARMA(4,1) and not an AR(1,4) (I think actually this notation does not even exist!). What I'd suggest is to use the usual Box and Jenkins approach and look at the ACF and PACF. It should give you two possible scenarios:

  1. significant lags in the ACF up to the fourth one (included): if this is the case (and assuming, as in reality, that you do not know the true data generating process) you may want to try to fit an ARMA(4,1) or an AR(4). OLS estimation does not suffer from redundant regressors added in your specification. Indeed you still get unbiased estimates both for the non-null coeffcients $\rho_{i}=0$, $i \in \{1,4\}$ and the null coefficients $\rho_{i}=0$, $i \in \{2,3\}$. Of course you will lose efficiency, i.e. higher standard errors, due to the loss of degrees of freedom;
  2. first and fourth lag in the ACF significant, while the second and the third one not: in this case, according to the principle of parsimony, one would be tempted to go directly for an AR(1) or an ARMA(1,1), assuming that the fourth lag being significant is just a first-type error. One has to be careful here (and in general when one has to deal with non-consecutive lags), as the fourth lag being significant could be a signal for the presence of seasonality. Accordingly, first look at your data and at the unit of time and then, if it is reasonable try to fit an appropriate specification.

To be precise, your model should not present seasonality, as in general model with seasonality could appear in the following form (just an example, think about quarterly data for sales) \begin{equation} y_{t} = \rho_{4}y_{t-4} + \rho_{8}y_{t-8} + \varepsilon_{t} \end{equation}

Coming up to R, if you want to go for the first approach what you have done is correct. If instead you want to try something more "hand-crafted" as solution 2, you can try to regress it manually, using OLS formulas or the function lm.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.