Funkcja transferu interwencji ARIMA - Jak wizualizować efekt

11

Mam comiesięczny szereg czasowy z interwencją i chciałbym oszacować wpływ tej interwencji na wynik. Zdaję sobie sprawę, że seria jest raczej krótka, a efekt nie jest jeszcze zakończony.

Dane

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

wprowadź opis zdjęcia tutaj

Metodologia

1) W przypadku tej auto.arimafunkcji zastosowano serię przedinterwencyjną (do października 2013 r.) . Sugerowanym modelem był ARIMA (1,0,0) ze średnią niezerową. Fabuła ACF wyglądała dobrze.

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) Biorąc pod uwagę wykres pełnej serii, odpowiedź impulsowa została wybrana poniżej, przy T = październik 2013,

wprowadź opis zdjęcia tutaj

które według cryer i chan można dopasować w następujący sposób do funkcji 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

Resztki z tego wydawały się OK:

wprowadź opis zdjęcia tutaj

Fabuła wyposażenia i aktualności:

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

wprowadź opis zdjęcia tutaj

Pytania

1) Czy ta metodologia jest odpowiednia do analizy interwencji?

2) Czy mogę spojrzeć na oszacowanie / SE dla składników funkcji przeniesienia i stwierdzić, że efekt interwencji był znaczący?

3) Jak wizualizować efekt funkcji przenoszenia (wykreślić?)

4) Czy istnieje sposób oszacowania, o ile interwencja zwiększyła produkcję po miesiącach „x”? Chyba do tego (a może # 3) pytam, jak pracować z równaniem modelu - gdyby to była prosta regresja liniowa ze zmiennymi fikcyjnymi (na przykład), mogłabym uruchamiać scenariusze z interwencją i bez i mierzyć wpływ - ale nie jestem pewien, jak pracować z tego typu modelem.

DODAJ

Na żądanie przedstawiamy resztki z dwóch parametryzacji.

Najpierw od dopasowania:

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

wprowadź opis zdjęcia tutaj

Następnie z tego dopasowania

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

wprowadź opis zdjęcia tutaj

B_Miner
źródło
Czy byłoby dobrze, jeśli dostarczę Ci rozwiązanie za pomocą oprogramowania SAS?
przepowiednia
Jasne, byłbym ciekawy, jeśli wpadniesz na lepszy model.
B_Miner
Ok, model jest trochę lepszy niż pierwotnie proponowany, ale podobny do @javlacalle.
przepowiednia

Odpowiedzi:

12

Model AR (1) z interwencją zdefiniowaną w równaniu podanym w pytaniu można dopasować, jak pokazano poniżej. Zauważ, jak transferzdefiniowano argument ; potrzebujesz również jednej zmiennej wskaźnikowej xtransfdla każdej interwencji (puls i zmiana przejściowa):

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%

Efekt interwencji można określić ilościowo w następujący sposób:

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)

Możesz wykreślić efekt interwencji w następujący sposób:

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

Całkowity efekt interwencji

ω21ω21

Liczbowo są to szacunkowe wzrosty określone ilościowo w każdym punkcie czasowym spowodowane interwencją w październiku 2013 r .:

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

Interwencja zwiększa wartość obserwowanej zmiennej w październiku 2013 r. O około 75%

Możemy również tworzyć interwencje ręcznie i przekazywać je stats::arimajako zewnętrzne regresory. Interwencje to puls plus przejściowa zmiana z parametrem0.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

Interwencje te są równoważne z wartościami odstającymi dodatków (AO) i zmianami przejściowymi (TC) zdefiniowanymi w opakowaniu tsoutliers. Możesz użyć tego pakietu do wykrycia tych efektów, jak pokazano w odpowiedzi @forecaster lub do zbudowania używanych wcześniej regresorów. Na przykład w tym przypadku:

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

Edytuj 1

Widziałem, że podane równanie można przepisać jako:

(ω0+ω1)ω0ω2B1ω2BPt

i można to określić tak jak za pomocą transfer=list(c(1, 1)).

Jak pokazano poniżej, ta parametryzacja prowadzi w tym przypadku do oszacowania parametrów, które wymagają innego efektu niż poprzednia parametryzacja. Przypomina mi to efekt innowacyjnej wartości odstającej, a nie impulsu plus przejściowej zmiany.

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

Nie znam się na oznaczeniu opakowania, TSAale myślę, że efekt interwencji można teraz określić ilościowo w następujący sposób:

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

parametryzacja efektu interwencji 2

Efekt ten można teraz opisać jako gwałtowny wzrost w październiku 2013 r., Po którym następuje spadek w przeciwnym kierunku; wtedy efekt interwencji zanika szybko na przemian pozytywne i negatywne skutki niszczenia masy.

Ten efekt jest nieco szczególny, ale może być możliwy w rzeczywistych danych. W tym momencie przyjrzałbym się kontekstowi twoich danych i wydarzeniom, które mogły mieć wpływ na dane. Na przykład, czy nastąpiła zmiana zasad, kampania marketingowa, odkrycie ..., które mogą wyjaśnić interwencję w październiku 2013 r. Jeśli tak, czy bardziej sensowne jest, aby to wydarzenie miało wpływ na dane, jak opisano wcześniej lub jak stwierdziliśmy z początkową parametryzacją?

18.9415.42 ). Wykres oryginalnej serii nie sugeruje wyraźnego dopasowania do ostrych zmian związanych z pomiarem drugiej zmiennej interwencyjnej.

0.9

Edytuj 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 dla różnych wartości omega

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

Prognozy można uzyskać i wyświetlić w następujący sposób:

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

zaobserwowane i prognozowane wartości

Pierwsze prognozy odpowiadają stosunkowo dobrze obserwowanym wartościom (szara linia przerywana). Pozostałe prognozy pokazują, w jaki sposób seria będzie kontynuować ścieżkę do pierwotnej średniej. Przedziały ufności są jednak duże, co odzwierciedla niepewność. Dlatego powinniśmy być ostrożni i zmieniać model w miarę rejestrowania nowych danych.

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
źródło
To wspaniale, dziękuję! Miałem kilka dalszych obserwacji, jeśli nie masz nic przeciwko. 1) Czy zastosowany przeze mnie proces jest prawidłowy? 2) Czy uznałbyś dopasowanie modelu za „wystarczająco dobre”, aby wykorzystać szacunki do oszacowania efektu interwencji? 3) Czy nie powinienem być w stanie korzystać z mojej parametryzacji, tj. Transfer = lista (c (1,1)) jako ekwiwalent i uzyskiwać dość bliskie wyniki? Przykład, który
podążyłem
@B_Miner Masz rację, zredagowałem moją odpowiedź.
javlacalle,
Zgadzam się z tobą, że w przypadku dwóch modeli pierwsza parametryzacja (być może przy usuniętym nieistotnym impulsie) byłaby najlepsza. Dlaczego dwie parametryzacje nie dają tego samego modelu, kiedy uważam, że powinny, jest tajemnicą. Wyślę pocztą e-mail do twórcy pakietu (który również napisał książkę, w której wspomniano o ich równoważności).
B_Miner
Dane to liczba certyfikatów depozytów otwieranych miesięcznie. Interwencją był wzrost średniej stopy procentowej, która wzrosła od 13 października. Poziom stopy procentowej pozostał względnie stały od 13 października. Wydawało mi się, że po gwałtownym wzroście popyt na produkt zaczął się zmniejszać - Nie jestem pewien, czy powróci do poprzedniej średniej, czy też ustabilizuje się na jakimś podwyższonym (z poprzedniego) poziomie.
B_Miner
B_miner, na podstawie danych, których naprawdę nie możemy stwierdzić, czy popyt ustabilizuje się na nowy sposób.
prezenter
4

Czasami mniej znaczy więcej. Mając 30 obserwacji w ręku, przesłałem dane do AUTOBOX, oprogramowania, które pomogłem opracować. Przesyłam następującą analizę w nadziei na zdobycie nagrody +200 (tylko żartuję!). Naszkicowałem rzeczywiste i oczyszczone wartości, sugerując wizualnie wpływ „ostatniej aktywności”. wprowadź opis zdjęcia tutaj. Model, który został opracowany automatycznie, pokazano tutaj. wprowadź opis zdjęcia tutaji tutaj wprowadź opis zdjęcia tutaj. Resztki tej raczej prostej serii z przesuniętym poziomem są tutaj przedstawione wprowadź opis zdjęcia tutaj. Statystyki modelu są tutaj wprowadź opis zdjęcia tutaj. Podsumowując, były interwencje, które można zidentyfikować empirycznie, czyniąc proces ARIMA; dwa impulsy i 1 zmiana poziomu wprowadź opis zdjęcia tutaj. Wykres Rzeczywiste / Dopasowanie i Prognoza dodatkowo podkreśla analizę.wprowadź opis zdjęcia tutaj

Chciałbym na przykład zobaczyć wykres reszt z wcześniej określonych i, moim zdaniem, potencjalnie nadmiernie określonych modeli.

IrishStat
źródło
Nie znam Autobox, ale czy część szumu modelu jest taka sama jak pierwotnie: średnia niezerowa i AR (1)?
B_Miner
Czy to wyjście mówi, że jedyną „interwencją” w okresie od 13 października do bieżącego okresu jest pojedynczy impuls dla 13 października, a następnie seria wraca do normalnego średniego poziomu?
B_Miner
Dodałem resztki z obu parametryzacji. Moim zdaniem wydaje się, że pierwszy wymieniony przeze mnie (ten, który pierwotnie pasuje do javlacalle) jest lepszy. Zgodzić się?
B_Miner
1) Część hałasu to AR (1) z niezerową średnią
IrishStat
1) Część szumowa to AR (1) ze niezerową średnią; 2) Istnieją 2 okresy interwencji 22 i okres 3, a po 13 października powraca do nowego poziomu, który rozpoczął się we wrześniu 13; 3) Biorąc pod uwagę wybór między dwoma, o których wspomniałeś, zgadzam się, ALE wolę model AUTOBOX ze względu na jego prostotę i wydajność. Możesz dowiedzieć się więcej o AUTOBOX na autobox.com/cms
IrishStat
3

R

Poniżej znajduje się kod:

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

Poniżej znajduje się szacunek, w październiku 2013 r. Nastąpił wzrost o ~ 236,3 jednostki przy standardowym błędzie ~ 481,8, a następnie ma on rozkład niszczący. Funkcja automatycznie identyfikuje AR (1). Musiałem wykonać kilka iteracji i wprowadzić różnicowanie zarówno sezonowe, jak i niesezonowe do 0, co znajduje odzwierciedlenie w args.tsmethod w funkcji 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

Poniżej znajduje się fabuła, tsoutlier to jedyny znany mi pakiet, który może ładnie wydrukować tymczasowe zmiany na wykresie.

wprowadź opis zdjęcia tutaj

Ta analiza, miejmy nadzieję, dostarczyła odpowiedzi na twoje 2, 3 i 4 pytania, aczkolwiek przy użyciu innej metodologii. Zwłaszcza wykres i współczynniki zapewniły efekt tej interwencji oraz to, co by się stało, gdyby nie interwencja.

Mam również nadzieję, że ktoś inny może powielić ten wykres / analizę przy użyciu modelowania funkcji przenoszenia w R. Nie jestem pewien, czy można to zrobić w R, może ktoś inny może mnie sprawdzić.

Synoptyk
źródło
Dzięki. Tak, ten wykres jest tym, co chciałbym zrobić z modelu arimax - spójrz na interwencję i bez interwencji (i odejmij). Myślę, że funkcja filtru w R może być użyta do wygenerowania wartości funkcji przenoszenia dla każdego miesiąca (a następnie po prostu wykreśl ją w celu wizualizacji), ale nie mogę wymyślić, jak to zrobić dla funkcji interwencji dowolnego impulsu.
B_Miner