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