ARIMA Intervention Transfer Function - Visualisierung des Effekts

11

Ich habe eine monatliche Zeitreihe mit einer Intervention und möchte die Auswirkung dieser Intervention auf das Ergebnis quantifizieren. Mir ist klar, dass die Serie ziemlich kurz ist und der Effekt noch nicht abgeschlossen ist.

Die Daten

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

Geben Sie hier die Bildbeschreibung ein

Die Methodik

1) Die Vorinterventionsreihe (bis Oktober 2013) wurde mit der auto.arimaFunktion verwendet. Das vorgeschlagene Modell war ARIMA (1,0,0) mit einem Mittelwert ungleich Null. Die ACF-Handlung sah gut aus.

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) In Anbetracht des Diagramms der vollständigen Reihe wurde die Impulsantwort unten mit T = Okt 2013 gewählt.

Geben Sie hier die Bildbeschreibung ein

die laut cryer und chan mit der arimax-funktion wie folgt angepasst werden können:

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

Die Residuen davon schienen in Ordnung zu sein:

Geben Sie hier die Bildbeschreibung ein

Die Handlung von angepassten und tatsächlichen:

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

Geben Sie hier die Bildbeschreibung ein

Die Fragen

1) Ist diese Methode für die Interventionsanalyse korrekt?

2) Kann ich die Schätzung / SE für die Komponenten der Übertragungsfunktion betrachten und sagen, dass der Effekt der Intervention signifikant war?

3) Wie kann man den Effekt der Übertragungsfunktion visualisieren (zeichnen?)

4) Gibt es eine Möglichkeit zu schätzen, um wie viel die Intervention die Leistung nach 'x' Monaten erhöht hat? Ich denke, dafür (und vielleicht # 3) frage ich mich, wie man mit einer Gleichung des Modells arbeitet - wenn dies eine einfache lineare Regression mit Dummy-Variablen wäre (zum Beispiel), könnte ich Szenarien mit und ohne Intervention ausführen und die Auswirkungen messen - aber ich bin mir nur nicht sicher, wie ich diese Art von Modell arbeiten soll.

HINZUFÜGEN

Auf Anfrage sind hier die Residuen der beiden Parametrisierungen aufgeführt.

Zuerst aus der Passform:

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

Geben Sie hier die Bildbeschreibung ein

Dann aus dieser Passform

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

Geben Sie hier die Bildbeschreibung ein

B_Miner
quelle
Wäre es in Ordnung, wenn ich Ihnen eine Lösung mit SAS-Software anbieten würde?
Prognostiker
Klar, ich wäre gespannt, ob Sie ein besseres Modell finden.
B_Miner
Ok, das Modell ist wenig besser als ursprünglich vorgeschlagen, aber ähnlich wie @javlacalle.
Prognostiker

Antworten:

11

Ein AR (1) -Modell mit dem in der in der Frage angegebenen Gleichung definierten Eingriff kann wie unten gezeigt angepasst werden. Beachten Sie, wie das Argument transferdefiniert ist. Sie benötigen außerdem eine Indikatorvariable xtransffür jede der Interventionen (den Puls und die vorübergehende Änderung):

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

ω0ω1coeftest

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%.

Der Interventionseffekt kann wie folgt quantifiziert werden:

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)

Sie können den Effekt der Intervention wie folgt darstellen:

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

Totaler Interventionseffekt

ω21ω21

Numerisch sind dies die geschätzten Zuwächse, die zu jedem Zeitpunkt quantifiziert wurden, der durch die Intervention im Oktober 2013 verursacht wurde:

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::arima0,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

ω20,9231xregω2

Diese Eingriffe entsprechen einem im Paket definierten additiven Ausreißer (AO) und einer vorübergehenden Veränderung (TC) tsoutliers. Mit diesem Paket können Sie diese Effekte erkennen, wie in der Antwort von @forecaster gezeigt, oder die zuvor verwendeten Regressoren erstellen. Zum Beispiel in diesem Fall:

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

Bearbeiten 1

Ich habe gesehen, dass die von Ihnen angegebene Gleichung wie folgt umgeschrieben werden kann:

(ω0+ω1)- -ω0ω2B.1- -ω2B.P.t

und es kann so angegeben werden, wie Sie es verwendet haben transfer=list(c(1, 1)).

Wie unten gezeigt, führt diese Parametrisierung in diesem Fall zu Parameterschätzungen, die einen anderen Effekt als die vorherige Parametrisierung beinhalten. Es erinnert mich eher an die Wirkung eines innovativen Ausreißers als an einen Puls und eine vorübergehende Veränderung.

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

Ich bin mit der Notation des Pakets nicht sehr vertraut, TSAaber ich denke, dass die Wirkung der Intervention jetzt wie folgt quantifiziert werden kann:

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

Parametrisierung des Interventionseffekts 2

Der Effekt kann nun als starker Anstieg im Oktober 2013 beschrieben werden, gefolgt von einem Rückgang in die entgegengesetzte Richtung; dann verschwindet die Wirkung des Eingriffs schnell und wechselt abwechselnd die positiven und negativen Auswirkungen des abnehmenden Gewichts.

Dieser Effekt ist etwas eigenartig, kann aber in realen Daten möglich sein. An dieser Stelle würde ich den Kontext Ihrer Daten und die Ereignisse betrachten, die die Daten beeinflusst haben könnten. Hat es beispielsweise eine Richtlinienänderung, eine Marketingkampagne, eine Entdeckung usw. gegeben, die die Intervention im Oktober 2013 erklären könnte? Wenn ja, ist es sinnvoller, dass sich dieses Ereignis auf die zuvor beschriebenen oder von uns gefundenen Daten auswirkt mit der anfänglichen Parametrierung?

- -18.94- -15.42

0,9

Bearbeiten 2

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

AIC für verschiedene Omega-Werte

ω2=0,880,9ω2=1

ω2=0,9

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

Die Vorhersagen können wie folgt abgerufen und angezeigt werden:

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

beobachtete und prognostizierte Werte

Die ersten Vorhersagen stimmen relativ gut mit den beobachteten Werten überein (grau gepunktete Linie). Die verbleibenden Prognosen zeigen, wie die Serie den Weg zum ursprünglichen Mittelwert fortsetzen wird. Die Konfidenzintervalle sind dennoch groß und spiegeln die Unsicherheit wider. Wir sollten daher vorsichtig sein und das Modell überarbeiten, wenn neue Daten aufgezeichnet werden.

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")
javlacalle
quelle
Das ist großartig, danke! Ich hatte ein paar Follow-ups, wenn es Ihnen nichts ausmacht. 1) Ist der von mir verfolgte Prozess korrekt? 2) Würden Sie die Anpassung des Modells als "gut genug" betrachten, um die Schätzungen zur Quantifizierung der Wirkung der Intervention zu verwenden? 3) Sollte ich meine Parametrisierung, dh Übertragung = Liste (c (1,1)), nicht als Äquivalent verwenden und ziemlich genaue Ergebnisse erzielen können? Das Beispiel, dem ich aus einem Lehrbuch folgte, schlug vor, dass ich dazu in der Lage sein sollte, aber in diesem Beispiel sind die Ergebnisse nicht nahe ...
B_Miner
@B_Miner Du hast recht Ich habe meine Antwort bearbeitet.
Javlacalle
Ich stimme Ihnen zu, dass von den beiden Modellen die erste Parametrisierung (möglicherweise ohne entfernten nicht signifikanten Impuls) am besten passt. Warum die beiden Parametrisierungen nicht dasselbe Modell ergeben, wenn ich glaube, dass sie es sollten, ist ein Rätsel. Ich werde dem Paketentwickler eine E-Mail senden (der auch das Buch geschrieben hat, in dem ihre Äquivalenz erwähnt wird).
B_Miner
Die Daten sind die Anzahl der pro Monat eröffneten Einlagenzertifikate. Die Intervention war ein Anstieg des durchschnittlichen Zinssatzes, der ab dem 13. Oktober anstieg. Das Zinsniveau ist seit dem 13. Oktober relativ konstant geblieben. Es schien mir, dass nach dem Anstieg die Nachfrage nach dem Produkt allmählich nachließ - Ich bin mir nicht sicher, ob es zum vorherigen Mittelwert zurückkehren oder sich auf einem erhöhten (vom vorherigen) Niveau niederlassen wird.
B_Miner
B_miner, basierend auf den Daten, die wir nicht wirklich schließen können, wenn sich die Nachfrage auf einen neuen Mittelwert einpendelt.
Prognostiker
4

Manchmal ist weniger mehr. Mit 30 Beobachtungen übermittelte ich die Daten an AUTOBOX, eine Software, die ich mitentwickelt habe. Ich reiche die folgende Analyse in der Hoffnung ein, die +200 Belohnung zu erhalten (nur ein Scherz!). Ich habe die tatsächlichen und gereinigten Werte grafisch dargestellt, um die Auswirkungen der "jüngsten Aktivitäten" visuell darzustellen. Geben Sie hier die Bildbeschreibung ein. Hier wird das automatisch entwickelte Modell gezeigt. Geben Sie hier die Bildbeschreibung einund hier Geben Sie hier die Bildbeschreibung ein. Die Residuen dieser eher einfachen pegelverschobenen Reihe werden hier vorgestellt Geben Sie hier die Bildbeschreibung ein. Die Modellstatistik finden Sie hier Geben Sie hier die Bildbeschreibung ein. Zusammenfassend gab es Interventionen, die empirisch identifiziert werden konnten, um einen ARIMA-Prozess zu erstellen. zwei Impulse und 1 Pegelverschiebung Geben Sie hier die Bildbeschreibung ein. Das Diagramm Ist / Anpassung und Prognose hebt die Analyse weiter hervor.Geben Sie hier die Bildbeschreibung ein

Ich jedenfalls möchte die Darstellung der Residuen aus den zuvor spezifizierten und meiner Meinung nach möglicherweise überspezifizierten Modellen sehen.

IrishStat
quelle
Ich bin mit Autobox nicht vertraut, aber ist der Geräuschanteil des Modells derselbe wie ursprünglich: ein Mittelwert ungleich Null und ein AR (1)?
B_Miner
Sagt diese Ausgabe, dass der einzige "Eingriff" im Zeitraum vom 13. Oktober bis zum aktuellen Zeitraum ein einzelner Impuls für den 13. Oktober ist und die Serie dann zu ihrem normalen Mittelwert zurückkehrt?
B_Miner
Ich habe die Residuen aus beiden Parametrisierungen hinzugefügt. Meines Erachtens scheint der erste, den ich aufgelistet habe (der ursprünglich von javlacalle passt), besser zu sein. Zustimmen?
B_Miner
1) Der Rauschteil ist ein AR (1) mit einem Mittelwert ungleich Null
IrishStat
1) Der Rauschanteil ist ein AR (1) mit einem Mittelwert ungleich Null; 2) Es gibt 2 Interventionsperioden 22 und 3 und nach dem 13. Oktober kehrt es auf ein neues Niveau zurück, das am 13. September begann. 3) Angesichts der Wahl zwischen den beiden, die Sie erwähnt haben, stimme ich zu, ABER ich bevorzuge das AUTOBOX-Modell wegen seiner Einfachheit und Effizienz. Weitere Informationen zu AUTOBOX
IrishStat
3

Basierend auf meinem ähnlichen Beitrag zu Ihrer früheren Frage habe ich die Funktion tso im Paket tsoutliers in verwendetR.und es wurde automatisch eine vorübergehende Änderung im Oktober 2013 erkannt. Bitte beachten Sie, dass sich die vorübergehende Änderung von der Rampenverschiebung in der Übertragungsfunktion unterscheidet, nach der Sie suchen. Ich glaube nicht, dass es ein Paket / eine Funktion gibt, von der ich weiß, dass sie die Übertragungsfunktion visualisieren kann. Hoffentlich würde dies einige Einblicke liefern. Ich habe keine Protokolltransformation verwendet, sondern sie direkt modelliert. Das Paket tsoutliers kann als automatische Interventionserkennung angesehen werden.

Unten ist der 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

Nachstehend die Schätzung: Im Oktober 2013 gab es einen Anstieg um ~ 2356,3 Einheiten mit einem Standardfehler von ~ 481,8 und hat danach einen abnehmenden Effekt. Die Funktion identifizierte AR (1) automatisch. Ich musste einige Iterationen durchführen und sowohl saisonale als auch nicht saisonale Differenzierungen auf 0 vornehmen, was sich in der args.tsmethod in der tso-Funktion widerspiegelt.

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

Unten ist die Handlung, tsoutlier ist das einzige Paket, von dem ich weiß, dass es temporäre Änderungen so schön in einer Handlung drucken kann.

Geben Sie hier die Bildbeschreibung ein

Diese Analyse lieferte hoffentlich eine Antwort auf Ihre 2, 3 und 4 Fragen, wenn auch mit einer anderen Methode. Insbesondere die Darstellung und die Koeffizienten lieferten den Effekt dieser Intervention und was wäre passiert, wenn Sie diese Intervention nicht gehabt hätten.

Ich hoffe auch, dass jemand anderes dieses Diagramm / diese Analyse mithilfe der Übertragungsfunktionsmodellierung in R replizieren kann. Ich bin mir nicht sicher, ob dies in R möglich ist. Vielleicht kann mich jemand anderes tatsächlich überprüfen.

Prognostiker
quelle
Vielen Dank. Ja, diese Handlung möchte ich mit dem Arimax-Modell machen - mit und ohne Intervention betrachten (und subtrahieren). Ich denke, die Filterfunktion in R kann verwendet werden, um den Übertragungsfunktionswert für jeden Monat zu generieren (und ihn dann nur zur Visualisierung zu zeichnen), aber ich kann nicht herausfinden, wie dies für eine beliebige Impulsinterventionsfunktion zu tun ist.
B_Miner