Konsekwencje modelowania niestacjonarnego procesu przy użyciu ARiMR?

Rozumiem, że powinniśmy używać ARIMA do modelowania niestacjonarnych szeregów czasowych. Ponadto wszystko, co czytam, mówi, że ARMA powinna być używana tylko do stacjonarnych szeregów czasowych.

Próbuję zrozumieć, co dzieje się w praktyce, kiedy błędnie klasyfikujesz model i zakładasz, d = 0że szereg czasowy jest niestacjonarny? Na przykład:

controlData <- arima.sim(list(order = c(1,1,1), ar = .5, ma = .5), n = 44)

dane kontrolne wyglądają następująco:

 [1]   0.0000000   0.1240838  -1.4544087  -3.1943094  -5.6205257
 [6]  -8.5636126 -10.1573548  -9.2822666 -10.0174493 -11.0105225
[11] -11.4726127 -13.8827001 -16.6040541 -19.1966633 -22.0543414
[16] -24.8542959 -25.2883155 -23.6519271 -21.8270981 -21.4351267
[21] -22.6155812 -21.9189036 -20.2064343 -18.2516852 -15.5822178
[26] -13.2248230 -13.4220158 -13.8823855 -14.6122867 -16.4143756
[31] -16.8726071 -15.8499558 -14.0805114 -11.4016515  -9.3330560
[36]  -7.5676563  -6.3691600  -6.8471371  -7.5982880  -8.9692152
[41] -10.6733419 -11.6865440 -12.2503202 -13.5314306 -13.4654890

Zakładając, że nie wiedziałem, że dane są ARIMA(1,1,1), mógłbym rzucić okiem pacf(controlData).

pacf (controlData)

Następnie używam Dickey-Fuller, aby sprawdzić, czy dane nie są stacjonarne:

require('tseries')
adf.test(controlData)

# Augmented Dickey-Fuller Test
#
# data:  controlData
# Dickey-Fuller = -2.4133, Lag order = 3, p-value = 0.4099
# alternative hypothesis: stationary

adf.test(controlData, k = 1)

# Augmented Dickey-Fuller Test
#
#data:  controlData
# Dickey-Fuller = -3.1469, Lag order = 1, p-value = 0.1188
# alternative hypothesis: stationary

Mogę więc założyć, że dane to ARIMA (2,0, *). Następnie użyć, auto.arima(controlData)aby spróbować jak najlepiej dopasować?

require('forecast')
naiveFit <- auto.arima(controlData)
naiveFit
# Series: controlData 
# ARIMA(2,0,1) with non-zero mean 
# 
# Coefficients:
#          ar1      ar2     ma1  intercept
#      1.4985  -0.5637  0.6427   -11.8690
# s.e.  0.1508   0.1546  0.1912     3.2647
#
# sigma^2 estimated as 0.8936:  log likelihood=-64.01
# AIC=138.02   AICc=139.56   BIC=147.05

Tak więc, mimo że przeszłe i przyszłe dane to ARIMA (1,1,1), mogę ulec pokusie sklasyfikowania ich jako ARIMA (2,0,1). tsdata(auto.arima(controlData))też wygląda dobrze.

Oto, co znalazłby poinformowany modelarz:

informedFit <- arima(controlData, order = c(1,1,1))
# informedFit
# Series: controlData 
# ARIMA(1,1,1)                    
#
# Coefficients:
#          ar1     ma1
#       0.4936  0.6859
# s.e.  0.1564  0.1764
#
# sigma^2 estimated as 0.9571:  log likelihood=-62.22
# AIC=130.44   AICc=131.04   BIC=135.79

1) Dlaczego te kryteria informacyjne są lepsze niż wybrany model auto.arima(controlData)?

Teraz po prostu graficznie porównuję rzeczywiste dane i dwa modele:

plot(controlData)
lines(fitted(naiveFit), col = "red")
lines(fitted(informedFit), col = "blue")

tsPlots

2) Jako adwokat diabła, jakie konsekwencje zapłaciłbym, stosując model ARIMA (2, 0, 1) jako model? Jakie jest ryzyko tego błędu?

3) Najbardziej niepokoją mnie jakiekolwiek implikacje dla prognoz na wiele okresów. Zakładam, że byłyby mniej dokładne? Szukam tylko dowodu.

4) Czy zaproponowałbyś alternatywną metodę wyboru modelu? Czy są jakieś problemy z moim rozumowaniem jako „niedoinformowanego” modelarza?

Jestem naprawdę ciekawy, jakie są inne konsekwencje tego rodzaju błędnej klasyfikacji. Szukałem niektórych źródeł i nic nie znalazłem. Cała literatura, którą mogłem znaleźć, dotyczy tylko tego tematu, zamiast tego stwierdzenie, że dane powinny być nieruchome przed wykonaniem ARiMR, a jeśli są niestacjonarne, to należy je różnicować d razy.

Dzięki!

r time-series arima stationarity Clark Henry
źródło

Mam wrażenie, że jest to analogiczne do założenia „błędów ortogonalnych” w regresji przekrojowej (tzn. Wpływa na błędy standardowe, ale nie na współczynniki), ale naprawdę interesuje mnie odpowiedź.

shadowtalker

Odpowiedzi:

Mam wrażenie, że na to pytanie nie ma unikalnej, w pełni ogólnej odpowiedzi, więc zbadam tylko najprostszy przypadek i to w nieco nieformalny sposób.

Załóżmy, że prawdziwym mechanizmem generowania danych jest z zwykle składową białego szumu o średniej zerowej iid, . Powyższe implikuje również to

\begin{matrix} (1) & y_{t} = y_{t - 1} + u_{t}, t = 1, . . ., T., y_{0} = 0 \end{matrix}

$y_t = y_{t-1} + u_t,\;\; t=1,...,T,\;\; y_0 =0 \tag{1}$

u_{t}

$u_t$

E (u_{t}^{2}) = σ_{u}^{2}

$E(u_t^2)= \sigma^2_u$

\begin{matrix} (2) & y_{t} = \sum_{ja = 1}^{t} u_{ja} \end{matrix}

$y_t = \sum_{i=1}^tu_i \tag{2}$

Możemy określić model, nazywają go modelować $A$

\begin{matrix} (3) & y_{t} = β y_{t - 1} + u_{t}, t = 1, . . ., T., y_{0} = 0 \end{matrix}

$y_t = \beta y_{t-1} + u_t,\;\; t=1,...,T,\;\; y_0 =0 \tag{3}$

i otrzymujemy oszacowanie dla postulowanego (omówimy metodę szacowania tylko w razie potrzeby). $\hat \beta$ $\beta$

Tak więc będzie przewidywanie -wyprzedzeń $k$

\begin{matrix} (4) & {\hat{y}}_{T. + k} = {\hat{β}}^{k} y_{T.} \end{matrix}

$\hat y_{T+k} = \hat \beta^k y_T \tag{4}$

i jego MSE będzie

M. S. {mi}_{ZA} [{\hat{y}}_{T. + k}] = mi {({\hat{β}}^{k} y_{T.} - y_{T. + k})}^{2)}

$MSE_A[\hat y_{T+k}] = E\left(\hat \beta^k y_T-y_{T+k}\right)^2$

\begin{matrix} (5) & = mi {[({\hat{β}}^{k} - 1) y_{T.} - \sum_{ja = T. + 1}^{T. + k} u_{ja}]}^{2)} = mi [({\hat{β}}^{k} - 1)^{2)} y_{T.}^{2)}] + k σ_{u}^{2)} \end{matrix}

$=E\left[(\hat \beta^k-1) y_T -\sum_{i=T+1}^{T+k}u_i \right]^2 = E\big[(\hat\beta^k-1)^2 y_T^2\big]+ k\sigma^2_u \tag{5}$

(zanika środkowa część kwadratu, a także iloczyny przyszłych błędów).

Powiedzmy teraz, że różnicowaliśmy nasze dane i określiliśmy model $B$

\begin{matrix} (6) & Δ y_{t} = γ Δ y_{t - 1} + u_{t} \end{matrix}

$\Delta y_t = \gamma \Delta y_{t-1} + u_t \tag{6}$

i uzyskałem oszacowanie . Nasz zróżnicowany model można napisać $\hat \gamma$

\begin{matrix} (7) & y_{t} = y_{t - 1} + γ (y_{t - 1} - y_{t - 2)}) + u_{t} \end{matrix}

$y_t = y_{t-1} + \gamma (y_{t-1}-y_{t-2}) + u_t \tag{7}$

więc prognozując poziom tego procesu, będziemy mieli

{\hat{y}}_{T. + 1} = y_{T.} + \hat{γ} (y_{T.} - y_{T. - 1})

$\hat y_{T+1} = y_{T} + \hat \gamma (y_{T}-y_{T-1})$

co w rzeczywistości, biorąc pod uwagę prawdziwy MZD będzie

\begin{matrix} (8) & {\hat{y}}_{T. + 1} = y_{T.} + \hat{γ} u_{T.} \end{matrix}

$\hat y_{T+1} = y_{T} + \hat \gamma u_T \tag {8}$

Jest to łatwe do sprawdzenia, a następnie, że dla modelu , $B$

{\hat{y}}_{T. + k} = y_{T.} + (\hat{γ} + {\hat{γ}}^{2)} + . . . + {\hat{γ}}^{k}) u_{T.}

$\hat y_{T+k} = y_{T} + \big(\hat \gamma + \hat \gamma^2+...+\hat \gamma^k)u_T$

Oczekujemy teraz, że biorąc pod uwagę każdą „sprawdzoną i wypróbowaną” procedurę szacunkową, uzyskamy ponieważ jej prawdziwa wartość wynosi , z wyjątkiem sytuacji, gdy mamy za mało danych lub w bardzo „złym” kształcie . Możemy więc powiedzieć, że w większości przypadków będziemy mieli $|\hat \gamma|<1$ $0$

\begin{matrix} (9) & {\hat{y}}_{T. + k} = y_{T.} + \frac{\hat{γ} - {\hat{γ}}^{k + 1}}{1 - \hat{γ}} u_{T.} \end{matrix}

$\hat y_{T+k} = y_{T} + \frac {\hat \gamma - \hat \gamma ^{k+1}}{1-\hat \gamma}u_T \tag{9}$

a więc

\begin{matrix} (10) & M. S. {mi}_{b} [{\hat{y}}_{T. + k}] = mi [{(\frac{\hat{γ} - {\hat{γ}}^{k + 1}}{1 - \hat{γ}})}^{2)} u_{T.}^{2)}] + k σ_{u}^{2)} \end{matrix}

$MSE_B[\hat y_{T+k}] = E\left[\left(\frac {\hat \gamma - \hat \gamma ^{k+1}}{1-\hat \gamma}\right)^2u_T^2\right] + k\sigma^2_u \tag{10}$

podczas gdy powtarzam dla wygody

\begin{matrix} (5) & M. S. {mi}_{ZA} [{\hat{y}}_{T. + k}] = mi [({\hat{β}}^{k} - 1)^{2)} y_{T.}^{2)}] + k σ_{u}^{2)} \end{matrix}

$MSE_A[\hat y_{T+k} ] = E\big[(\hat\beta^k-1)^2 y_T^2\big]+ k\sigma^2_u \tag{5}$

Tak więc, aby zróżnicowany model działał lepiej pod względem przewidywania MSE, chcemy

M. S. {mi}_{b} [{\hat{y}}_{T. + k}] \leq M. S. {mi}_{ZA} [{\hat{y}}_{T. + k}]

$MSE_B[\hat y_{T+k}] \leq MSE_A[\hat y_{T+k}]$

\Rightarrow mi [{(\frac{\hat{γ} - {\hat{γ}}^{k + 1}}{1 - \hat{γ}})}^{2)} u_{T.}^{2)}] \leq mi [({\hat{β}}^{k} - 1)^{2)} y_{T.}^{2)}]

$\Rightarrow E\left[\left(\frac {\hat \gamma - \hat \gamma ^{k+1}}{1-\hat \gamma}\right)^2u_T^2\right] \leq E\big[(\hat\beta^k-1)^2 y_T^2\big]$

Podobnie jak w przypadku estymatora w modelu , rozszerzamy tę samą uprzejmość na estymator w modelu : słusznie oczekujemy, że będzie „bliski jedności”. $B$ $A$ $\hat \beta$

Oczywiste jest, że jeśli tak się stanie, że , ilość po prawej stronie nierówności będzie miała tendencję do zwiększania się bez ograniczeń, ponieważ , liczba kroków do prognozy, wzrośnie. Z drugiej strony ilość po lewej stronie pożądanej nierówności może wzrosnąć wraz ze wzrostem , ale ma górną granicę . Więc w tym scenariuszu oczekujemy differenced modelu Fair lepsze pod względem przewidywania MSE w porównaniu do modelu . $\hat \beta >1$ $k$ $k$ $B$ $A$

Załóżmy jednak, że bardziej korzystny jest przypadek modelu , gdzie . Wówczas ilość po prawej stronie również jest związana. Następnie jako musimy sprawdzić, czy $A$ $\hat \beta <1$ $k \rightarrow \infty$

mi [{(\frac{\hat{γ}}{1 - \hat{γ}})}^{2)} u_{T.}^{2)}] \leq mi [y_{T.}^{2)}] = T. σ_{u}^{2)} ? ?

$E\left[\left(\frac {\hat \gamma}{1-\hat \gamma}\right)^2u_T^2\right] \leq E\big[y_T^2\big]= T\sigma^2_u\;\; ??$

( jest wygodą - w rzeczywistości obie wielkości będą już bliskie swojej supremie już dla małych wartości ). $k \rightarrow \infty$ $k$

Zauważ, że termin powinien być „raczej bliski” , więc model ma przewagę nad tym aspektem. $\left(\frac {\hat \gamma }{1-\hat \gamma}\right)^2$ $0$ $B$

Nie możemy oddzielić pozostałej oczekiwanej wartości, ponieważ estymator nie jest niezależny od . Ale możemy przekształcić nierówność w $\hat \gamma$ $u_T$

Cov [{(\frac{\hat{γ}}{1 - \hat{γ}})}^{2)}, u_{T.}^{2)}] + mi [{(\frac{\hat{γ}}{1 - \hat{γ}})}^{2)}] \cdot σ_{u}^{2)} \leq T. σ_{u}^{2)} ? ?

$\operatorname{Cov}\left[\left(\frac {\hat \gamma}{1-\hat \gamma}\right)^2,\,u_T^2\right] + E\left[\left(\frac {\hat \gamma}{1-\hat \gamma}\right)^2\right]\cdot \sigma^2_u \leq T\sigma^2_u\;\; ??$

\Rightarrow Cov [{(\frac{\hat{γ}}{1 - \hat{γ}})}^{2)}, u_{T.}^{2)}] \leq (T. - mi [{(\frac{\hat{γ}}{1 - \hat{γ}})}^{2)}]) \cdot σ_{u}^{2)} ? ?

$\Rightarrow \operatorname{Cov}\left[\left(\frac {\hat \gamma}{1-\hat \gamma}\right)^2,\,u_T^2\right] \leq \left (T-E\left[\left(\frac {\hat \gamma}{1-\hat \gamma}\right)^2\right]\right)\cdot \sigma^2_u \;\; ??$

Oczekuje się, że kowariancja po lewej stronie będzie niewielka, ponieważ estymator zależy od wszystkich błędówPo drugiej stronie nierówności pochodzi ze stacjonarnego zestawu danych, więc oczekuje się, że oczekiwana wartość powyższej funkcji będzie znacznie mniejsza niż wielkość próbki (ponieważ więcej w tej funkcji będzie się wahać w ). $\hat \gamma$ $T$ $\hat \gamma$ $(0,1)$

Podsumowując, nie dyskutując o żadnej konkretnej metodzie szacowania, uważam, że byliśmy w stanie nieformalnie wykazać, że zróżnicowany model powinien działać lepiej pod względem przewidywania MSE.

Alecos Papadopoulos
źródło

To jest dobre pytanie.

Jak rozumiem, właśnie zastanawiałeś się nad pacf, ale to nie wystarczy. Zarówno ACF, jak i PACF są niezbędne, aby wybrać najlepszy model.

Z drugiej strony testy stacjonarne są słabe i wrażliwe i wymagają przetestowania dużej liczby opóźnień.

Ponadto zaleca się, aby szeregi czasowe były nieruchome przed zastosowaniem dowolnego modelu. Z grubsza mówiąc, modele ARIMA biorą pod uwagę szczególny przypadek niestacjonarności (najlepiej w modzie).

Co do twoich pytań, nie jestem pewien co do funkcji auto.arima, ale jestem pewien, że liczba punktów danych w twoim przykładzie jest niewielka. Symulacja modelu przy użyciu dużej liczby punktów danych dobrze odpowie na twoje pytania. Radzę również rozważyć ACF szeregów czasowych, a także PACF. Przy wyborze modelu ogólną zasadą jest wybór najprostszego modelu (zwróć uwagę, że najprostszy model po unieruchomieniu szeregów czasowych).

Odsyłam cię do tego odniesienia. Ta książka nie odpowiada na wszystkie twoje pytania, ale daje pewne wskazówki.

----- sekcja uzupełniająca ------- @nsw biorąc pod uwagę trend w twoich danych. Jeśli weźmiesz pod uwagę model stacjonarny, spowoduje to prognozę w górę / w dół, ale w rzeczywistości modele ARMA są zaprojektowane do przewidywania płaskich danych. Zmieniłem twój kod, aby odzwierciedlić tę różnicę:

wymagają („prognozy”)

wymagają („serów”)

controlData <- arima.sim (lista (kolejność = c (1,1,1), ar = .5, ma = .5), n = 1000 )

acf (controlData)

ts.plot (controlData)

naiveFit <- arima (controlData, order = c (2,0,1))

trueFit <- arima (controlData, order = c (1,1,1))

PrnaiveFit <-forecast.Arima (naiveFit, 10)

PrtrueFit <- prognozy.Arima (trueFit, 10)

matplot (cbind (PrnaiveFit $ mean, PrtrueFit $ mean), type = 'b', col = c ('red', 'green'), ylab = c ('przewidzieć jon'), pch = c ('n', „t”))

TPArrow
źródło

Pytanie dotyczy tego , dlaczego preferuje się „uszeregowanie szeregów czasowych”. To tak naprawdę nie odpowiada na to pytanie.

shadowtalker

@ssdecontrol Masz ogólnie rację. Naprawdę bardziej martwię się dorozumianymi konsekwencjami dla przewidywania po błędnym sprecyzowaniu. Ale nie chcę za bardzo walczyć z Hamedem. Nadal odpowiedział na moje ostatnie pytanie: „Czy to właściwy sposób na wybór modelu?”. Ale żeby to jeszcze raz powtórzyć, to mój najmniejszy problem.

Clark Henry