I have a monthly time series with an intervention and I would like to quantify the effect of this intervention on the outcome. I realize the series is rather short and the effect is not yet concluded.

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")

enter image description here

The methodology

1) The pre-intervention series (up until October 2013) was used with the auto.arima function. The model suggested was ARIMA(1,0,0) with non-zero mean. The ACF plot looked good.

pre <- window(cds, start=c(2012, 01), end=c(2013, 09))

mod.pre <- auto.arima(log(pre))

# Coefficients:
#          ar1  intercept
#       0.5821     7.9652
# s.e.  0.1763     0.0810
# 
# sigma^2 estimated as 0.02709:  log likelihood=7.89
# AIC=-9.77   AICc=-8.36   BIC=-6.64

2) Given the plot of the full series, the pulse response was chosen below, with T = Oct 2013,

enter image description here

which according to cryer and chan can be fit as follows with the arimax function:

mod.arimax <- arimax(log(cds), order=c(1, 0, 0),
                     seasonal=list(order=c(0, 0, 0), frequency=12),
                     include.mean=TRUE,
                     xtransf=data.frame(Oct13=1 * (seq(cds) == 22)),
                     transfer=list(c(1, 1)))
mod.arimax

# Series: log(cds) 
# ARIMA(1,0,0) with non-zero mean 
# 
# Coefficients:
#          ar1  intercept  Oct13-AR1  Oct13-MA0  Oct13-MA1
#       0.7619     8.0345    -0.4429     0.4261     0.3567
# s.e.  0.1206     0.1090     0.3993     0.1340     0.1557
# 
# sigma^2 estimated as 0.02289:  log likelihood=12.71
# AIC=-15.42   AICc=-11.61   BIC=-7.22

The residuals from this appeared OK:

enter image description here

The plot of fitted and actuals:

plot(fitted(mod.arimax), col="red", type="b")
lines(window(log(cds), start=c(2012, 02)), type="b")

enter image description here

The Questions

1) Is this methodology correct for intervention analysis?

2) Can I look at estimate/SE for the components of the transfer function and say that the effect of the intervention was significant?

3) How can one visualize the transfer function effect (plot it?)

4) Is there a way to estimate how much the intervention increased the output after 'x' months? I guess for this (and maybe #3) I am asking how to work with an equation of the model - if this were simple linear regression with dummy variables (for example) I could run scenarios with and without the intervention and measure the impact - but I am just unsure how to work this this type of model.

ADD

Per request, here are the residuals from the two parametrizations.

First from the fit:

fit <- arimax(log(cds), order=c(1, 0, 0),
              xtransf=
              data.frame(Oct13a=1 * (seq_along(cds) == 22),
                         Oct13b=1 * (seq_along(cds) == 22)),
              transfer=list(c(0, 0), c(1, 0)))

plot(resid(fit), type="b")

enter image description here

Then, from this fit

mod.arimax <- arimax(log(cds), order=c(1, 0, 0),
                     seasonal=list(order=c(0, 0, 0), frequency=12),
                     include.mean=TRUE,
                     xtransf=data.frame(Oct13=1 * (seq(cds) == 22)),
                     transfer=list(c(1, 1))) 

mod.arimax
plot(resid(mod.arimax), type="b")

enter image description here

  • Would it be ok if I provide you solution using SAS software ? – forecaster Jul 22 '14 at 14:29
  • Sure, I would be curious if you come up with a better model. – B_Miner Jul 23 '14 at 12:26
  • Ok, the model is little better than originally proposed, but similar to @javlacalle. – forecaster Jul 23 '14 at 21:17
up vote 9 down vote accepted
+200

An AR(1) model with the intervention defined in the equation given in the question can be fitted as shown below. Notice how the argument transfer is defined; you also need one indicator variable in xtransf for each one of the interventions (the pulse and the transitory change):

require(TSA)
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")

fit <- arimax(log(cds), order = c(1, 0, 0), 
              xtransf = data.frame(Oct13a = 1 * (seq_along(cds) == 22), 
                                   Oct13b = 1 * (seq_along(cds) == 22)),
              transfer = list(c(0, 0), c(1, 0)))
fit
# Coefficients:
#          ar1  intercept  Oct13a-MA0  Oct13b-AR1  Oct13b-MA0
#       0.5599     7.9643      0.1251      0.9231      0.4332
# s.e.  0.1563     0.0684      0.1911      0.1146      0.2168
# sigma^2 estimated as 0.02131:  log likelihood = 14.47,  aic = -18.94

You can test the significance of each intervention by looking at the t-statistic of the coefficients $\omega_0$ and $\omega_1$. For convenience, you can use the function coeftest.

require(lmtest)
coeftest(fit)
#            Estimate Std. Error  z value  Pr(>|z|)    
# ar1        0.559855   0.156334   3.5811 0.0003421 ***
# intercept  7.964324   0.068369 116.4896 < 2.2e-16 ***
# Oct13a-MA0 0.125059   0.191067   0.6545 0.5127720    
# Oct13b-AR1 0.923112   0.114581   8.0564 7.858e-16 ***
# Oct13b-MA0 0.433213   0.216835   1.9979 0.0457281 *  
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

In this case the pulse is not significant at the $5\%$ significance level. Its effect may be already captured by the transitory change.

The intervention effect can be quantified as follows:

intv.effect <- 1 * (seq_along(cds) == 22)
intv.effect <- ts(
  intv.effect * 0.1251 + 
  filter(intv.effect, filter = 0.9231, method = "rec", sides = 1) * 0.4332)
intv.effect <- exp(intv.effect)
tsp(intv.effect) <- tsp(cds)

You can plot the effect of the intervention as follows:

plot(100 * (intv.effect - 1), type = "h", main = "Total intervention effect")

Total intervention effect

The effect is relatively persistent because $\omega_2$ is close to $1$ (if $\omega_2$ were equal to $1$ we would observe a permanent level shift).

Numerically, these are the estimated increases quantified at each time point caused by the the intervention in October 2013:

window(100 * (intv.effect - 1), start = c(2013, 10))
#           Jan      Feb      Mar      Apr      May Jun Jul Aug Sep      Oct
# 2013                                                              74.76989
# 2014 40.60004 36.96366 33.69046 30.73844 28.07132                         
#           Nov      Dec
# 2013 49.16560 44.64838

The intervention increases the value of the observed variable in October 2013 by around a $75\%$. In subsequent periods the effect remains but with a decreasing weight.

We could also create the interventions by hand and pass them to stats::arima as external regressors. The interventions are a pulse plus a transitory change with parameter $0.9231$ and can be built as follows.

xreg <- cbind(
  I1 = 1 * (seq_along(cds) == 22), 
  I2 = filter(1 * (seq_along(cds) == 22), filter = 0.9231, method = "rec", 
              sides = 1))
arima(log(cds), order = c(1, 0, 0), xreg = xreg)
# Coefficients:
#          ar1  intercept      I1      I2
#       0.5598     7.9643  0.1251  0.4332
# s.e.  0.1562     0.0671  0.1563  0.1620
# sigma^2 estimated as 0.02131:  log likelihood = 14.47,  aic = -20.94

The same estimates of the coefficients as above are obtained. Here we fixed $\omega_2$ to $0.9231$. The matrix xreg is the kind of dummy variable that you may need to try different scenarios. You could also set different values for $\omega_2$ and compare its effect.

These interventions are equivalent to an additive outlier (AO) and a transitory change (TC) defined in the package tsoutliers. You can use this package to detect these effects as shown in the answer by @forecaster or to build the regressors used before. For example, in this case:

require(tsoutliers)
mo <- outliers(c("AO", "TC"), c(22, 22))
oe <- outliers.effects(mo, length(cds), delta = 0.9231)
arima(log(cds), order = c(1, 0, 0), xreg = oe)
# Coefficients:
#          ar1  intercept    AO22    TC22
#       0.5598     7.9643  0.1251  0.4332
# s.e.  0.1562     0.0671  0.1563  0.1620
# sigma^2 estimated as 0.02131:  log likelihood=14.47
# AIC=-20.94   AICc=-18.33   BIC=-14.1

Edit 1

I've seen that the equation that you gave can be rewritten as:

$$ \frac{(\omega_0 + \omega_1) - \omega_0 \omega_2 B}{1 - \omega_2 B} P_t $$

and it can be specified as you did using transfer=list(c(1, 1)).

As shown below, this parameterization leads, in this case, to parameter estimates that involve a different effect compared to the previous parameterization. It reminds me the effect of an innovational outlier rather than a pulse plus a transitory change.

fit2 <- arimax(log(cds), order=c(1, 0, 0), include.mean = TRUE, 
  xtransf=data.frame(Oct13 = 1 * (seq(cds) == 22)), transfer = list(c(1, 1)))
fit2
# ARIMA(1,0,0) with non-zero mean 
# Coefficients:
#          ar1  intercept  Oct13-AR1  Oct13-MA0  Oct13-MA1
#       0.7619     8.0345    -0.4429     0.4261     0.3567
# s.e.  0.1206     0.1090     0.3993     0.1340     0.1557
# sigma^2 estimated as 0.02289:  log likelihood=12.71
# AIC=-15.42   AICc=-11.61   BIC=-7.22

I'm not very familiar with the notation of package TSA but I think that the effect of the intervention can now be quantified as follows:

intv.effect <- 1 * (seq_along(cds) == 22)
intv.effect <- ts(intv.effect * 0.4261 + 
  filter(intv.effect, filter = -0.4429, method = "rec", sides = 1) * 0.3567)
tsp(intv.effect) <- tsp(cds)
window(100 * (exp(intv.effect) - 1), start = c(2013, 10))
#              Jan         Feb         Mar         Apr         May Jun Jul Aug
# 2014  -3.0514633   1.3820052  -0.6060551   0.2696013  -0.1191747            
#      Sep         Oct         Nov         Dec
# 2013     118.7588947 -14.6135216   7.2476455

plot(100 * (exp(intv.effect) - 1), type = "h", 
  main = "Intervention effect (parameterization 2)")

intervention effect parameterization 2

The effect can be described now as a sharp increase in October 2013 followed by a decrease in the opposite direction; then the effect of the intervention vanishes quickly alternating positive and negative effects of decaying weight.

This effect is somewhat peculiar but may be possible in real data. At this point I would look at the context of your data and the events that may have affected the data. For example, has there been a policy change, marketing campaign, discovery,... that may explain the intervention in October 2013. If so, is it more sensible that this event has an effect on the data as described before or as we found with the initial parameterization?

According to the AIC, the initial model would be preferred because it is lower ($-18.94$ against $-15.42$). The plot of the original series does not suggest a clear match with the sharp changes involved in the measurement of the second intervention variable.

Without knowing the context of the data, I would say that an AR(1) model with a transitory change with parameter $0.9$ would be appropriate to model the data and measure the intervention.

Edit 2

The value of $\omega_2$ determines how fast the effect of the intervention decays to zero, so that's the key parameter in the model. We can inspect this by fitting the model for a range of values of $\omega_2$. Below, the AIC is stored for each of these models.

omegas <- seq(0.5, 1, by = 0.01)
aics <- rep(NA, length(omegas))
for (i in seq(along = omegas)) {
  tc <- filter(1 * (seq_along(cds) == 22), filter = omegas[i], method = "rec", 
               sides = 1)
  tc <- ts(tc, start = start(cds), frequency = frequency(cds))
  fit <- arima(log(cds), order = c(1, 0, 0), xreg = tc)
  aics[i] <- AIC(fit)
}
omegas[which.min(aics)]
# [1] 0.88

plot(omegas, aics, main = "AIC for different values of the TC parameter")

AIC for different values of omega

The lowest AIC is found for $\omega_2 = 0.88$ (in agreement with the value estimated before). This parameter involves a relatively persistent but transitory effect. We can conclude that the effect is temporary since with values higher than $0.9$ the AIC increases (remember that in the limit, $\omega_2=1$, the intervention becomes a permanent level shift).

The intervention should be included in the forecasts. Obtaining forecasts for periods that have already been observed is a helpful exercise to assess the performance of the forecasts. The code below assumes that the series ends in October 2013. Forecasts are then obtained including the intervention with parameter $\omega_2=0.9$.

First we fit the AR(1) model with the intervention as a regressor (with parameter $\omega_2=0.9$):

tc <- filter(1 * (seq.int(length(cds) + 12) == 22), filter = 0.9, method = "rec", 
             sides = 1)
tc <- ts(tc, start = start(cds), frequency = frequency(cds))
fit <- arima(window(log(cds), end = c(2013, 10)), order = c(1, 0, 0), 
             xreg = window(tc, end = c(2013, 10)))

The forecasts can be obtained and displayed as follows:

p <- predict(fit, n.ahead = 19, newxreg = window(tc, start = c(2013, 11)))

plot(cbind(window(cds, end = c(2013, 10)), exp(p$pred)), plot.type = "single", 
     ylab = "", type = "n")
lines(window(cds, end = c(2013, 10)), type = "b")
lines(window(cds, start = c(2013, 10)), col = "gray", lty = 2, type = "b")
lines(exp(p$pred), type = "b", col = "blue")
legend("topleft",
       legend = c("observed before the intervention",
           "observed after the intervention", "forecasts"),
       lty = rep(1, 3), col = c("black", "gray", "blue"), bty = "n")

observed and forecasted values

The first forecasts match relatively well the observed values (gray dotted line). The remaining forecasts show how the series will continue the path to the original mean. The confidence intervals are nonetheless large, reflecting the uncertainty. We should therefore be cautions and revise the model as new data are recorded.

$95\%$ confidence intervals can be added to the previous plot as follows:

lines(exp(p$pred + 1.96 * p$se), lty = 2, col = "red")
lines(exp(p$pred - 1.96 * p$se), lty = 2, col = "red")
  • This is great, thank you! I had a couple follow-ups if you dont mind. 1) Is the process I followed correct? 2) Would you consider the fit of the model "good enough" in order to use the estimates to quantify the effect of the intervention? 3) Should I not be able to use my parametrization, i.e. transfer=list(c(1,1)) as an equivalent and get pretty close results? The example I was following from a textbook suggested I should be able to, but in this example, the results are not close... – B_Miner Jul 23 '14 at 15:15
  • @B_Miner You are right I have edited my answer. – javlacalle Jul 23 '14 at 20:08
  • I agree with you that of the two models, the first parametrization ( maybe with the non-significant pulse removed) would be the best fit. Why the two parametrization do not yield the same model, when I believe they should, is a mystery. I will email the package developer (who also wrote the book that mentions their equivalence). – B_Miner Jul 24 '14 at 0:10
  • The data is the number of certificates of deposits opened per month. The intervention was an increase in average interest rate, that spiked up starting on Oct 13. The level of the interest rate has remained relatively constant since Oct 13. It seemed to me that after the spike, the demand for the product started to taper - I am not sure if it will return to the previous mean or settle at some elevated (from prior) level. – B_Miner Jul 24 '14 at 0:12
  • B_miner, based on the data we cant really conclude, if the demand will settle down to a new mean. – forecaster Jul 24 '14 at 3:22

Sometimes less is more. With 30 observations in hand I submitted the data to AUTOBOX , a piece of software that I have helped develop. I submit the following analysis in the hope of acquiring the +200 reward (just kidding !) . I have plotted the Actual and Cleansed Values a visually suggesting the impact of "recent activity". enter image description here . The model that was automatically developed is shown here. enter image description here and here enter image description here . The residuals from this rather simple level-shifted series are presented here enter image description here . The model statistics are here enter image description here . In summary there were an interventions that could be empirically identified rendering an ARIMA process ; two pulses and 1 level shiftenter image description here . The Actual/Fit and Forecast graph further highlights the analysis.enter image description here

I for one would like to see the plot of the residuals from the previously specified and in my opinion potentially overly-specified models.

  • I am not familiar with Autobox, but is the noise part of the model the same as originally I had : a non-zero mean and an AR(1)? – B_Miner Jul 24 '14 at 0:17
  • Is this output saying that the only "intervention" in the Oct 13 to current time period is a single pulse for Oct 13 and then the series returns to its normal mean level? – B_Miner Jul 24 '14 at 0:21
  • I added the residuals from both parametrizations. To my eye, it seems like the first one I listed (the one originally fit by javlacalle) is better. Agree? – B_Miner Jul 24 '14 at 0:46
  • 1) The noise part is an AR(1) with a non-zero mean – IrishStat Jul 24 '14 at 10:24
  • 1) The noise part is an AR(1) with a non-zero mean ; 2) There are 2 interventions period 22 and period 3 and after oct 13 it returns to a new level that started at sept 13 ; 3 ) Given the choice between the two that you mentioned, I agree BUT I prefer the AUTOBOX model for it's simplicity and efficiency. You can find out more about AUTOBOX from autobox.com/cms – IrishStat Jul 24 '14 at 10:32

Based on my similar post to your earlier question, I used tso function in tsoutliers package in $R$ and it automatically detected a temporary change at October 2013. Please note that temporary change is different from ramp shift in transfer function which is what you are after. I don't think there is a package/function that I'm aware of that would be able to visualize transfer function. Hopefully this would provide some insights. I did not use log transformation, I modeled it directly. tsoutliers package can be thought of as an automatic intervention detection.

Below is the code:

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")
arimatr <- tsoutliers::tso(cds,args.tsmethod=list(d=0,D=0))
plot(arimatr)
arimatr

Below is the estimate, there was a ~2356.3 unit increase in October 2013 with a standard error of ~481.8 and has a decaying effect thereafter. The function automatically identified AR(1). I had to do couple of iteration and make both seasonal and non seasonal differencing to 0, which is reflected in the args.tsmethod in the tso function.

Series: cds 
ARIMA(1,0,0) with non-zero mean 

Coefficients:
         ar1  intercept       TC22
      0.5969  3034.6560  2356.2914
s.e.  0.1495   206.5202   481.7981

sigma^2 estimated as 209494:  log likelihood=-219.03
AIC=446.06   AICc=447.73   BIC=451.53

Outliers:
  type ind    time coefhat tstat
1   TC  22 2013:10    2356 4.891

Below is the plot, tsoutlier is the only package that I know of that can print temporary changes this nicely in a plot.

enter image description here

This analysis hopefully provided answer to your 2, 3 and 4 questions albeit using a different methdeology. Especially the plot and the coefficients provided the effect of this intervention and what would have happened if you did not have this intervention.

Also hoping someone else can replicate this plot/analysis using transfer function modeling in R. I'm not sure if this could be done in R, may be someone else can fact check me on this.

  • Thanks. Yeah, this plot is what I would like to do from the arimax model - look at with and without the intervention (and subtract). I think the filter function in R can be used to generate the transfer function value for each month (and then just plot it to visualize) but I cant figure out how to do it for an arbitrary pulse intervention function. – B_Miner Jul 18 '14 at 13:18

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

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