Función de transferencia de intervención ARIMA: cómo visualizar el efecto

Tengo una serie de tiempo mensual con una intervención y me gustaría cuantificar el efecto de esta intervención en el resultado. Me doy cuenta de que la serie es bastante corta y el efecto aún no está concluido.

Los datos

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

La metodología

1) La serie previa a la intervención (hasta octubre de 2013) se utilizó con la auto.arimafunción. El modelo sugerido fue ARIMA (1,0,0) con una media distinta de cero. La trama de ACF se veía bien.

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) Dada la trama de la serie completa, la respuesta de pulso se eligió a continuación, con T = Oct 2013,

que según cryer y chan se pueden ajustar de la siguiente manera con la función arimax:

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

Los residuos de esto parecían estar bien:

La trama de ajustados y reales:

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

Las preguntas

1) ¿Es esta metodología correcta para el análisis de intervención?

2) ¿Puedo mirar la estimación / SE para los componentes de la función de transferencia y decir que el efecto de la intervención fue significativo?

3) ¿Cómo se puede visualizar el efecto de la función de transferencia (graficarlo?)

4) ¿Hay alguna manera de estimar cuánto aumentó la intervención la producción después de 'x' meses? Supongo que para esto (y tal vez # 3), estoy preguntando cómo trabajar con una ecuación del modelo; si se tratara de una regresión lineal simple con variables ficticias (por ejemplo), podría ejecutar escenarios con y sin la intervención y medir el impacto. pero no estoy seguro de cómo trabajar este tipo de modelo.

AÑADIR

Por solicitud, aquí están los residuos de las dos parametrizaciones.

Primero del ajuste:

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

Entonces, de este ajuste

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

Un modelo AR (1) con la intervención definida en la ecuación dada en la pregunta puede ajustarse como se muestra a continuación. Observe cómo transferse define el argumento ; También necesita una variable indicadora xtransfpara cada una de las intervenciones (el pulso y el cambio transitorio):

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

$\omega_0$$\omega_1$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

$5\%$

El efecto de la intervención se puede cuantificar de la siguiente manera:

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)

Puede trazar el efecto de la intervención de la siguiente manera:

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

$\omega_2$$1$$\omega_2$$1$

Numéricamente, estos son los aumentos estimados cuantificados en cada punto de tiempo causado por la intervención en octubre de 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

$75\%$

stats::arima$0.9231$

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

$\omega_2$$0.9231$xreg$\omega_2$

Estas intervenciones son equivalentes a un valor atípico aditivo (AO) y un cambio transitorio (TC) definido en el paquete tsoutliers. Puede usar este paquete para detectar estos efectos como se muestra en la respuesta de @forecaster o para construir los regresores utilizados anteriormente. Por ejemplo, en este caso:

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

Editar 1

He visto que la ecuación que diste puede reescribirse como:

y se puede especificar como lo hizo usando transfer=list(c(1, 1)).

Como se muestra a continuación, esta parametrización conduce, en este caso, a estimaciones de parámetros que implican un efecto diferente en comparación con la parametrización previa. Me recuerda el efecto de un valor atípico innovador en lugar de un pulso más un cambio transitorio.

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

No estoy muy familiarizado con la notación de paquete, TSApero creo que el efecto de la intervención ahora se puede cuantificar de la siguiente manera:

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

El efecto se puede describir ahora como un fuerte aumento en octubre de 2013 seguido de una disminución en la dirección opuesta; entonces el efecto de la intervención se desvanece rápidamente alternando los efectos positivos y negativos del peso en descomposición.

Este efecto es algo peculiar, pero puede ser posible en datos reales. En este punto, miraría el contexto de sus datos y los eventos que pueden haber afectado los datos. Por ejemplo, ha habido un cambio de política, campaña de marketing, descubrimiento, ... que pueda explicar la intervención en octubre de 2013. Si es así, ¿es más sensato que este evento tenga un efecto en los datos como se describió anteriormente o como encontramos? con la parametrización inicial?

$-18.94$$-15.42$

$0.9$

Editar 2

$\omega_2$$\omega_2$

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

$\omega_2 = 0.88$$0.9$$\omega_2=1$

$\omega_2=0.9$

$\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)))

Los pronósticos se pueden obtener y mostrar de la siguiente manera:

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

Los primeros pronósticos coinciden relativamente bien con los valores observados (línea punteada gris). Los pronósticos restantes muestran cómo la serie continuará el camino hacia la media original. Sin embargo, los intervalos de confianza son grandes, lo que refleja la incertidumbre. Por lo tanto, debemos ser precavidos y revisar el modelo a medida que se graben nuevos datos.

$95\%$

lines(exp(p$pred + 1.96 * p$se), lty = 2, col = "red")
lines(exp(p$pred - 1.96 * p$se), lty = 2, col = "red")

A veces menos es más. Con 30 observaciones en mano, envié los datos a AUTOBOX, un software que he ayudado a desarrollar. Presento el siguiente análisis con la esperanza de adquirir la recompensa +200 (¡es broma!). He trazado los valores reales y purificados visualmente sugiriendo el impacto de la "actividad reciente". . El modelo que se desarrolló automáticamente se muestra aquí. y aquí . Aquí se presentan los residuos de esta serie de cambio de nivel bastante simple . Las estadísticas del modelo están aquí . En resumen, hubo intervenciones que podrían identificarse empíricamente para generar un proceso ARIMA; dos pulsos y 1 cambio de nivel . El gráfico Actual / Ajuste y Pronóstico resalta aún más el análisis.

Por mi parte, me gustaría ver la gráfica de los residuos de los modelos previamente especificados y, en mi opinión, potencialmente demasiado especificados.

$R$

Debajo está el código:

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

A continuación se muestra la estimación, hubo un aumento de ~ 2356.3 unidades en octubre de 2013 con un error estándar de ~ 481.8 y tiene un efecto de descomposición a partir de entonces. La función identificó automáticamente AR (1). Tuve que hacer un par de iteraciones y hacer una diferenciación estacional y no estacional a 0, que se refleja en el método args.ts en la función tso.

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

A continuación se muestra la trama, tsoutlier es el único paquete que conozco que puede imprimir cambios temporales muy bien en una trama.

Con suerte, este análisis proporcionó respuesta a sus 2, 3 y 4 preguntas, aunque utilizando una metodología diferente. Especialmente la gráfica y los coeficientes proporcionaron el efecto de esta intervención y lo que habría sucedido si no hubiera tenido esta intervención.

También espero que alguien más pueda replicar este gráfico / análisis utilizando el modelado de la función de transferencia en R. No estoy seguro de si esto podría hacerse en R, puede ser que alguien más pueda verificarme al respecto.