Nieliniowa metoda najmniejszych kwadratów dla modelu trendu

W wielu pakietach ekonometrycznych można użyć Nieliniowej Metody Najmniejszych Kwadratów (NMNK), aby uzyskać bezpośrednio poszukiwany parametr. NMNK ma podobnie jak MNK podstawy teoretyczne - o ile słuszność MNK wynika z twierdzenia Gaussa-Markowa, o tyle NMNK jest słuszna z powodu uogólnionych twierdzeń Gaussa-Markowa [1, 2, 3], a także dowiedzionej asymptotycznej zgodności jej estymatora [4].

NMNK nie jest jednolitą metodą. W Gretlu stosowany jest algorytm Levenberga–Marquardt. Algorytm ten jest kombinacją dwóch znanych metod numerycznych: metody gradientu prostego oraz algorytmu Gaussa-Newtona. W metodzie gradientu prostego suma kwadratów błędów jest korygowana w kierunku spadku nachylenia funkcji (gradientu), bo tam znajduje się minimum. W metodzie Gaussa-Newtona wykorzystuje się z kolei przybliżenie funkcji za pomocą wzoru Taylora; kwadraty różnic pomiędzy prawdziwą funkcją a funkcją Taylora powinny być jak najmniejsze [5].

Na chwilę zostawmy to. Wcześniej (tu oraz tu) pokazałem, jak poprawnie retransformować model trendu, aby uzyskać warunkową oczekiwaną stopę zwrotu. Mamy następujący model trendu:

(1)

Powinniśmy zdefiniować składnik losowy. Powiedzmy, że jest on stacjonarny, posiada wartość oczekiwaną = 0 oraz stałą wariancję. Jeżeli zmienna losowa P(T) ma rozkład log-normalny, to jej wartość oczekiwana wynosi:

(2)

Dzięki NMNK możemy oszacować parametry bezpośrednio z modelu (1), a więc bez potrzeby transformacji logarytmicznej. Przykładowo, biorąc te same dane dla mbanku co w wymienionych na początku artykułach (rocznie 1994-2015; 22 obserwacje), możemy stworzyć w Gretlu model (1). Po zaimportowaniu danych, trzeba kliknąć: Model -> Nieliniowa Metoda Najmniejszych Kwadratów. Pojawi się okno i tam wpisujemy:
.................
scalar a = 1
scalar b = 1
Zamkniecie = exp(a+b*time)
deriv a = exp(a+b*time)
deriv b = exp(a+b*time)*time
.................

gdzie:
Zamkniecie - to szereg czasowy, w naszym przypadku kurs mbanku (zmienna P(T));
time - zmienna czasowa (wygenerowana w opcji Dodawanie zmiennych -> time - zmienna czasowa.
a, b - współczynniki regresji, które trzeba od początku zadeklarować. Są to skalary w naszej funkcji, które przyjmą jakąś początkową wartość, np. 1. Standardowo deklarujemy więc scalar a = 1, scalar b = 1. Mogą być one dowolne, ale może się zdarzyć, że nie wystarczy iteracji, aby zbiegły do prawdziwej wartości wg algorytmu. Jeśli tak się stanie, zmieniamy metodą prób błędów.

deriv - pochodna funkcji po danym współczynniku. Np. pochodna d(Zamkniecie)/d(a) = exp(a+b*time), dlatego jest zapis:
deriv a = exp(a+b*time).
I tak samo dla b: pochodna d(Zamkniecie)/d(b) = exp(a+b*time)*time. Dlatego jest zapis:
deriv b = exp(a+b*time)*time.

Mimo że a i b są to skalary modelu, to są one w funkcji szacowane, czyli pierwotnie są zmiennymi.

Niezbyt przyjemne. Można się obyć bez pochodnej, ale w helpie Gretl rekomenduje jej obliczanie.

Przed kliknięciem OK, możemy zaznaczyć jeszcze opcję "Odporne błędy standardowe". Wtedy OK.

Jak widać na tablicy zbieżność nastąpiła po 27 iteracjach. Gdybyśmy wpisali scalar a = 10, scalar b = 10, zbieżność nastąpiłaby dopiero po 240 iteracjach.

Uzyskany parametr b = 0,0996, natomiast wolny parametr a możemy zignorować, bo interesują nas stopy zmian. W porównaniu do logarytmicznego modelu, gdzie uzyskaliśmy logarytmiczną stopę ponad 11%, może się wydać dziwne, że teraz jest niższa wartość. Przecież otrzymaliśmy od razu stopę, która powinna być właśnie wyższa od tej z modelu logarytmicznego MNK. Co jest tego przyczyną? Stało się tak, bo nasz model jest niepełny. Założyliśmy błędnie, że wygląda jak na poniższym rysunku:

W tym modelu składnik losowy jest stacjonarny, może więc to być zmienna losowa o rozkładzie normalnym. Ta zmienna po prostu jest argumentem funkcji wykładniczej. Ale w rzeczywistości powinna być funkcją czasu. Gdy za tę zmienną losową wstawimy proces Browna, czyli skumulowaną zmienną losową o rozkładzie normalnym, to otrzymamy geometryczny proces ruchu Browna. Wtedy nasz model wygląda mniej więcej tak:

Ogólnie geometryczny proces ruchu Browna ma postać (zob. np. [6]):

(3)

gdzie B(T) jest zwykłym procesem ruchu Browna. B(0) = 0.

Proces (3) w danym punkcie czasu może przyjąć różne wartości zgodnie z pewnym rozkładem prawdopodobieństwa. Przestaje być więc rozumiany jako proces, a staje się "ruchem" - zmienną stacjonarną. Ten ruch posiada zawsze rozkład log-normalny. Oznacza to, że możemy wyprowadzić analogicznie do (1) wartość oczekiwaną tej zmiennej.

Zauważamy, że składnik losowy w modelu (1) to jest to samo co (odchylenie standardowe razy proces ruchu Browna) w modelu (3), co oznacza, że musi się tak samo przekształcać jak w modelu (2), gdzie występuje wartość oczekiwana procesu. I stąd podstawiając do (3) powyższe wyprowadzenie wariancji, dostaniemy wartość oczekiwaną geometrycznego (procesu) ruchu Browna:

(4)

W tym miejscu rodzi się pytanie czy rzeczywiście znamy wartość początkową P(0). To zależy od tego czy zaczynamy od T = 1 czy T = 0. Musimy zdecydować czy pierwsza wartość już decyduje o funkcji zależnej od czasu czy dopiero kolejna. Jeśli założymy, że P(0) jest znane, to podstawimy pierwszą cenę zamknięcia mbanku z 1994 r., czyli 30, 68. W przeciwnym wypadku P(0) musi zostać oszacowane. Wybór sposobu ma istotne znaczenie dla uzyskanych wyników: dodanie wolnego wyrazu zmienia postać modelu. Ale jeżeli wybierzemy model z wyrazem wolnym, to estymacja NMNK przyniesie ten sam rezultat co dla modelu (2) - czyli powracamy do modelu (1), który przecież uznaliśmy za niepełny. Wstawienie P(0) nie rozwiąże tego problemu, ale możemy sprawdzić do czego doprowadzi.

Wstawiamy pierwszą cenę zamknięcia P(0) = 30,68. Możemy wtedy powtórzyć całą procedurę w Gretlu, używając modelu (4):

.................
scalar b = 1
Zamkniecie = 30.68*exp(b*time)
deriv b = 30.68*exp(b*time)*time
.................

Uzyskane parametry to:

Dostaniemy wykres trendu:

Zatem u = 12,5, czyli efektywna oczekiwana stopa exp(12,5) - 1 = 13,3%.

Mimo iż obliczyliśmy parametr u, to taki estymator jest niezgodny w tym sensie, że nie będzie dążył do prawdziwej wartości parametru u. Żeby była zgodność, składnik losowy musi być stacjonarny (zob. twierdzenia w [4]). Najprawdopodobniej potrzeba użyć uogólnionej nieliniowej metody najmniejszych kwadratów, która pozwoliłaby uwzględnić rosnącą wariancję i autokorelację składnika losowego.

Ostatecznie można ten problem rozwiązać. Łatwo zauważyć, że tylko w punkcie dla T = 1 geometryczny proces ruchu Browna będzie tożsamy z modelem (1), dla którego składnik losowy jest stacjonarny. W tym punkcie obydwie krzywe się przetną. Czyli w tym punkcie model (2) równa się modelowi (4):

Wyrazy wolne możemy do siebie dostosować, aby zachować identyczną interpretację punktu startu. Dlatego przyjmiemy lnP(0) = a. Stąd:

(5)

Ze względu na stacjonarność składnika losowego dla T = 1 estymator u równy b plus połowa wariancji może stać się w końcu zgodny. Wcześniej oszacowany współczynnik b był estymatorem zgodnym (zakładaliśmy początkowo stacjonarność składnika losowego), chociaż występował w modelu niepełnym. Połowa wariancji też jest estymatorem zgodnym. Suma estymatorów zgodnych daje estymator zgodny (wynika to z własności granicy ciągów). Stąd parametr u równy b plus połowa wariancji stanie się estymatorem zgodnym. Parametr b został wcześniej oszacowany, ale wariancję również znamy, bo jest to przecież ten sam składnik, który był obliczany dla MNK. Estymowaliśmy go wcześniej za pomocą błędu standardowego reszt (tu) i wyniósł 0,346. Możemy podstawić:

u = b + 0,5var
u = 0,0996 + 0,5*0,346^2 = 0,1596

Aby uzyskać oczekiwaną (efektywną) stopę zwrotu wstawimy exp(0,1596) - 1 = 17,3%. Dla porównania w Transformacja lognormalnego modelu z nieznanym parametrem wyszła wartość zbliżona, bo 18,7%.

Literatura:
[1] Louton, T., The Gauss-Markov Theorem for Nonlinear Models, Dec. 1982,
[2] Kariya, T., A Nonlinear Version of the Gauss-Markov Theorem, Jun. 1985,
[3] Kariya, T., A Maximal Extension of the Gauss–Markov Theorem and Its Nonlinear Version, 2002,
[4] Wu, C-F, Asymptotic Theory of Nonlinear Least Squares Estimation, May 1981,
[5] Marquardt, D. W., An Algorithm for Least-Squares Estimation of Nonlinear Parameters, Jun. 1963.
[6] https://en.wikipedia.org/wiki/Geometric_Brownian_motion