Istota metody największej wiarygodności

Idea metody największej wiarygodności zaproponowanej przez R. F 'skera (1912) jest następująca. Niech Y będzie badaną zmienną w populacji o roz a zie zależnym od nieznanego parametru 0 (lub wektora parametrów 0) oraz ( i, i, Yt) oznacza próbę prostą wylosowaną z tej populacji. Ponacto, niec^

f { y \ 6 ) oznacza funkcję prawdopodobieństwa lub gęstość zmiennej losowej zależną od parametru 0. Estymatorem uzyskanym metodą największej wiarygodności parametru 0 jest funkcja elementów próby maksymalizująca un cję postaci:

Funkcja (1) nosi nazwę funkcji wiarygodności. W praktyce często łatwiej wyznacza się maksimum logarytmu funkcji wiarygodności.

L ( y l , y 2, - , y T',0 ) = Y j lnf ^ y > ’^ (2 ) (=1

niż samej funkcji (1), przy czym należy zauważyć, iż obie te funkcje osiągają maksimum w tym samym punkcie. Jeżeli funkcję we wzorze (1) zastąpi się przez funkcję prawdopodobieństwa lub gęstości rozkładu należącego do liniowo-wykładniczej rodziny rozkładów (np. rozkładu gamma czy normalnego), to poprzez maksymalizację tak zdefiniowanej funkcji otrzymuje się MNPW - estymator parametru 0.2

Szczególny przypadek procesów generowanych przez przełącznikowy model Markowa stanowią mieszaniny niezależnych i identycznych rozkładów (i . i . d mixture distributions). Niech reżim, w którym znajduje się proces w chwili t będzie indeksowany przez nieobserwowalną zmienną losową st, która przyjmuje N możliwych wartości (st = 1, 2, ..., N). Innymi słowy proces stochastyczny generujący badane zjawisko jest mieszaniną kilku rozkładów normalnych o różnych średnich i wariancjach:

y , / O), ~ , er 2) (3)

gdzie:

■ nieobserwowalna zmienna st , będąca jednorodnym łańcuchem Markowa o N stanach i macierzy prawdopodobieństw przejścia P = \ p tJf]NxA, określa reżim, w jakim znajduje się zmienna y, w chwili t

■ zbiór <X>t zawiera wszystkie dostępne informacje dotyczące obserwowalnego procesu y, do chwili t włącznie.

Macierz prawdopodobieństw przejścia dla łańcucha M arkowa posiada następujące własności macierzy stochastycznej:

Pu P u ••• Pn\ P u P u P n 2

P =

_P\n P2N ••• P m N

Pij > 0 d l a i , j e { l , 2 , . . . , N } o raz ^ p ;j = 1 (4) j=i

gdzie:

Pij = Pr(s, = j / s,.i = i) oznacza prawdopodobieństwo przejścia ze stanu i w chwili t- 1 do stanu j w chwili t.

Macierz P wykorzystywana jest do określenia zależności pomiędzy rozkładami bezwarunkowymi nieobserwowalnej zmiennej reżimowej st w

por. Domański Cz., Pruska K., Nieklasyczne metody statystyczne PWE, Warszawa 2000, s. 126

90

kolejnych okresach:

N

P rO ,+, = j Z ® , ; 0 ) = '£ jP iJ -P r(sl = / / < D m ; 0 ) dl aj = 1, 2 ,..., N (5) i=i

Następnie definiuje się warunkową funkcję gęstości dla zmiennej obserwowalnej yt:

1 ~(y, ~ M )2

= — exP ( ^ } dl aj = 1,2, ...,N (6)

^ ItzcTj 2 a )

gdzie wektor 0 =(piI, p2, oj, o2, pn, p22) zawiera parametry modelu przełącznikowego postaci (3), które należy oszacować.

Z kolei warunkowa funkcja gęstości (6) określa funkcję gęstości rozkładu bezwarunkowego zmiennej y,:

/ O , / ' i 0) = X / C V i »ii = * / ; 9 )'=

i=i

E / O ' , / * / = i7 ® ,_ ,; 0 ) (7) i=i

Warunkowy rozkład zmiennej sterującej zmianami reżimu st można przedstawić w postaci ilorazu łącznej funkcji gęstości y, i s, przez funkcję gęstości rozkładu brzegowego:

P r ( i , = 27 r „ O , _ , ; « ) = ■' r ' ' -/ O , -/ i, = -/ ; 0) • Pr(j, = 7 / o , . , ; 0)

(8)

= / , O , . 1; 0 ) - P r ( i i = i7O>,_,;0) /=1

dlaj = 1 , 2 , . . . , N. Mając do dyspozycji obserwowalne dane dotyczące zmiennej y,, można wskazać z określonym przez powyższy szacunek prawdopodobieństwem, iż obserwacja w chwili t została wygenerowana w nieobserwowalnym reżimie j. W zależności od zakresu zbioru informacji O, w oparciu, o który przeprowadza się powyższe wnioskowania można wyróżnić dwa typy prawdopodobieństw określonych przez relację (8):

7 prawdopodobieństwa filtrowane (filter probability) wyznaczane na podstawie zbioru informacji Ot dla t < T

^ prawdopodobieństwa wygładzone (smoothed probability) wyznaczane na podstawie wszystkich informacji zawartych w próbie statystycznej 0 T.

Aby móc wnioskować na temat, który z reżimów generuje daną wartość zmiennej objaśnianej y, niezbędne jest uzyskanie oszacowań wektora parametrów 0.

Dysponując obserwacjami zmiennej y, dla t = 1, 2, ..., T można zbudować

funkcję wiarygodności:

L{ yx , y 2,...,yT;6) = Y J log / O , 1 '> &) (9) i=i

gdzie funkcję gęstości rozkładu bezwarunkowego zmiennej y, określa relacja (7).

Oceny wektora parametrów 0 można uzyskać m.in. m etodą największej wiarygodności w wyniku maksymalizacji wyrażenia (9), uwzględniając następujące ograniczenia dotyczące elementów tego wektora:

¿ P r (j, = j l <D,_,;0) = 1, Pr(i, = j l OM;^ )> 0 dlaj = 1 , 2 , N.

Udowodniono, iż ocena wektora parametrów 0 uzyskana metodą największej wiarygodności stanowi rozwiązanie następującego systemu nieliniowych równań:4

Y,y< -p rf a

Mj = J ń j — dla j = 1, 2, N (10)

¿ P r ( s , = j / y „ Ó ^ l;0) /=i

Z C y , - f i j ) 2 -VKs, = J / y l

c r ) = ^ F ---- dl aj = 1 , 2 , N (11)

¿ P r ( j f = j / y „ O M ;0 ) i=i

Zadana przez równość (3) mieszanina funkcji gęstości posiada tę własność, iż funkcja wiarygodności dla logarytmu tej gęstości określona wzorem (9) nie posiada globalnego maksimum. Osobliwość powstaje kiedy tylko na jeden z rozkładów nałoży się warunek, iż jego średnia jest dokładnie równa jednej z obserwacji (gi = yj, powiedzmy) a wariancja

zmierza do zera (cr,2 —» 0). W takim punkcie logarytmiczna funkcja wiarygodności staje się nieograniczona.

4 *

w przypadku gdy Pr(s, = j / y„ <J>M; 0 ) przyjmuje wartości pomiędzy zerem i jedynką dla pewnych obserwacji, ocena ¡J.- jest ważoną średnią wszystkich obserwacji pochodzących z danej próby, gdzie waga odpowiadająca obserwacji y, jest proporcjonalna do prawdopodobieństwa oznaczającego iż obserwacja w momencie t została wygenerowana przez reżim j. Ściślej, obserwacji, która musi pochodzić z reżimu j przypisuje się większą wagę w procesie estymacji parametru g,. Podobnie, cr2 jest ważoną średnią kwadratów odchyleń obserwacji yt od średniej jUj, podczas gdy p j jest w zasadzie frakcją obserwacji, które wydaje się, że pochodzą z reżimu j; por. Hamilton J. D , Time Series Analysis, Princeton University Press, Princeton, New Jersey 1994,s. 699 -700

92

p j = Pr(sf = y / O l. 1) = r ' , | ] P r ( j ( = j I y„<&,_p,6) dl aj = 1, 2 , N (12) (=i

Rezultat ten można osiągnąć poprzez zastosowanie standardowych metod numerycznych bądź algorytmu EM.

2. Algorytm EM (Expectations M aximization)

Jednym ze sposobów maksymalizacji skonstruowanej dla potrzeb przełącznikowego modelu Markowa funkcji wiarygodności (9) jest zastosowanie odpowiedniej wersji algorytmu Expectations Maximization, zaproponowanego pierwotnie przez Dempstera, Lairda, Rubina (1977).5 Zastosowanie algorytmu EM do klasy modeli przełącznikowych wymaga wyprowadzenia zależności pomiędzy oszacowanymi na podstawie dostępnych informacji prawdopodobieństwami filtrowanymi, że obserwacja yt została wygenerowana prze reżim j (8), a warunkami narzuconymi na parametry poprzez układ równań wynikający z potrzeby maksymalizacji funkcji wiarygodności (9):

ł - - i y , y T,0) = o S/A

- ^ T ( y „ ...,y T,d ) = 0 (13)

5<j,.

^ -(.y

, yr ,0) = o S p ,

dla i = 1, 2 , N.

Zależność ta pozwala na znalezienie w kolejnych iteracjach algorytmu coraz lepszych (w sensie kryterium związanego z funkcją wiarygodności) oszacowań wektora parametrów modelu. Po ustaleniu wartości początkowych dla wszystkich parametrów modelu, w każdej iteracji algorytmu EM wykonywane są dwa kroki. Pierwszy krok (expectations) polega na wyznaczeniu prawdopodobieństw filtrowanych P r(i, = j / y , , O , , ; # ) zgodnie z wzorem (8).

Krok drugi (maximization) prowadzi do wyznaczenia za pom ocą układu równań (13) wektora parametrów maksymalizującego funkcję wiarygodności, co stanowi rozwiązanie systemu nieliniowych równań (10) - (12). Procedurę wyznaczania ocen wektora parametrów 0 przedstawiono schematycznie na poniższym rysunku:

Rys. 1 Algorytm EM

Źródło: Diebold F. X., Lee J-H., Weinbach G. C., Regime switching with time-varying transition probabilities, w Hargreaves C. P. (ed.), Nonstationary Time Series Analysis and

Cointegration, Oxford University Press, New York 1994, s. 289

Procedurę rozpoczyna podanie dowolnej wartości początkowej wektora 0, oznaczonej jako 9(0), która służy do oszacowania wygładzonego prawdopodobieństwa Pr(st = j / y,; 0<O)) z równości (8). Z kolei oszacowane prawdopodobieństwo wykorzystuje się do estymacji parametrów wektora 0 opisanych równaniami (10) - (12). Oszacowane w ten sposób parametry f t j , <J2j , p j utworzą nowy wektor 0(1). Ta nowa ocena wektora parametrów 0(1) może zostać wykorzystana do ponownego wyznaczenia wartości Pr(st = j / yt; 0(1)) oraz wyrażeń zadanych przez relacje (10) - (12). Oszacowania parametrów

otrzymane dla wektora 0(1), utworzą z kolei nowy wektor 0(2).

Estymacja ustaje w momencie spełnienia zadanego kryterium zbieżności, np. gdy różnica pomiędzy 0(m+1) i 0(m) będzie mniejsza od zadanej z góry wartości.

Udowodniono, iż każda kolejna iteracja w tym algorytmie powiększa wartość funkcji wiarygodności. Oczywiście, jeśli w wyniku przeprowadzanych iteracji 94

zostanie osiągnięty punkt równowagi taki, że 0(m) = 0(m+1), oznacza to, iż algorytm znalazł taką ocenę 0 , która maksymalizuje funkcję wiarygodności. Hamilton (1990) udowodnił, że ciąg otrzymanych w ten sposób oszacowań jest zbieżny do lokalnego maksimum funkcji wiarygodności.6

3. Podejście bayesowskie do estymacji parametrów przełącznikowego modelu