Analiza Szeregów Czasowych/Techniki analizy szeregów czasowych

Z Skrypty dla studentów Ekonofizyki UPGOW

(Różnice między wersjami)
(→Metoda Yule-Walkera)
m (→Algorytm innowacyjny)
 
(Nie pokazano 33 wersji pomiędzy niniejszymi.)
Linia 55: Linia 55:
Jeżeli tylko <math> \hat{\gamma}(0) > 0 </math> to możemy podzielić obie strony przez <math> \hat{\gamma}(0) </math> i otrzymamy wtedy
Jeżeli tylko <math> \hat{\gamma}(0) > 0 </math> to możemy podzielić obie strony przez <math> \hat{\gamma}(0) </math> i otrzymamy wtedy
 +
: <math>
 +
\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}
 +
</math>
 +
gdzie
 +
: <math> \hat{\rho}_p = (\hat{\rho}(1), \hat{\rho}(2), \dots, \hat{\rho}(p) )' = \hat{\gamma}_p/\hat{\gamma}(0)</math> - wektor autokorelacji próby.
 +
Z tak określonym <math> \hat{\bar{\varphi}} </math> można pokazać, że powyższy model jest losowy. Autokowariancje <math> \gamma_F(h) </math>  fitowanego modelu powinny spełniać
 +
:<math> \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. </math>
 +
Każdej niesingularnej macierzy <math> \Gamma_p </math> 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 <math> \bar{\varphi}_m = \bar{\varphi}_{m1}, \bar{\varphi}_{m2}, \dots, \bar{\varphi}_{mm}) </math> wykaże dla indeksów większych od ''p'' niskie wartości <math> \bar{\varphi}_{mm}. </math>
====Algorytm Durbina - Levinsona====
====Algorytm Durbina - Levinsona====
-
====Algorytm Burga====
+
Wiemy jak wstępnie oszacować parametry metodą Yule-Walkera
 +
: <math> X_t - \varphi_{m1} X_{t-1} - \dots - \varphi_{mp} X_{t-m} = Z_t, ~~Z_t \sim BS(0, \sigma^2). </math>
 +
z odpowiednimi
 +
: <math> \hat{\varphi}_m = (\hat{\varphi}_{m1}, \dots, \hat{\varphi}_{mm}) = \hat{R}^{-1}_m \hat{\rho}_m </math> - predyktory <math> \varphi_i, </math>
 +
: <math> \hat{v}_m = \hat{\gamma}(0) [1 - \hat{\rho'}_p \hat{R}^{-1}_p \hat{\rho}_p]</math> - predyktor <math> \sigma^2. </math>
 +
 
 +
; Algorytm DL : pozwala na sukcesywne obliczanie kolejnych rzędów procesu AR.
 +
Jeżeli <math> \hat{\gamma}(0) > 0 </math> to dla ''m = 1, 2, ..., n-1'' możemy rekursywnie dopasować model AR
 +
: <math> \hat{\varphi}_{11} = \hat{\rho}(1),\ </math>
 +
: <math> \hat{v}_1 = \hat{\gamma}(0) [1 - \hat{\rho}^2(1)], \ </math>
 +
: <math> \hat{\varphi}_{mm} = [\hat{\gamma}(m) - \sum_{j=1}^{m-1} \hat{\varphi}_{m-1,j} \hat{\gamma}(m-j)] / \hat{v}_{m-1}, </math>
 +
: <math> \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}, </math>
 +
: <math> \hat{v}_m = \hat{v}_{m-1} (1 - \hat{\varphi}^2_{mm}). </math>
 +
 
 +
<source lang="matlab">
 +
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
 +
</source>
 +
 
 +
Przy okazji, korzystając z powyższego algorytmu dostajemy funkcję autokorelacji częściowej próby <math> \hat{\alpha}(m) = \hat{\varphi}_{mm} </math>. Wiemy, że dla modelu AR(p) dla wszystkich ''m > p'' mamy <math> \alpha (m) = \varphi_{mm} = 0 </math>. Dodatkowo estymator <math> \varphi_{mm} </math> 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 <math> \varphi_{mm} </math> 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 <math> \varphi_{mm} </math> dla ''m > p'' powinno zawierać się w granicy 95% zgodności, czyli leżeć pomiędzy granicami <math> \pm 1.96/\sqrt{n} </math> 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 <math> \varphi_1 = 0.9, \varphi = -0.3 </math> i wariancją <math> \sigma^2 = 1 </math>. Następnie proces ten poddajemy działaniu algorytmu DL aż do ''p=5''. Widzimy, że oszacowanie parametrów wielomianu <math> \varphi </math> dla pierwszych dwóch wartości różni się o kilka procent dla wszystkich ''p > 1'', natomiast oszacowane wielkości <math> \varphi_3, ...,  \varphi_5</math> są znikomo małe (maksymalnie 4% <math> \varphi_1 </math> lub 12% <math> \varphi_2 </math>) 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, <math> \varphi_1 </math> na poziomie 98% a <math> \varphi_2 </math> na poziomie mniejszym niż 1%, dając. Reasumując
 +
* oryginalny model
 +
: <math> X_t - 0.9 X_{t-1} - 0.3 X_{t-2} = Z_t \ </math>
 +
* dopasowany model [dla asc_durbin_levinson(x, 2)]
 +
: <math> \hat{X}_t - 0.91710 \hat{X}_{t-1} - 0.30064 \hat{X}_{t-2} = \hat{Z}_t \ </math>
 +
 
 +
<source lang="Matlab">
 +
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
 +
</source>
 +
 
 +
<!-- ====Algorytm Burga==== -->
===MA===
===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'''.
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
 +
: <math> X_t = Z_t + \hat{\theta_{m1}} Z_{t-1} + \dots + \hat{\theta_{mm}} Z_{t-m}, \{ Z_t \} = BS(0,\sigma^2), \,</math>
 +
dla ''m = 1,2,...'' używając '''algorytmu innowacyjnego'''.
====Algorytm innowacyjny====
====Algorytm innowacyjny====
 +
Będziemy szacować wektor MA, czyli wektor <math> \hat{\mathbf{\theta}}_m = (\hat{\theta}_{m1}, \dots, \hat{\theta}_{mm})', \,</math> oraz wariancję procesu ''BS'' <math> \hat{v}_m. \,</math>
 +
 +
; Algorytm :
 +
Jeżeli <math> \{ X_t \} \,</math> ma średnią zerową oraz wartość oczekiwana korelatora <math> E(X_iX_j) = \kappa(i,j) </math> jest zdefiniowana przez niesingularną macierz kowariancji <math>[\kappa(i,j)]_{i,j=1}^n \,</math> dla każdego ''n=1,2,...'' wtedy wartości <math> \hat{X}_{n+1}, n \ge 0 \,</math> oraz odpowiednie średnio-kwadratowe błędy <math>v_n, n \ge 1 \,</math> dane są poprzez równanie
 +
: <math> \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. </math>
 +
a estymatory modelu liczymy rekurencyjnie z równań
 +
: <math>
 +
\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,n+1) - \sum_{j=0}^{n-1} \theta^2_{n,n-j} v_j
 +
\end{align}
 +
</math>
 +
 +
Poniżej prezentujemy skrypt Matlab / GNU Octave realizujący algorytm innowacyjny.
 +
 +
<source lang="matlab">
 +
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
 +
</source>
 +
 +
; Przykład :
 +
Poniżej prezentujemy przykład użycia powyższego algorytmu do dopasowania sztucznych danych. Generujemy proces MA(2), z <math> \theta_1 = 0.5, \theta_2 = -0.3 </math> i wariancją <math> \sigma^2 = 0.3 </math>. Następnie proces ten poddajemy działaniu algorytmu innowacyjnego aż do ''q=5''. Widzimy, że oszacowanie parametrów wielomianu <math> \theta </math> dla pierwszych dwóch wartości różni się o kilka procent dla wszystkich ''q > 1'', natomiast oszacowane wielkości <math> \theta_3, ...,  \theta_5</math> są znikomo małe (maksymalnie 3% <math> \theta_1 </math> lub 4% <math> \theta_2 </math>) 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, <math> \theta_1 </math> na poziomie 62% a <math> \theta_2 </math> na poziomie 81%, dając
 +
* oryginalny model
 +
: <math> X_t = Z_t + 0.5 Z_{t-1} - 0.3 Z_{t-2}, \sigma^2 = 0.3 \ </math>
 +
* dopasowany model [dla asc_innovation(x, 2)]
 +
: <math> \hat{X}_t = \hat{Z}_t + 0.31128 \hat{Z}_{t-1} - 0.24571 \hat{Z}_{t-2}, \hat{\sigma}^2 = 0.36761 \ </math>
 +
 +
<source lang="Matlab">
 +
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
 +
</source>
===ARMA===
===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====
====Rekurencyjny algorytm dopasowania ARMA====
 +
W tym przypadku zakładamy, że mamy do czynienia z przyczynowym procesem ARMA(p,q)
 +
: <math> 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). \,</math>
 +
Przyczynowość zakłada, że
 +
: <math> X_t = \sum_{j=0}^\infty \Psi_j Z_{t-j}, \,</math>
 +
czyli, że
 +
: <math>
 +
\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}
 +
</math>
 +
Oczywiście dla <math>j > q, \theta_j = 0</math> oraz dla <math> j>p, \varphi_j = 0 </math>. Aby obliczyć <math>\Psi_j </math> możemy użyć np: [[Analiza_Szeregów_Czasowych/Techniki_analizy_szeregów_czasowych#Algorytm_innowacyjny|algorytmu innowacyjnego]] a następnie podstawić <math>\Psi_i = \hat{\theta}_{mi}. , i=1,\dots,p+q</math>. Dostaniemy dzięki temu równanie
 +
: <math> \hat{\theta}_{mj} = \theta_j + \sum_{i=1}^{\min(j,p)} \varphi_i \hat{\theta}_{mj-i}, ~~~j=1,\dots,p+q. </math>
 +
Dzięki tej transformacji znaleźliśmy wstępnie oszacowane parametry <math> \hat{\bar{\varphi}} </math> i <math> \hat{\bar{\theta}} </math> (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 <math> \hat{\bar{\varphi}} </math>
 +
: <math>
 +
\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}.
 +
</math>
-
====Estymacja modelem ''maximum Likelihood''====
+
Rozwiązując powyższe równanie macierzowe dostaniemy cały wektor <math> \hat{\bar{\varphi}} </math>, który następnie użyjemy do obliczenia
 +
: <math> \theta_j = \hat{\theta}_{mj} - \sum_{i=1}^{\min(j,p)} \varphi_i \hat{\theta}_{mj-i}, ~~~\mathbf{j=1,\dots,q}. </math>
 +
Na samym końcu znajdujemy oszacowanie wariancji procesu
 +
: <math> \hat{\sigma}^2 = \hat{v}_m. \,</math>
 +
 +
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 [[Analiza_Szeregów_Czasowych/Techniki_analizy_szeregów_czasowych#Algorytm_innowacyjny|algorytmu innowacyjnego]].
 +
 +
<source lang="matlab">
 +
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
 +
</source>
 +
 +
; Przykład :
 +
Poniżej prezentujemy przykład użycia powyższego algorytmu do dopasowania sztucznych danych. Generujemy proces ARMA(2,1):
 +
: <math> X_t - 0.2 X_{t-1} - 0.1 X_{t-2} = Z_t + 0.3 Z_{t-1}, ~~~\sigma^2 = 0.1. \ </math>
 +
Następnie proces ten poddajemy działaniu powyższego algorytmu ''armafit''.
 +
Wartość wariancji obliczona została na poziomie 99% wartości prawdziwej, <math> \theta_1 </math> na poziomie 96%, <math> \varphi_1 </math> na poziomie 32% a <math> \varphi_2 </math>na poziomie 19% dając model
 +
* dopasowany model [dla asc_armafit(x, 2, 1)]
 +
: <math> X_t - 0.33490 X_{t-1} - 0.28179 X_{t-2} = Z_t + 0.26164 Z_{t-1}, ~~~\sigma^2 = 0.10983. \ </math>
 +
 +
<source lang="Matlab">
 +
octave:20> x = arma_rnd([0.2,0.1], [0.3], 0.1, 10000, 10000);
 +
octave:21> [phi, theta, v] = asc_armafit(x, 2, 1)
 +
armafit(x,p,q)
 +
p =  2
 +
q =  1
 +
--------------
 +
phi =
 +
 +
  0.33490
 +
  0.28179
 +
 +
theta =  0.26164
 +
v =  0.10983
 +
</source>
 +
 +
====Estymacja modelem maksimum wiarygodności====
 +
Zakładamy, że mamy do czynienia z przyczynowym procesem ARMA(p,q)
 +
: <math> 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), \theta_0 = 1. \,</math>
 +
Przyczynowość zakłada, że dla <math> |u| \le 1 </math> część równania związana z procesem AR będzie niezerowa (<math> 1 + \varphi_{1} u - \dots - \varphi_{p} u^p \ne 0 </math>). Podobnie zakładamy z częścią odpowiedzialną za losowość w układzie dla <math> |u| < 1 </math>, <math> 1 + \theta_{1} u - \dots - \theta_{q} u^q \ne 0 </math> (czyli zakładamy zaszumienie układu).
 +
 +
Zadaniem naszym jest oszacować wektory <math> \hat{\bar{\varphi}} </math> i <math> \hat{\bar{\theta}} </math> oraz wariancję szumu metodą największej wiarygodności.
 +
Musimy zatem zmaksymalizować funkcję wiarygodności Gaussa daną poprzez
 +
 +
: <math> L(\hat{\bar{\varphi}}, \hat{\bar{\theta}}, \hat{\sigma}) = \frac{1}{\sqrt{(2 \pi \hat{\sigma}^2)^n r_0 \cdots r_{n-1}}} \exp \Big[ - \frac{\hat{\sigma}^2}{2} \sum_{j=1}^n  \frac{(X_j - \hat{X}_j)^2}{r_{j-1}}\Big].</math>
 +
 +
W powyższym wzorze
 +
* <math> \hat{X}_i </math> to jednokrokowe predyktory <math> X_i </math> obliczane po wstępnym szacowaniu parametrów modelu,
 +
* <math> r_i </math> - to ich odpowiednie średniokwadratowe błędy (szacowania).
 +
Predyktory te obliczamy z wzorów rekurencyjnych (patrz program ''asc_one_step_predictors'')
 +
: <math> \hat{X}_1 = 0, </math>
 +
: <math> \hat{X}_i = \sum_{j=1}^i \theta_{ij} (X_{i-j} - \hat{X}_{i-j}) , ~~~2 \le i < \max(p,q)</math>
 +
: <math> \hat{X}_i = \sum_{j=1}^p\varphi_{j} X_{i-j} + \sum_{j=1}^q \theta_{ij} (X_{i-j} - \hat{X}_{i-j}) , ~~~i \ge \max(p,q)</math>
 +
 +
Błędy ''r'' szacowane są z kolei za pomocą wzorów rekurencyjnych (pamiętając, że <math>\theta_{kj} = 0</math> dla ''j > q'' i danego ''k'')
 +
: <math> r_1 = 1, </math>
 +
: <math> r_k = 1 - \sum_{j=0}^{k-1} \theta^2_{k,k-j} r_j, k = 2, \dots, n.</math>
 +
 +
; Algorytm :
 +
# Wstępnie szacujemy wektory <math> \hat{\bar{\varphi}} </math> i <math> \hat{\bar{\theta}} </math> oraz wariancję szumu  dowolnym algorytmem (rekurencyjnym dla ARMA, czy innowacyjnym dla modelu losowego).
 +
# Obliczamy predyktory jednokrokowe i ich błędy dla konkretnych wartości <math> \hat{\bar{\varphi}}, \hat{\bar{\theta}}.</math>
 +
# Liczymy funkcję <math> S = \sum_{j=1}^n  \frac{(X_j - \hat{X}_j)^2}{r_{j-1}}. </math> Jest ona podstawą szacowania maksimum wiarygodności. Możemy teraz postąpić dwojako
 +
## Znajdujemy minimum funkcji kwadratów błędów S.
 +
## Znajdujemy minimum zredukowanej funkcji wiarygodności <math> l = \ln(S/n) + \sum_{j=1}^n \ln r_{j-1} / n. </math>
 +
# Mając parametry minimalizujące jedną z powyższych funkcji kończymy obliczając wariancję szumu z wzoru <math> \hat{\sigma}^2 = S / n </math> (dla minimalnej wartości danej statystyki.
 +
 +
Poniżej znajdziecie program obliczający predyktory jednokrokowe oraz realizujący powyższy algorytm znajdywania parametrów modelu w oparciu o minimalizację wybranej statystyki (''l'' lub ''S'').
 +
 +
<source lang="matlab">
 +
function [phi, theta, v, stat] = asc_armafit_likelihood(x, P, Q, method='likeli')
 +
 +
    if ( isempty (x) || (!isvector(x)) )
 +
    error ("armafit: pierwszy agrument musi byc niepustym wektorem");
 +
    return;
 +
  endif
 +
 +
  %dlugosc wektora
 +
  lx = length(x);
 +
 +
  if ( !isnumeric(P) || P < 0 || !isnumeric(Q) || Q < 0)
 +
    error ("armafit: drugi i trzeci argument musi byc liczba naturalna p,q >= 0\n");
 +
    return;
 +
  endif
 +
 +
  v    = [];
 +
  phi  = [];
 +
  theta = [];
 +
 
 +
  S_params = zeros(P * Q, P + Q);
 +
  S        = zeros(P*Q,1);
 +
  indexes  = zeros(P*Q,2);
 +
  for p = 1 : P
 +
    for q = 1 : Q
 +
      i = (p - 1) * Q + q;
 +
     
 +
      indexes(i,1:2) = [p,q];
 +
     
 +
      [phi_tmp, theta_tmp, v_tmp] = asc_armafit(x, p, q);
 +
      S_params(i, 1 : p + q) = [phi_tmp', theta_tmp'];
 +
     
 +
      [xbar, r] = asc_one_step_predictors(x, phi_tmp, theta_tmp);
 +
      S(i) = sum(power(x(2:end)' - xbar(2:end),2)./r(1:end-1));     
 +
    endfor
 +
  endfor
 +
 
 +
  if method == "square" 
 +
    stat = S;
 +
  else
 +
    stat = log(S / lx) + sum(log(r)) / lx;   
 +
  endif
 +
 +
  [S_value, S_index] = min(stat); 
 +
  v    = min(S) / lx;
 +
  phi  = S_params(S_index, 1 : indexes(S_index, 1));
 +
  theta = S_params(S_index, indexes(S_index,1) + 1 : indexes(S_index, 1) + indexes(S_index, 2));
 +
 +
endfunction
 +
</source>
 +
 +
<source lang="matlab">
 +
function [xbar, r] = asc_one_step_predictors(x, phi, theta)
 +
 +
    if ( isempty (x) || (!isvector(x)) )
 +
    error ("OSP: pierwszy agrument musi byc niepustym wektorem");
 +
    return;
 +
  endif
 +
 +
  %dlugosc wektora
 +
  lx = length(x);
 +
  q  = length(theta);
 +
  p  = length(phi);
 +
 
 +
  if ( !isvector(phi) || isempty(phi) || !isvector(theta) || isempty(theta))
 +
    error ("OSP: drugi i trzeci argument musi byc niepustym wektorem\n");
 +
    return;
 +
  endif
 +
 
 +
  xbar = zeros(1, lx);
 +
  r    = zeros(1, lx);
 +
  %m    = q;
 +
  m    = max(q,p); %TO DO & RETHINK!
 +
 
 +
  if q == m
 +
    vec = theta;
 +
  else
 +
    vec = phi;
 +
  endif
 +
 
 +
  % m pierwszych predyktorów
 +
  for i = 2 : m
 +
    suma = 0;
 +
    for j = 1 : i
 +
      suma += vec(j) * (x(i+1-j) - xbar(i+1-j));
 +
    endfor
 +
    xbar(i) = suma;
 +
  endfor
 +
 
 +
  % lx - m - 1 kolejnych predyktorów
 +
  for i = m + 1: lx
 +
    suma_phi  = 0;
 +
    for k = 1 : p
 +
      suma_phi += phi(k) * x(i + 1 - k);
 +
    endfor
 +
    suma_theta = 0;
 +
    for j = 1 : q
 +
      suma_theta += theta(j) * (x(i+1-j) - xbar(i+1-j));
 +
    endfor
 +
    xbar(i) = suma_phi + suma_theta;
 +
  endfor
 +
 
 +
  % wektor błędów predyktorów
 +
  r            = [];
 +
  theta_r      = zeros(1, lx);
 +
  theta_r(1:q) = theta;
 +
 
 +
  r(1) = 1;
 +
  for n = 2 : lx
 +
    suma_r = 0;
 +
    for j = 1 : n - 1
 +
      suma_r += power(theta_r(n-j),2) * r(j);
 +
    endfor
 +
    r(n) = 1 - suma_r;
 +
  endfor
 +
 
 +
endfunction
 +
</source>
 +
 +
<!--
==Testy==
==Testy==
===AICC===
===AICC===
-
 
==Prognoza==
==Prognoza==
 +
-->

Aktualna wersja na dzień 18:19, 27 maj 2013

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}). \)
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


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,n+1) - \sum_{j=0}^{n-1} \theta^2_{n,n-j} v_j \end{align} \)

Poniżej prezentujemy skrypt Matlab / GNU Octave realizujący algorytm innowacyjny.

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.

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
Przykład 

Poniżej prezentujemy przykład użycia powyższego algorytmu do dopasowania sztucznych danych. Generujemy proces ARMA(2,1):

\( X_t - 0.2 X_{t-1} - 0.1 X_{t-2} = Z_t + 0.3 Z_{t-1}, ~~~\sigma^2 = 0.1. \ \)

Następnie proces ten poddajemy działaniu powyższego algorytmu armafit. Wartość wariancji obliczona została na poziomie 99% wartości prawdziwej, \( \theta_1 \) na poziomie 96%, \( \varphi_1 \) na poziomie 32% a \( \varphi_2 \)na poziomie 19% dając model

  • dopasowany model [dla asc_armafit(x, 2, 1)]
\( X_t - 0.33490 X_{t-1} - 0.28179 X_{t-2} = Z_t + 0.26164 Z_{t-1}, ~~~\sigma^2 = 0.10983. \ \)
octave:20> x = arma_rnd([0.2,0.1], [0.3], 0.1, 10000, 10000);
octave:21> [phi, theta, v] = asc_armafit(x, 2, 1)
armafit(x,p,q)
p =  2
q =  1
--------------
phi =
 
   0.33490
   0.28179
 
theta =  0.26164
v =  0.10983

Estymacja modelem maksimum wiarygodności

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), \theta_0 = 1. \,\)

Przyczynowość zakłada, że dla \( |u| \le 1 \) część równania związana z procesem AR będzie niezerowa (\( 1 + \varphi_{1} u - \dots - \varphi_{p} u^p \ne 0 \)). Podobnie zakładamy z częścią odpowiedzialną za losowość w układzie dla \( |u| < 1 \), \( 1 + \theta_{1} u - \dots - \theta_{q} u^q \ne 0 \) (czyli zakładamy zaszumienie układu).

Zadaniem naszym jest oszacować wektory \( \hat{\bar{\varphi}} \) i \( \hat{\bar{\theta}} \) oraz wariancję szumu metodą największej wiarygodności. Musimy zatem zmaksymalizować funkcję wiarygodności Gaussa daną poprzez

\( L(\hat{\bar{\varphi}}, \hat{\bar{\theta}}, \hat{\sigma}) = \frac{1}{\sqrt{(2 \pi \hat{\sigma}^2)^n r_0 \cdots r_{n-1}}} \exp \Big[ - \frac{\hat{\sigma}^2}{2} \sum_{j=1}^n \frac{(X_j - \hat{X}_j)^2}{r_{j-1}}\Big].\)

W powyższym wzorze

  • \( \hat{X}_i \) to jednokrokowe predyktory \( X_i \) obliczane po wstępnym szacowaniu parametrów modelu,
  • \( r_i \) - to ich odpowiednie średniokwadratowe błędy (szacowania).

Predyktory te obliczamy z wzorów rekurencyjnych (patrz program asc_one_step_predictors)

\( \hat{X}_1 = 0, \)
\( \hat{X}_i = \sum_{j=1}^i \theta_{ij} (X_{i-j} - \hat{X}_{i-j}) , ~~~2 \le i < \max(p,q)\)
\( \hat{X}_i = \sum_{j=1}^p\varphi_{j} X_{i-j} + \sum_{j=1}^q \theta_{ij} (X_{i-j} - \hat{X}_{i-j}) , ~~~i \ge \max(p,q)\)

Błędy r szacowane są z kolei za pomocą wzorów rekurencyjnych (pamiętając, że \(\theta_{kj} = 0\) dla j > q i danego k)

\( r_1 = 1, \)
\( r_k = 1 - \sum_{j=0}^{k-1} \theta^2_{k,k-j} r_j, k = 2, \dots, n.\)
Algorytm 
  1. Wstępnie szacujemy wektory \( \hat{\bar{\varphi}} \) i \( \hat{\bar{\theta}} \) oraz wariancję szumu dowolnym algorytmem (rekurencyjnym dla ARMA, czy innowacyjnym dla modelu losowego).
  2. Obliczamy predyktory jednokrokowe i ich błędy dla konkretnych wartości \( \hat{\bar{\varphi}}, \hat{\bar{\theta}}.\)
  3. Liczymy funkcję \( S = \sum_{j=1}^n \frac{(X_j - \hat{X}_j)^2}{r_{j-1}}. \) Jest ona podstawą szacowania maksimum wiarygodności. Możemy teraz postąpić dwojako
    1. Znajdujemy minimum funkcji kwadratów błędów S.
    2. Znajdujemy minimum zredukowanej funkcji wiarygodności \( l = \ln(S/n) + \sum_{j=1}^n \ln r_{j-1} / n. \)
  4. Mając parametry minimalizujące jedną z powyższych funkcji kończymy obliczając wariancję szumu z wzoru \( \hat{\sigma}^2 = S / n \) (dla minimalnej wartości danej statystyki.

Poniżej znajdziecie program obliczający predyktory jednokrokowe oraz realizujący powyższy algorytm znajdywania parametrów modelu w oparciu o minimalizację wybranej statystyki (l lub S).

function [phi, theta, v, stat] = asc_armafit_likelihood(x, P, Q, method='likeli')
 
    if ( isempty (x) || (!isvector(x)) ) 
    error ("armafit: pierwszy agrument musi byc niepustym wektorem");
    return;
  endif
 
  %dlugosc wektora
  lx = length(x);
 
  if ( !isnumeric(P) || P < 0 || !isnumeric(Q) || Q < 0)
    error ("armafit: drugi i trzeci argument musi byc liczba naturalna p,q >= 0\n");
    return;
  endif
 
  v     = [];
  phi   = [];
  theta = [];
 
  S_params = zeros(P * Q, P + Q);
  S        = zeros(P*Q,1);
  indexes  = zeros(P*Q,2);
  for p = 1 : P
    for q = 1 : Q
      i = (p - 1) * Q + q;
 
      indexes(i,1:2) = [p,q];
 
      [phi_tmp, theta_tmp, v_tmp] = asc_armafit(x, p, q);
      S_params(i, 1 : p + q) = [phi_tmp', theta_tmp'];
 
      [xbar, r] = asc_one_step_predictors(x, phi_tmp, theta_tmp);
      S(i) = sum(power(x(2:end)' - xbar(2:end),2)./r(1:end-1));      
    endfor
  endfor
 
  if method == "square"   
    stat = S;
  else
    stat = log(S / lx) + sum(log(r)) / lx;    
  endif
 
  [S_value, S_index] = min(stat);  
  v     = min(S) / lx;
  phi   = S_params(S_index, 1 : indexes(S_index, 1));
  theta = S_params(S_index, indexes(S_index,1) + 1 : indexes(S_index, 1) + indexes(S_index, 2));
 
endfunction

function [xbar, r] = asc_one_step_predictors(x, phi, theta)
 
    if ( isempty (x) || (!isvector(x)) ) 
    error ("OSP: pierwszy agrument musi byc niepustym wektorem");
    return;
  endif
 
  %dlugosc wektora
  lx = length(x);
  q  = length(theta);
  p  = length(phi);
 
  if ( !isvector(phi) || isempty(phi) || !isvector(theta) || isempty(theta))
    error ("OSP: drugi i trzeci argument musi byc niepustym wektorem\n");
    return;
  endif
 
  xbar = zeros(1, lx);
  r    = zeros(1, lx);
  %m    = q;
  m    = max(q,p); %TO DO & RETHINK!
 
  if q == m
    vec = theta;
  else
    vec = phi;
  endif
 
  % m pierwszych predyktorów
  for i = 2 : m
    suma = 0;
    for j = 1 : i
      suma += vec(j) * (x(i+1-j) - xbar(i+1-j));
    endfor
    xbar(i) = suma;
  endfor
 
  % lx - m - 1 kolejnych predyktorów
  for i = m + 1: lx
    suma_phi   = 0;
    for k = 1 : p
      suma_phi += phi(k) * x(i + 1 - k);
    endfor
    suma_theta = 0;
    for j = 1 : q
      suma_theta += theta(j) * (x(i+1-j) - xbar(i+1-j));
    endfor 
    xbar(i) = suma_phi + suma_theta;
  endfor
 
  % wektor błędów predyktorów
  r            = [];
  theta_r      = zeros(1, lx);
  theta_r(1:q) = theta;
 
  r(1) = 1;
  for n = 2 : lx
    suma_r = 0;
    for j = 1 : n - 1
      suma_r += power(theta_r(n-j),2) * r(j);
    endfor
    r(n) = 1 - suma_r;
  endfor
 
endfunction


Źródło „http://172.16.3.1/ekonofizyka/index.php/Analiza_Szereg%C3%B3w_Czasowych/Techniki_analizy_szereg%C3%B3w_czasowych”