Symulacja i automatyczne modelowanie ARIMA w języku R

W języku R można symulować i modelować procesy ARIMA(p, d, q) za pomocą wielu pakietów, ale zaczniemy od funkcji wbudowanej w sam R. 

Symulacja

Symulacje wykonamy za pomocą funkcji arima.sim. 

Wybór ziarna

Do odtwarzania zawsze tych samych wyników na liczbach pseudolosowych służy tzw. ziarno, czyli wartość startowa, która determinuje cały ciąg liczb udających losowe. Możemy więc wybrać dowolną liczbę całkowitą i uzyskamy nowe ziarno, a tym samym inny przebieg losowy, który jednak możemy za każdym razem odtwarzać. Podstawowa funkcja to set.seed(x), gdzie x to dowolna liczba całkowita.  Czy liczba ziarn jest nieskończona? Nie. W pomocy do funkcji set.seed dowiemy się, że domyślnie używany jest stosunkowo nowy algorytm z 2002 r., którego okres wynosi 2^19937 - 1, czyli tyle jest unikatowych ziarn, a potem powtarzają się od nowa. Tyle w teorii. A w praktyce taka liczba ziarn jest niemożliwa dla programu, gdyż to oznaczałoby liczbę składającą się z 6002 cyfr. Z moich testów wynika, że obecna wersja programu (64bit) jest ograniczona do paru miliardów. W wersji której używam dokładnie jest 4 294 967 295 ziaren. Ponieważ argumenty funkcji seed() nie muszą być dodatnie, to tę liczbę należy rozbić:  2147483647*2 + 1, z czego 2147483647 dotyczy liczb dodatnich oraz ujemnych, stąd mnożymy razy 2, a dodanie 1 dotyczy argumentu 0. Oznacza to ni mniej ni więcej, że pierwsze ziarno otrzymamy za pomocą funkcji seed(-2147483647), a największe seed(2147483647). Możemy to sprawdzić sami. Stwórzmy próbkę z rozkładu normalnego o średniej 0 i wariancji 1:

set.seed(2147483646)

sim1 = rnorm(100, 0, 1)

plot(sim1, type="l")

To samo, ale zwiększmy seed o 1:

set.seed(2147483647)

sim1 = rnorm(100, 0, 1)

plot(sim1, type="l")

Wykres zmienił się. A teraz:

set.seed(2147483648)

sim1 = rnorm(100, 0, 1)

plot(sim1, type="l")

Otrzymamy informację, że wystąpił błąd. W zależności od wersji, którą posiadamy, program może się inaczej zachować i albo w ogóle nie zmienić wykresu, albo zresetować ziarno. To samo dotyczy liczb ujemnych. We wszystkich symulacjach musimy więc ograniczyć argument ziarna od -2147483647 do 2147483647. Jednak w poniższych przykładach będę używał zawsze tego samego ziarna - seed(1) - które przy zmienianych parametrach modelu i tak będzie przynosić zupełnie różne szeregi. 

Przykład 1: brak autokorelacji

AR(0) to proces czysto losowy o rozkładzie normalnym. Aby go zasymulować, stosowałem wcześniej funkcję rnorm. Dzięki ziarnu zobaczymy, że arima i rnorm dadzą dokładnie to samo. Przykładowo bierzemy pierwsze ziarno i generujemy 100 danych za pomocą arima.sim

set.seed(1)

sim1 = arima.sim(list(order=c(0,0,0)), n=100)

plot(sim1, type="l")

Otrzymamy dokładnie ten sam wykres, gdy wykonamy kod z rnorm:

set.seed(1)

sim1 = rnorm(100,0,1)

plot(sim1, type="l")

Za pomocą funkcji acf i pacf możemy sprawdzić autokorelacje próby. O różnicy między ACF i PACF pisałem kiedyś w art. Jak sprawdzić stabilność modelu ekonometrycznego w czasie?. Zatem wiemy, że ACF służy do oceny stacjonarności, a PACF do oceny autokorelacji danego rzędu. W tym przypadku wpiszmy najpierw

acf(sim1)

Samo takie wpisanie wywoła wykres:

Zerowy rząd, który zawsze będzie równy 1 pomijamy, więc widać, że proces jest stacjonarny.

Potem 

pacf(sim1)

Brak autokorelacji - tak jak być powinno.

Przykład 2: autoregresja rzędu 1: AR(1). Powiedzmy, że współczynnik AR(1) - phi - wynosi 0,7. Stosujemy kod:

set.seed(1)

sim2 = arima.sim(list(order=c(1,0,0), ar=0.7), n=100)

plot(sim2, type="l")

Wykres:

Aby pokazać ACF i PACF obok siebie zastosujemy funkcję par(). Funkcja ustawia parametry grafiki w programie. Ponieważ robi to na stałe w sesji, to użyjemy jej dwukrotnie: raz do ustawienia dwóch wykresów, dwa do powrotu do ustawień domyślnych. Dodatkowo ustawimy tylko 5 rzędów acf i pacf, bo więcej nie potrzeba. To się robi za pomocą parametru lag.max:  

par(mfrow=c(2,1), mar=c(2,5,2,2))

acf(sim2, lag.max=5)

pacf(sim2,lag.max=5)

par(mfrow=c(1,1), mar=c(5, 4, 4, 2) + 0.1)

Jak widać funkcja acf jest niespójna z pacf, bo w tej drugiej prawidłowo od razu zaczynamy od rzędu 1, a nie bez sensu od 0. W każdym razie ACF wskazuje na występowanie dwóch rzędów autokorelacji, co jednak nie przekłada się na liczbę rzędów AR, która zostaje wykryta dopiero przez pacf. Im silniejszy współczynnik AR, tym większa będzie liczba rzędów na ACF, nie wpływając na liczbę rzędów PACF. Gdy współczynnik przekroczy 1, proces staje się niestacjonarny. Popatrzmy na kolejny przykład.

Przykład 3. AR(1) z siłą phi = 1,3. Najpierw wpiszmy kod:

sim3 = arima.sim(list(order=c(1,0,0), ar=1.3), n=100)

Wynik:

Otrzymaliśmy błąd, ponieważ AR nie może być niestacjonarny w modelu ARIMA. Aby uzyskać AR(1), musimy go sprowadzić do stacjonarności i do tego właśnie służy drugi parametr, d. Otrzymamy model autoregresyjny rzędu 1 z integracją stopnia 1, tj. ARI(1, 1) lub ogólniej ARIMA(1, 1, 0). Ze względu na integrację zamiast siły 1,3 wpisujemy różnicę o 1, czyli 0,3:

sim3 = arima.sim(list(order=c(1,1,0), ar=0.3), n=100)

plot(sim3, type="l")

I wykres będzie typowo rynkowy:

Po wykonaniu jeszcze raz kodu:

par(mfrow=c(2,1), mar=c(2,5,2,2))

acf(sim3, lag.max=5)

pacf(sim3,lag.max=5)

par(mfrow=c(1,1), mar=c(5, 4, 4, 2) + 0.1)

 dostaniemy coś ciekawego:

ACF wskazuje na niestacjonarność (powolne spadki autokorelacji), podczas gdy PACF informuje nadal tylko o 1 rzędzie. Można by zapytać dlaczego, skoro mamy tu do czynienia z procesem niestacjonarnym? No właśnie dzieje się tak, bo to co tutaj mamy to AR(1) z p >= 1, więc całkowita autokorelacja wskazuje na bardzo dużą liczbę rzędów (cena z dawnej przeszłości koreluje z ceną teraźniejszą), ale żeby dostać autokorelację tylko dla jednego rzędu, trzeba usunąć wszystkie inne, które na nią wpływają. Tak działa PACF. Czyli zauważmy, że mimo iż przeszliśmy z procesu stacjonarnego do niestacjonarnego tylko ACF się bardzo zmieniło, a PACF zmieniło jedynie poziom samej autokorelacji rzędu 1.

Przykład 4: średnia ruchoma rzędu 1: MA(1). Powiedzmy, że współczynnik MA(1) - theta  - wynosi 0,7. Stosujemy kod:

set.seed(1)

sim4 = arima.sim(list(order=c(0,0,1), ma=0.7), n=100)

plot(sim4, type="l")

par(mfrow=c(2,1), mar=c(2,5,2,2))

acf(sim4, lag.max=5)

pacf(sim4,lag.max=5)

par(mfrow=c(1,1), mar=c(5, 4, 4, 2) + 0.1)

Po pierwsze - mimo że ciągle używamy seed(1), to wynik znowu jest inny dla MA(1) niż dla AR(1). Po drugie zmieniła się struktura PACF - korelacja 1-go rzędu jest nadal wysoka, ale pojawiła się drugiego i czwartego ujemna. To świadczy o tym, że MA jest dużo mniej intuicyjna niż AR, ale z drugiej strony obejmuje bardziej skomplikowane procesy. Fakt, że mamy tu ujemną autokorelację 2-go rzędu nie jest przypadkiem. Dla lepszej ilustracji oraz dalszej analizy zróbmy to samo, ale zwiększmy próbę do 1000 oraz liczbę rzędów autokorelacji do 15. Za n wstawiamy do kodu 1000, reszta taka sama. Efekt:

Jak się okazuje PACF realizuje się w postaci na przemian dodatniej i ujemnej autokorelacji. Z czego to wynika? Popatrzmy, MA(1) przy zerowym wyrazie wolnym ma postać:

(1)

W okresie 0 składnik losowy z okresu -1 jest równy zero i wtedy:

Dla okresu 1 podstawiamy do (1) i wyznaczamy składnik losowy:

W okresie 2 i 3 analogicznie:

W końcu w okresie 4 na podstawie poprzednich danych już tylko wyznaczamy y(4):

Zauważmy, że dostaliśmy właśnie na przemian znaki plus i minus dla kolejnych rzędów zależnych od y z okresów poprzednich. Ten wzór będzie się powtarzał do nieskończoności, a to znaczy, że MA(1) można zapisać w postaci procesu AR, w którym następuje jakby częściowe "wygaszanie" wpływu poprzednich obserwacji. 

MA jest procesem ergodycznym i stacjonarnym (w przeciwieństwie do AR, który być nie musi). Zarówno ergodyczność jak i stacjonarność wiążą się prawem regresji do średniej. Sama stacjonarność nie wystarcza dla niego, bo może być proces stacjonarny bez żadnych zmian, jakaś stała [1]. Dopiero ergodyczność wymusza zmienność, gdyż wiąże się z rozkładem częstości. Ergodyczność zapewnia, że prędzej czy później zmienna zapełni cały swój rozkład, tak że zawsze uzyskamy np. kształt dzwonu. W każdym razie zawsze dostaniemy taki obraz PACF dla MA(1) z powodu powrotu do średniej, która zostaje jakby wzmocniona zarówno w stosunku do AR jak i procesu czysto losowego.

Można by jeszcze zapytać dlaczego MA nazywana jest średnią ruchomą albo średnią kroczącą. Pomijając względy historyczne, bo nie jest to rzeczywiście średnia, to zauważmy, że:

* dostaliśmy wyżej model zależny liniowo od wartości poprzednich, a współczynniki można porównać z wagami średniej ważonej (chociaż nie muszą się sumować do 1, dlatego to nie jest prawdziwa średnia),

* mimo że występuje zamiana znaków plus na minus, co wydaje się psuć ideę średniej kroczącej, to zauważmy, że MA dotyczy zmiennej stacjonarnej, a więc zawsze przyspiesza prawo regresji do średniej. Co więcej, możemy utworzyć MA kolejnych rzędów. To znaczy, że dla MA(2) drugi rząd będzie się zachowywał tak samo jak pierwszy z przesunięciem o 1, a to oznacza, że dla okresu 2 znak minus z MA(1) będzie jakby redukowany przez plus z MA(2). 

Przykład 5. ARMA(1, 0, 1). phi(1) = 0,7 i theta(1) = 0,7. Kod:

set.seed(1)

sim5 = arima.sim(list(order=c(1,0,1), ar=0.7, ma=0.7), n=1000)

plot(sim5, type="l")

par(mfrow=c(2,1), mar=c(2,5,2,2))

acf(sim5, lag.max=15)

pacf(sim5,lag.max=15)

par(mfrow=c(1,1), mar=c(5, 4, 4, 2) + 0.1) 

Znowu możemy dokładnie obserwować zmiany procesu, bo ziarno nie zmieniło się. Dostaliśmy mix AR i MA, dzięki czemu proces ma własności jednego i drugiego.

Przykład 6. ARIMA(1, 1, 1). Siła AR(1) = 0,7 i siła MA(1) = 0,7 oraz d = 1. Kod:

set.seed(1)

sim6 = arima.sim(list(order=c(1,1,1), ar=0.7, ma=0.7), n=1000)

plot(sim6, type="l")

par(mfrow=c(2,1), mar=c(2,5,2,2))

acf(sim6, lag.max=15)

pacf(sim6,lag.max=15)

par(mfrow=c(1,1), mar=c(5, 4, 4, 2) + 0.1)

Wprowadzenie jednoczesne AR i MA wywołało autokorelację rzędu 1 praktycznie równą jedności. I co więcej, nawet jeśli zmniejszymy parametr MA niewiele się zmieni. W ogóle wykres się praktycznie nie zmieni. Zobaczymy to w kolejnym przykładzie.

Przykład 7. ARIMA(1, 1, 1): phi(1) = 0,7 i theta(1) = 0,1 oraz d = 1. Kod:

set.seed(1)

sim7 = arima.sim(list(order=c(1,1,1), ar=0.7, ma=0.1), n=1000)

plot(sim7, type="l")

par(mfrow=c(2,1), mar=c(2,5,2,2))

acf(sim7, lag.max=15)

pacf(sim7,lag.max=15, plot=FALSE)

par(mfrow=c(1,1), mar=c(5, 4, 4, 2) + 0.1)

Przykład 8. Ciekawszy wynik dostaniemy, jeśli ustawimy phi na 0.1 thetę na -0.95. Funkcja PACF występuje dwukrotnie, pierwsza pokazuje wykres, druga - same współczynniki korelacji:

set.seed(1)

sim8 = arima.sim(list(order=c(1,1,1), ar=0.1, ma=-0.95), n=1000)

plot(sim8, type="l")

par(mfrow=c(2,1), mar=c(2,5,2,2))

acf(sim8, lag.max=15)

pacf(sim8,lag.max=15)

par(mfrow=c(1,1), mar=c(5, 4, 4, 2) + 0.1)

Mimo iż wykres stał się bardzo poszarpany, to nadal występuje dodatnia korelacja  i to wielu rzędów. Nawet bardzo ujemna część MA nie przykrywa małej dodatniej części AR.

Optymalne modelowanie

Kolejnym etapem jest fitting, czyli dopasowanie modelu. Najłatwiej użyć optymalizacji za pomocą auto.arima(). Jest to funkcja zawarta w pakiecie forecast, który należy zainstalować:

install.packages("forecast")

i włączyć

library(forecast)

Sprawdzimy, jak auto.arima radzi sobie ze wszystkimi podanymi wyżej przykładami. 

Ad 1: brak autokorelacji. Dla sim1:

dopas=auto.arima(sim1)

dopas

Otrzymamy:

Series: sim1
ARIMA(0,0,0) with zero mean

Prawidłowo.

Ad 2: phi(1) = 0,7. Za sim1 wstawiamy sim2, więc nie ma sensu powtarzać kodu. Wyniki:

Series: sim2
ARIMA(1,0,0) with zero mean

Coefficients:
ar1
0.6956
s.e. 0.0712

Prawidłowo.

Ad 3. phi = 1,3.

Series: sim3
ARIMA(1,1,0)


Coefficients:
ar1
0.3116
s.e. 0.0958

Prawidłowo.

Ad 4. theta = 0.7

Series: sim4 

ARIMA(0,0,2) with zero mean 

Coefficients:

         ma1      ma2

      0.6609  -0.0724

s.e.  0.0318   0.0345

Pierwszy raz błąd. Powinno być MA(1), a dostaliśmy MA(2), pomimo, iż było 1000 danych. Gdy wrócimy do PACF dla przykładu 4, to zobaczymy, że występują istotne ujemne autokorelacje. Jest to niewątpliwie przyczyna błędu. Mimo wszystko powinniśmy uzyskać optimum dla MA(1). Co w takiej sytuacji zrobić? Warto przyjrzeć się p-value. Niestety pakiet tego nie oferuje. Żeby je uzyskać możemy użyć kodu:

(1-pnorm(abs(dopas$coef)/sqrt(diag(dopas$var.coef))))*2

Wynik:

       ma1        ma2 

0.00000000 0.03583408  

p-value dla ma1 = 0 i ma2 = 0.0358. Nie można więc odrzucić modelu MA(2). 

Tak więc już napotykamy pierwsze problemy. Pakiet auto.arima domyślnie głównie opiera się na AIC, ale wg moich testów BIC lepiej się sprawdza dla ARIMA. Algorytm można jednak ustawić tak, aby kierował się tylko BIC. Kod musimy zmodyfikować. Domyślny kod miał ustawienia następujące:

dopas=auto.arima(sim1, ic = c("aicc", "aic", "bic"))

Zastąpimy go:

dopas=auto.arima(sim1, ic = c("bic"))

dopas

Wynik:

Series: sim4 

ARIMA(0,0,1) with zero mean 

Coefficients:

         ma1

      0.7039

s.e.  0.0258

Okazuje się, że kryterium BIC prawidłowo pokazało, że najlepszym modelem jest MA(1).

Ad 5. ARIMA(1, 0, 1). phi(1) = 0,7 i theta(1) = 0,7. Ponowne zastosowanie z BIC daje wyniki

ARIMA(1,0,1) with zero mean 

Coefficients:

         ar1     ma1

      0.6404  0.7329

s.e.  0.0266  0.0268

Czyli prawidłowo. Dla domyślnych ustawień ponownie algorytm nie działał dobrze, dając ARIMA(2, 0, 1).

Ad 6. ARIMA(1, 1, 1). phi(1) = 0,7 i theta(1) = 0,7. To samo, ok.

Ad 7. ARIMA(1, 1, 1). Siła AR(1) = 0,7 i siła MA(1) = 0,1. Dla BIC:

Series: sim7 

ARIMA(1,1,0) 

Coefficients:

         ar1

      0.7030

s.e.  0.0225

Tym razem BIC nie poradził sobie, najwyraźniej theta = 0,1 była za małą wartością. Ale phi jest właściwa. Sprawdźmy jeszcze domyślne kryteria:

Series: sim7 

ARIMA(2,1,2) 

Coefficients:

          ar1     ar2     ma1     ma2

      -0.2678  0.6085  1.0397  0.0686

s.e.   0.0485  0.0346  0.0554  0.0470

Domyślny algorytm poradził sobie jeszcze gorzej.

Ad 8. ARIMA(1, 1, 1) z phi 0.1, theta -0.95. Zarówno domyślny jak i kryterium BIC dało to samo:

Series: sim8 

ARIMA(1,1,2) 

Coefficients:

          ar1      ma1      ma2

      -0.9124  -0.0011  -0.9222

s.e.   0.0401   0.0304   0.0285

Obydwa kryteria błędnie wskazały zarówno liczbę rzędów, jak i wartości współczynników. Proces okazał się dla nich zbyt złożony. Zwiększmy więc liczbę danych z 1000 na 10000:

set.seed(1)

sim8 = arima.sim(list(order=c(1,1,1), ar=0.1, ma=-0.95), n=10000)

plot(sim8, type="l")

Co się okazuje, obydwa kryteria pokazują to samo i radzą sobie dość dobrze:

Series: sim8 

ARIMA(2,1,1) 

Coefficients:

         ar1     ar2      ma1

      0.1205  0.0157  -0.9582

s.e.  0.0105  0.0105   0.0032

Chociaż model błędnie dodaje phi(2), to jego p-value okazuje się nieistotne (powtarzamy kod na p-value):

(1-pnorm(abs(dopas$coef)/sqrt(diag(dopas$var.coef))))*2

      ar1       ar2       ma1 

0.0000000 0.1331158 0.0000000 

Wobec tego odrzucamy ar2, tak że automat znajduje prawidłowe rozwiązanie.

Podsumowując część optymalnego modelowania, dostajemy dowód, że w ARIMA należy posługiwać się BIC - przynajmniej stosując auto.arima tak jak wyżej. Oczywiście są jeszcze inne możliwe kryteria do sprawdzenia, ale pakiet oferuje jedynie AIC i BIC. 

Literatura:

[1] Hansen, B., Econometrics, 2022, link: https://www.ssc.wisc.edu/~bhansen/econometrics/

P. S. Bardzo dobre, zrozumiałe opracowanie.

Jaka jest różnica między testem t-Studenta i Chowa? Kiedy nie powinno się ich stosować i czym je zastąpić?

 Załóżmy dla uproszczenia, że stopy zwrotu z jakichś aktywów, np. funduszy akcji, mają rozkład normalny i nie są skorelowane w czasie. Dalej zastanawiamy się, w który fundusz zainwestować. Ponieważ rozkład jest gaussowski, są tylko dwa parametry do analizy: średnia i wariancja. Im większa średnia stopa zwrotu, tym lepiej oraz im mniejsza wariancja (i jej pierwiastek, czyli odchylenie std) tym lepiej. W sumie wystarczy użyć współczynnika Sharpe'a będącego w najprostszym ujęciu ilorazem średniej stopy zwrotu do odchylenia standardowego. Wydawałoby się więc, że wystarczy wyciągnąć stopy zwrotu z kilku ostatnich lat, np. z notowań mbanku lub stooq.pl, obliczyć wskaźnik Sharpe'a i wybrać fundusz o najlepszym wyniku. Jednak takie proste podejście jest niepoważne z naukowego punktu widzenia. Przykładowo, jak można traktować poważnie taką metodę, skoro przyjmujemy za pewnik, że średnia stopa zwrotu pozostaje stabilna w czasie? Niektórzy autorzy uzasadniają użycie stałych parametrów, odwołując się do ekonofizyki: M. Dobija i D. Dobija [1] stwierdzają, że premia za ryzyko jest wielkością wynikającą z praw natury, a przez to stabilną, wynoszącą średnio ok. 8% rocznie (w zależności od sposobu liczenia może wynieść nieco więcej). Chociaż analogie do fizyki są bardzo kuszące, to mimo wszystko nie muszą one odzwierciedlać w pełni rzeczywistości, bo nie wywodzą się z dedukcji ekonomicznej. Teoria CAPM wprost mówi, że im większe ryzyko rynkowe, tym większa premia za ryzyko. Owe 8% może istotnie dotyczyć pewnej średniej, czyli np. portfela rynkowego (bo wraz z premią za czas równą ok. 3% koszt kapitału wyniesie 11%), ale fundusze akcji budują własne portfele. W dynamice giełdowej fundusze zmieniają skład portfela, a wtedy też mogą zmieniać się ich zyski. 

Dochodzimy do wniosku, że potrzebne są nam narzędzia do oceny stabilności parametrów rozkładu stopy zwrotu. Z czego można więc skorzystać na szybko?

Test t

Jest to tak standardowe narzędzie, że znajdziemy je w każdej wersji Excela. Na dodatek testu t możemy użyć na dwa sposoby: poprzez formułę (T.TEST) lub poprzez specjalną opcję przeznaczoną do analizowania danych statystycznych (dodatek Analiza danych, w ang. Analysis Toolpak). Użyjemy formuły.

Jest kilka rodzajów testu t. Wg pomocy Excela wszystkie one stanowią test t-Studenta, choć trzeba wiedzieć, że oryginalny test t-Studenta dotyczy sytuacji, gdy wariancje w obu porównywanych próbach są równe. Gdy wariancje nie są równe, stosuje się tzw. test Welcha. W Excelu nazwa nie ma znaczenia. Generalnie, zasada jest taka, że dla zmiennych ekonomicznych lub giełdowych zawsze stosujemy test Welcha. Nawet jeśli wariancja jest taka sama, to test Welcha też będzie dobrze działać. Niektórzy autorzy uważają, że w ogóle zamiast testu Studenta należy stosować test Welcha, ponieważ jego moc jest zbliżona do tego pierwszego, gdy wariancje są jednorodne [2].

Czyli mamy dwie próby A i B. Nie jest konieczne, żeby miały taką samą liczebność, chociaż sami czujemy, że żeby prawidłowo porównać dwie próby, to jednak liczba danych powinna być w miarę podobna, jeśli oczywiście to możliwe.

Formuła jest banalnie prosta:

=T.TEST(A; B; 2; 3)

Nazwa może być myląca, bo  nie zwraca ona wartości statystyki t, tylko prawdopodobieństwo, że średnie są równe. Trzeci argument wynosi 2, bo chcemy sprawdzić zachowanie ogonów rozkładu z dwóch stron - dodatniej i ujemnej. Tak naprawdę ten argument mógłby w ogóle nie występować, bo ze względu na założenie symetryczności, zawsze ogon jednej strony będzie równy połowie dwóch ogonów w sensie prawdopodobieństwa - więc samemu można sobie już podzielić na dwa czy podwoić.

Istotny jest natomiast czwarty argument - typ testu. Będzie równy 1 dla testu sparowanego, tzn. gdy pozycja każdego elementu w jednej próbie ma znaczenie dla elementu o tej samej pozycji w drugiej próbie. To znaczy, każdy element w jednej próbie jest sparowany z elementem na takiej samej pozycji w drugiej próbie. Przykładowo mamy grupę pacjentów A. Każdy pacjent ma jakiś wynik bez stosowania leku. Następnie stosujemy na każdym pacjencie lek. Jest oczywiste, że kolejność pacjentów w grupie A ma znaczenie, jeśli chcemy porównać ich wyniki po użyciu leku. Dlatego, tworząc grupę B pacjentów, którzy przyjmą lek, musimy mieć pewność, że dany pacjent ma tę samą pozycję w obu grupach. Wynik w jednej grupie może się różnić w drugiej grupie, ale czy statystycznie wyniki będą różne? Zauważmy, że nawet jeśli dla kilku pacjentów lek będzie działał doskonale, ale dla reszty nie przyniesie żadnej poprawy, to badacz musi uznać, że lek nie działa. Widać więc, że ważne jest, aby grupa pacjentów była reprezentatywna w badaniach i samo grupowanie ma duże znaczenie. W ekonomii użycie sparowanego testu w niektórych sytuacjach ma sens, jak choćby porównanie zmian PKB okres do okresu. W przypadku stóp zwrotu na rynkach nie ma to już znaczenia.

Jeżeli dane nie muszą być sparowane, a wariancja w obydwu próbach jest równa, czyli zmienna jest  homoskedastyczna, to typ będzie równy 2.

W końcu, jeśli wariancje w obu próbach się różnią (zmienna jest heteroskedastyczna), to dajemy 3. To jest właśnie test Welcha. 

Jak wygenerować losowe obserwacje pochodzące z rozkładu normalnego?

Aby sprawdzić jak test działa, możemy wygenerować dwie próby o rozkładzie normalnym. W polskiej wersji używamy formuły ROZKŁAD.NORMALNY.ODW (w ang. wersji NORMINV) - tj. rozkład normalny odwrócony, czyli po prostu nasza zmienna o rozkładzie normalnym (formuła ROZKŁAD.NORMALNY, w ang. NORMDIST, zwraca prawdopodobieństwo tej zmiennej). Formuła zawiera 3 argumenty: prawdopodobieństwo zwracanej wartości, średnia oraz odchylenie std. Powiedzmy, że chcemy wygenerować 50 danych. Dwa ostatnie są oczywiste - założymy średnią równą 0 i odch. std 1, Gorzej z pierwszym argumentem. Przede wszystkim należy wiedzieć, że zwracana wartość pochodzi z rozkładu skumulowanego (dystrybuanty), a nie gęstości. To z kolei umożliwia zastosowanie pewnej sztuczki: za prawdopodobieństwo wstawiamy formułę LOS(), w ang. RAND(). Generuje ona wartość losową z przedziału od 0 do 1, ale - uwaga - z rozkładu jednostajnego. Czyli nasza formuła ma postać:

=ROZKŁAD.NORMALNY.ODW(LOS(); 0; 1)

Formułę tę kopiujemy 50 razy i za każdym razem uzyskamy nowy wynik, który będzie pochodził już z rozkładu normalnego.

-----------------------------------------------------------------------------------------

Dla zainteresowanych: jak działa ta formuła? Intuicyjnie wydawałoby się, że powinna ona przynieść zmienną także z rozkładu jednostajnego: przecież dystrybuanta to jest funkcja rosnąca, a więc każdej wylosowanej obserwacji z przedziału od 0 do 1 przyporządkowujemy zawsze jedną liczbę. Skoro zmienna niezależna pochodzi z rozkładu jednostajnego, to jak to się dzieje, że przekształca się w rozkład normalny? Problem leży w istnieniu dodatkowych warunków, jak średnia i wariancja, a szczególnie chodzi tu o wariancję. W art. Od zasady maksymalnej entropii do regresji logistycznej pokazałem, że maksymalizacja entropii Shannona bez dodatkowych warunków prowadzi do rozkładu jednostajnego. Po dodaniu warunku średniej, otrzymamy rozkład logistyczny. Gdyby jeszcze tam dodać warunek wariancji, dostalibyśmy rozkład normalny.

To jednak wyjaśnienie jest bardziej pośrednie, bo nie dotyczy samej odwrotnej dystrybuanty (kumulanty). Gdyby bezpośrednio to wyjaśnić, to trzeba odwołać się do jej własności. Jedną z nich jest fakt taki, że jeśli zmienna p ma rozkład jednostajny od 0 do 1, to zmienna 1/F(p) posiada dystrybuantę F (zob. art. w Wikipedii). Własność tę wykorzystują różne generatory liczb losowych. W szczególności dla rozkładu normalnego będzie to transformacja Boxa-Mullera.

-----------------------------------------------------------------------------------------

Generujemy więc dwie losowe próby, A i B, po 50 danych za pomocą przedstawionej formuły: 

=T.TEST(A; B; 2; 2) lub 

=T.TEST(A; B; 2; 3) 

Nie tu ma większego znaczenia którą wybierzemy. Klawiszem F9 będziemy zmieniać wylosowane wartości. Zauważymy, że niemal zawsze wartość testu da wartość powyżej 0,05. To znaczy, że hipotezy zerowej mówiącej, że średnie są równe, nie powinniśmy odrzucać. Czyli test działa prawidłowo.

W kolejnym kroku próbę B nieco zmieniamy: podwyższamy średnią do 0,2. Czyli tworzymy formułę =ROZKŁAD.NORMALNY.ODW(LOS(); 0,2; 1). 

Patrzymy na test t. Co się okaże? Najczęściej dostaniemy wartość testu powyżej 0,05, a to znaczy, że test źle działa. Wiemy przecież, że pierwsza próba ma ciągle średnią 0, a druga 0,2, więc powinniśmy odrzucić hipotezę zerową. Różnica 0,2 okazuje się po prostu za mała. Zwiększmy średnią B do 0,5. Tym razem test radzi sobie lepiej, bo najczęściej daje wartości poniżej 0,05, choć nadal ma czasem problemy. Zwiększmy więc średnią do 0,7. Tym razem niemal zawsze dostaniemy wartości testu poniżej 0,05, czyli odrzucamy hipotezę zerową, że średnie są równe. Dostaliśmy wynik prawidłowy. 

Czyli już widzimy, że przy 50 danych różnica między średnimi nie może być minimalna, aby test dawał wiarygodne wyniki. Ale sytuacja jeszcze się pogorszy, gdy zwiększymy odchylenie std w B. Powiedzmy, że zwiększamy je dwukrotnie, tak że dostaniemy formułę =ROZKŁAD.NORMALNY.ODW(LOS(); 0,7; 2). Wtedy już stosujemy tylko test Welcha: 

=T.TEST(A; B; 2; 3). 

Okazuje się, że test działa słabo dla 50 danych. Dopiero po zwiększeniu liczebności obu prób do 100, statystyka testu wygląda lepiej: błędów dużo mniej.

Chociaż Excel nadaje się do robienia podstawowych i niezbyt złożonych operacji statystycznych, to jeżeli chcemy zrobić coś więcej, musimy ściągnąć jakieś dodatki albo samemu stworzyć formułę lub makro. Mimo wszystko najlepiej chyba posłużyć się specjalistycznym oprogramowaniem. 

Test Chowa

Kiedyś poruszyłem temat stabilności modelu ekonometrycznego. Jednym z nich był test Chowa, który nieco przypomina test t-Studenta z tego powodu, że także porównuje dwie próby. Różnica polega na tym, że test Chowa bada stabilność wszystkich parametrów w regresji liniowej zbudowanej na podstawie połączonych obu prób. Można by pomyśleć, że jeśli w regresji ustalimy zerowe nachylenie, to powinniśmy dostać te same rezultaty dla Chowa co dla t. Jednak tak nie jest, bo test t porównuje różnicę między jedną średnią a drugą, natomiast test Chowa (inaczej test F) stosunek średniej z sumy połączonych próbek do średniej pierwszej próby. Okazuje się jednak, że w jednym przypadku oba testy dają praktycznie te same rezultaty w sensie prawdopodobieństwa błędu (p), a dzieje się to w sytuacji najczęściej stosowanym, czyli gdy obie próbki mają taką samą liczebność (tj. w przypadku Chowa potencjalna zmiana średniej następuje w połowie badanego zakresu).

Tym razem użyję języka R do testu Chowa. W zasadzie to od R w ogóle mógłbym zacząć i nie bawić się z Excelem, ale z drugiej strony Excela zna każdy, więc łatwiej w ten sposób znaleźć punkt zaczepienia. Chociaż nie jest to konieczne, to polecane jest stosowanie R-Studio. I jeszcze nota dot. Gretla. Możemy także go użyć do analiz, chociaż R jest bardziej wszechstronny, dlatego od teraz na nim się skupię. Gretl ma w tym temacie tylko jedną przewagę: analiza zmian strukturalnych modelu jest już w programie zainstalowana, a w R musimy sami zainstalować pakiet. Dodatkową niedogodnością w R jest wymóg jeszcze załadowania pakietu.

W konsoli albo skrypcie wpisujemy komendę install.packages(strucchange) lub  install.packages("strucchange") i zatwierdzamy ENTER. W R-Studio można też wybrać Tools -> Install packages -> wpisujemy nazwę strucchange. 

Sporą niedogodnością, być może w przyszłości rozwiązaną, jest to, że aby użyć danego pakietu, musimy użyć aż trzech komend. Pierwsze to właśnie instalacja, drugie to załadowanie pakietu przez komendę library(strucchange) lub library("strucchange"). Tu jeszcze mała uwaga. Ponieważ będę na blogu używał potem różnych pakietów do R, warto od razu mieć kod sprawdzający czy dany pakiet jest zainstalowany i jeśli nie, to należy go zainstalować. Następnie włączyć. W sumie początkowy kod byłby następujący:

#jeżeli to wymagane, instalacja pakietu (tu strucchange) i włączenie go

if (require(strucchange)==FALSE) {

  install.packages(strucchange)

  library(strucchange)

}

Trzecią komendą będzie użycie samego pakietu. Ponieważ składa się on z wielu funkcji, nie przywołujemy już samego strucchange, który jest klasą funkcji, ale funkcję, np. sctest(x, ...) , gdzie x to jakiś obiekt. Na tej stronie  znajdziemy odnośniki do dokumentacji w formacie pdf.  Samemu ciężko się przebić przez to wszystko, więc trzeba po kolei zacząć od prostych komend.

Podobnie jak w Excelu musimy najpierw stworzyć dwie próby do testowania. Każda ma mieć 50 danych o rozkładzie normalnym. Służy do tego funkcja rnorm(). 

(1) Średnia 0 i odch std 1. W skrypcie wpisujemy kod:

n1 = 50

n2 = 50

sim1 = rnorm(n1, 0, 1)

sim2 = rnorm(n2, 0, 1)

gdzie n1 i n2 to liczebność próbek, sim1 i sim2 to nazwane zmienne symulujące rozkład Gaussa. 

Musimy teraz połączyć obie próbki. Najprostsze to użycie funkcji c:

sim12 = c(sim1,sim2)

gdzie sim12 to połączenie próbki 1 i 2 w jedną.

W końcu wykonujemy test Chowa. To właśnie teraz przywołujemy funkcję sctest:

sctest(sim12 ~ 1, type="Chow", point=n1)

Falka oznacza tworzenie regresji. Chociaż jest zrozumiałe, dlaczego jest ona użyta (notacja matematyczna), to jej wpisanie na klawiaturze nie jest zbyt wygodne (shift+tylda+spacja). W każdym razie sim12 to zmienna zależna, a sama cyfra 1 po tyldzie oznacza po prostu, że tworzymy regresję liniową bez nachylenia, jedynie ze stałą, której parametr chcemy wyznaczyć. Formuła sctest umożliwia analizę pod różnym kątem i jednym z nim z nich jest test Chowa, stąd type = "Chow". Ostatni element wprost nawiązuje do tego testu, bo wymaga wprowadzenia punktu, w którym sprawdzamy zmianę struktury modelu. Chcemy sprawdzić czy dzieje się to w połowie próby, więc będzie to n1 lub n1 + 1.

Wykonujemy kod. W R-Studio za pomocą CTRL+ENTER albo klikamy Run w oknie skryptu.

Hipoteza zerowa to brak zmiany struktury. Wobec tego gdy p > 0,05 lub jeśli mamy niewielką próbę, to gdy p > 0,1, uznajemy, że brak podstaw do odrzucenia hipotezy. Kod wykonujemy wielokrotnie. Test działa prawidłowo - praktycznie za każdym razem dostaniemy p > 0,05. 

(2) Średnia powyżej 0 i odch std 1.

Zmieńmy więc średnią drugiej próby, sim2, na 0,2. Tzn. sim2 = rnorm(n1, 0.2, 1)

Tym razem p najczęściej znajduje się powyżej 0,1, co wskazywałoby błędnie na niezmienność średniej. Przy 50 danych zmiana średniej o 0,2 jest po prostu za niska. Podnieśmy ją do 0,5: sim2 = rnorm(n1, 0.5, 1).

Po wielokrotnym użyciu testu wnioskujemy, że jest lepiej, ale nadal dość słabo. Dopiero podniesienie średniej sim2 do 0,7 będzie dawało bardzo dobre rezultaty, najczęściej p < 0,05. Co ważne, nie zawsze graficznie widoczna ta różnica. 3 przykłady poniżej (sim2 = 0,7), za pomocą komendy plot(sim12, type="l")

 

Porównanie testu Chowa z testem t

Wcześniej napisałem, że oba testy dają identyczne rezultaty jeśli obie próby są równoliczne. Możemy to teraz sami sprawdzić. Dotychczas test t robiłem w Excelu, ale skoro przechodzimy na R, trzeba to samo zrobić w tym języku. W R test t wykonujemy za pomocą komendy t.test. W naszym przykładzie będzie to 

t.test(sim1,sim2). 

I tak np. dostaniemy dla niego

t = -2.3179, df = 97.889, p-value = 0.02254

a dla testu Chowa 

F = 5.3726, p-value = 0.02254

Z czego ta równość wynika? Podnieśmy statystykę t, czyli -2,3179, do kwadratu. Otrzymamy 5,3726, czyli F. 

Jednak ta równowaga obu testów powstaje tylko w przypadku, gdy punkt zmiany struktury wynosi n1, czyli od połowy całej próby. Powiedzmy, że chcemy sprawdzić zmianę średniej nie od połowy zakresu, ale nieco wcześniej, np. od 40-tej obserwacji. Z punktu widzenia testu Welcha próbkę 1 zmniejszymy o 10, a do próbki 2 dokładamy te odjęte wartości. Wykonajmy teraz kod:

n1 = 40

n2 = 60

sim1 = rnorm(n1, 0, 1)

sim2 = rnorm(n2, 0.7, 1)

sim12 = c(sim1,sim2)

t.test(sim1,sim2)

sctest(sim12 ~ 1, type="Chow", point=n1)

Dostałem dla testu Welcha:

t = -2.0099, df = 79.143, p-value = 0.04785

i Chowa:

F = 4.1691, p-value = 0.04385

P-value już nie są takie same (t^2 będzie różne od F). Kolejne przesunięcie punktu zmiany spowoduje jeszcze większą różnicę.

Problem okresu zmiany parametru

W rzeczywistości rzadko wiemy z góry, kiedy moment zmiany może nastąpić. Oryginalny test Chowa [3] zakładał sytuację, gdy badacz ma przed sobą taką oto sytuację. Najpierw jakaś zmienna wykazywała trend przez wiele lat, więc można było znaleźć parametry jej regresji. Mijają kolejne lata i badacz ma zadanie sprawdzić czy trend się zmienił od poprzedniego badania. Czyli do starej regresji dodawane są kolejne obserwacje i można teraz sprawdzić czy nowa regresja (stara powiększona o nowe obserwacje) zmieniła parametry. Widać więc, że jest to odwrotna analiza do naszej. W naszej wybieramy punkt załamania, a Chow nie wybierał żadnego punktu - był on dany z góry. 

Mimo tej różnicy, przeprowadzana przez nas analiza jest prawidłowa. Logicznie biorąc nie ma znaczenia, czy dostajemy nową próbkę czy też dzielimy całą próbę na dwie części. W tym drugim podejściu powstaje jednak problem, polegający na arbitralnym wyborze okresu zmiany struktury modelu.

W powyższym przykładzie normalnie założylibyśmy połowę zakresu dla testu Chowa. Czyli sctest miałby postać:

sctest(sim12 ~ 1, type="Chow", point=(n1+n2)/2)

Przykładowo otrzymałem  dla Welcha p = 0,008 i Chowa p = 0,021. Innym razem odpowiednio 0,02 i 0,03. Test Welcha zawsze wskaże mniejsze p, czyli da lepsze wyniki, ale generalnie test Chowa nadal sobie dobrze radzi mimo iż błędnie wskazaliśmy moment zmiany. Dzieje się tak, bo odchylenie 10 okresów nie powoduje dużego zamieszania: nadal prrzed połową średnia będzie mniejsza niż po niej.

Niestety jeśli zmniejszymy średnią sim2 do 0,5 sytuacja się pogorszy na niekorzyść Chowa. Np. otrzymamy dla t p-value = 0,014 i dla Chowa 0,092. Problem ten można "rozwiązać" zwiększając liczebność próby. Np. wykonując kod:

n1 = 100

n2 = 200

sim1 = rnorm(n1, 0, 1)

sim2 = rnorm(n2, 0.5, 1)

sim12 = c(sim1,sim2)

t.test(sim1,sim2)

sctest(sim12 ~ 1, type="Chow", point=(n1+n2)/2)

test Chowa będzie już dawał bardzo dobre wyniki mimo odchylenia 50 okresów.

Problem małej próby

Oczywiście nie jest to żadne prawdziwe rozwiązanie, bo w przypadku szeregów czasowych to właśnie liczebność stanowi zawsze ograniczenie. Chociaż zazwyczaj w przypadku zbyt małej liczby danych rocznych czy kwartalnych, zwiększa się częstość do miesięcznej albo tygodniowej, to trzeba pamiętać, że to spowoduje, że także średnia spadnie i uchwycenie zmiany może być trudniejsze. Dodatkową przeszkodą jest zwiększona zmienność dla większych częstości.

Jeżeli fundusz funkcjonuje od 10 lat, to mamy zaledwie 40 kwartałów. Np. jeżeli w próbie 50-elementowej od 35 elementu nastąpi zmiana średniej z 0 na 0,7, to test Chowa:

n1 = 35

n2 = 15

sim1 = rnorm(n1, 0, 1)

sim2 = rnorm(n2, 0.7, 1)

sim12 = c(sim1,sim2)

sctest(sim12 ~ 1, type="Chow", point=(n1+n2)/2)

będzie bardzo chwiejny. Istotność należy już ustawić na 0,1, a i tak co druga próba będzie dawać p > 0,1. 

Najgorsze jest to, że jeżeli zmiana średniej nastąpi tylko 5 obserwacji więcej niż w połowie zakresu (czyli n1 = 30 i n2 = 20), to dodatkowo test będzie często zawodził. Nie można wtedy polegać na nim.

Test QLR

W samej idei test QLR (Quandt Likelihood Ratio) rozwiązuje powyższy kłopot bardzo prosto - wykonuje test Chowa po wszystkich punktach zakresu (poza skrajnymi punktami początku i końca sumarycznej próby), a następnie znajduje maksimum otrzymanych statystyk F, które powinno przekroczyć poziom uznany za przypadkowy. Pozornie wydaje się, że można by napisać kod z pętlą, która to robi. To by z kolei znaczyło, że wystarczy wyznaczyć max F i sprawdzić czy p-value > poziom istotności. Tylko że sprawdzając każdy punkt potencjalnej zmiany struktury, za każdym razem wyznaczamy jakieś p-value, a więc dla wielu punktów p staje się zmienną o pewnym rozkładzie częstości. Jeżeli wybieramy sobie z próbki najmniejszą wartość tego rozkładu (któremu odpowiada największe F), to przy dużej próbie zawsze trafimy na p mniejsze od poziomu istotności. Przy mniejszej próbie też będzie na to duża szansa.  Będzie to przypadkowe odchylenie zupełnie jak w grze w totolotka, gdy przy wielu losowaniach w końcu trafimy na liczbę wygrywającą. W naszym przypadku każde losowanie to kolejna obserwacja w szeregu czasowym. Testy statystyczne polegają na oszacowaniu prawdopodobieństwa, że takie odchylenie to losowy przypadek, a nie realna cecha rozkładu zmiennej. Widać, że potrzebna jest metoda oceniająca czy przebicie progu istotności nie jest efektem wyłącznie samej budowy rozkładu prawdopodobieństwa, w tym wypadku rozkładu F. Quandt [4] posłużył się ilorazem wiarygodności (likelihood ratio). 

W pakiecie strucchange do wykonania testu QLR służy funkcja Fstats. W sumie będą 4 nowe linijki kodu.  Pierwsza to funkcja Fstats, którą dla skojarzenia nazwiemy qlr. Druga  nie jest potrzebna, ale przyda się grafika - plot. Jest to funkcja wewnętrznie związana ze strucchange, dlatego zawiera także argument alpha. Ustawimy go na 0,1. Trzecia jest najważniejsza, to sctest, ale tym razem głównym argumentem nie będzie model liniowy, tylko sam qlr. Drugi tutaj argument "supF" oznacza supremum statystyki F, czyli właśnie jego max. Czwarta linia - breakpoints - to sam moment zmiany, w naszym przypadku test powinien wskazać 35. 

n1 = 35

n2 = 15

sim1 = rnorm(n1, 0, 1)

sim2 = rnorm(n2, 0.7, 1)

sim12 = c(sim1,sim2)

qlr = Fstats(sim12 ~ 1)

plot(qlr, alpha=0.1)

sctest(qlr, type="supF")

breakpoints(qlr)

Przykładowo dostaniemy takie wyniki:

Fstats prawidłowo pokazał, że zmiana struktury nastąpiła w 35 punkcie, ponieważ jest to najwyższy punkt statystyki F. Dodatkowo podano, że odpowiada to 68% całkowitego zakresu.

Dostaliśmy więc bardzo ciekawe narzędzie do oceny stabilności średniej, które można wykorzystać do analizy portfela, zarówno akcji, jak i funduszy. Zróbmy przykład dla obu.

----------------------------------------------------------------------------------------

Aspekty techniczne kodowania w R:

Standardowym katalogiem roboczym R jest folder Pobrane / Downloaded , więc tam możemy ściągać dane np. ze stooq.pl. Zaletą portalu jest możliwość ściągnięcia miesięcznych, kwartalnych i rocznych danych, ale wadą jest to, że można nie zauważyć, że ostatni element nie należy do grupy o oczekiwanej częstości, tu kwartalnej, ale jest z ostatniego dnia. Dlatego trzeba pamiętać, że ostatni element usuwamy, jeżeli jego data wyprzedza datę bieżącą. Zrobimy to też za pomocą kodu.

Po ściągnięciu pliku csv, od razu tworzymy nową zmienną o nazwie plik, którą definiujemy jako tabelę wczytaną za pomocą komendy read.csv. Istotną opcją jest wybór separatora, bo czasami będzie to średnik, czasami przecinek. Pliki stooq mają zawsze więcej niż 1 kolumnę, bo 1 to zawsze data. Stąd problem separatora rozwiążemy za pomocą funkcji if, która przypisze przecinek, jeśli liczba kolumn zmiennej plik nadal będzie równa 1. 

Powszechnym problemem w kodowaniu jest format daty. Jedną z najważniejszych rzeczy jest przyswojenie sobie, że to co wygląda tak samo, takie być nie musi. Data może być tekstem albo liczbą. Przyzwyczailiśmy się, że datę traktujemy zawsze jak liczbę nie do końca zdając sobie z tego sprawy. Dlatego zawsze musimy zamienić datę-tekst na datę-liczbę. W naszym przypadku jest to proste, bo wystarczy użyć funkcji as.Date(). 

Pierwsza zmienna to daty, druga - cena. Obie są danymi zewnętrznymi. Stopa to trzecia zmienna- stopa zwrotu. Zmienną tę możemy utworzyć na dwa sposoby. Pierwszy to klasyczna pętla, jednak jest ona mniej wydajna (przy dużej liczbie danych wolniejsza) niż druga metoda zwana wektoryzacją. Ta druga metoda jest także stosowana w Excelu poprzez formuły tablicowe. 

----------------------------------------------------------------------------------------

CD Projekt

Dla akcji przykładem niech będzie CDR od 2017 r. - miesięcznie (63 obserwacji). Cena i stopa zwrotu kształtowały się następująco:

Wykonujemy kod:

#jeżeli to wymagane, instalacja pakietu (tu strucchange) i włączenie go

if (require(strucchange)==FALSE) {

  install.packages(strucchange)

  library(strucchange)

}

#czytamy plik csv

plik = read.csv('cdr_m.csv', sep=";")

if (ncol(plik)==1) {

  plik = read.csv('cdr_m.csv', sep=",")

}

#jeżeli ostatnia data przekracza wczorajszą datę, to usuwamy ostatnią obserwację

daty = plik[,1]

daty = as.Date(daty)

if (daty[length(daty)]>Sys.Date()) {

  plik = plik[-nrow(plik),]

  daty = plik[,1]

  daty = as.Date(daty)

}

cena = plik[,2]

#1 sposób obliczenia stopy zwrotu - klasyczna pętla, niewydajny

for (t in 1:length(cena)) {

  if (t>1) {

    stopa[t-1]=cena[t]/cena[t-1] - 1

  }

}

#2 sposób - wektoryzacja

stopa = diff(cena)/cena[-length(cena)]

#dalej już qlr

qlr = Fstats(stopa ~ 1)

sctest(qlr, type="supF")

breakpoints(qlr)

#tworzymy wykresy

par(mfrow=c(2,1), mar=c(2,5,2,2))

plot(qlr, alpha=0.1)

plot(x=daty, y=cena, type="l")

Dostajemy statystyki:

 wykres:

Jedna uwaga: ponieważ mamy tu tylko 63 obserwacje, wybrałem istotność 10%. 

Wygląda na to, że średnia stopa zwrotu CDR zmieniła się od połowy 2020 r. Koresponduje to z widocznym załamaniem trendu rosnącego. Uzyskano p-value = 0,053, czyli hipoteza mówiąca, że nastąpiła zmiana średniej, jest wiarygodna.

Allianz FIO Akcji Małych i Średnich Spółek

Kolejnym przykładem miał być fundusz, więc wybrałem Allianz Akcji MiŚ spółek. Dane także pobrałem ze stooq.pl. Symbol na portalu to 1394_N. Aby zachować porównywalność z poprzednim testem, tak samo przyjąłem początek 2017 i 63 danych miesięcznych. 

Kod:

#jeżeli to wymagane, instalacja pakietu (tu strucchange) i włączenie go

if (require(strucchange)==FALSE) {

  install.packages(strucchange)

  library(strucchange)

}

#czytamy plik csv

plik = read.csv('1394_n_m.csv', sep=";")

if (ncol(plik)==1) {

  plik = read.csv('1394_n_m.csv', sep=",")

}

#jeżeli ostatnia data przekracza wczorajszą datę, to usuwamy ostatnią obserwację

daty = plik[,1]

daty = as.Date(daty)

if (daty[length(daty)]>Sys.Date()) {

  plik = plik[-nrow(plik),]

  daty = plik[,1]

  daty = as.Date(daty)

}

cena = plik[,2]

#stopa zwrotu (wektoryzacja)

stopa = diff(cena)/cena[-length(cena)]

#dalej już qlr

qlr = Fstats(stopa ~ 1)

sctest(qlr, type="supF")

breakpoints(qlr)

#tworzymy wykresy

par(mfrow=c(2,1), mar=c(2,5,2,2))

plot(qlr, alpha=0.1)

plot(x=daty, y=cena, type="l")

Wyniki:

I wykresy:

W pierwszej połowie 2020 r. notowania funduszu jakby zmieniły kierunek i wielu intuicyjnie uznałoby, że to zmiana trendu ze spadkowego na zwyżkujący. Test QLR rzeczywiście pokazał, że coś się zadziało na stopach zwrotu w tym okresie, bo statystyka F osiągnęła wysoki szczyt. Niemniej obiektywnie oceniając, nie ma powodu, by wierzyć, że to realna zmiana trendu: p-value wynosi 0,17, a więc sporo powyżej 0,1. Uznajemy, że średnia stopa zwrotu funduszu nie zmieniła się.

Teoretycznie są dwie konsekwencje testów na stabilność średniej stopy zwrotu. Po pierwsze, jeżeli otrzymamy przesłankę, że nastąpiła zmiana tej średniej, to dostajemy argument, by inwestować lub nie inwestować w dane aktywo. Po drugie, jak dla mnie ważniejsze i bardziej wiarygodne, otrzymujemy informację, czy przy budowie portfela można porównywać średnie stopy zwrotu z różnych aktywów. Np. wskaźnik Sharpe'a może przekonywać, że dane akcje są lepsze, nie biorąc pod uwagę, że mogła nastąpić zmiana struktury i od jakiegoś czasu są tak naprawdę gorsze.

To co tu pokazałem, to oczywiście tylko niewielki wycinek problemu stabilności parametrów stóp zwrotu. Nie poruszyłem sprawy wpływu potencjalnej zmiany wariancji na skuteczność testu Chowa i QLR. Dalej, co jeśli występuje choćby lekka autokorelacja stóp zwrotu lub ich warunkowej wariancji? W końcu, jak na wyniki wpłynie fakt, że stopy zwrotu na giełdzie nie mają rozkładu normalnego? Przecież od początku zakładaliśmy tu klasyczną regresję liniową. Na te wszystkie wątpliwości należy jeszcze odpowiedzieć.

Literatura:

[1] Dobija, M., Dobija, D.,  Premia za ryzyko jako stała ekonomiczna, Ekonomika Organizacja Przedsiębiorstwa, 6/2006.

[2] Rasch, D., Kubinger, K.D., Moder, K.,  The two-sample t-test: pre-testing its assumptions does not pay off, Statistical Papers volume 52, p. 219–231 (2011).

[3] Chow, G., Tests of Equality Between Sets of Coefficients in Two Linear Regressions, Econometrica, Vol. 28, No. 3 (Jul., 1960), pp. 591-605.

[4] Quandt, R. E., Tests of the Hypothesis that a Linear Regression System Obeys Two Separate Regimes, Journal of the American Statistical Association, Vol. 55, No. 290 (Jun., 1960), pp. 324-330.

Walka z fake-newsami? Może niech zaczną od siebie?

Od ok. 2 lat trwa na szeroką skalę walka z tzw. fake-newsami. Niestety wg mnie jest to ryzykowne, bo jej ofiarą mogą paść sami fact-checkerzy. Problem polega na asymetrii między fact-checkerem a dezinformującym medium: ten pierwszy musi z samej swojej natury wykazywać się najwyższą starannością i rzetelnością. Ten drugi rzecz jasna nie musi. Jednak staranność i rzetelność wymaga wysiłku i czasu na analizę, co ogranicza jego możliwości. Ponieważ te drugie nie dbają o prawdę, ale szybkość, to będą powielać znacznie więcej fake-newsów niż ci pierwsi będą w stanie je wyłapać. Może to prowadzić do sytuacji, że ci pierwsi tracąc czujność, chcąc coś obalić, wpadają we własne sidła - same podają fałszywe informacje.

I tak jest Dziennik Gazeta Prawna (DGP) - powszechnie szanowana i uważana za medium wiarygodne, z ogromną ilością porad prawnych i opinii. I właśnie trafiłem na krótki artykuł, który kompromituje dziennik. Mianowicie, chodzi o Czy NBP „drukuje puste pieniądze”? z 16 grudnia 2021.

Już pierwsze zdanie to fałsz:

NBP nie „drukuje pustych pieniędzy” jak twierdzą autorzy fake newsów publikowanych w mediach społecznościowych. 

Zgodnie z art. 227 ust. 1 Konstytucji centralnym bankiem państwa jest Narodowy Bank Polski, któremu przysługuje wyłączne prawo emisji pieniądza oraz ustalania i realizowania polityki pieniężnej, odpowiada on także za wartość polskiego pieniądza. Dalej, zgodnie z art. 33 ust. 1 ustawy z dnia 29 sierpnia 1997 r. o Narodowym Banku Polskim wprowadzanie do obiegu banknotów lub monet odbywa się na podstawie zarządzenia Prezesa Narodowego Banku Polskiego ogłoszonego w Monitorze Polskim, które ustala termin wprowadzenia do obiegu i wielkość emisji oraz wzór i wartość nominalną banknotów lub monet, a w przypadku monet także stop, próbę i masę.

Drukowanie pieniędzy czyli ich emisja, zależy od polityki pieniężnej, a więc czy będzie to pusty pieniądz. Już z Konstytucji wynika więc, że NBP nie tylko drukuje pieniądze, ale że drukuje puste pieniądze - jeśli chce zwiększyć inflację. Z ustawy o NBP wynika, że ustala on wielkość nowej podaży pieniądza.

W dalszej części artykułu autor stara się przekonać czytelnika, że to banki komercyjne "drukują" pieniądze:

We współczesnych gospodarkach rynkowych pieniądz kreowany jest przez banki komercyjne poprzez udzielanie kredytów. Banki, udzielając kredytu, tworzą nowy pieniądz, którego nie było przed udzieleniem tego kredytu.

Mamy rozdwojenie jaźni? O co tu chodzi? Autor po prostu nie rozumie w ogóle tematu, który porusza. Drukowanie i emisja pieniądza to jedno. Kreacja pieniądza to drugie. To są są dwie różne sprawy.

W uproszczeniu można powiedzieć, że są dwa rodzaje "drukowania pieniędzy":

(1) emisja pieniądza, która nie musi oznaczać materialnej gotówki w kasach, ale w formie cyfrowej na rachunkach bankowych;

(2) kreowanie pieniądza poprzez pożyczanie pieniędzy innym podmiotom, wykorzystując maszynerię rachunkowości w postaci aktywów i pasywów: jeśli pożyczę od banku pieniądze na inwestycję, to bank ma należności, a ja mam zobowiązania. Jeżeli jednak ja te pieniądze przeznaczę na zakup obligacji firmy, to zyskuję teraz należności w firmie, a firma ma zobowiązania u mnie. Trick polega na tym, że choć matematycznie zobowiązania redukują się z należnościami, to globalnie księgowo następuje podwojenie aktywów - są to należności banku (początkowa pożyczka) i moje należności (moje obligacje). 

Jak widać pkt 1 to prawdziwa kreacja pieniądza, a druga to po prostu kreatywna księgowość. Łatwo udowodnić, że ta druga to nie jest żadne drukowanie nowych pieniędzy: gdyby nagle bank wymagał ode mnie spłaty całego długu, to musiałbym ja wymagać spłaty całego długu od firmy. Firma musiałaby znowu od kogoś innego. Cała gospodarka by się zawaliła! Jak widać ta kreacja pieniądza - kreatywna księgowość - w całości opiera się na zaufaniu podmiotów gosp. do siebie. Można powiedzieć, że cała współczesna gospodarka rynkowa na tym się opiera. To z tego m.in. powodu nastąpił taki kryzys w 2007 r., kiedy się okazało się, że to zaufanie było za duże - bo właśnie całe to pompowanie opierało się na kreatywnej księgowości. Niektórzy zwalają winę kryzysu na wolny rynek, wskazując na jego niedoskonałości w postaci przesadnego zaufania. No nie - gdyby nie przepisy prawne, zrzucające całą odpowiedzialność na banki lub państwo, to nie byłoby tego zaufania. 

Wróćmy jeszcze na moment do tego artykułu DGP. Autor - niepodpisany zresztą - stwierdza:

Bank centralny nie zwiększa podaży pieniądza, kupując obligacje. Te zakupy są dokonywane głównie od banków komercyjnych, które – w zamian za sprzedane papiery wartościowe – otrzymują płatność na swoje rachunki w NBP.

Podkreśliłem, że otrzymują płatność. A skąd niby NBP ma mieć na nią pieniądze? Z kosmosu? Trzyma w sejfie prezesa? Widzimy, że autor pisze o rzeczach, o których nie ma pojęcia.

Oczywiście jest jeszcze możliwość, że to tekst kogoś pracującego w NBP, który celowo wprowadza czytelników w błąd. Biorąc jednak pod uwagę niekompetencję osób pracujących w tej instytucji od kilku lat, nie ma to znaczenia. 

P. S. Dla zainteresowanych: dokładniejsze wyjaśnienie związku między gospodarką a pieniądzem przedstawiłem tutaj.

GARCH zamiast wstęgi Bollingera

A więc prognoza ropy za pomocą AT nie sprawdziła się. Wygrała ekonomia - blokada podaży rosyjskiej ropy wywołała eksplozję cenową ropy na całym świecie. Można mówić o odwróconym czarnym łabędziu, którego też brałem pod uwagę dwa tygodnie temu. Gdy Rosja rozpoczęła inwazję na Ukrainę, cena oscylowała wokół 100$, teraz jest to 125, tygodniowo wstęga Bollingera mocno przekroczona, RSI(14) także

W krótkim okresie na statystyce nie można nigdy polegać, ale w długim już tak. Dlatego nie ma wątpliwości, że nastąpi powrót. Zazwyczaj pada teraz pytanie: kiedy? Ale to złe pytanie, bo wiadomo, że to nastąpi prędzej niż później. Ważniejsze, jak mocny będzie powrót? Bo zamiast oczekiwanego powrotu do średniej, możemy dostać prawie poziomy kierunek w kolejnych tygodniach. Gdyby tak było, to znaczyłoby, że wariancja stóp zwrotu spadła niemal do zera.

Zamiast Bollingera, którego zresztą nie jestem zwolennikiem, znacznie lepiej używać modelu klasy ARCH. Oprócz nieporównywalnej mocy statystycznej, pozwala on także prognozować odchylenia w kolejnym okresie.

Najprostszym sposobem, nie wymagającym żadnej wiedzy programistycznej, jest używanie Gretla. Program ten, jak na darmowy, całkiem nieźle się rozwija, obecna wersja pochodzi z lutego 2022. Przy instalacji możemy zaznaczyć język polski, choć musimy się nastawić przeskakiwanie z polskiego na angielski, bo program tutaj ma pewne braki.

Wstępne ustawienia: Tygodniowe notowania cen (na zamknięcia) ropy pochodzą z https://stooq.pl/q/d/?s=cl.f . Dane od 07.01.2018 do 6.03.2022. W gretlu otwieramy plik z danymi, np. .csv. Problem zazwyczaj dotyczy ustawienia dat, które trzeba ręcznie ustawić, bo program często nie potrafi czytać prawidłowo dat. Jeśli pyta o typ danych, wybieramy szeregi czasowe. Kiedy dane są gotowe, zaznaczamy zmienną - cenę ropy. Wybieramy Dodawanie zmiennych -> procentowe zmiany dla wybranych zmiennych -> wpisujemy nazwę dla stopy zwrotu, np. 'ropa_zwrot'. 

Model GARCH: Wybieramy Model -> Modele jednowymiarowych szeregów czasowych -> GARCH variants. Upewniamy się, że zmienna zależna to ropa_zwrot. Resztę ustawień na razie zostawiamy. Dajemy OK. Dostałem model

Wszystkie współczynniki istotne stat., czyli jesteśmy na dobrej drodze. Jednak to początek. W otrzymanym oknie klikamy na ikonę wykresu, wybieramy Density:

Jedna uwaga: dlaczego podpisano residuals (reszty), a nie returns (stopy zwrotu)? Zauważmy, że rozkład teoretyczny i empiryczny mają średnią zero. Różnica między resztami / błędami a stopami zwrotu jest taka, że te pierwsze mają średnią zawsze równą zero, a te drugie niekoniecznie. Ich wariancje są już sobie równe (dowód - ostatni wzór tutaj). Zazwyczaj tam gdzie są modele klasy ARCH, określenie residuals / errors, czytamy to jako zmienna (tu stopa zwrotu) pomniejszona o swoją średnią. Trzeba jednak dodać, że w tym modelu pojawiają się też standaryzowane reszty - wtedy dodatkowo taka reszta dzielona jest przez odchylenie standardowe, aby sprowadzić do porównywalności ze standardowym rozkładem normalnym. To jest właśnie ten przypadek (opis brzmi std residuals vs std normal). 

Rozkład empiryczny jest leptokurtyczny, czyli stopy zwrotu nie mają rozkładu normalnego. Dla pewności można zrobić test normalności. Zmienna -> Testy normalności:

................................................................

Test na normalność rozkładu ropa_zwrot:

 Test Doornika-Hansena = 121,681, z wartością p 3,77906e-027

 Test Shapiro-Wilka = 0,901681, z wartością p 9,5057e-011

 Test Lillieforsa = 0,0966947, z wartością p ~= 0

 Test Jarque'a-Bera = 277,187, z wartością p 6,45201e-061

................................................................

Wszędzie p < 0,01 , co daje pewność, że to nie jest rozkład Gaussa.

Trzeba zmienić w opcjach GARCH variants -> Rozkład -> zmieniamy na t (czyli t-Studenta).  

 

Parametry są istotne stat - efekt ARCH i GARCH ewidentnie występuje.

Sprawdzamy znów Density:

Wygląda to dużo lepiej - rozkład t pozwala imitować efekty leptokurtozy.

Ok, mamy model, to możemy teraz zobaczyć graficznie, jak zmienia się warunkowe odchylenie standardowe w porównaniu z resztami. Klikamy ikonę wykresu -> Time series:

To co tu dostaliśmy, to nie są wartości otoczone przedziałem ufności - ten zresztą ma zakres ok. 2 odchyleń standardowych. A tutaj mamy jedno odchylenie - warunkowe. Interesujące jest, że niemal zawsze dane mieszczą się w tym zakresie. Z 2-3 razy było duże wychylenie poza ten obszar (rok 2020), ale tak wielkiego jak teraz jeszcze nie było.

Powiedzmy, że chcemy porównać dodatkowo warunkowe odchylenie z niewarunkowym - na jednym wykresie. Tego już tak łatwo nie dostaniemy. Wg opisu wykresu mamy tu reszty - już niestandaryzowane - otoczone niestandaryzowanymi odchyleniami warunkowymi. Klikamy na oknie modelu "Zapisz bundle jako" -> wybieramy uhat (series), zapisujemy zmienną, a następnie to samo dla h (series). W głównym oknie gretl dostaliśmy dwie nowe zmienne w postaci niestandaryzowanych reszt (uhat) oraz warunkowej wariancji (h). Teraz klikamy Dodawanie zmiennych -> Definiowanie nowej zmiennej. Wpisujemy csd_plus = h^0.5. Następnie Dodawanie zmiennych -> Definiowanie nowej zmiennej -> wpisujemy csd_minus = -h^0.5. Otrzymaliśmy warunkowe odchylenia std na plus i minus. I w końcu Dodawanie zmiennych -> Definiowanie nowej zmiennej -> wpisujemy: 

series usd_plus = sd(uhat). 

Analogicznie Dodawanie zmiennych -> Definiowanie nowej zmiennej:

series usd_minus = -sd(uhat). 

Uzyskaliśmy niewarunkowe odchylenia std (na plus i minus).

Wybieramy Widok -> Wykres szeregów czasowych. Zaznaczamy 5 zmiennych: uhat, csd_plus, csd_minus, usd_plus i usd_minus. Dostaniemy wykres:

Co warto zauważyć, niewarunkowe odchylenie std jest stosunkowo wysoko umiejscowione - jest to efekt roku 2020, który zaburzył ogólne statystyki. Jest to dobry przykład dowodzący, że błędnym podejściem jest klasyczne opieranie się na zwykłym odchyleniu standardowym. 

Niemniej w ostatnim tygodniu tempo cen przekroczyło znacznie powyżej 3 odchyleń std.  Zdarzenie tego rodzaju w ciągu 2 lat jest mało prawdopodobne. To znaczy, że musi teraz nastąpić powrót do średniej.

W sumie powstają 2 wnioski:

(1) nie powinniśmy się sugerować statycznym odchyleniem std, tylko warunkowym, dynamicznym;

(2) nie tylko ceny, ale też wariancje warunkowe są nieprzewidywalne. Można pobawić się w prognozy, ale lepiej posłużyć się GARCH do oceny powrotu ceny do średniej. Czyli otrzymujemy ulepszoną wersję wstęgi Bollingera. 

Napad Rosji na Ukrainę nie powinien być szokiem dla inwestorów

Krótki komentarz do wpisu z 26 lutego bloggera Podtwórca:

Życie jest przewrotne. Do samego końca żyłem iluzją, że Rosja nie zacznie pełnoskalowej inwazji na Ukrainę. Że jakoś to będzie: uznają "niepodległość" Donbasu, obstawią armią i będą jak w Gruzji używać "republiki" do destabilizowania sąsiada. Jeszcze w środę wieczorem napisałem, że na giełdzie wchodzimy w fazę kapitulacji, więc teoretycznie powinienem spodziewać się czegoś grubszego. A jednak rano na wieść o agresji zamurowało mnie, jak większość świata

Rzeczywiście większość świata była zszokowana, jedna z polskich chyba dziennikarek stwierdziła, że nie wierzy w to, że ktoś przewidywał, że dojdzie do takiej agresji. Jest to zdumiewająca pewność siebie - dokładnie ta sama, która nie pozwalała na myśl, że Rosjanie zaatakują. Przecież dzień w dzień byliśmy bombardowani wręcz komunikatami z USA, że za chwilę może dojść do tej inwazji. Mało tego, jeszcze w grudniu 2021 prezydent Ukrainy ostrzegał, że Rosja zaatakuje w styczniu lub lutym jego kraj. Dwa dni przed atakiem USA wprost zapowiedziały, że w ciągu następnych kilku dni Rosja przejmie - wg swoich planów - Kijów  i wprowadzi siłą swoje rządy. Wszystko było przewidywane prze wiele tygodni! Po czymś takim RACJONALNY człowiek powinien chociaż przygotować się mentalnie, nastawić się, że może naprawdę dojść w końcu do zbrojnej napaści. Nie oznaczało to bycia pewnym - bo nikt oprócz Rosjan nie był. Ale twierdzenie, że nikt nie przewidywał, a nawet doznanie szoku jest po prostu zwykłym lekceważeniem rzeczywistości i ignorancją na niepokojące informacje.

Może właśnie dlatego, że nie byłem zszokowany, napisałem w dniu rozpoczęcia inwazji 2 wpisy czysto giełdowe: Ropa blisko szczytu? oraz ... a WIG blisko dołka? (niektórym mogły one nawet wydać się nie na miejscu w takim dniu, no trudno). 

Dlatego lepiej nie przyznawać się, że jest / było się zszokowanym czy nawet zaskoczonym. NIE MA TU ŻADNEGO ZASKOCZENIA. Wystarczy przeanalizować na sucho informacje pojawiające się przed wojną. 

Jedynie co dobre, to to, że dostaliśmy w sumie kolejny dowód, że ludzie nie są racjonalni, nie potrafią prawidłowo oceniać sytuacji - i ci bardziej racjonalni mogą to wykorzystywać, oczywiście także na rynkach. Powtarzam jednak, że nie chodzi o to by być nadmiernie pewnym siebie, bo nikt nie miał 100% pewności, że to się stanie. Chodzi tylko o branie tego pod uwagę, tak by uniknąć właśnie szoku, w którym obecnie pogrążony jest cały świat. I oczywiście nie można mylić szoku ze zmartwieniem, poczuciem grozy, niepokoju czy strachu, bo to zupełnie coś innego (to uczucie jest akurat racjonalne na pewnym poziomie, również je odczuwam).

Nie wiem czy tylko mnie, ale owa napaść, która nastąpiła pod koniec pandemii przypomina 1 wojnę św., która jednak miała miejsce przed grypą hiszpanką. To skojarzenie rzecz jasna budzi niepokój, bo zapowiada 3 wojnę św., ale z drugiej strony nie można przesadzać z analogiami. To co napawa optymizmem w tej kwestii to fakt, że świat - szczególnie Europa - się obudził. Właściwie polityka międzynarodowa zmienia się na naszych oczach o 180 stopni, jeśli chodzi o nastawienie do Rosji. Oprócz nakładanych na ten kraj sankcji przygotowania wojskowe na wypadek 3 wojny samo w sobie będą ochroną przed nią. 

Ogromny błąd Rosji polega na tym, że zaatakowała znienacka i bardzo agresywnie, w tym ludność cywilną. W 2014 r. działała o wiele bardziej sprytnie i niemal kawałkami salami anektowała kolejne terytoria. Metodą gotowanej żaby usypiała czujność Europy i przyzwyczajała do swoich zbrodni. Ten więc "szokujący" atak doprowadzi Rosję do katastrofy ekonomicznej. Nie wiadomo czy kraje sąsiedzkie, np. skonfliktowana z nią Czeczenia, nie będzie chciała tego wykorzystać i przejąć część ros. terytoriów. Nie wiadomo czy Kazachstan nie dokona jakichś rozliczeń z sąsiadem. Jeśli miałoby dojść do 3 wojny św., to myślę, że nie takiej jaką sobie wyobrażamy. Byłoby to zbyt oczywiste. 

Wątpię też, by Putin chciał ryzykować wojną nuklearną. Jest prawdopodobne, że będzie straszył nią, ale już teraz wie, że jakikolwiek realny taki ruch doprowadzi do całkowitej anihilacji jego kraju.

... a WIG blisko dołka?

Kontynuacja poprzedniego wpisu. Skoro już zauważyliśmy odwrócenie korelacji ropy z wigiem, to przyjrzyjmy się teraz tym samym wskaźnikom, ale od strony właśnie WIG. Okazuje się, że dostajemy jakby lustrzane odbicie ceny ropy(!) Fibo 61,8% dziś osiągnięty, lekko nagięty (ropa też lekko nagięła), dzienny Bollinger wyprzedany, dzienny RSI wpadł pod wyprzedaż:

Tygodniowe tak samo:

Różnica jedynie leży w miesięcznych wskaźnikach - tutaj jesteśmy na środku:

Niepełne wyprzedanie miesięcznego RSI nie oznacza jednak, że nie może nastąpić odwrócenie trendu.* Przykładowo, sytuacja przypomina tę z roku 2011:

Jak pamiętamy, ropa nawet miesięcznie jest wykupiona. Czyli wydawałoby się, że nie jest to żadne lustrzane odbicie. Myślenie to się zmieni, gdy uświadomimy sobie, że przecież ogólna korelacja powinna pozostać dodatnia, a miesięczne częstości dotyczą średnio i długoterminowych zależności. Teraz wszystko się układa w całość: wkrótce nastąpi powrót dodatniej korelacji wig-ropa. Można zarysować następujący scenariusz: średnio lub długoterminowo WIG będzie dalej spadał - już razem z ropą (a także zagranicznymi indeksami), natomiast krótkoterminowo lub średnioterminowo odbije się. Być może będzie to trzecia fala:

* W ogóle istnieją pewne obiektywne przesłanki, dla których miesięczny RSI nie musi wcale spadać w te poziomy poniżej 30 i rosnąć powyżej 70. Przypomnijmy, że RSI wskazuje, jaką część kanału min-max zmian ceny (w uśrednieniu), stanowi średnia max zmiany ceny (nie mówimy o wartościach bezwzględnych). Czyli pokazuje jak w uśrednieniu oscyluje między min a max. Problem leży w ustaleniu długości okresu, w jakim uśredniać. Standardem jest 14 - i to ma sens dla zmian dziennych, nawet tygodniowych, bo jest spora szansa, że średnia obejmie ostatni "cykl". Jednak w ciągu 14 miesięcy sporo się zmienia. Powiedzmy, że przez 14 miesięcy  miesięcznie byłyby same wzrosty (hossa). Wtedy licznik będzie równy mianownikowi, bo min w mianowniku będzie zerowy - RSI równa się 1. Kiedy więc teraz następuje zmiana trendu i seria spadków, to nawet użycie średniej wykładniczej (która nadaje większą wagę niedawnej przeszłości) nie spowoduje, że RSI wpadnie do 30 - bo nadal duży wpływ ma seria nieprzerwanych wzrostów z 14 miesięcy. Stosowanie RSI ma sens, jeśli średnia obejmuje cały cykl - i w tym sensie dobrze połączyć ten wskaźnik z analizą spektralną. Kiedy ten warunek zostanie spełniony, wtedy stosowanie RSI ma racjonalne uzasadnienie w postaci twierdzenia Pukelsheima-Vysochanskij'a-Petuniego (dokładniej o tym można poczytać tutaj).

Ropa blisko szczytu?

Na pewno wielu zauważyło, że ostatnio giełda koreluje ujemnie z ceną ropy. Zobaczmy ostatni miesiąc:

Normalnie korelacja jest lekko dodatnia (ok. 30% na tygodniowych stopach zwrotach od 2018). To odwrócenie schematu jest oczywiście spowodowane "czynnikiem rosyjskim". Po agresji Rosji na Ukrainę, która dziś nabrała pełnoskalowego wymiaru (można mówić oficjalnie o wojnie, którą wywołał Putin) podaż ropy z Rosji będzie oczywiście mocno zredukowana, jeśli nie wstrzymana. Popyt się nie zmienia, więc będzie szukał innych, prawdopodobnie droższych źródeł - dlatego cena wystrzeliła. Jednocześnie wojna ma zawsze negatywne znaczenie dla gospodarek, stąd giełdy spadają. 

Jednakże z punktu widzena AT mocno wspartej przez teorię optymalizacji (zob. ten wpis) cena ropy właśnie osiągnęła max. Stosując Zasięg Fibonacciego 161,8%, otrzymamy poziom ok. 100:

W dodatku Bollinger (ustawiony standarderm) też wskazuje wykupienie - na każdej częstości: dziennej:

tygodniowej:

miesięcznej:

Dodatkowo RSI wchodzi w poziomy wykupienia - znowu dzienne, tygodniowe i miesięczne.

Można postawić tezę, że nastąpi wkrótce odwrócenie. Oczywiście wojna może okazać się "czarnym łabędziem" i cena przebije dalsze poziomy. 

W tym momencie ropa już przebiła poziom 100 i wygląda, że jest to już powyżej Fibo. 

Jeżeli nastąpi odwrócenie, to czy utrzyma się ujemna korelacja z WIG? Odwrócenie trendu ropy oznaczałoby, że nastąpiło szybkie zdyskontowanie wojny. Z ekonomicznego punktu widzenia może to polegać na tym, że pod wpływem dramatycznych wydarzeń jakie ostatnio mają miejsce na Ukrainie, kraje UE wraz z USA zaczną mocno współpracować ze sobą w celu zachowania bezpieczeństwa energetycznego i zapewnić dostarczenia paliw. W normalnych czasach zmiany odbywają się wolno, ale pandemia koronawirusa udowodniła, że procedury można znacznie przyspieszyć, jeśli sytuacja tego wymaga. Jeżeli więc Fibo się sprawdzi, to możemy oczekiwać tego typu informacji  na dniach.

Nie nabierajmy się na nazwy funduszy

W czasach, gdy stopy procentowe znajdują się poniżej naturalnego poziomu (obecnie poniżej stopy inflacji) osoby oszczędzające powinny zapomnieć o lokatach bankowych i skierować uwagę na bardziej ryzykowne formy oszczędzania - są zmuszeni zostać inwestorami. Wiadomo, że ktoś kto chce w miarę regularnie oszczędzać, nie będzie wszystkiego wrzucał w akcje, więc wybierze obligacje, surowce albo fundusze. Najłatwiej zainwestować w te ostatnie, chociaż liczba różnych funduszy może przytłoczyć. Wiadomo, że ich analiza i ocena jest czasochłonna, stąd posługując się subiektywnymi heurystykami wielu będzie dokonywać wyboru przez pryzmat samej nazwy oferty.

Np. na stronie mbanku można na ten moment znaleźć listę aż 298 funduszy: https://www.mbank.pl/serwis-ekonomiczny/notowania-funduszy/#

Jeśli szukamy alternatywy dla lokat czy obligacji skarbowych, wybierzemy jakieś dłużne, najwyżej mieszane. I tak przeglądając listę dłużnych natrafiłem na Allianz Obligacji Inflacyjnych, który mnie zainteresował. Ale właśnie zaintrygowała mnie nazwa "obligacja inflacyjna". W opisie czytamy:

Celem inwestycyjnym Allianz Obligacji Inflacyjnych jest osiąganie przychodów z lokat netto tego Subfunduszu. Subfundusz inwestuje co najmniej 70% aktywów netto w instrumenty dłużne, których cena jest powiązana z kształtowaniem się wskaźnika inflacji. Instrumenty te są emitowane przede wszystkim przez Skarb Państwa lub są poręczane lub gwarantowane przez Skarb Państwa jak również są emitowane przez państwa członkowskie i państwa należące do OECD. Emitentami instrumentów dłużnych, których cena jest powiązana z kształtowaniem się wskaźnika inflacji, mogą być również inne podmioty, w tym przedsiębiorstwa. Akcje i podobne instrumenty mogą stanowić od 0 do 10% aktywów netto. Inne fundusze mogą stanowić od 0 do 30% aktywów netto.

Następnie popatrzyłem na wykres, który powinien wedle tego opisu korelować z inflacją, a zamiast tego widzę spore spadki od paru miesięcy:

 

Czyli jeżeli przewidywałem, że inflacja mocno skoczy i chciałem się zabezpieczyć przed nią kupując jednostki tego funduszy, mocno bym się rozczarował. W powyższym opisie dwukrotnie mówi się o powiązaniu handlowanych instrumentów z inflacją. Z logicznego punktu widzenia mają rację - związek istnieje, szkoda tylko że przeciwny.

Zrobiłem porównanie rocznych stóp zwrotu tego funduszu ze stopami inflacji na podst. GUS. Notowania na koniec każdego roku od 2016 do 2021 pobrałem na dwa sposoby: ze stooq.pl i mbanku. Wyniki się nieco różnią. Stopy zwrotu ze stooq.pl w por. z inflacją prezentuje poniższa tabela:

Tylko pierwszy rok fundusz przyniósł znakomite wyniki, kilka razy wyżej od inflacji. Potem jednak od 2018 osiągał gorsze wyniki - średnio o 3,1 pkt proc.

Stopy zwrotu z mbanku w por. z inflacją prezentuje poniższa tabela:

Wyniki niewiele się różnią, chociaż różnica zwiększyła się na niekorzyść funduszu od 2018 - średnio 3,2 pkt proc.

Zwróćmy uwagę na 10,2% inflacji w 2021 i 3,5-3,6% stopy zwrotu "obligacji inflacyjnych". O ile obliczenia stopy zwrotu są oczywiste, o tyle inflacji już mniej, dlatego załączam link z analizą i obliczeniami rocznej inflacji (na podstawie miesięcznych zmian):

link

W sytuacji, gdy mamy najwyższą od ponad 20 lat inflację, wygląda to wyjątkowo źle dla Allianza. Bo tak naprawdę nie chodzi o to, że ten fundusz jest źle zarządzany czy coś takiego, bo gdy porównamy sobie go z innymi funduszami dłużnymi, to sytuacja wygląda identycznie. Jest tylko jedna różnica: żaden inny fundusz (nie tylko dłużny, ale w ogóle) nie ma w swojej nazwie "inflacyjny".

Dlaczego fundusze dłużne tak spadają? Powody są co najmniej dwa: rosnące rzeczywiste stopy procentowe NBP oraz oczekiwania, że będą one dalej podnoszone (bo muszą). Rosnąca stopa to rosnąca rentowność obligacji, czyli ich stopa dyskontowa. Gdy stopa dyskontowa rośnie, cena spada - ale pod warunkiem, że przepływy pieniężne nie rosną. Jeśli cash flow rośnie tak samo jak stopa, to cena się nie zmieni. Można by więc rozumować, że cena obligacji zmiennokuponowych (skorelowanych z WIBOR) nie będzie spadać, a tylko te o stałym kuponie czy dyskontowe. Zapomina się przy tym, że ten lepszy cash w liczniku - choć istotnie kompensuje wzrost rentowności w mianowniku - bierze się stąd, że firma więcej zapłaci inwestorom, a to zawsze jakieś ryzyko utraty płynności.

W sumie po Nationale-Nederlanden, Allianz to kolejny fundusz, który traci w moich oczach.

Wpis chwilowy - "eksperci" nagle zdziwieni skalą inflacji

Kolejne rekordy inflacji "szokują". Jak żenujące są słowa tzw. ekspertów, którzy twierdzą, że nie dało się przewidzieć skali wzrostu cen. Za co oni biorą pieniądze? To ja zrobiłem darmowe analizy, które przewidywały, że do tego dojdzie, którymi nikt jak widać się nie zainteresował. Już w 2020 r. wskazałem na ryzyko 20% inflacji, najpierw na skutek zablokowania podaży towarów: Czy powinniśmy przygotować się na inflację rzędu 20%?, potem wskazałem na niebezpieczeństwo zbyt niskiej stopy (Inflacja a luka PKB) i w końcu ostrzegałem przed tym scenariuszem, który już praktycznie ma miejsce, a który wówczas wydawał się absurdalny: Widmo 20-procentowej inflacji powraca. I co, i dziś okazuje się, że inflacja PPI skoczyła w grudniu do 14,2%.  Brak realnej odpowiedzialności za sytuację dowodzi, że należy albo ograniczyć rolę NBP w niektórych obszarach, albo wprowadzić sankcje personalne, bo nie może być tak, że prezesi i ich doradcy biorą grube pieniądze za bezmyślne naśladowanie innych krajów (bo do tego nie są potrzebne żadne kwalifikacje) albo za zaklinanie rzeczywistości i udawanie, że w ogóle coś się robi.