Z Skrypty dla studentów Ekonofizyki UPGOW
Analiza Szeregów Czasowych <<< Modelowanie szeregów czasowych | Matlab / GNU Octave >>>
Spis treści |
Estymacja parametrów modeli
Szacowanie parametrów modeli rządzących szeregami czasowymi to niełatwe zagadnienie. Jest to również przedostatni krok w analizie szeregów czasowych. Ostatnim krokiem jest predykcja przyszłych wartości szeregu w oparciu o dane posiadane. ,
Aby oszacować jaki to model ARMA(p,q) stoi za analizowanym szeregiem czasowym musimy wykonać kilka kroków
- jakie p i q należy wybrać
- oszacować średnią oraz współczynniki AR \( \varphi_i \ \) oraz MA \( \theta_j \ \), i=1,...,p, j=1,...,q,
- oszacować wariancję szumu \(\sigma^2 \ \) dla wybranych parametrów,
- sprawdzić poprawność wybranego modelu (najlepiej dla różnych zestawów parametrów p i q).
Ostateczna decyzja, czy dany model dobrze reprezentuje dane zależy od kilku możliwych testów.
Zakładamy obecnie, że fitować będziemy model ARMA do danych których średnia wynosi 0
- \( \langle X_t \rangle = EX_t = 0. \)
Jeżeli \( \{ Y_t \} \) oznacza oryginalne dane, to \( X_t = Y_t - EY_t \).
Będziemy dopasowywać dane do modelu
- \( X_t - \varphi_1X_{t-1} - \dots - \varphi_p X_{t-p} = Z_t + \theta_1Z_{t-1} - \dots - \varphi_q Z_{t-q}, \{ Z_t \} = BS(0,\sigma^2) \ \)
Czyli dla wybranych przez na p i q celem będzie znalezienie wektorów \( \bar{\varphi} = (\varphi_1, \dots, \varphi_p) \) oraz \( \bar{\theta} = (\theta_1, \dots, \theta_q) \).
AR
W przypadku, gdy posiadane dane mogą być przybliżone poprzez model autoregresji rzędu p (tj: q = 0), dość dobrym estymatorem wektora \( \bar{\varphi} \) okazuje się być prosty algorytm porównujący autokowariancję próby oraz teoretyczną wyliczoną z modelu AR(p). Metoda ta nosi nazwę Yule-Walkera.
Metoda Yule-Walkera
Niech \( \{ X_t \} \ \) będzie procesem AR(p) o średniej zerowej
- \( X_t - \varphi_1 X_{t-1} - \dots - \varphi_p X_{t-p} = Z_t, ~~Z_t \sim BS(0, \sigma^2). \)
Postaramy się teraz oszacować \( \bar{\varphi} \) oraz \( \sigma^2 \).
Korzystając z założenia, że \( \{ X_t \} \ \) jest losowy \( \big( X_t = \sum \Psi_j Z_{t-j} \ \big),\) otrzymujemy równania Yole - Walkera
- \( \begin{align} \Gamma_p \bar{\varphi} &= \gamma_p \\ \sigma^2 &= \gamma(0) - \bar{\varphi}' \gamma_p \end{align} \)
gdzie
- \( \Gamma_p = [\gamma(i-j)]_{i,j=1}^p \) - macierz kowariancji,
- \( \gamma_p = (\gamma(1), \gamma(2), \dots, \gamma(p) )'. \)
Z równań tych możemy łatwo oszacować \( \gamma(0), \dots, \gamma(p) \) znając \( \sigma^2, \bar{\varphi} \). Z drugiej strony, jeżeli zamienimy występujące tutaj funkcje autokowariancji na odpowiednie funkcje autokowariancji próby otrzymamy estymatory Yule-Walkera
- \( \begin{align} \hat{\Gamma}_p \hat{\varphi} &= \hat{\gamma}_p, \\ \hat{\sigma}^2 &= \hat{\gamma}(0) - \hat{\varphi}' \hat{\gamma}_p, \end{align} \)
gdzie odpowiednie
- \( \hat{\Gamma}_p = [\hat{\gamma}(i-j)]_{i,j=1}^p \) - macierz kowariancji próby,
- \( \hat{\gamma}_p = (\hat{\gamma}(1), \hat{\gamma}(2), \dots, \hat{\gamma}(p) )' \) - wektor kowariancji próby.
Jeżeli tylko \( \hat{\gamma}(0) > 0 \) to możemy podzielić obie strony przez \( \hat{\gamma}(0) \) i otrzymamy wtedy
- \( \begin{align} \hat{\bar{\varphi}} &= \hat{R}^{-1}_p \hat{\rho}_p, \\ \hat{\sigma}^2 &= \hat{\gamma}(0) [1 - \hat{\rho'}_p \hat{R}^{-1}_p \hat{\rho}_p], \end{align} \)
gdzie
- \( \hat{\rho}_p = (\hat{\rho}(1), \hat{\rho}(2), \dots, \hat{\rho}(p) )' = \hat{\gamma}_p/\hat{\gamma}(0)\) - wektor autokorelacji próby.
Z tak określonym \( \hat{\bar{\varphi}} \) można pokazać, że powyższy model jest losowy. Autokowariancje \( \gamma_F(h) \) fitowanego modelu powinny spełniać \[ \gamma_F(h) - \hat{\varphi}_1 \gamma_F(h-1) - \dots - \hat{\varphi}_p \gamma_F(h-p) = \left \{ {{\sigma^2, h = 0} \atop {0, h = 1,\dots, p}} \right. \] Każdej niesingularnej macierzy \( \Gamma_p \) odpowiada proces AR(p). Niestety nie możemy stwierdzić w ogólności tego samego dla procesu MA(q).
Generalnie (czyli teoretycznie) rząd p nie jest nam znany, gdy zaczynamy dopasowywać dane do modelu AR(p). Jeżeli założymy, że prawdziwym rzędem modelu jest rząd p a zamierzamy dopasować model rzędu m to powinniśmy oczekiwać, że oszacowany estymator \( \bar{\varphi}_m = \bar{\varphi}_{m1}, \bar{\varphi}_{m2}, \dots, \bar{\varphi}_{mm}) \) wykaże dla indeksów większych od p niskie wartości \( \bar{\varphi}_{mm}. \)
Algorytm Durbina - Levinsona
Wiemy jak wstępnie oszacować parametry metodą Yule-Walkera
- \( X_t - \varphi_{m1} X_{t-1} - \dots - \varphi_{mp} X_{t-m} = Z_t, ~~Z_t \sim BS(0, \sigma^2). \)
z odpowiednimi
- \( \hat{\varphi}_m = (\hat{\varphi}_{m1}, \dots, \hat{\varphi}_{mm}) = \hat{R}^{-1}_m \hat{\rho}_m \) - predyktory \( \varphi_i, \)
- \( \hat{v}_m = \hat{\gamma}(0) [1 - \hat{\rho'}_p \hat{R}^{-1}_p \hat{\rho}_p]\) - predyktor \( \sigma^2. \)
- Algorytm DL
- pozwala na sukcesywne obliczanie kolejnych rzędów procesu AR.
Jeżeli \( \hat{\gamma}(0) > 0 \) to dla m = 1, 2, ..., n-1 możemy rekursywnie dopasować model AR
- \( \hat{\varphi}_{11} = \hat{\rho}(1),\ \)
- \( \hat{v}_1 = \hat{\gamma}(0) [1 - \hat{\rho}^2(1)], \ \)
- \( \hat{\varphi}_{mm} = [\hat{\gamma}(m) - \sum_{j=1}^{m-1} \hat{\varphi}_{m-1,j} \hat{\gamma}(m-j)] / \hat{v}_{m-1}, \)
- \( \begin{bmatrix} \hat{\varphi}_{m1} \\ \vdots \\ \hat{\varphi}_{m,m-1}\end{bmatrix} = \hat{\bar{\varphi}}_{m-1} - \hat{\varphi}_{m,m} \begin{bmatrix} \hat{\varphi}_{m-1,m-1} \\ \vdots \\ \hat{\varphi}_{m-1,1}\end{bmatrix}, \)
- \( \hat{v}_m = \hat{v}_{m-1} (1 - \hat{\varphi}^2_{mm}). \)
| Kod MS3 (Matlab / Octave) |
|---|
function [phi, v] = asc_durbin_levinson(x, p) if ( isempty (x) || (!isvector(x)) ) error ("DL: pierwszy agrument musi byc niepustym wektorem"); return; endif %dlugosc wektora lx = length(x); if ( !isnumeric(p) || p < 1) error ("DL: drugi argument musi byc liczba naturalna 0 < |h| <= dim(wektor)\n"); return; endif %%% % algorytm Durbina - Levinsona fi = []; v = []; % krok 0 if (p == 1) p; phi = asc_kor(x,1); v = asc_kow(x,0) * (1.0 - power(phi,2.0)); else [fi, w] = asc_durbin_levinson(x, p - 1); p; suma = 0.0; for j=1:p-1 suma += fi(p-1,j) * asc_kow(x,p-j); endfor fi(p, p) = (asc_kow(x,p) - suma) / w(p-1); for m=1:p-1 fi(p,m) = fi(p-1,m) - fi(p,p)*fi(p-1,p-m); endfor w(p) = w(p-1) * (1.0 - power(fi(p,p),2.0)); phi = fi; v = w; end endfunction |
Przy okazji, korzystając z powyższego algorytmu dostajemy funkcję autokorelacji częściowej próby \( \hat{\alpha}(m) = \hat{\varphi}_{mm} \). Wiemy, że dla modelu AR(p) dla wszystkich m > p mamy \( \alpha (m) = \varphi_{mm} = 0 \). Dodatkowo estymator \( \varphi_{mm} \) dla m > p ma z dobrym przybliżeniem rozkład normalny N(0,1/n), czyli - jeżeli model AR(p) jest dobrze dobrany to dla dostatecznie dużych m wartości \( \varphi_{mm} \) powinny się zgadzać z rozkładem normalnym N(0,1/n). W szczególności, jeżeli rząd AR jest równy p to \( \varphi_{mm} \) dla m > p powinno zawierać się w granicy 95% zgodności, czyli leżeć pomiędzy granicami \( \pm 1.96/\sqrt{n} \) z prawdopodobieństwem bliskim 0.95 (tylko 5% wartości może się nie zgadzać, czyli wypaść poza podane granice). Jeżeli tak jest w rzeczywistości to obliczony model AR(p) będzie dobrym oszacowaniem procesu losowego stojącego za danym szeregiem czasowym.
- Przykład
- Poniżej prezentujemy przykład użycia powyższego algorytmu do dopasowania sztucznych danych. Generujemy proces AR(2), z \( \varphi_1 = 0.9, \varphi = -0.3 \) i wariancją \( \sigma^2 = 1 \). Następnie proces ten poddajemy działaniu algorytmu DL aż do p=5. Widzimy, że oszacowanie parametrów wielomianu \( \varphi \) dla pierwszych dwóch wartości różni się o kilka procent dla wszystkich p > 1, natomiast oszacowane wielkości \( \varphi_3, ..., \varphi_5\) są znikomo małe (maksymalnie 4% \( \varphi_1 \) lub 12% \( \varphi_2 \)) i dodanie ich do modelu nie zmieni w zasadzie niczego, dlatego można je odrzucić jako mało istotne. Wartość wariancji obliczona została na poziomie 99% wartości prawdziwej, \( \varphi_1 \) na poziomie 98% a \( \varphi_2 \) na poziomie mniejszym niż 1%, dając. Reasumując
- oryginalny model
- \( X_t - 0.9 X_{t-1} - 0.3 X_{t-2} = Z_t \ \)
- dopasowany model [dla asc_durbin_levinson(x, 2)]
- \( \hat{X}_t - 0.91710 \hat{X}_{t-1} - 0.30064 \hat{X}_{t-2} = \hat{Z}_t \ \)
octave:69> x = arma_rnd([0.9,-0.3], [0], 1, 10000, 10000); octave:70> [phi, v] = asc_durbin_levinson(x, 5) phi = 0.70512 0.00000 0.00000 0.00000 0.00000 0.91710 -0.30064 0.00000 0.00000 0.00000 0.92346 -0.32003 0.02114 0.00000 0.00000 0.92403 -0.32871 0.04618 -0.02712 0.00000 0.92404 -0.32872 0.04625 -0.02732 0.00022 v = 1.09071 0.99213 0.99169 0.99096 0.99096
Algorytm Burga
TBW
MA
Dla przypadków, gdy q > 0 metoda zaprezentowana wcześniej nie do końca zdaje egzamin. W tym przypadku posługujemy się algorytmem innowacyjnym.
Tak samo, jak możemy dopasować modele AR do posiadanych szeregów stosując algorytm Durbina-Levinsona do sutokorelacji częściowych, w podobny sposób, możemy rekursywnie dopasowywać dane do procesów średniej kroczącej
- \( X_t = Z_t + \hat{\theta_{m1}} Z_{t-1} + \dots + \hat{\theta_{mm}} Z_{t-m}, \{ Z_t \} = BS(0,\sigma^2), \,\)
dla m = 1,2,... używając algorytmu innowacyjnego.
Algorytm innowacyjny
Będziemy szacować wektor MA, czyli wektor \( \hat{\mathbf{\theta}}_m = (\hat{\theta}_{m1}, \dots, \hat{\theta}_{mm})', \,\) oraz wariancję procesu BS \( \hat{v}_m. \,\)
- Algorytm
Jeżeli \( \{ X_t \} \,\) ma średnią zerową oraz wartość oczekiwana korelatora \( E(X_iX_j) = \kappa(i,j) \) jest zdefiniowana przez niesingularną macierz kowariancji \([\kappa(i,j)]_{i,j=1}^n \,\) dla każdego n=1,2,... wtedy wartości \( \hat{X}_{n+1}, n \ge 0 \,\) oraz odpowiednie średnio-kwadratowe błędy \(v_n, n \ge 1 \,\) dane są poprzez równanie
- \( \hat{X}_{n+1} = \left \{ {{0, n = 0} \atop {\sum_{j=0}^n \theta_{nj} (X_{n+1-j} - \hat{X}_{n+1-j}), n \ge 1}} \right. \)
a estymatory modelu liczymy rekurencyjnie z równań
- \( \begin{align} v_0 &= \kappa(1,1), \\ \theta_{n,n-k} &= v_k^{-1} \Big( \kappa(n+1,k+1) - \sum_{j=0}^{k-1} \theta_{k,k-j} \theta_{n,n-j} v_j, ~~k=0,\dots,n-1\\ v_n &= \kappa(n+1,k+1) - \sum_{j=0}^{k-1} \theta^2_{n,n-j} v_j \end{align} \)
Poniżej prezentujemy skrypt Matlab / GNU Octave realizujący algorytm innowacyjny.
| Matlab / Octave |
|---|
function [theta,v] = asc_innovation(x, qIn) if ( isempty (x) || (!isvector(x)) ) error ("I: pierwszy agrument musi byc niepustym wektorem"); return; endif %dlugosc wektora lx = length(x); if ( !isnumeric(qIn) || qIn < 1) error ("I: drugi argument musi byc liczba naturalna q > 0\n"); return; endif %%% % algorytm innowacyjny %%% % pracujemy na wektorach w = []; teta = []; %%% % wartości startowe % dla v(1) i theta(1,1) v_zero = asc_kow(x,0); w = v_zero; teta = asc_kor(x,1); %%% % właściwy algorytm innowacyjny for q = 1 : qIn % dla rzedu = MA(1)...MA(q) %%% % kolejne wartości \theta_{q,q-k} for k=0:q-1 suma = 0; for j=0:k-1 if k > j suma += teta(q,q-j) * teta(k,k-j) * w(j+1); endif endfor teta(q,q-k) = (asc_kow(x,q-k) - suma) / w(k+1); %%% % kolejne wartości v_q suma = 0; for j=0:q-1 if j == 0 fau = v_zero; else fau = w(j); endif if q > j suma = power(teta(q,q-j),2) * fau; endif endfor w(q) = asc_kow(x,0) - suma; endfor endfor theta = teta; v = w; endfunction |
- Przykład
Poniżej prezentujemy przykład użycia powyższego algorytmu do dopasowania sztucznych danych. Generujemy proces MA(2), z \( \theta_1 = 0.5, \theta_2 = -0.3 \) i wariancją \( \sigma^2 = 0.3 \). Następnie proces ten poddajemy działaniu algorytmu innowacyjnego aż do q=5. Widzimy, że oszacowanie parametrów wielomianu \( \theta \) dla pierwszych dwóch wartości różni się o kilka procent dla wszystkich q > 1, natomiast oszacowane wielkości \( \theta_3, ..., \theta_5\) są znikomo małe (maksymalnie 3% \( \theta_1 \) lub 4% \( \theta_2 \)) i dodanie ich do modelu nie zmieni w zasadzie niczego, dlatego można je odrzucić jako mało istotne. Wartość wariancji obliczona została na poziomie 81% wartości prawdziwej, \( \theta_1 \) na poziomie 62% a \( \theta_2 \) na poziomie 81%, dając
- oryginalny model
- \( X_t = Z_t + 0.5 Z_{t-1} - 0.3 Z_{t-2}, \sigma^2 = 0.3 \ \)
- dopasowany model [dla asc_innovation(x, 2)]
- \( \hat{X}_t = \hat{Z}_t + 0.31128 \hat{Z}_{t-1} - 0.24571 \hat{Z}_{t-2}, \hat{\sigma}^2 = 0.36761 \ \)
octave:90> x = arma_rnd([0], [0.5,-0.3], 0.3, 10000, 10000); octave:91> [theta, v] = asc_innovation(x,5) theta = 0.25308 0.00000 0.00000 0.00000 0.00000 0.31128 -0.24571 0.00000 0.00000 0.00000 0.32776 -0.25557 0.01019 0.00000 0.00000 0.33206 -0.25853 0.01077 -0.00110 0.00000 0.33087 -0.25614 0.00534 0.00319 -0.01660 v = 0.37838 0.36761 0.36478 0.36405 0.36442
ARMA
Tak więc wiemy już jak szacować parametry procesów AR(p) i MA(q) w przypadkach gdy szeregi czasowe mogą być nimi opisane. Najbardziej jednak ogólną metodą szacowania parametrów modeli będzie znalezienie modelu ARMA(p,q) jako dopasowania do analizowanych danych. W przypadku gdy będziemy analizować np: dane AR(3) metodami dopasowywania do ARMA(p,q), powinniśmy z dobrym przybliżeniem dostać model ARMA(3,0).
Rekurencyjny algorytm dopasowania ARMA
W tym przypadku zakładamy, że mamy do czynienia z przyczynowym procesem ARMA(p,q)
- \( X_t - \varphi_{1} X_{t-1} - \dots - \varphi_{p} X_{t-p} = Z_t + \hat{\theta_{1}} Z_{t-1} + \dots + \hat{\theta_{q}} Z_{t-q}, ~~\{ Z_t \} = BS(0,\sigma^2). \,\)
Przyczynowość zakłada, że
- \( X_t = \sum_{j=0}^\infty \Psi_j Z_{t-j}, \,\)
czyli, że
- \( \begin{align} \Psi_0 &= 1,\\ \Psi_j &= \theta_j + \sum_{i=1}^{\min(j,p)} \varphi_i \Psi_j-i, ~~~j \ne 0. \end{align} \)
Oczywiście dla \(j > q, \theta_j = 0\) oraz dla \( j>p, \varphi_j = 0 \). Aby obliczyć \(\Psi_j \) możemy użyć np: algorytmu innowacyjnego a następnie podstawić \(\Psi_i = \hat{\theta}_{mi}. , i=1,\dots,p+q\). Dostaniemy dzięki temu równanie
- \( \hat{\theta}_{mj} = \theta_j + \sum_{i=1}^{\min(j,p)} \varphi_i \hat{\theta}_{mj-i}, ~~~j=1,\dots,p+q. \)
Dzięki tej transformacji znaleźliśmy wstępnie oszacowane parametry \( \hat{\bar{\varphi}} \) i \( \hat{\bar{\theta}} \) (wektory parametrów modeli AR i MA). Jeżeli z powyższego równania odrzucimy część związaną z losową częścią średniej kroczącej i zostawimy parametry dla j = q + 1, ..., q + p, to możemy łatwo zobaczyć, że spełnione będzie poniższe równanie na parametry wektora \( \hat{\bar{\varphi}} \)
- \( \begin{bmatrix} \hat{\theta}_{m,q} & \hat{\theta}_{m,q-1} & \hat{\theta}_{m,q-2} & \cdots & \hat{\theta}_{m,q+1-p} \\ \hat{\theta}_{m,q+1} & \hat{\theta}_{m,q} & \hat{\theta}_{m,q-1} & \cdots & \hat{\theta}_{m,q+2-p} \\ \hat{\theta}_{m,q+2} & \hat{\theta}_{m,q+1} & \hat{\theta}_{m,q} & \cdots & \hat{\theta}_{m,q+3-p} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \hat{\theta}_{m,q+p-1} & \hat{\theta}_{m,q+p-2} & \hat{\theta}_{m,q+p-3} & \cdots & \hat{\theta}_{m,q} \end{bmatrix} \begin{bmatrix} \hat{\varphi}_{1} \\ \hat{\varphi}_{2} \\ \hat{\varphi}_{3} \\ \vdots \\ \hat{\varphi}_{p}\end{bmatrix} = \begin{bmatrix} \hat{\theta}_{m,q+1} \\ \hat{\theta}_{m,q+2} \\ \hat{\theta}_{m,q+3} \\ \vdots \\ \hat{\theta}_{m,q+p}\end{bmatrix}. \)
Rozwiązując powyższe równanie macierzowe dostaniemy cały wektor \( \hat{\bar{\varphi}} \), który następnie użyjemy do obliczenia
- \( \theta_j = \hat{\theta}_{mj} - \sum_{i=1}^{\min(j,p)} \varphi_i \hat{\theta}_{mj-i}, ~~~\mathbf{j=1,\dots,q}. \)
Na samym końcu znajdujemy oszacowanie wariancji procesu
- \( \hat{\sigma}^2 = \hat{v}_m. \,\)
Jak zwykle prezentujemy program w języku MATLAB / GNU Octave implementujący powyższy algorytm. Poniższy algorytm odwołuje się do powyższej zaprezentowanego algorytmu innowacyjnego.
| Kod MS3 (Matlab / Octave) |
|---|
function [phi, theta, v] = asc_armafit(x, p, q) if ( isempty (x) || (!isvector(x)) ) error ("armafit: pierwszy agrument musi byc niepustym wektorem"); return; endif lx = length(x); if ( !isnumeric(p) || p < 0 || !isnumeric(p) || q < 0) error ("armafit: drugi i trzeci argument musi byc liczba naturalna p,q >= 0\n"); return; endif disp("armafit(x,p,q)"); p = p q = q disp("--------------"); v = []; phi = []; theta = []; [preTheta, v] = asc_innovation(x,p+q); preTheta = preTheta(p+q,:); MpreTheta = []; for i=1:p for j=1:p index = q + j - i; if index < 1 localPreTheta = 0; else localPreTheta = preTheta(index); endif MpreTheta(i,j) = localPreTheta; end end phi = inv(MpreTheta) * preTheta(q+1:q+p)'; suma = []; for j = 1:q suma(j) = 0; for i = 1:min(j,p) suma(j) += phi(i) * preTheta(j-i+1); end end theta = preTheta(1:q) - suma; % oszacowanie phi = phi; theta = theta'; v = v(p+q); endfunction |