Intervention analysis in time-series regression with seasonal ARIMA errors

Intervention analysis in Box-Jenkins framework crosspoinds to time-series regression with arma errors if the noise is stationary or arima errors if the noise is non-stationary.

For a seasonal time series data with increasing trend, the noise model can be express as

$$N_t = \frac{\Theta(B)}{(1-B)(1-B^{12})\Phi(B)} \eta_t$$

If there is a step $S_t$ (0 before intervention and 1 after intervention) and a pulse $P_t$ (1 at intervention and 0 elsewhere) interventions, the model then can be expressed as

$$Y_t=\beta_1S_t+\beta_2P_t+\frac{\Theta(B)}{(1-B)(1-B^{12})\Phi(B)} \eta_t$$

Also because there may different responses to the interventions, say graduate change in level is by $\frac{\omega S_t}{1-\delta B}$ or decayed responses $\frac{\omega P_t}{1-\delta B}$.

$$Y_t=\frac{\omega S_t}{1-\delta B}+\frac{\omega P_t}{1-\delta B}+\frac{\Theta(B)}{(1-B)(1-B^{12})\Phi(B)} \eta_t$$

Therefore my question is:

if the data is seasonal time series, then in the practice, does it mean we need to perform difference $(1-B)(1-B^{12})S_t$ and $(1-B)(1-B^{12})P_t$ along with $(1-B)(1-B^{12})Y_t$ anyways when consider those interventions?

Thanks and Regards

-

3 Answers

The differencing implied by the denominator of your error term must be applied to $Y_t$, $S_t$ and $P_t$. That is, your model is equivalent to $$\nabla\nabla_{12}Y_t=\frac{\omega \nabla\nabla_{12}S_t}{1-\delta B}+\frac{\omega \nabla\nabla_{12}P_t}{1-\delta B}+\frac{\Theta(B)}{\Phi(B)} \eta_t,$$ where $\nabla\nabla_{12} = (1-B)(1-B^{12})$. This is a transfer function model with ARMA errors which is how it would actually be estimated.

If you intended that the pulse and step apply to the differenced $Y$ series, then you need to doubly integrate $S$ and $P$ in the model (as suggested by @IrishStat). That is $$Y_t=\frac{\omega S_t}{\nabla\nabla_{12}(1-\delta B)}+\frac{\omega P_t}{\nabla\nabla_{12}(1-\delta B)}+\frac{\Theta(B)}{\nabla\nabla_{12}\Phi(B)} \eta_t.$$

-
Thanks Rob. One follow up question about such a model in R. In the article of "The ARIMAX model muddle", you mentioned that "The arima() function in R (and Arima() and auto.arima() from the forecast package) fits a regression with ARIMA errors. Unfortunately there is a bug in the code so that if a non-stationary error is specified, the data are not differenced before estimation of the regression coefficients." But does R do the estimation of the regression as (1-B)(1-B^12)[Y(t)-X(t)]= (Theta/Phi)*Eta(t) – Fred Nov 9 '11 at 21:04
For example set.seed(1) v = rnorm(100,1,1) x = cumsum(v) t <- 1:length(x) lm(x~t-1)returns the coefficient 1.117. if in the arima model, arima(x, order=c(1,1,0),xreg=1:length(x)), the regression coefficient is 1.117 too. From the model specification above, it actually can work like that, so in R differences by our own seems unecessary. – Fred Nov 9 '11 at 21:18
First, arima() does not fit a transfer function model of the kind in this question. You would need the arimax function from the TSA package for that. The example you give should give an error because the differenced x variable is a constant, and the arima model already has a constant. The fact that R happily estimates the model is a problem. In any case, the coefficients are only equal to the first three decimal places. – Rob Hyndman Nov 10 '11 at 4:30
In arima(), if there is a difference specified in c(p,1,q), then arima() will suppress the constant term. However, if your arima model is drifting, like in the example above, then arima(x, order=c(1,1,0)) will return the wrong answer which forces the series goes through zero, unless it is performed as arima(diff(x),order=c(1,0,0)) which returns a constant term, or arima(x,order=c(1,1,0),xreg=1:length(x)), but the interpretaton of constant term is different. – Fred Nov 10 '11 at 4:50
However if we consider xreg=1:length(x) as an input series for the transfer function model, then xreg=1:length(x) is no longer necessary in arima(), because it is a constant, where arima(diff(x),order=c(1,0,0)) will return a constant for us. Therefore doesn't it mean arima(diff(x),order=c(1,0,0)) and arima(x,order=c(1,1,0),xreg=1:length(x)) are identical methods? – Fred Nov 10 '11 at 4:55
show 1 more comment

If you wish to estimate the model that you specified you would specify regular and seasonal differencing on Y and provide the two doubly integrated intervention series. It appears you are doing intervention modelling and not intervention detection prior to intervention modelling. The differencing operator in the noise component will essentially return your two indicators to the required status.

-
 Yes, it is intervention modelling. Because the modelling procedure in the software will be performed on differenced Y series, and differenced S(t) and P(t) series. Does it mean even the intervention is gradual change or decayed responses like in the second equation above, the differenced series will still capture the features? – Fred Nov 8 '11 at 20:56 Fred: Specify the doubly integrated series SSt= St*{[1-B][1-B**12]}**-1 and similarly create PPt from Pt. AS you said the coefficients will capture the instantaneous effect and the gradual decay effect – IrishStat Nov 8 '11 at 21:32

I don't know enough to totally parse your question, but I believe the usual practice is to first do a regression with indicator variables $S_t$ and $P_t$, then do your ARIMA on the residuals.

-
 If the errors are non-stationary, that won't give consistent results. If the errors are stationary, it won't give the most efficient estimates. – Rob Hyndman Nov 9 '11 at 9:49 @Rob: good comment! I'll leave my answer (even though it had a negative rating) because I think the comment would be quite useful to other more-naive folks such as myself. – Wayne Nov 9 '11 at 15:47