# Intervention With Differencing

When conducting an intervention analysis with time series data (aka Interrupted Time series) as discussed here for example one requirement I have is to estimate the total gain (or loss) due to the intervention - i.e. number of units gained or lost (the Y variable).

Not entirely understanding how to estimate the intervention function using a filter function within R, I went about it in a brute force manner, hoping this is general enough to work in any situation.

Lets say that given the data

 cds<- structure(c(2580L, 2263L, 3679L, 3461L, 3645L, 3716L, 3955L,
3362L, 2637L, 2524L, 2084L, 2031L, 2256L, 2401L, 3253L, 2881L,
2555L, 2585L, 3015L, 2608L, 3676L, 5763L, 4626L, 3848L, 4523L,
4186L, 4070L, 4000L, 3498L), .Dim = c(29L, 1L), .Dimnames = list(
NULL, "CD"), .Tsp = c(2012, 2014.33333333333, 12), class = "ts")


we decide that the best fitting model is as follows, with the intervention function as

$m_t= \frac{\omega_0}{(1-\delta B)}X_t$ where $X_t$ is a pulse at October 2013.

fit4 <- arimax(log(cds), order = c(1,1,0),include.mean=FALSE,
xtransf = data.frame(Oct13 = 1*(seq_along(cds)==22)),
transfer = list(c(1,0))
,xreg=1*(seq_along(cds)==3))
fit4

#    ARIMA(1,1,0)

#    Coefficients:
#              ar1    xreg  Oct13-AR1  Oct13-MA0
#          -0.0184  0.2718     0.4295     0.4392
#    s.e.   0.2124  0.1072     0.3589     0.1485

#    sigma^2 estimated as 0.02176:  log likelihood=13.85
#    AIC=-19.71   AICc=-16.98   BIC=-13.05


I have two questions:

1) Even though we have differenced the ARIMA errors , to assess the intervention function which was then technically fit using the differenced series $\bigtriangledown X_t$ is there anything we need to do in order to "change back" the estimate of $\omega_0$ or $\delta$ from using $\bigtriangledown X_t$ to $X_t$?

2) Is this correct: In order to determine the gain of the intervention, I constructed the intervention $m_t$ from the parameters. Once I have $m_t$ then I compare the fitted values from the model fit4 (exp() to reverse the log) to exp( fitted values minus $m_t$ ) and determine that over the observed period, the intervention resulted in 3342.37 extra units.

Is this process the correct one to determine the gain generally from an intervention analysis?

    int_vect1<-1*(seq_along(cds)==22)
wo<- 0.4392
delta<-0.4295

mt<-rep(0,length(int_vect1))

for (i in 1:length(int_vect1))
{

if (i>1)
{
mt[i]<-wo*int_vect1[i]+delta*mt[i-1]
}

}

mt

sum(exp(fitted(fit4)) - (exp(fitted(fit4) - mt)))

• Wondering if anyone can shed light on the proper way to estimate the impact of an intervention - generally, if the procedure I demonstrated is in fact correct? Aug 12, 2014 at 21:32
• This is a great question. I suspect that the choice of indicator function might not be the best choice depending on the nature of the intervention. Perhaps an exponential decay function, $m_t = \exp(\alpha(t - i)) \text{if} i \ge t \text{else} 0$ for some scaling parameter $\alpha$. Aug 17, 2014 at 15:46

Assuming this is toy example:

1) Even though we have differenced the ARIMA errors , to assess the intervention function which was then technically fit using the differenced series ▽Xt is there anything we need to do in order to "change back" the estimate of ω0 or δ from using ▽Xt to Xt?

When you difference the data, you should difference the response/intervention variables. When you back difference (transform) after you model then this would automatically take care of differencing** I know this is very easy when you use SAS Proc ARIMA. I dont know how to do this R.

Second Question:

2) Is this correct: In order to determine the gain of the intervention, I constructed the intervention mt from the parameters. Once I have mt then I compare the fitted values from the model fit4 (exp() to reverse the log) to exp( fitted values minus mt ) and determine that over the observed period, the intervention resulted in 3342.37 extra units.

To determine, gain in intervention, you need to take exponent and then subtract -1, this would give the proportion or incremental effect. To demonstrate this in your case, see below. For the first month, the impact was 55% of original sales and rapidly decays. Cumulativelt you have 4580 units of incremental effect (Oct 13 thru Feb 2014. (I referred to Forecasting Principle and Applications by Delurgio P: 518. There is an excellent voluminous chapter on intervention analysis).

Someone please correct if this methodology is correct ?

Pulse intervention + decay is clearly not sufficient in this case, I would do a pulse + permanent level shift as shown in the diagram (e) below which is from the classic paper by Box and Tiao.

• Hi @forecaster. How did you get 3170 as the effect? Here is what I did, I looked at the fitted values of the model, which was 8.64245833 (still on the log scale). Then, exp(8.64245833) = 5667.244674. Then, I took 8.64245833 - 0.4392 = 8.20325833. Since exp(8.64245833) - exp(8.20325833) = 2014.411599 this is the effect. exp(8.64245833)/ exp(8.20325833) = 1.55 which looked to me like support for this. Nov 27, 2014 at 23:25
• You used the actuals and the modeled effect I guess, versus my approach which used the model for both. I used the idea of what does the model say with and without the effect. Which is right? Nov 27, 2014 at 23:28
• Hi @B_miner, with logarithmic transformation scale we need to look at rate of change. The approach I outlined is direct approach per the textbook, I quoted. However, your approach is reasonable as well. I'll screenshot the textbook pages in the near future. Nov 27, 2014 at 23:56
• The rate of change was 0.55 which is also the rate of change in the model approach I took. I wonder which approach is more right? I lean towards mine since the approach is based on the model (fitted veruss actuals). If the model is very close to actual, the two approaches will be the sample. I would like to see the pages. I see the book looks out of print? Nov 28, 2014 at 0:26
• yes, the book is out of print. The book example is a permanent change vs. pulse intervention in your example. I think your approach is straight forward and accurate. Nov 28, 2014 at 0:50

@forecaster After allowing AUTOBOX to identify 3 outliers using 29 values (not inappropriate in y experience) a useful model was found and here . The residual acf plot does not suggest an under-specified model . The Actual/Fit/Forecast plot is with Fit/Forecast here . Forecaster had (correctly) previously mentioned how a pulse variable can morph into a level/step variable when a denominator coefficient of nearly 1.0 is introduced. In finding two level shifts (the most recent one starting in 9/2013) and a pulse at 10/2013 , the model presents a clearer picture. In terms of the impact of the pulse at 10/13 it is simply the value of the coefficient. HTH

• Which of the two questions were you answering? Nov 27, 2014 at 19:59
• The first question premised a model that assumed a log transform which I believe is not warranted. The pulse at 10/2013 = 1710 which is the estimate of the effect at 10/2013 Nov 28, 2014 at 2:08
• @B_Miner you could say the level shift at 9/2013 raised things by 1480 thus the net lift at 10/2013 would be 1710+1480=3190 Nov 28, 2014 at 21:36
• This post seems to be more of an extended comment on the question than it is an answer to any part of the question. Perhaps it could be augmented to address the questions directly?
– whuber
Dec 1, 2014 at 22:10
• The false premise of the first question is the crux of my answer : taking logs and incorporating unnecessary differencing is questionable/incorrect in my opinion ... thus my "answer" is in part to correct the premise and to suggest the effect of the impact at 10/2013 is simply the sum of a temporary and permanent change. Even though the OP has accepted another answer I haven't. Dec 1, 2014 at 22:23