Un AR(1) con el modelo de intervención definido en la ecuación dada en la pregunta se puede montar como se muestra a continuación. Observe cómo el argumento transfer está definido; también necesita una variable de indicador en xtransf para cada una de las intervenciones (el pulso y la transitoriedad de cambio):
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
Usted puede probar la significación de cada intervención mirando la estadística t de los coeficientes $\omega_0$$\omega_1$. Para su conveniencia, usted puede usar la función 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
En este caso, el pulso no es significativa en el $5\%$ de nivel de significación. Su efecto puede ser capturadas por la transitoriedad de cambio.
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)
Usted puede trazar el efecto de la intervención de la siguiente manera:
plot(100*(intv.effect - 1), type = "h", main = "Total intervention effect")
![Total intervention effect]()
El efecto es relativamente persistente, porque $\omega_2$ está cerca de a $1$ (si $\omega_2$ eran iguales a $1$ observaríamos un permanente cambio de nivel).
Numéricamente, estos son el aumento previsto de la cuantificados en cada punto de tiempo causada 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
La intervención aumenta el valor de la variable observada en octubre de 2013
por alrededor de un $75\%$. En los períodos posteriores el efecto se mantiene, pero con una disminución de peso.
También podríamos crear las intervenciones de la mano y paso a stats::arima externo como regresores. Las intervenciones son de pulso, además de un cambio transitorio con el parámetro $0.9231$ y puede ser construido de la siguiente manera.
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
Las mismas estimaciones de los coeficientes, como anteriormente se obtienen. Aquí se fija $\omega_2$$0.9231$. La matriz xreg es el tipo de variable ficticia que puede que tenga que probar diferentes escenarios. También se podría establecer diferentes valores de $\omega_2$ y comparar su efecto.
Estas intervenciones son equivalentes a un valor atípico aditivo (AO) y un transitorio de cambio (TC) definidos en el paquete tsoutliers. Usted puede usar este paquete para detectar estos efectos, como se muestra en la respuesta por @pronosticador o para construir los regresores utilizados antes. 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
Edición 1
He visto que la ecuación que dio puede ser reescrita como:
$$
\frac{(\omega_0 + \omega_1) - \omega_0 \omega_2 B}{1 - \omega_2 B} P_t
$$
y puede ser especificado como lo hizo usando transfer=list(c(1,1)).
Como se muestra a continuación, esta parametrización conduce, en este caso, para las estimaciones de los parámetros que implican un efecto diferente en comparación con el anterior de la parametrización. Me recuerda el efecto de un innovador outlier en lugar de un pulso, además de 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 TSA pero creo que el efecto de la intervención puede ahora ser cuantificados 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)")
![intervention effect parameterization 2]()
El efecto puede ser descrito ahora como un fuerte aumento en octubre de 2013, seguido por
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 de la descomposición de peso.
Este efecto es un tanto peculiar, pero que puede ser posible en datos reales. En este punto me gustaría ver en el contexto de los datos y los eventos que pueden haber afectado a los datos. Por ejemplo, ha sido un cambio de política, campaña de marketing, descubrimiento,... que puede explicar la intervención en octubre de 2013. Si es así, es más sensato que este evento tiene un efecto en los datos como se describe antes o como nos encontramos con la parametrización inicial?
De acuerdo a la AIC, el modelo inicial sería preferible porque es más baja ($-18.94$ contra $-15.42$). La trama de la serie original no sugieren una clara coincidencia con los cambios bruscos que intervienen en la medición de la segunda intervención de la variable.
Sin conocer el contexto de los datos, yo diría que un AR(1) modelo con un cambio transitorio con el parámetro $0.9$ sería apropiado para un modelo de los datos y la medida de la intervención.
Edit 2
El valor de $\omega_2$ determina la rapidez del efecto de la intervención decae a cero, por lo que el parámetro clave en el modelo. Podemos inspeccionar este al ajustar el modelo para un rango de valores de $\omega_2$. A continuación, el AIC se almacena para cada uno de estos modelos.
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]()
El menor AIC se encuentra por $\omega_2 = 0.88$ (de acuerdo con el valor estimado de antes). Este parámetro implica una relativamente persistente pero el efecto transitorio. Podemos concluir que el efecto es temporal ya que con valores superiores a $0.9$ la AIC aumenta (recuerde que en el límite, $\omega_2=1$, la intervención se convierte en un permanente cambio de nivel).
La intervención debe ser incluido en las previsiones. La obtención de pronósticos para los períodos que ya se ha observado es un ejercicio útil para evaluar el desempeño de los pronósticos. El código siguiente se supone que la serie termina en octubre de 2013. Los pronósticos se obtienen entonces, incluyendo la intervención con el parámetro $\omega_2=0.9$.
Primero nos ajuste el AR(1) con el modelo de la intervención como un regresor (con el parámetro $\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 pueden ser obtenidos y se muestra 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")
![observed and forecasted values]()
La primera previsiones relativamente bien los valores observados (gris línea de puntos). El resto de las previsiones muestran cómo la serie seguirá la ruta de acceso a la original que decir. Los intervalos de confianza son, no obstante, grandes, reflejando la incertidumbre. Por lo tanto, debemos tener cuidado y revisar el modelo como los nuevos datos que se registran.
$95\%$ intervalos de confianza puede ser añadido a la anterior parcela de la siguiente manera:
lines(exp(p$pred + 1.96 * p$se), lty = 2, col = "red")
lines(exp(p$pred - 1.96 * p$se), lty = 2, col = "red")
. El modelo que fue automáticamente desarrollado se muestra aquí. 
y aquí 
. Los residuos de este nivel simple desplazado a la serie que aquí se presenta 
. El modelo estadísticas están aquí 
. En resumen, hubo una de las intervenciones que podrían ser identificado empíricamente la representación de un proceso ARIMA ; dos pulsos y 1 cambio de nivel
. El Real/Ajuste y Previsión gráfica resalta aún más el análisis.
