Model autoregresyjny (AR)

AS/ Model autoregresyjny (AR)

Model autoregresyjny (rzędu [math]M[/math]) opisuje procesy dyskretne, w których wartość sygnału w danej chwili jest sumą liniowej kombinacji [math]M[/math] wartości poprzednich i nieskorelowanego szumu [math]\epsilon[/math]

[math] s[n] = \sum_{i=1}^M a_i s[n-i] + \epsilon[n] [/math]

W każdej realizacji tego samego procesu (dla tych samych współczynników [math]a_i[/math] i wartości początkowych sygnału), [math]\epsilon_t[/math] są niezależnymi liczbami losowymi, więc o wartości [math]s(t)[/math] w konkretnej chwili [math]t[/math] możemy mówić tylko językiem prawdopodobieństwa.

Trzy przykładowe realizacje procesu AR 3-go rzędu ([math]M=3[/math]) o tych samych współczynnikach i wartościach początkowych.

Mimo tego, na podstawie współczynników AR możemy określić wiele ogólnych własności sygnału, np. wartość oczekiwaną [math]\bar{s}[/math] (w praktyce estymowaną przez wartość średnią) i wariancję (jej estymatorem jest suma kwadratów odchyleń wartości sygnału od wartości oczekiwanej), a nawet widmo mocy. Można również rozważać szersze klasy modeli tego typu, jak np. model MA (ruchomej średniej, ang. moving average), gdzie uśredniamy [math]\epsilon_t[/math] zamiast [math]s(t)[/math], czy proces mieszany ARMA, opisany między innymi w klasycznych pozycjach „Analizie szeregów czasowych”, autorstwa Boxa i Jenkinsa oraz w „Metodach analizy szeregów czasowych” autorstwa Piersola i Bendata.

AR(1)

Najprostszym przykładem jest proces AR pierwszego rzędu (nazywany liniowym procesem Markowa), w którym wartość w danej chwili zależy wyłacznie od wartości w chcili poprzedniej i szumu: [math] s[n] = a s[n-1] + \epsilon_n [/math]

podstawiając trzy kolejne wyrazy

[math]s[n] = \epsilon_n + a s[n-1] [/math]

[math]s[n-1] = \epsilon_{n-1} + a s[n-2] [/math]

[math]s[n-2] = \epsilon_{n-2} + a s[n-3][/math]

dostaniemy

[math] s[n] = [/math]

[math] \epsilon_n + a s[n-1] = [/math]

[math] \epsilon_n + a \left( \epsilon_{n-1} + a s[n-2] \right) = [/math]

[math] \epsilon_n + a \left( \epsilon_{n-1} + a (\epsilon_{n-2} + a s[n-3]) \right) = [/math]

[math] \epsilon_n + a \epsilon_{n-1} + a^2 \epsilon_{n-2} + a^3 s[n-3] [/math]

W ogólnym przypadku [math]N[/math] wyrazów będzie to suma

[math] s[n] = \sum_{i=0}^{N-1} a^i \epsilon_{n-i} + a^N s[n-N] [/math]

Dla [math]N \rightarrow \infty[/math] zależność od pierwszego elementu [math]s[n-N][/math] zanika i dostejemy asymptotyczną reprezentację

[math] s[n] = \epsilon_n + a\epsilon_{n-1} + a^2\epsilon_{n-2} +\ldots = \sum_{i=0}^{\infty} a^i \epsilon_{n-i} [/math]

Jeśli wartość oczekiwana [math]\epsilon_i[/math] wynosi 0 ([math]E(\epsilon_i)=0[/math]) a wariancja [math]\sigma^2(\epsilon_i)=\sigma_\epsilon^2[/math], to wariancja w punkcie [math]n[/math]

[math]\begin{matrix} \sigma^2_{s[n]} = E\left( (\epsilon_n + a\epsilon_{n-1} + a^2\epsilon_{n-2}+\ldots+a^{n-1}\epsilon_1)^2\right) =\\ = \sigma_\epsilon^2 \left(1+a^2+a^4+\ldots+a^{2n-2} \right) = \left\{ \begin{matrix} \sigma_\epsilon^2 \left(\frac{1-a^{2n}}{1-a^2} \right) & |a|\ne 1\\ n \sigma_\epsilon^2 & |a|=1 \end{matrix} \right. \end{matrix}[/math]

Autokowariancja [math]E(s[n] s[n+\tau])[/math]

[math]\begin{matrix} E\left( (\epsilon_n + a\epsilon_{n-1} + a^2\epsilon_{n-2}+\ldots+a^{n-1}\epsilon_1) (\epsilon_{n+\tau} + a\epsilon_{n+\tau-1} +\ldots+a^{n+\tau-1}\epsilon_1)\right) =\\ = \sigma_\epsilon^2 \left(a^\tau+a^{\tau+2}+\ldots+a^{\tau+2(n-1)} \right) = \left\{ \begin{matrix} \sigma_\epsilon^2 a^\tau \left(\frac{1-a^{2n}}{1-a^2} \right) & |a|\ne 1\\ n \sigma_\epsilon^2 & |a|=1 \end{matrix} \right. \end{matrix}[/math]

Dla [math]|a|\ne 1[/math] przy [math]n\rightarrow\infty[/math] [math] \sigma^2_{x[n]} \stackrel{n\rightarrow\infty}{\longrightarrow} \frac{\sigma^2_\epsilon}{1-a^2} \;\;\; ; \;\;\; \sigma_{x[n], x[n+\tau]} \stackrel{n\rightarrow\infty}{\longrightarrow} \frac{\sigma^2_\epsilon a^\tau}{1-a^2} [/math]

Autokowariancja [math] \rho(\tau) = \frac{ \sigma_{x[n], x[n+\tau]} }{ \sigma^2_{x[n]} } \stackrel{n\rightarrow\infty}{\longrightarrow} a^{|\tau|} [/math]

Proces jest asymptotycznie stacjonarny do rzędu 2, czyli wariancja i średnia nie zależą od czasu.

Dla [math]a=1[/math] proces ten obrazuje tzw. błądzenie przypadkowe.

Na podstawie znajomości samego współczynnika [math]a[/math] modelu AR(1) policzyliśmy np. funkcję autokorelacji modelu, co daje już znajomość widma procesu (z przytoczonego poniżej twierdzenia Wienera-Chinczyna). Podobnie w procesach wyższych rzędów (1) znajomość współczynników [math]\{a_i\}_{i=1..M}[/math] daje nam dokładną wiedzę o własnościach generowanych przez nie procesów, bez znajomości sygnału [math]s[n][/math], którego wartości mogą różnić się w kolejnych realizacjach ze względu na element stochastyczny — szum [math]\epsilon[/math].

W praktyce analizy sygnału postępujemy odwrotnie — do konkretnej realizacji dopasowujemy model AR. Głównym problemem jest wybór rzędu modelu, estymacja współczynników [math]a_i[/math] najlepiej pasujących do danego sygnału posiada stabilne rozwiązania.

Jeśli dozwolimy, aby sygnał zależał również bezpośrednio od poprzednich wartości szumu [math]\epsilon[/math], dostajemy pełną postać procesu ARMA(L,M) (ang. auto-regressive moving average):

[math] \sum_{i=1}^L b_i\epsilon_{n-i} = \sum_{j=1}^M a_j s[n-j] [/math]

Twierdzenie Wienera-Chinczyna

Transformata Fouriera funkcji autokorelacji jest równa kwadratowi modułu transformaty Fouriera.

Dowód Kładąc [math]f = g[/math] we wzorze na funkcję korelacji sygnałów f i g, dostajemy

[math] \mathcal{F} \left( \int_{-\infty}^{\infty} f(t) f(t+\tau) dt \right) = [/math] [math] \int_{-\infty}^{\infty} e^{-i\omega \tau} \left( \int_{-\infty}^{\infty} f(t) f(t+\tau) dt \right) d\tau = [/math] [math] \int_{-\infty}^{\infty} e^{-i\omega(t+\tau)} e^{i\omega t} \int_{-\infty}^{\infty} f(t) f(t+\tau) dt d\tau = [/math] [math] \int_{-\infty}^{\infty} e^{-i\omega(t+\tau)} f(t+\tau) d\tau \int_{-\infty}^{\infty} e^{i\omega t} f(t) dt = [/math] [math] \hat{f}(\omega) \overline{\hat{f}(\omega)} = |\hat{f}(\omega)|^2 [/math]

Parametryczna estymacja widma mocy sygnałów

Pokazaliśmy powyżej (na przykładzie błądzenia przypadkowego), że znając współczynniki (parametry) modelu AR możemy z nich wyliczyć funkcję autokorelacji odpowiadającego im procesu, bez znajomości konkretnej realizacji sygnału. Z kolei z funkcji autokorelacji możemy z pomocą powyższego twierdzenia obliczyć widmo. To widmo będzie wyliczone a nie estymowane, ale nie odnosi się bezpośrednio do sygnału, od którego zaczynaliśmy, tylko do procesu opisanego wyestymowanymi parametrami modelu AR.

Wielozmienny model AR

Model AR opisuje wartość sygnału w chwili [math]t[/math] jako kombinację liniową jego wartości w chwilach poprzednich (oraz szumu). W przypadku wielowymiarowym możemy włączyć do tego opisu wartości wszystkich sygnałów [math]s_i[/math], czyli wektora [math]\vec{s}(t)[/math]. Wielozmienny model AR (MVAR, multivariate autoregressive ) można wówczas opisać wzorem:

[math] \vec{s}(t)=\sum_{i=1}^p A(i) \vec{s}(t-i) + \vec{\epsilon}(t) , [/math]

gdzie [math]\vec{\epsilon}(t)[/math] będzie wektorem szumów, zaś [math]A(i)[/math] będą macierzami współczynników modelu. Przechodząc do przestrzeni częstości otrzymamy:

[math] \vec{s}(\omega)=A^{-1}(\omega)\vec{\epsilon}(\omega)=H(\omega)\vec{\epsilon}(\omega), [/math]

gdzie [math]H(\omega)[/math] jest macierzą przejścia. MVAR jest modelem typu "czarna skrzynka", gdzie na wejściu występują szumy, na wyjściu sygnały, a system jest opisany przez macierz przejścia. Zawiera on informacje o własnościach widmowych sygnałów i związkach między nimi.

Na podstawie macierzy [math]H(\omega)[/math] można obliczyć macierz gęstości widmowej zawierającą widma mocy dla pojedynczych kanałów jak również funkcje wzajemnej gęstości mocy pomiędzy kanałami. Stosując tego typu podejście, w którym wszystkie sygnały generowane przez pewien proces są rozpatrywane jednocześnie, można policzyć z macierzy spektralnej nie tylko koherencje zwykłe pomiędzy dwoma kanałami, ale również koherencje wielorakie opisujące związek danego kanału z pozostałymi i koherencje cząstkowe opisujące bezpośrednie związki między dwoma kanałami po usunięciu wpływu pozostałych kanałów. W przypadku gdy pewien kanał 1 będzie wpływał na kanały 2 i 3, obliczając koherencję zwykłą znajdziemy związek między 2 oraz 3, chociaż nie są one ze sobą bezpośrednio powiązane, natomiast koherencja cząstkowa nie wykaże związku między nimi.

Macierz [math]H(\omega)[/math] jest niesymetryczna, a jej wyrazy pozadiagonalne mają sens przyczynowości Grangera, co oznacza, że uwzględnienie wcześniejszej informacji zawartej w jednym z sygnałów zmniejsza błąd predykcji drugiego sygnału. Opierając się na tej własności zdefiniowano Kierunkową Funkcję Przejścia (DTF, directed transfer function) jako znormalizowany element pozadiagonalny [math]H(\omega)[/math]. DTF opisuje kierunek propagacji i skład widmowy rozchodzących się sygnałów.

Otrzymamy w ten sposób całościowy opis zmian wszystkich sygnałów jednocześnie. Co ciekawe, obliczona na tej podstawie funkcja charakteryzująca zależności między sygnałami [math]s_i[/math] (funkcja przejścia) nie jest symetryczna, w przeciwieństwie do np. korelacji. Dzięki temu może służyć wnioskowaniu nie tylko o sile zależności między poszczególnymi sygnałami składowymi, ale też o kierunku przepływu informacji między nimi. W przybliżeniu odpowiada to informacji, w którym z sygnałów struktury odpowiadające danej częstości pojawiają się wcześniej.