Z Skrypty dla studentów Ekonofizyki UPGOW

Wersja z 17:29, 29 gru 2010

Analiza Szeregów Czasowych
<<< Modelowanie szeregów czasowych | Matlab / GNU Octave >>>

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 Durbina - Levinsona 

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.

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.

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;
    end %for k=0:q-1
 
  end %for q=1:qIn
 
  theta = teta;
  v     = w;
 
endfunction

ARMA

Rekurencyjny algorytm dopasowania ARMA

Estymacja modelem maximum Likelihood

Testy

AICC

Prognoza