Using a function of intervention time as a covariate in arimax

Question

I am trying to estimate the impact of an intervention on a time series. The problem is that the intervention effect looks curve-linear (i.e. increasing, plateauing, and then decreasing), and I am unsure how to model it.

Here is a sample of my data with 75 time periods on either side of the intervention (a shorter sample won't show the problem):

        series <- structure(c(28.57, 31.71, 35.57, 33.57, 33.71, 33.43, 35.57, 
29.57, 33.71, 32.86, 32.43, 30.43, 32, 31.14, 34, 34.29, 33.29, 
32.71, 32.43, 33.29, 27.71, 28, 32.71, 31.71, 29.71, 33.29, 33.14, 
33.71, 33.57, 31, 31.43, 32.71, 30.29, 31.86, 33.14, 37, 31.43, 
32.14, 35.29, 42.43, 35.57, 31.57, 27.86, 31.43, 24, 23.86, 28.57, 
38, 33.71, 33.86, 34.29, 37.86, 35, 37.43, 31.86, 35.86, 35.86, 
33.86, 33.71, 34, 40.86, 37.43, 28, 38, 35.43, 34.43, 38.29, 
32.57, 34.29, 38.86, 35.29, 37.86, 31.14, 35.43, 30.86, 28.57, 
34.14, 38.29, 36, 33.14, 37.29, 39.71, 40.14, 45.43, 39.14, 36.86, 
44.86, 38.14, 40.14, 46.43, 46.57, 44.71, 46.71, 42.86, 42.29, 
43.71, 45.86, 34.71, 41.14, 52.43, 47.14, 44.57, 43.57, 37, 46.14, 
49.29, 49.86, 46.71, 45.86, 49.14, 43.86, 46.71, 35.29, 43.43, 
41.29, 40.86, 42.14, 44.43, 36.57, 40.29, 40.71, 40, 39.57, 37.43, 
36.86, 37.43, 39.86, 33.71, 33.57, 36, 38.43, 43.14, 47.14, 38.71, 
42.57, 41.29, 40.43, 37, 43.71, 42.43, 38.14, 41.86, 36.14, 42.43, 
47.43, 46.57, 42.43, 38.86, 35.14, 38.57, 31.43), .Tsp = c(-75, 
75, 1), class = "ts")

The intervention occurs at $t=0$. If we plot the series, it looks like this:

To the naked eye, there seems to be a a gradual increase in the series that begins right after the intervention and starts dropping off after about 40 periods. After 75 periods it appears to be back to normal.

I've previously estimated effects like this using transfer functions (see e.g. Cryer and Chan (2008): "Time Series Analysis"), but I can't find any form that fits this curve-linear pattern. The series ends 75 periods after the intervention, so a transfer function with a step input fits quite well (the ARIMA(1,1,1)-specification has been determined previously):

    intervention <- structure(list(step = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)), .Names = "step", row.names = c(NA, 
-151L), class = "data.frame")

require(TSA)

arimax(series,order=c(1,1,1),
       xtransf=intervention, transfer=list(c(1,0)))

Call:
arimax(x = series, order = c(1, 1, 1), xtransf = intervention, transfer = list(c(1, 
    0)))

Coefficients:
         ar1      ma1  step-AR1  step-MA0
      0.3087  -0.9272    0.8526    1.4303
s.e.  0.0965   0.0461    0.0689    0.6254

sigma^2 estimated as 12.8:  log likelihood = -404.74,  aic = 817.47

The problem is that a permanent level shift is very implausible/impossible for this series (the contents of which are confidential, unfortunately). Estimating that it never returns to its pre-intervention level doesn't make any sense.

So I've tried to model the parabolic shape of the intervention effect by including a linear and a quadratic post-intervention term as exogenous regressors. Like so:

xreg_var <- structure(list(linear = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 
18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 
34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 
50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 
66, 67, 68, 69, 70, 71, 72, 73, 74, 75), quadratic = c(0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 4, 9, 16, 25, 36, 49, 64, 
81, 100, 121, 144, 169, 196, 225, 256, 289, 324, 361, 400, 441, 
484, 529, 576, 625, 676, 729, 784, 841, 900, 961, 1024, 1089, 
1156, 1225, 1296, 1369, 1444, 1521, 1600, 1681, 1764, 1849, 1936, 
2025, 2116, 2209, 2304, 2401, 2500, 2601, 2704, 2809, 2916, 3025, 
3136, 3249, 3364, 3481, 3600, 3721, 3844, 3969, 4096, 4225, 4356, 
4489, 4624, 4761, 4900, 5041, 5184, 5329, 5476, 5625)), .Names = c("linear", 
"quadratic"), row.names = c(NA, -151L), class = "data.frame")

require(forecast)

Arima(series,order=c(1,1,1),
      xreg = xreg_var)

Series: series 
ARIMA(1,1,1)                    

Coefficients:
         ar1      ma1  linear  quadratic
      0.3181  -0.9055  0.4361    -0.0055
s.e.  0.1179   0.0701  0.1470     0.0018

sigma^2 estimated as 13.45:  log likelihood=-406.29
AIC=822.57   AICc=822.99   BIC=837.63

The fit is similar and this parabolic model makes way more sense theoretically. However, I've never seen a function of the intervention time used as a covariate in an arimax model like this before.

Is it ok to do?

Best,

Bertel

Answer 1

[...] this parabolic model makes way more sense theoretically. However, I've never seen a function of the intervention time used as a covariate in an arimax model like this before. Is it ok to do?

I think it is OK to do. Time series modelling is often about approximating complicated processes in parsimonious ways. If you want to describe a pattern (a reaction to an impulse) in a simple way and you have in mind a functional form that seems to match to what you see, why not use it.

Once you have fitted the model, you can assess whether it captures the pattern well. If model residuals still contain substantial patterns in them, you might reconsider the functional form. But that will become clear once you have fitted the model; it should not discourage you from fitting it in the first place, though.