Pojęcie procesu stacjonarnego

Pojęcie procesu stacjonarnego

Definicja 1. Procesem stochastycznym nazywamy rodzinę zmiennych (wektorów) losowych {Xt}t∈T okre- ślonych na tej samej przestrzeni probabilistycznej.

Jeśli T ⊆ R, to indeks numerujący zmienne losowe często interpretuje się jako czas. Tak będzie w naszych rozważaniach. Dalej ograniczymy się do przypadku, gdy T = Z.

Definicja 2. Proces stochastyczny {X_t}_t∈Z nazywamy stacjonarnym (w węższym sensie), jeśli

∀ t₁, t₂, . . . , t_k∈ Z ∀ h ∈ Z (X_t1+h, X_t2+h, . . . , X_tk+h)= (X^D _t1, X_t2, . . . , X_tk).

Przykładami procesów stacjonarnych są stacjonarne łańcuchy Markowa tzn. łańcuchy Markowa ma- jące rozkład stacjonarny, dla których rozkładem początkowym jest ich rozkład stacjonarny.

Definicja 3. Proces stochastyczny {Xt}_t∈Z nazywamy stacjonarnym w szerszym sensie, gdy (i) dla każdego t ∈ Z istnieje EXt² (tzn. gdy EX_t²< ∞),

(ii) dla dowolnych s, t ∈ Z zachodzi EX^s= EXt (tzn. wartość oczekiwana procesu jest stała w czasie), (iii) dla dowolnych s, t, h ∈ Z zachodzi Cov(X^s, Xs+h) = Cov(Xt, Xt+h) (tzn. kowariancja procesu w

dwóch punktach czasu zależy jedynie od odległości między nimi).

Fakt 1. Jeśli proces {Xt}_t∈Z jest stacjonarny w szerszym sensie, to dla dowolnych s, t ∈ Z mamy V ar(Xs) = V ar(Xt) (tzn. wariancja procesu jest stała w czasie).

Dowód. Wystarczy w punkcie (iii) definicji procesu stacjonarnego w szerszym sensie przyjąć h = 0.

Fakt 2.

(i) Proces {Xt} stacjonarny w węższym sensie, dla którego istnieje EX_t², jest procesem stacjonarnym w szerszym sensie.

(ii) Proces {X_t} stacjonarny w szerszym sensie, który jest gaussowski (tzn. taki że dla dowolnych t₁, t₂, . . . , t_k ∈ Z wektor losowy (Xt1, X_t2, . . . , X_tk) ma [wielowymiarowy] rozkład normalny), jest procesem stacjonarnym w węższym sensie.

Punkt (ii) jest konsekwencją faktu, iż wielowymiarowy rozkład normalny jest w całości opisywany za pomocą wektora wartości oczekiwanych i macierzy kowariancji.

Od tego momentu określenie proces stacjonarny zawsze będzie się odnosiło do procesu stacjonarnego w szerszym sensie.

Skoro dla procesu stacjonarnego kowariancja procesu w dwóch punktach czasu zależy jedynie od odległości między nimi, to dla badania takiego procesu wprowadza się następującą wielkość:

Definicja 4. Niech {Xt}_t∈Z będzie procesem stacjonarnym. Funkcję γ : Z → R daną wzorem γ(h) = Cov(Xt, Xt+h) = Cov(X0, Xh)

nazywamy funkcją autokowariancji procesu {X_t}_t∈Z. γ(h) będziemy nazywali autokowariancją rzędu h.

Własności funkcji autokowariancji i autokorelacji

Twierdzenie 1 (Własności funkcji autokowariancji). Funkcja autokowariancji γ procesu stacjonarnego ma następujące własności:

(i) γ(0) ≥ 0,

(ii) |γ(h)| ≤ γ(0) dla każdego h ∈ Z, (iii) γ(−h) = γ(h) dla każdego h ∈ Z.

Dowód. (i) γ(0) = Cov(X0, X0) = V ar(X0) ≥ 0.

(ii) wynika z nierówności Schwarza: dla dowolnych zmiennych losowych X i Y takich, że EX² i EY² istnieją, mamy

|Cov(X, Y )| ≤p

V ar(X) · V ar(Y ).

(iii) γ(−h) = Cov(X0, X_−h) = Cov(X_−h, X0) = Cov(X_−h+h, X0+h) = Cov(X0, Xh) = γ(h).

Definicja 5. Niech {Xt}_t∈Z będzie procesem stacjonarnym. Funkcję ρ : Z → R daną wzorem ρ(h) = Corr(X₀, X_h) = Cov(X₀, X_h)

pV ar(X0) · V ar(Xh) =Cov(X₀, X_h)

V ar(X0) =γ(h) γ(0) nazywamy funkcją autokorelacji procesu {Xt}_t∈Z.

ρ(h) będziemy nazywali autokorelacją rzędu h.

Twierdzenie 2 (Własności funkcji autokorelacji). Funkcja autokorelacji ρ procesu stacjonarnego ma następujące własności:

(i) ρ(0) = 1,

(ii) |ρ(h)| ≤ 1 dla każdego h ∈ Z, (iii) ρ(−h) = ρ(h) dla każdego h ∈ Z.

Procesy stacjonarne będziemy utożsamiali z ich wartością oczekiwaną i funkcją autokowiariancji. In- nymi słowy: dwa procesy stacjonarne z taką samą wartością oczekiwaną i taką samą funkcją autokowa- riancji będą dla nas nierozróżnialne.

Definicja 6. Funkcję κ : Z → R nazywamy nieujemnie określoną, jeżeli dla każdego n ∈ N⁺, dla dowol- nych a1, a2, . . . , an∈ R i dowolnych t1, t2, . . . , tn∈ Z zachodzi

n

X

i=1 n

X

j=1

aiκ(ti− tj)aj ≥ 0

tzn. gdy dla każdego n ∈ N⁺ i dowolnych t1, t2, . . . , tn ∈ Z macierz M = (mij)ⁿ_i,j=1, gdzie mij = κ(ti−tj), jest nieujemnie określona.

Twierdzenie 3 (Charakteryzacja funkcji autokowariancji). Parzysta funkcja γ : Z → R jest funkcją autokowariancji pewnego procesu stacjonarnego wtedy i tylko wtedy, gdy jest nieujemnie określona.

Dowód. Załóżmy, że funkcja γ : Z → R jest funkcją autokowariancji pewnego procesu stacjonarnego {X_t}_t∈Z. Pokażemy, że jest nieujemnie określona. Weźmy dowolne t₁, t₂, . . . , t_n ∈ Z. Wówczas γ(ti − t_j) = Cov(X₀, X_ti−tj) = Cov(X_0+tj, X_(ti−tj)+tj) = Cov(X_tj, X_ti). Z tego wynika, że macierz Γ = (γ(t_i− t_j))ⁿ_i.j=1 jest macierzą kowariancji wektora losowego (X_t1, X_t2, . . . , X_tn)⁰, a zatem jest nieujemnie określona. W takim razie funkcja γ jest nieujemnie określona.

Z kolei jeśli mamy parzystą i nieujemnie określoną funkcję γ : Z → R, to zawsze można skonstruować proces gaussowski o wartości oczekiwanej 0 takiej funkcji autokowiariancji, gdyż wielowymiarowy rozkład normalny jest w całości opisywany za pomocą wektora wartości oczekiwanych i macierzy kowariancji.

Aby pokazać, że jakaś funkcja jest funkcją autokowariancji pewnego procesu, często zamiast poka- zywać, że jest nieujemnie określona, łatwiej jest skonstruować proces stacjonarny, którego jest funkcją autokowiariancji.

Przykład 1

Funkcja γ(h) = σ²cos θh jest funkcją autokowariancji procesu stacjonarnego {X_t} postaci:

X_t= A cos θt + B sin θt,

gdzie A i B są nieskorelowanymi zmiennymi losowymi takimi że EA = EB = 0 i V ar(A) = V ar(B) = σ², gdzie σ²jest ustaloną liczbą dodatnią zaś θ ustaloną liczbą rzeczywistą.

Aby się o tym przekonać, upewnijmy się najpierw, że proces {X_t} zdefiniowany j.w. jest stacjonarny.

Ponieważ EA = EB = 0, więc EA² = V ar(A) + (EA)² = V ar(A) = σ² i EB² = V ar(B) + (EB)² = V ar(B) = σ². Mamy kolejno:

EX_t² = E(A cos θt + B sin θt)²= E(A²cos²θt + AB sin 2θt + B²sin²θt) =

= cos²θtEA²+ sin 2θtEAEB + sin²θtEB²=

= cos²θt · σ²+ sin²θt · σ²= (cos²θt + sin²θt)σ²= σ²< ∞,

EXt= E(A cos θt + B sin θt) = cos θtEA + sin θtEB = cos θt · 0 + sin θt · 0 = 0 = const, Cov(Xt, Xt+h) =

= E(Xt− EXt)(Xt+h− EXt+h) = EXtXt+h=

= E(A cos θt + B sin θt)(A cos θ(t + h) + B sin θ(t + h)) =

= E[A²cos θt cos θ(t + h) + AB(cos θt cos θ(t + h) + sin θt cos θ(t + h)) + B²sin θ(t + h) sin θt] =

= cos θt cos θ(t + h)EA²+ EAEB(cos θt cos θ(t + h) + sin θt cos θ(t + h)) + sin θ(t + h) sin θtEB²=

= cos θt cos θ(t + h) · σ²+ sin θ(t + h) sin θt · σ²= [cos θt cos θ(t + h) + sin θ(t + h) sin θt]σ²=

= cos(θ(t + h) − θt) · σ²= σ²cos θh

Przykłady procesów stacjonarnych

Przykład 2

Ciąg nieskorelowanych zmiennych losowych {Zt}_t∈Z o tej samej wartości oczekiwanej µ i tej samej wa- riancji σ² nazywamy białym szumem i oznaczamy: {Zt} ∼ W N (µ, σ²). Funkcje autokowariancji i autokorelacji białego szumu W N (µ, σ²) mają postać:

γ(h) =







σ², h = 0 0, h 6= 0

, ρ(h) =







1, h = 0 0, h 6= 0

.

Biały szum to cegiełki, z których zbudowanych jest wiele procesów stacjonarnych.

Przykład 3

Ciąg niezależnych zmiennych losowych {Z_t}_t∈Z o jednakowym rozkładzie o wartości oczekiwanej µ i wariancji σ² nazywamy IID-szumem i oznaczamy: {Zt} ∼ IID(µ, σ²). Funkcje autokowariancji i autokorelacji IID-szumu są takie same jak funkcje autokowariancji i autokorelacji białego szumu.

IID-szum jest białym szumem, ale niekoniecznie biały szum jest IID-szumem.

Jeśli {Zt} ∼ IID(µ, σ²) i Zt∼ N (µ, σ²), t ∈ Z, to {Z^t}_t∈Znazywamy gaussowskim białym szumem.

Przykład 4

Proces {Xt}_t∈Zpostaci:

X_t= Z_t+ θZ_t−1, gdzie {Zt} ∼ W N (0, σ²),

przy czym θ ∈ R i σ²> 0 są ustalonymi liczbami, nazywamy procesem średniej ruchomej rzędu 1 i oznaczamy: {X_t} ∼ MA(1). Ponieważ EZ_t = 0, więc EZ_t² = V ar(Z_t) + (EZ_t)² = V ar(Z_t) = σ². W takim razie

EX_t² = E(Zt+ θZt−1)²= E(Z_t²+ 2θZtZt−1+ θ²Z_t−1² ) = EZ_t²+ 2θEZtEZt−1+ θ²EZ_t−1² =

= σ²+ θ²σ²= (1 + θ²)σ²< ∞,

EX_t= E(Z_t+ θZ_t−1) = EZ_t+ θEZ_t−1= 0,

Cov(X_t, X_t+h) = E(X_t− EXt)(X_t+h− EXt+h) = EX_tX_t+h= E(Z_t+ θZ_t−1)(Z_t+h+ θZ_t+h−1) =

= EZtZt+h+ θEZtZt+h−1+ θEZt−1Zt+h+ θ²EZt−1Zt+h−1. Dla h = 0 mamy:

Cov(Xt, Xt) = EZ_t²+ 2θEZtEZt−1+ θ²EZ_t−1² = (1 + θ²)σ². Dla h = 1 mamy:

Cov(X_t, X_t+1) = EZ_tEZ_t+1+ θEZ_t²+ θEZ_t−1EZ_t+h+ θ²EZ_t−1EZ_t= θσ². Dla h > 1 mamy:

Cov(Xt, Xt+h) = EZtEZt+h+ θEZtEZt+h−1+ θEZt−1EZt+h+ θ²EZt−1EZt+h−1= 0.

Istnieje drugi moment procesu {Xt}, wartość oczekiwana jest stała w czasie a kowariancja procesu w dwóch punktach czasu zależy jedynie od odległości miedzy nimi, a zatem proces {Xt} jest stacjonarny.

Jego funkcje autokowariancji i autokorelacji wyrażają się za pomocą następujących wzorów:

γ(h) =











(1 + θ²)σ², h = 0 θ²σ², |h| = 1

0, |h| > 1

, ρ(h) =











1, h = 0 θ²

1 + θ², |h| = 1 0, |h| > 1

(Uzupełnienia dla h < 0 dokonujemy w oparciu o parzystość funkcji autokowariancji i autokorelacji).

Przykład 5

Proces {Xt}_t∈Zpostaci:

Xt= Zt+ θ1Zt−1+ θ2Zt−2+ . . . + θqZt−q, gdzie {Zt} ∼ W N (0, σ²),

przy czym θ1, θ2, . . . , θq ∈ R i σ²> 0 są ustalonymi liczbami, nazywamy procesem średniej ruchomej rzędu q i oznaczamy: {Xt} ∼ MA(q). Jest on procesem stacjonarnym a jego funkcje autokowariancji i autokorelacji wyrażają się następująco:

γ(h) =











σ²Pq

j=0θ_j², h = 0 σ²Pq−|h|

j=0 θjθj+h, 1 ≤ |h| ≤ q 0, |h| > q

, ρ(h) =















1, h = 0

Pq−|h|

i=0 θ_iθ_i+h Pq

j=0θ²_j , 1 ≤ |h| ≤ q 0, |h| > q

Przykład 6

Jeśli proces stacjonarny {Xt}_t∈Zjest postaci:

Xt= φXt−1+ Zt, gdzie{Zt} ∼ W N (0, σ²),

przy czym φ ∈ R i σ² > 0 są ustalonymi liczbami, to nazywamy go procesem autoregresji rzędu 1 i oznaczamy: {Xt} ∼ AR(1). Można pokazać, że jeśli |φ| < 1, to proces spełniający powyższą zależność rekurencyjną jest stacjonarny. Wówczas jego funkcje autokowariancji i autokorelacji wyrażają się jako:

γ(h) = σ²

1 − φ² φ^|h|, ρ(h) = φ^|h|.

Funkcja autokorelacji cząstkowej

Definicja 7. Funkcją autokorelacji cząstkowej procesu stacjonarnego {Xt}_t∈Z nazywamy funkcję α : Z → R zdefiniowaną w następujący sposób:

α(h) =







0, h = 0 φh,h, h 6= 0

,

gdzie φh= (φh,1, φh,2, . . . , φh,h)⁰ jest rozwiązaniem równania Γhφh= γh, w którym Γh= (γ(i − j))^h_i,j=1 i γh= (γ(1), γ(2), . . . , γ(h))⁰.

α(h) będziemy nazywali autokorelacją cząstkową rzędu h.

Intuicyjnie autokorelacja cząstkowa rzędu h mierzy stopień liniowej zależności między X0 i Xh po wyeliminowaniu liniowego wpływu zmiennych losowych X1, X2, . . . , Xh−1 dla h > 0 lub zmiennych loso- wych Xh+1, Xh+2, . . . , X₋₁ dla h < 0.

Procesy liniowe

Dużą klasę procesów stacjonarnych można zbudować z kombinacji liniowych białego szumu. Procesy te nazywają się procesami liniowymi.

Definicja 8. Proces stochastyczny {Xt}_t∈Z postaci

Xt=

∞

X

j=−∞

ψjZt−j,

gdzie {Zt} ∼ W N (0, σ²) iP∞

j=−∞|ψj| < ∞, nazywamy procesem liniowym. Zbieżność w definicji (jeśli faktycznie mamy do czynienia z nieskończoną liczbą składników) będziemy interpretowali jako zbieżność z prawdopodobieństwem 1.

W tym miejscu należałoby wykazać, że warunek P∞

j=−∞|ψj| < ∞ jest wystarczający do tego, aby zbieżnośćP∞

j=−∞ψjZt−j zachodziła z prawdopodobieństwem 1, jednak pominiemy ten dowód.

Twierdzenie 4. Proces liniowy przy oznaczeniach jak w powyższej definicji jest procesem stacjonarnym o wartości oczekiwanej 0 i funkcji autokowariancji postaci:

γ(h) = σ²

∞

X

j=−∞

ψ_jψ_j+h.

Dowód. Zanim przystąpimy do właściwego dowodu, zwróćmy uwagę, że E|Z_t| < ∞ i EZ_t²= V ar(Z_t) + (EZ_t)²= V ar(Z_t) + 0 = V ar(Z_t) = σ². Co więcej, ponieważ 2ab ≤ a²+ b², więc dla dowolnych k, l ∈ Z mamy 2E|Zk||Zl| ≤ E(Z_k²+ Z_l²) = EZ_k²+ EZ_l²= σ²+ σ²= 2σ², a zatem E|Zk||Zl| ≤ σ².

Najpierw pokażemy, że wartość oczekiwana Xtistnieje:

E|X_t| = E      

∞

X

j=−∞

ψ_jZ_t−j      

≤ E

∞

X

j=−∞

|ψ_j| · |Z_t−j|≤^?

∞

X

j=−∞

|ψ_j| · E|Z_t−j| = E|Z₀| ·

∞

X

j=−∞

|ψ_j| < ∞.

W miejscu oznaczonym gwiazdką skorzystaliśmy z lematu Fatou. W takim razie możemy skorzystać z twierdzenia Lebesgue’a o zbieżności ograniczonej (zmajoryzowanej):

EX_t= E

∞

X

j=−∞

ψ_jZ_t−j=

∞

X

j=−∞

ψ_j· EZ_t−j =

∞

X

j=−∞

ψ_j· 0 = 0.

W następnej kolejności wykażemy, że EX_t²istnieje:

EX_t²= E|X_t²| = E      

∞

X

i=−∞

ψiZt−i

!

·

∞

X

j=−∞

ψjZt−j

!     

= E      

∞

X

i=−∞

∞

X

j=−∞

ψiψjZt−iZt−j

≤

≤ E

∞

X

i=−∞

∞

X

j=−∞

|ψi||ψj||Zt−i||Zt−j|≤^?

∞

X

i=−∞

∞

X

j=−∞

|ψi||ψj| · E|Zt−i||Zt−j| ≤

≤

∞

X

i=−∞

∞

X

j=−∞

|ψi||ψj| · σ²= σ²

∞

X

i=−∞

|ψi|

!

·

∞

X

j=−∞

|ψj|

!

< ∞.

W miejscu oznaczonym gwiazdką ponownie skorzystaliśmy z lematu Fatou.

Skoro istnieje drugi moment zmiennej losowej Xt, to w poniższym rachunku możemy zamienić kolej- ność wartości oczekiwanej i sumy, ponownie korzystając z twierdzenia Lebesgue’a o zbieżności ograniczo- nej (zmajoryzowanej):

Cov(X_t, X_t+h) = E(X_t− EX_t

| {z }

0

)(X_t+h− EX_t+h

| {z }

0

) = EX_tX_t+h=

= E

∞

X

i=−∞

ψiZt−i

!

·

∞

X

j=−∞

ψjZt+h−j

!

= E

∞

X

i=−∞

∞

X

j=−∞

ψiψjZt−iZt+h−j=

=

∞

X

i=−∞

∞

X

j=−∞

ψiψj· EZt−iZt+h−j=X X

−i=h−j

ψiψj· EZt−iZt+h−j+X X

−i6=h−j

ψiψj· EZt−iZt+h−j

| {z }

=EZ_t−iEZ_t+h−j=0

=

=

∞

X

i=−∞

ψ_iψ_i+h· EZ_t−i² =

∞

X

i=−∞

ψ_iψ_i+h· σ²= σ²

∞

X

i=−∞

ψ_iψ_i+h.

Powyższe twierdzenie jest szczególnym przypadkiem nieco ogólniejszego faktu, który można udowodnić w sposób analogiczny do powyższego.

Twierdzenie 5. Niech {X_t}_t∈Z będzie procesem stacjonarnym o wartości oczekiwanej µ_X i funkcji au- tokowariancji γ_X. Wówczas proces {Y_t}_t∈Z zdefiniowany jako

Y_t=

∞

X

j=−∞

ψ_jX_t−j

jest procesem stacjonarnym o wartości oczekiwanej µY i funkcji autokowariancji γY zadymi jako

µY = µX

∞

X

j=−∞

ψj i γY(h) =

∞

X

i=−∞

∞

X

j=−∞

ψiψjγX(i − j + h).

Procesy ARMA

Definicja 9. Procesem autoregresji rzędu p i średniej ruchomej rzędu q (w skrócie: ARMA(p,q)) nazywamy proces stacjonarny {Xt}_t∈Z będący rozwiązaniem równania

X_t= φ₁X_t−1+φ₂X_t−2+. . .+φ_pX_t−p+Z_t+θ₁Z_t−1+θ₂Z_t−2+. . .+θ_qZ_t−q, gdzie {Z_t} ∼ W N (0, σ²), (1) przy czym φ1, φ₂, . . . , φ_p, θ₁, θ₂, . . . , θ_q ∈ R i σ² > 0 są ustalonymi liczbami (co oznaczamy: {X_t} ∼ ARM A(p, q)).

Równanie (1) można zapisać jako:

Xt−φ1Xt−1−φ2Xt−2−. . .−φpXt−p= Zt+θ1Zt−1+θ2Zt−2+. . .+θqZt−q, gdzie {Zt} ∼ W N (0, σ²). (2) Niech B oznacza operator cofnięcia o 1 tzn. Bxt= xt−1, B²xt= xt−2, B³xt= xt−3 itd. Niech

Φ(s) = 1 − φ₁s − φ₂s²− . . . − φ_ps^p i Θ(s) = 1 + θ₁s + θ₂s²+ . . . + θ_ps^q. Korzystając z tych oznaczeń, możemy zapisać równanie (2) jako:

Φ(B)Xt= Θ(B)Zt, gdzie {Zt} ∼ W N (0, σ²). (3) Wielomiany Φ i Θ noszą nazwy odpowiednio wielomianu autoregresji i wielomianu średniej ruchomej.

Stopień wielomianu autoregresji wynosi p zaś stopień wielomianu średniej ruchomej wynosi q.

Szczególnymi przykładami procesów ARMA są:

ˆ proces autoregresji (ang. autoregression) rzędu p zdefiniowany jako proces stacjonarny będący roz- wiązaniem równania:

X_t= φ₁X_t−1+ φ₂X_t−2+ . . . + φ_pX_t−p+ Z_t czyli Φ(B)X_t= Z_t, gdzie {Z_t} ∼ W N (0, σ²), tzn. proces ARMA(p,0), który oznaczamy krótko: AR(p),

ˆ proces średniej ruchomej (ang. moving average) rzędu q zdefiniowany jako proces będący rozwiąza- niem równania:

Xt= Zt+ θ1Zt−1+ θ2Zt−2+ . . . + θqZt−q czyli Xt= Θ(B)Zt, gdzie {Zt} ∼ W N (0, σ²), (4) tzn. proces ARMA(0,q), który oznaczamy krótko: MA(q).

Należy przypomnieć, że proces stanowiący rozwiązanie równania (4) zawsze jest procesem stacjonar- nym. Ponadto z samej postaci równania (4) wynika, że proces średniej ruchomej jest procesem liniowym.

Przypomnijmy również, że dla procesu MA(q) mamy ρ(h) = 0 dla |h| > q.

Z kolei dla procesu AR(p) mamy α(h) = 0 dla |h| > p.

W tym momencie należy postawić trzy ważne pytania:

ˆ Problem stacjonarności: Czy dowolny proces będący rozwiązaniem równania (1) jest procesem sta- cjonarnym?

ˆ Problem wynikowości: Czy rozwiązanie {Xt}_t∈Z równania (1) jest skorelowane tylko z przeszłością białego szumu a nie z jego przyszłością tzn. czy Cov(Xt, Zs) = 0 dla s > t?

ˆ Problem odwracalności: Czy jeśli {Xt}_t∈Zjest rozwiązaniem równania (1), to Ztmożna przedstawić jako liniową funkcję Xt, Xt−1, Xt−2, . . . ?

Proces spełniający drugi postulat nazywamy procesem wynikowym. Proces spełniający trzeci postulat nazywamy procesem odwracalnym.

Twierdzenie 6.

(i) Jeżeli wszystkie pierwiastki wielomianu Φ leżą poza okręgiem |s| = 1, to pewien proces liniowy {X_t}_t∈Z postaci

X_t=

∞

X

j=−∞

ψ_jZ_t−j, gdzie {Z_t} ∼ W N (0, σ²), (5)

jest jedynym stacjonarnym rozwiązaniem równania (1).

(ii) Jeżeli wszystkie pierwiastki wielomianu Φ leżą poza kołem |s| ≤ 1, to stacjonarne rozwiązanie równania (1) jest procesem wynikowym.

(iii) Jeżeli wszystkie pierwiastki wielomianu Θ leżą poza kołem |s| ≤ 1, to proces {Xt}_t∈Z zdefiniowany w (5) jest odwracalny.

Charakterystyki próbkowe procesów stacjonarnych

Zakładamy, że proces stacjonarny {Xt}_t∈Z obserwujemy w chwilach t = 1, 2, . . . , n tzn. mamy dane zmienne losowe X1, X2, . . . , Xn. Realizację procesu stochastycznego indeksowanego dyskretnym czasem nazywamy szeregiem czasowym. Możemy więc powiedzieć, że mamy dany stacjonarny szereg czasowy X1, X2, . . . , Xn.

Średnią próbkową procesu będziemy nazywali statystykę

X = 1 n

n

X

j=1

Xj.

Fakt 3. X jest nieobciążonym estymatorem wartości oczekiwanej procesu stacjonarnego.

Autokowariancją próbkową rzędu h będziemy nazywali statystykę

ˆ

γ(h) = 1 n

n−h

X

j=1

(X_j− X)(X_j+h− X), 0 ≤ h ≤ n,

ˆ

γ(h) = ˆγ(−h), −n ≤ h ≤ 0.

Autokorejacją próbkową rzędu h będziemy nazywali statystykę

ˆ

ρ(h) = ˆγ(h) ˆ

γ(0), −n ≤ h ≤ n.

Losowe funkcje ˆγ : {−n, . . . , −2, −1, 0, 1, 2, . . . , n} → R i ˆρ : Z → R będziemy nazywali odpowiednio próbkową funkcją autokowariancji i próbkową funkcją autokorelacji.

Zwróćmy uwagę, że jeśli h jest bliskie n lub −n, to suma występująca w definicji ˆγ(h) składa się z niewielu składników, więc można się spodziewać słabej jakości estymacji γ(h) i ρ(h) poprzez ˆρ(h) i ˆρ(h), dlatego też nie wykorzystuje się ˆγ(h) i ˆρ(h) dla h bliskich n lub −n, choć formalnie można je obliczyć.

W literaturze można się także spotkać z dzieleniem przez n − h w miejsce dzielenia przez n przy obliczaniu autokowariancji próbkowej rzędu h, jednak dzielenie przez n zapewnia nieujemną określoność funkcji ˆγ.

Definicje autokowariancji próbkowej i autokorelacji próbkowej zachowują sens także w przypadku procesów niestacjonarnych. Wielkości te mogą być użyte do wykrywania niestacjonarności oraz możliwych jej przyczyn.

Twierdzenie 7. Niech {X_t}_t∈Z będzie procesem stacjonarnym takim że X_t = P∞

j=−∞ψ_jZ_t−j, gdzie {Zt} ∼ IID(0, σ²) orazP∞

j=−∞|ψj| < ∞ i EZ_j⁴< ∞ (lub P∞

j=−∞ψ_j²|j| < ∞). Ustalmy k ≥ 1. Niech ρ_k = (ρ(1), ρ(2), . . . , ρ(k)) i ˆρ_k = ( ˆρ(1), ˆρ(2), . . . , ˆρ(k)). Wówczas

√n( ˆρk− ρk)−−−−→^D

n→∞ N (0, W ), gdzie W = (w_ij)^k_i,j=1 i w_ij =P∞

l=1[ρ(l + i) + ρ(l − i) − 2ρ(l)ρ(i)] · [ρ(l + j) + ρ(l − j) − 2ρ(l)ρ(j)].

Wniosek 1. Niech {Xt}_∼IID(0, σ²). Ustalmy k ≥ 1. Niech ρk= (ρ(1), ρ(2), . . . , ρ(k)) i ˆρk= ( ˆρ(1), ˆρ(2), . . . , ˆρ(k)).

Wówczas √

n( ˆρk− ρk)−−−−→^D

n→∞ N (0, I).

Wniosek 2. Niech {Xt} ∼ IID(0, σ²). Wówczas dla każdego h ≥ 1

√n ˆρ(h)−−−−→^D

n→∞ N (0, 1).

W szczególności dla każdego α ∈ (0,¹₂) mamy

n→∞lim P√

n | ˆρ(h)| ≥ Φ⁻¹ 1 −α

2



= lim

n→∞P



|ˆρ(h)| ≥ 1

√nΦ⁻¹ 1 −α

2



= α.

Powyższy fakt jest wykorzystywany przy sporządzaniu standardowych wykresów próbkowej funkcji autokorelacji. Za pomocą słupków są na nich zaznaczone liczby ˆρ(h), natomiast promień zaznaczonego pasa wokół zera wynosi ^√1_nΦ⁻¹ 1 − ^α₂
. Z pomocą takiego rysunku łatwo można przeprowadzić test hipotezy

H : ρ(h) = 0 przeciwko K : ρ(h) 6= 0.

Jeśli słupek wychodzi poza zaznaczony pas, to przyjmujemy K; jeśli nie, to stwierdzamy, że brak jest podstaw do odrzucenia H na rzecz K.

Cząstkową autokorelacją próbkową rzędu h będziemy nazywali statystykę zdefiniowaną w następujący sposób:

ˆ α(h) =







0, h = 0 φ_h,h, 1 ≤ |h| ≤ n

,

gdzie φh= (φh,1, φh,2, . . . , φh,h)⁰ jest rozwiązaniem równania ˆΓhφh= ˆγh, w którym ˆΓh = (ˆγ(i − j))^h_i,j=1 i ˆγh= (ˆγ(1), ˆγ(2), . . . , ˆγ(h))⁰.

Losową funkcję ˆα : {−n, . . . , −2, −1, 0, 1, 2, . . . , n} → R będziemy nazywali próbkową funkcją autoko- relacji cząstkowej.

Podobnie jak w przypadku próbkowej funkcji autokowariancji i próbkowej funkcji autokorelacji esty- macja α(h) poprzez ˆα(h) dla h bliskich n lub −n jest słabej jakości, więc w praktyce nie posługujemy się

ˆ

α(h) dla takich h.

O próbkowej funkcji autokorelacji cząstkowej można udowodnić twierdzenie graniczne podobne do tego dotyczącego próbkowej funkcji autokorelacji. W podobny sposób można je wykorzystać do konstrukcji testu dla testowania hipotezy o tym, że dla ustalonego h zachodzi α(h) = 0 przeciwko hipotezie głoszącej, że α(h) 6= 0.

Estymacja trendu i sezonowości

Będziemy zakładali, że obserwujemy proces stochastyczny {Xt}_t∈Zpostaci:

Xt= mt+ st+ et,

gdzie {mt} nazywamy trendem i jest to funkcja deterministyczna, {st} nazywamy składnikiem sezonowym lub krótko: sezonowością i jest to deterministyczna funkcja okresowa, która ma tę własność, że jeśli T jest jej okresem, toPT

j=1s_j= 0, zaś {e_t} jest procesem stacjonarnym. Model w takiej postaci nazywamy addytywnym modelem szeregu czasowego. Możliwe są modele, w których s_t ≡ 0 (wówczas mówimy o modelu bez składnika sezonowego) lub w których m_t≡ 0 (wówczas mówimy o modelu bez trendu).

Inną spotykaną postacią modelu jest następująca:

Xt= mt· st· et,

gdzie mt, st, et > 0 a ponadto jeśli T jest okresem czynnika sezonowego, to QT

j=1sj = 1. Model taki nazywamy multiplikatywnym modelem szeregu czasowego. Model multiplikatywny łatwo sprowadzić do modelu addytywnego przez zlogarytmowanie, stąd nie będziemy rozważać modelu multiplikatywnego.

Niekiedy bywają rozważane także inne modele szeregów czasowych, ale nie będziemy o nich mówić, bowiem model addytywny (i ewentualnie model multiplikatywnym) wydają się być najczęściej stosowane.

Analizę danych indeksowanych czasem warto rozpocząć od narysowania wykresu ich przebiegu. Jeśli trend jest rosnący i wraz z czasem rośnie też rozrzut danych wokół trendu, warto spróbować ustabilizować ten rozrzut, nakładając na dane przekształcenie z rodziny Boxa-Coxa. (To najczęściej spotykana sytuacja, w której przekształcenia Boxa-Coxa mogą się przydać, jednak nie wyczerpuje ona wszystkich możliwości).

Naszym zadaniem będzie wyznaczenie estymatorów funkcji {m_t} i {st} oraz znalezienie probabili- stycznej struktury procesu {e_t}. Zadanie to nazywany estymacją szeregu czasowego.

Filtry liniowe

Definicja 10. Przekształcenie postaci:

Y_t=

s

X

j=−r

a_jX_t−j, r, s ≥ 0 (6)

nazywamy filtrem liniowym.

W R do nakładania filtrów liniowych na dane służy funkcja filter.

Szczególną uwagę należy zwrócić na dwa filtry liniowe najczęściej używane w estymacji szeregów czasowych.

a) Średnią ruchomą rzędu d nazywamy filtr postaci:

Y_t= 1

2q + 1(X_t−q+ X_t−q+1+ . . . + X_t+q−1+ X_t+q), q < t ≤ n − q, jeśli d = 2q + 1 (czyli d jest nieparzyste) lub

Yt= 1 2q

 1

2Xt−q+ Xt−q+1+ . . . + Xt+q−1+1 2Xt+q



, q < t ≤ n − q, jeśli d = 2q (czyli d jest parzyste).

Średnia ruchoma rzędu d jest zatem filtrem liniowym postaci (6). Jeśli d = 2q + 1, to r = s = q i (a_−q, a_−q+1, . . . , a₀, . . . , a_q−1, a_q) =

 1

2q + 1, 1

2q + 1, . . . , 1

2q + 1, . . . , 1

2q + 1, 1 2q + 1



. Jeśli d = 2q, to r = s = q i (a−q, a−q+1, . . . , a0, . . . , aq−1, aq) =

 1 2 · 2q, 1

2q, . . . , 1

2q, . . . , 1 2q, 1

2 · 2q

 .

Średnia ruchoma nie zniekształca trendu liniowego natomiast średnia ruchoma o rzędzie równym długości okresu składnika sezonowego pozwala usunąć składnik sezonowy tzn. jeśli Xt = mt+ st+ et, mt= at + b i sT +t= st, to

ˆ dla T = 2q + 1 mamy 1

2q + 1(Xt−q+ Xt−q+1+ . . . + Xt+q−1+ Xt+q) =

= 1

2q + 1 h

a ((t − q) + (t − q + 1) + . . . + t + · · · + (t + q − 1) + (t + q))

| {z }

(2q+1)t

+(2q + 1)b+

+ (st−q+ st−q+1+ . . . + st+q−1+ st+q)

| {z }

0

+(et−q+ et−q+1+ . . . + et+q−1+ et+q)i

=

at + b + 1

2q + 1(et−q+ et−q+1+ . . . + et+q−1+ et+q), q < t ≤ n − q,

ˆ dla T = 2q mamy

1 2q

 1

2Xt−q+ Xt−q+1+ . . . + Xt+q−1+1 2Xt+q



=

= 1 2q

"

a 1

2(t − q) + (t − q + 1) + . . . + t + · · · + (t + q − 1) +1 2(t + q)



| {z }

2qt

+2qb+

+ (st−q+ st−q+1+ . . . + st+q−1+ st+q)

| {z }

0

+ 1

2et−q+ et−q+1+ . . . + et+q−1+1 2et+q

#

=

at + b + 1 2q

 1

2et−q+ et−q+1+ . . . + et+q−1+1 2et+q



, q < t ≤ n − q.

W R średnią ruchomą można zrealizować za pomocą wyżej wspomnianej funkcji filter. Jej argument filter zawiera wektor (a−q, a_−q+1, . . . , a₀, . . . , a_q−1, a_q) z wcześniejszego opisu.

b) Różnicowaniem z krokiem h > 0 nazywamy filtr liniowy postaci:

Y_t= X_t− X_t−h, t > h.

Różnicowanie z dowolnym krokiem h > 0 pozwala na usunięcie trendu liniowego tzn. jeśli X_t= m_t+e_t i mt= at + b, to

Xt− Xt−h= a(t − (t + h)) + (et− et−h) = ah + (et− et−h), h < t ≤ n.

Jeśli dodatkowo założymy, że w modelu występuje składnik sezonowy o okresie T , to różnicowanie z krokiem T pozwala na usunięcie składnika sezonowego tzn. jeśli Xt= mt+ et, mt= at + b i sT +t = st, to

Xt− Xt−T = a(t − (t − T )) + (st− st−T)

| {z }

0

+(et− et−T) = aT + (et− et−T), T < t ≤ n.

Różnicowanie krokiem 1 ma bardzo klarowną interpretację: pozwala obliczyć przyrost obserwowanego zjawiska między kolejnymi momentami pomiaru np. jeśli dane są gromadzone w odstępach jednego dnia, różnicowanie krokiem 1 pozwala obliczyć przyrost dzienny. Również różnicowanie danych o charakterze okresowym krokiem równym długości okresu ma klarowną interpretację, bowiem pozwala obliczyć przy- rost obserwowanego zjawiska w danym momencie okresu z jednego okresu na drugi np. jeśli gromadzone są dane miesięczne i zjawisko wykazuje roczną okresowość, różnicowanie krokiem 12 pozwala obliczyć przyrost obserwowanego zjawiska w danym miesiącu rok do roku.

W R różnicowanie można przeprowadzić za pomocą funkcji diff. Argument lag oznacza krok (do- myślnie 1), natomiast jeśli chcemy różnicowanie złożyć kilka razy ze sobą, liczbę iteracji różnicowania po- dajemy w argumencie differences (domyślnie 1). Operację odwrotną do różnicowania wykonuje funkcja diffinv o podobnych argumentach.

Stosując średnią ruchomą jak i różnicując proces stacjonarny {et}, otrzymujemy proces stacjonarny (patrz: twierdzenie 5).

Estymacja trendu

Zakładamy, że mamy do czynienia z danymi postaci:

Xt= mt+ st+ et lub Xt= mt+ et. Naszym celem będzie estymacja trendu mt.

T1. Trend możemy dopasować za pomocą metody najmniejszych kwadratów. Dla przykładu dopa- sowanie trendu wielomianowego stopnia k sprowadza się do estymacji za pomocą metody najmniejszych kwadratów parametrów modelu liniowego o wektorze zmiennych objaśnianych i macierzy planu postaci odpowiedni:















 X1

X2

... X_n















 i

















1 1 1² . . . 1^k 1 2 2² . . . 2^k . . . . 1 n n² . . . n^k















 .

Dzięki takiemu podejściu otrzymamy analityczną postać estymatora trendu tzn. wzór a nie tylko wartości. Metoda wymaga jednak wstępnego określenia, wśród jakiej parametrycznej rodziny krzywych poszukujemy estymatora trendu.

Powyższą macierz można łatwo wygenerować za pomocą funkcji poly. Jako argument x podajemy jej zakres czasu, jako degree stopień wielomianu. Oprócz tego przyjmujemy jeszcze raw=TRUE (domyślna wartość do FALSE) – w przeciwnym razie nie otrzymamy wielomianów postaci 1, x, x², x³itd., ale pewien układ wielomianów ortogonalnych.

T2. Jako estymator trendu możemy także przyjąć wynik działania na danych średniej ruchomej, przy czym jeśli rozważamy model bez składnika sezonowego, to możemy posłużyć się średnią ruchomą dowolnego rzędu, natomiast jeśli rozważamy model ze składnikiem sezonowym, to musimy się posłużyć średnią ruchomą o rzędzie równym okresowi składnika sezonowego.

Stosując średnią ruchomą w celu estymacji trendu, nie otrzymujemy analitycznej postaci estymatora, a jedynie jego wartości. Poza tym stosując średnią ruchomą rzędu h, estymacji możemy dokonać tylko dla h < t ≤ n − h.

W kontekście estymacji trendu działanie średniej ruchomej bywa określanie jako wygładzanie (mówimy o wygładzaniu danych za pomocą średniej ruchomej, inną popularną metodą jest wygładzanie wykładnicze – tutaj nie będziemy go omawiać). Oczekujemy zatem, aby estymator trendu uzyskany w ten sposób był odpowiednio gładki (oczywiście w sensie potocznym, gdyż dla danych dyskretnych nie może być mowy o gładkości w sensie analitycznym).

T3. Jeśli rozważamy model ze składnikiem sezonowym i możemy uznać, że trend jest stały na każdym okresie, to za estymator trendu na danym okresie możemy przyjąć średnią arytmetyczną obserwacji na tym okresie.

Przy takiej estymacji trendu nie otrzymujemy jego analitycznej postaci a jedynie wartości. Dodatkowo należy zauważyć, że dość naturalnym jest założenie, iż trend ma w miarę ciągły przebieg (nie możemy mówić w sensie ścisłym o ciągłości, ponieważ obserwacji dokonujemy w czasie dyskretnym), co powoduje, że jeśli założymy, że trend jest stały na każdym okresie, to zmienność trendu jest dość mała.

Po dokonaniu estymacji trendu możemy odjąć estymator trendu od obserwacji. O tym kroku popu- larnie mówi się, że usuwamy trend. Należy pamiętać, że jeśli posłużyliśmy się metodą T2, to uzyskamy w ten sposób wartości tylko dla h < t ≤ n − h.

Jeśli w danych obecny jest trend ale nie wykazują one sezonowej struktury, zamiast dokonywania estymacji trendu możemy postąpić jeszcze inaczej: możemy usunąć go poprzez jednokrotne (trend li- niowy) lub wielokrotne różnicowanie krokiem 1. Jak to już wcześniej powiedziano, różnicowanie krokiem

1 ma bardzo klarowną interpretcję: pozwala obliczyć przyrost obserwowanego zjawiska między kolejnymi momentami pomiaru np. jeśli dane są gromadzone w odstępach jednego dnia, różnicowanie krokiem 1 pozwala obliczyć przyrost dzienny. (O usuwaniu trendu za pomocą średniej ruchomej patrz dalej: procesy ARIMA).

Estymacja sezonowości

Oprócz oceny „na oko” na podstawie wykresu danych, że w danych występuje składnik sezonowy i, jeśli tak, ile wynosi jego okres, możemy też spojrzeć na wykres próbkowej funkcji autokorelacji (autoko- wariancji) – jeśli występuje składnik sezonowy, to funkcja ta zachowuje się okresowo.

Istnieją bardziej profesjonalne narzędzia służące do wykrywania obecności składnika sezonowego i rozpoznawania długości jego okresu. Są one częścią spektralnej teorii szeregów czasowych (nieobecnej na tym wykładzie).

Dalsze postępowanie zakłada, że wiemy, iż składnik sezonowy jest obecny i ile wynosi jego okres.

Jeśli mamy do czynienia z danymi, z których usunęliśmy trend bądź brak jest trendu, ale występuje w nich sezonowość, to estymacja składnika sezonowego może przebiegać na jeden z poniższych sposobów:

S1. Pamiętając, że każdą funkcję okresową można w sposób jednostajny przybliżyć wielomianem trygonometrycznym, możemy przyjąć, że składowa sezonowa o okresie T jest wielomianem trygonome- trycznym pewnego stopnia k postaci:

st= a0+

k

X

j=1



ajcos2πj

T t + bjsin2πj T t

 .

Współczynniki a0, a1, a2, . . . , ak, b1, b2, . . . , bk estymujemy za pomocą metody najmniejszych kwadratów tzn. posługujemy się macierzą planu postaci:

















1 sin^2π_T sin^2π·2_T . . . sin^2π·k_T cos^2π_T cos^2π·2_T . . . cos^2π·k_T 1 sin 2 ·^2π_T sin 2 ·^2π·2_T . . . sin 2 ·^2π·k_T cos 2 · ^2π_T cos 2 · ^2π·2_T . . . cos 2 · ^2π·k_T . . . .

1 sin n · ^2π_T sin n ·^2π·2_T . . . sin n ·^2π·k_T cos n ·^2π_T cos n ·^2π·2_T . . . cos n ·^2π·k_T

















W ten sposób otrzymamy analityczną postać estymatora sezonowości tzn. wzór a nie tylko wartości.

Metoda wymaga jednak podania stopnia wielomianu trygonometrycznego, za pomocą którego będziemy modelowali sezonowość.

W R do skonstruowania powyższej macierzy można użyć funkcji harmonic z pakietu tsModel. Jej kolejne argumenty oznaczają: x – zakres czasu, nfreq – stopień dopasowywanego wielomianu trygono- metrycznego (czyli k), period – okres (czyli T ), intercept – wartość logiczna mówiąca, czy w macierzy ma być umieszczona kolumna jedynek.

Trzeba być ostrożnym, aby nie zbudować macierzy planu o liniowo zależnych kolumnach. Jeśli taka sytuacja będzie miała miejsce, z macierzy postaci jak wyżej należy usunąć odpowiednie kolumny, tak aby nie zmieniając przestrzeni, na którą odbędzie się rzutowanie, pozostawić układ liniowo niezależnych kolumn.

S2. Obliczamy liczby w₁, w₂, . . . , w_T w takim sposób, że w_i jest średnią arytmetyczną obserwacji w dostępnych punktach postaci i + kT , k ∈ Z. Wówczas przyjmujemy

ˆ

si= wi− 1

T(w1+ w2+ . . . + wT), i = 1, 2, . . . , T i przedłużamy okresowo {si} na t = T + 1, T + 2, . . . , n.

Postępując w sposób powyższy, nie otrzymujemy analitycznej postaci estymatora składowej sezonowej, a jedynie jego wartości.

Po dokonaniu estymacji sezonowości możemy odjąć estymator składnika sezonowego od obserwacji (wcześniej pozbawionych trendu). O tym kroku popularnie mówi się, że usuwamy sezonowość.

Należy jeszcze dodać, że dla danych z liniowym trendem (lub bez trendu) i sezonowością możemy najpierw starać się usunąć składnik sezonowy bez jego estymacji. Dokonamy tego poprzez różnicowanie

danych krokiem równym długości okresu składnika sezonowego (por. rachunek dotyczący różnicowania).

Różnicowanie danych o charakterze okresowym krokiem równym długości okresu ma, jak wcześniej wspo- mniano, klarowną interpretację, bowiem pozwala obliczyć przyrost obserwowanego zjawiska w danym momencie okresu z jednego okresu na drugi np. jeśli gromadzone są dane miesięczne i zjawisko wykazuje roczną okresowość, różnicowanie krokiem 12 pozwala obliczyć przyrost obserwowanego zjawiska w danym miesiącu rok do roku. Należy przy tym pamiętać, że otrzymamy wówczas wartości w zakresie czasu o jeden okres krótszym od danych wejściowych. Potem możemy np. przejść do estymacji trendu za pomocą jednej z metod podanych wcześniej, dopasowania modelu stacjonarnego, jeśli w danych nie ma już trendu, lub też – jeśli jest trend – do dopasowania procesu ARIMA (patrz danej).

W wyniku usuwania trendu i sezonowości chcemy otrzymać dane, które mogą stanowić realizację pro- cesu stacjonarnego (niekoniecznie takiego jak wyjściowy {et}), w szczególności nie zawierają one trendu ani sezonowości. Następnie do tych danych dopasowujemy jakiś znany nam model procesu stacjonarnego.

Funkcja decompose w R dokonuje estymacji trendu za pomocą metody T2 i następnie po usunięciu trendu dokonuje estymacji sezonowości za pomocą metody S2.

Jako argument funkcji decompose musimy podać obiekt klasy "ts" służący w R do reprezentowania szeregu czasowego. Obiekt taki tworzymy za pomocą funkcji ts. Innymi jako argument funkcji decompose musimy podać wartość funkcji ts na danych. Argument frequency tej funkcji oznacza długość okresu.

Funkcji decompose działa jedynie wtedy, gdy jako frequency podano liczbę większą od 1 i w danych mie- szą się przynajmniej dwa okresy. W szczególności za pomocą funkcji decompose nie da się przeprowadzić estymacji i eliminacji trendu dla danych pozbawionych składnika sezonowego.

Wartość random funkcji decompose to dane po usunięciu trendu i sezonowości.

Metody T1 i S1 można połączyć, od razu przybliżając dane sumą wielomianu i wielomianu trygono- metrycznego, co odpowiada połączeniu obu macierzy planu.

Bardziej zaawansowane metody dopasowania trendu i sezonowości oparte na tzw. lokalnym wygładza- niu w R można zrealizować np. za pomocą funkcji stl.

Dopasowanie modelu ARMA do danych

Do danych, które potencjalnie przedstawiają się jako realizacja procesu stacjonarnego, możemy pró- bować dopasować proces ARMA.

Podstawowym narzędziem diagnostycznym będą tu wykresy próbkowej funkcji autokorelacji i próbko- wej funkcji autokorelacji cząstkowej wraz z zaznaczonymi na nich pasami wokół osi OX, które pozwalają na odczytywanie wyników testów hipotez mówiących o zerowaniu się autokorelacji i autokorelacji prób- kowej poszczególnych rzędów.

W R do narysowania próbkowej funkcji autokorelacji wraz z odpowiednim pasem wokół osi OX służy funkcja acf. (Jeśli w funkcji acf zrezygnujemy z domyślnej wartości type="correlation" na rzecz type="covariance", to otrzymamy wykres próbkowej funkcji autokowariancji, jednak bez jakiegokolwiek pasa wokół osi OX służącego do testowania). Do narysowania próbkowej funkcji autokorelacji wraz z odpowiednim pasem wokół osi OX służy funkcja pacf. Wykres danych, próbkowej funkcji autokorelacji i próbkowej funkcji autokorelacji cząstkowej można narysować jednocześnie za pomocą funkcji tsdisplay z pakietu forecast.

Przypomnijmy, że dla procesu MA(q) funkcja autokorelacji zeruje się dla argumentów co do wartości bezwzględnej przekraczających q, a zatem jeśli dla argumentów co do wartości bezwzględnej przekracza- jących q próbkowa funkcja autokorelacji mieści się w pasie wokół zera, to możemy próbować dopasować do danych proces średniej ruchomej (rzędu co najwyżej q).

Przypomnijmy też, że dla procesu AR(p) funkcja autokorelacji cząstkowej zeruje się dla argumentów co do wartości bezwzględnej przekraczających p, wobec czego jeśli dla argumentów co do wartości bez- względnej przekraczających p próbkowa funkcja autokorelacji cząstkowej mieści się w pasie wokół zera, to możemy próbować dopasować do danych proces autoregresji (rzędu co najwyżej p).

W wielu sytuacjach proces autoregresji okazuje się bardziej intuicyjny do zrozumienia niż proces średniej ruchomej.

Opisane poniżej funkcje potrafią dopasować model ARMA plus stała tzn. proces, który po odjęciu pewnej stałej (wycentrowaniu) jest procesem ARMA (gdyż sam proces ARMA ma średnią 0).

Proces autoregresji możemy dopasować za pomocą funkcji ar. Jeśli przyjmiemy aic=TRUE, najpierw za pomocą rachunku zbliżonego do kryterium Akaike (przy założeniu, że biały szum jest gaussowski) zostanie wybrany rząd autoregesji, a następnie – stosownie do wartości argumentu method – za pomocą jednego z wybranych algorytmów (Yule’a-Walkera, Burga, uogólnionej metody najmniejszych kwadratów, metody

największej wiarogodności) zostanie przeprowadzona estymacja współczynników wielomianu autoregresji.

Jeśli przyjmiemy aic=FALSE, za pomocą wybranego algorytmu zostanie dopasowany proces autoregresji takiego rzędu jak podany w argumencie max.order.

Proces autoregresji i średniej ruchomej możemy dopasować za pomocą funkcji arima. Korzysta ona z metody największej wiarogodności (zakładając, że biały szum jest gaussowski). Rząd autoregresji i śred- niej ruchomej podajemy jako pierwszą i trzecią współrzędną wektora stanowiącego argument order. Jeśli dopasowujemy model ARMA, druga współrzędna tego wektora musi pozostać równa 0 (inaczej będziemy mieli model ARIMA). Musimy pozostawić domyślną wartość argumentu seasonal (gdyż w przeciwnym razie otrzymamy model SARIMA); argument ten nie ma nic z wspólnego ze składową sezonową.

Możemy również wybrać rząd autoregresji i średniej ruchomej za pomocą kryterium informacyj- nego Akaike, skorygowanego kryterium informacyjnego Akaike (AICc) bądź bayesowskiego kryterium informacyjnego (przy założeniu, że biały szum jest gaussowski). Służy do tego funkcja auto.arima z pakietu forecast. Należy w niej przyjąć d=0 (w przeciwnym razie dopasujemy proces ARIMA) oraz seasonal=FALSE (lub alternatywnie max.P=0 i max.Q=0, w przeciwnym razie będziemy mieli do czy- nienia z modelem SARIMA). W zależności od tego, jaką wartość logiczną przypiszemy do argumentu stepwise, zostanie wykonana procedura krokowa (zachłanna) lub brutalne przeszukiwanie wszystkich modeli (drugi wariant jest bardziej czasochłonny, ale oczywiście może dać lepszy wynik, pierwszy wariant czasem może nie wystarczyć).

Uwaga! Zamiast wstępnego usuwania trendu i sezonowości, można jednocześnie dopasować trend, se- zonowość i proces ARMA. Wystarczy że w funkcji arima bądź auto.arima w argumencie xreg umieścimy macierz złożoną z kolumn macierzy planu (poza wyrazem wolnym), która rozpina przestrzeń liniową, na którą chcemy rzutować dane.

Po dokonaniu dopasowania powinniśmy sprawdzić, czy reszty (residuals) mogą być realizacją białego szumu (w szczególności czy mogą uchodzić za nieskorelowane) oraz – jeśli posłużyliśmy się metodami wykorzystującymi założenie, że biały szum jest gaussowski – czy reszty mogą być realizacją niezależnych zmiennych losowych o tym samym rozkładzie normalnym. Jeśli tak nie jest, to znaczy, że model został źle dopasowany... i zabawa trwa nadal. Być może do dobrego dopasowania modelu potrzebna jest znajomość szerszej klasy procesów stacjonarnych niż ARMA.

Do sprawdzenia nieskorelowania reszt ponownie posłużą nam wykresy próbkowej funkcji autokowa- riancji i próbkowej funkcji autokowariancji cząstkowej. Można też użyć testu Ljunga-Boxa lub testu Boxa- Pierce’a.

W obu tych testach zakładamy, że mamy do czynienia z procesem IID(µ, σ²). Ustalamy h ≥ 1.

Testujemy

H : ρ(1) = ρ(2) = . . . = ρ(h) = 0 K : ρ(1) 6= 0 ∨ ρ(2) 6= 0 ∨ . . . ∨ ρ(h) 6= 0.

W teście Ljunga-Boxa statystyka testowa ma postać:

T = n (n + 2)

h

X

k=1

ˆ ρ²(k) n − k.

Statystyka testowa testu Boxa-Pierce’a jest nieco uproszczoną wersją statystyki testowej testu Ljunga- Boxa i ma postać:

T = n

h

X

k=1

ˆ ρ²(k).

W obu przypadkach T −−−−→^D

n→∞ χ²_h. Odrzucamy H na rzecz K na poziomie istotności α ∈ (0, 1), gdy T > (χ²_h)⁻¹(1 − α) (przyjmujemy ˆp = 1 − χ²_h(T )).

Jeśli test wykonujemy dla reszt z dopasowania procesu ARMA(p,q), bardziej adekwatny w użyciu może się okazać rozkład χ²_h−(p+q).

W R test Ljunga-Boxa i test Boxa-Pierce’a realizuje funkcja Box.test. Jeden z dwóch testów wybie- ramy za pomocą argumentu type. Wykonując test dla reszt z dopasowania procesu ARMA(p,q), warto jest przyjąć fitdf=p + q (ze wspomnianych wyżej względów).

Wykres reszt, próbkowej funkcji autokorelacji dla nich oraz p-wartości testu Ljunga-Boxa (przy fitdf=0) dla różnych h wyświetla funkcja tsdiag wywołana na obiekcie klasy "Arima" (który jest wy- nikiem działania funkcji arima bądź auto.arima).

Procesy ARIMA

Definicja 11. Proces {X_t}_t∈Znazywamy scałkowanym procesem autoregresji i średniej ruchomej (ang. autoregresive integrated moving average) rzędu p, d i q, jeśli na skutek jego d-krotnego różnicowania krokiem 1 otrzymamy proces ARMA(p,q) tzn. jeśli (1 − B)^dXtjest procesem ARMA(p,q).

Równoważnie możemy powiedzieć, że proces ARIMA(p,d,q) jest rozwiązaniem {Xt}_t∈Z równania Φ(B)(1 − B)^dXt= Θ(B)Zt, gdzie Z ∼ W N (0, σ²),

stopień wielomianu Φ wynosi p zaś stopień wielomianu Θ wynosi q, przy czym jest takim rozwiązaniem, że {(1 − B)^dXt}_t∈Zjest procesem stacjonarnym.

Na uwagę zasługuje fakt, że proces ARIMA nie musi być procesem stacjonarnym, a jedynie po jego d-krotnym różnicowaniu krokiem 1 uzyskujemy proces stacjonarny. Z tego powodu procesy ARIMA na- dają się do modelowania zjawisk, które wykazują trend lecz nie wykazują sezonowości (gdyż, jak pa- miętamy, różnicowanie dowolnym krokiem pozwala usunąć trend liniowy; trend wielomianowy wyższego rzędu można niejednokrotnie usunąć przez wielokrotne różnicowanie). Wówczas przy estymacji możemy pominąć estymację trendu i od razu starać się dopasować do danych proces ARIMA. Za pomocą procesu ARIMA możemy jednak modelować również reszty po usunięciu trendu i sezonowości.

Jako narzędzie diagnostyczne użytecznie przed próbą dopasowania procesu ARIMA do danych mo- żemy wykorzystać funkcje autokorelacji i autokorelacji cząstkowej dla danych poddanych d-krotnemu różnicowaniu krokiem 1 (po uprzednim usunięciu trendu i sezonowości). Jeśli próbkowa funkcja autokore- lacji mieści się w pasie wokół zera dla argumentów przekraczających q, to możemy próbować dopasować do danych proces ARIMA(0,d,q⁰), gdzie q⁰≤ q. Z kolei próbkowa funkcja autokorelacji cząstkowej mieści się w pasie wokół zera dla argumentów przekraczających p, to możemy próbować dopasować do danych proces ARIMA(p⁰,d,0), gdzie p⁰≤ p.

Do estymacji procesu ARIMA możemy wykorzystać funkcję arima, jak to opisano wcześniej. Argument order składa się z liczb p, d i q.

Aby dopasować proces ARIMA do danych wraz z estymacją rzędów na podstawie kryterium infor- macyjnego (zakładając, że biały szum jest gaussowski), można posłużyć się funkcją auto.arima, jak to opisano wcześniej. Wówczas argument d musimy pozostawić z domyślną wartością NA. Przypisanie argu- mentowi d ustalonej wartości spowoduje ustalenie w wyniku odpowiedniego rzędu na równy tej wartości.

Należy wystrzegać się przyjmowania zbyt dużych wartości d. Jeśli dokonamy zbyt wielu różnicować, możemy utracić możliwość dobrego dopasowania procesu stacjonarnego do danych. Takie zjawisko na- zywa się przeróżnicowaniem. Można je rozpoznać m.in. po tym, że wykres próbkowej funkcji autokorelacji dla wielokrotnie zróżnicowanych danych ma strukturę wyraźnie okresową (dopuszczalnie z malejącą am- plitudą), mimo że oryginalne dane struktury okresowej nie posiadały.

Po dokonaniu estymacji procesu ARIMA, powinniśmy sprawdzić, czy reszty mogą być realizacją białego szumu oraz – jeśli posłużyliśmy się metodami wykorzystującymi założenie, że biały szum jest gaussowski – czy reszty mogą być realizacją niezależnych zmiennych losowych o tym samym rozkładzie normalnym.