Trzy wykresy do oceny

Na stronie OECD można znaleźć gotowe prognozy wzrostu gospodarczego na kolejny rok i to nawet z podziałem na kwartały - link. Porównałem te prognozy z faktycznymi zmianami PKB z sąsiedniej strony.  Poniżej uzyskane efekty dla Niemiec, Polski i USA:

Korelacja prognozy z faktycznymi zmianami - dla Polski i Niemiec wynosi ok. 0,4, dla USA nieco powyżej 0,3. To znaczy, że OECD dość dobrze prognozuje PKB. Kolejne kwartały zapowiadają się słabo. Dla Niemiec 1-1,3%, USA 1-2%, Polska 2-4%. P/E S&P 500 > 27, co jest nie do utrzymania. Chociaż można spotkać opinie, że czeka nas powtórka z bańki technologicznej z 2001 r., to w pewnym stopniu ta bańka już ma miejsce.

Link do analizy

Czy z nieruchomości da się coś jeszcze wycisnąć?

W art. Czy warto kupować akcje deweloperów? pisałem, że inwestowanie w firmy deweloperskie jest nieopłacalne. Ich naturalnym rynkiem jest rynek dłużny. Dodałem jednak, że nie chodzi mi o spekulację, która ma zawsze posmak psychologiczny. Trochę żałuję, że nie kupiłem paru spółek, bo przecież było niemal pewne, że kursy będą dalej rosły. Z drugiej strony powstrzymywała mnie myśl, że to zbyt oczywiste - przecież wszyscy wiedzą, że program z kredytem 2% czy teraz 0% doprowadzi do dużego wzrostu cen nieruchomości, a więc i zysków deweloperów, stąd ich kursy musiały wzrosnąć. Od sierpnia 2022 do marca 2023 WIG-NIERUCHOMOSCI wzrósł o 70%, a WIG o 50% (poziom zero odpowiada poziomowi 100 na wykresie). 

W rzeczywistości wiedza o tym, że ceny towarów (jakimi są też nieruchomości) będą rosły w niedługim okresie nie implikuje wzrostu cen akcji w tym samym okresie. Jeżeli rynek jest efektywny, to powinien zareagować natychmiast, gdy taki program dopłat do kredytów ogłoszono, a nie być dyskontowany stopniowo. Biorąc pod uwagę zapoznanie się z programem, w praktyce mogłoby to trwać parę dni. 

Są więc trzy możliwości. Albo rynek nie wierzył do końca, że taki nieracjonalny i populistyczny program będzie realizowany, albo wpływ na zyski nie będzie wcale taki wielki, albo rynek jest nieefektywny. Pierwsze odrzucam, ponieważ program zapoczątkował PIS, który znany jest z realizacji populistycznych programów. Sprawdźmy pobieżnie drugie w BiznesRadar.

Z komponentów WIG-NRUCH wybieramy tylko te, które już mają raport za cały rok (lub 4-ty kwartał) 2023. Nie wszystkie spółki mają, a musimy poprawnie porównywać.  Na dziś w BiznesRadarze tylko 12 z 23 ma ten raport. Suma ich udziałów w indeksie równa się 82. Aby skorygować każdy z tych 12-tu udziałów do 100 stosujemy równanie na proporcję albo od razu mnożymy udział * 100/82 (lub w ogólnym wzorze 100 przez sumę nowych udziałów). tak zmodyfikowany udział mnożymy teraz przez roczną zmianę zysku - ja wybrałem dla lepszej porównywalności EBIT. Uzyskujemy średnioroczną zmianę zysku w indeksie ważoną udziałami spółek w indeksie. Okazuje się, że zysk operacyjny wzrósł ledwo o 3,67%. Sugeruje to realną stratę deweloperów w 2023 r. Co więcej, BiznesRadar także podaje wartości dla całego sektora i wychodzi mu spadek EBIT o 15,7%. Nie wiem, na ile są to dane prawidłowe, ale słabość jest potwierdzona. Poniższa tabela pokazuje szczegółowe dane:

Dodatkowo sprawdziłem jeszcze to samo, ale dla zysku netto. Efekt? Niemal identyczny, jednak wg BiznesRadaru zysk sektora spadł o 7%, co by wskazywało, że jest trochę lepiej, podczas gdy moje obliczenia pokazują że jest +3,1%, czyli nieco gorzej niż EBIT. Tabela:

OK, ale branża deweloperska charakteryzuje się tym, że zyski nie muszą korespondować z przepływami z działalności operacyjnej - księgowy zysk powstaje dopiero z chwilą zarejestrowania aktu notarialnego. Mieszkania są kupowane w ratach na podstawie umów przedwstępnych, a więc będzie to widoczne w przepływach, ale nie w rachunku zysków i strat.

Wobec tego należy sprawdzić jeszcze przepływy operacyjne w 2023. Okazuje się, że spadły o ponad 187% (natomiast w sektorze zmiana wynosiła -23%):

Głównym ciężarem dla tego spadku było ECHO (spadek o 265%), bez którego byłby wzrost 78%. W ten sposób dałoby się wyjaśnić wzrost indeksu o 70%. Niemniej odejmownie spółki, żeby "pasowało" jest zwykłym nadużyciem. ECHO stanowi 7% całego indeksu i prawie 9% użytego do powyższych obliczeń.  

Całą analizę można pobrać stąd.

Widzimy, że wcale nie jest oczywisty tak duży wzrost. Rynek zdyskontował parę lat wprzód kredyt 2% i tak samo powinno być z kredytem 0%. Pytanie tylko czy to już nastąpiło? Aby to ocenić można oczywiście próbować wyceniać akcje, ale na szybko można zastosować rozwinięcie Fibonacciego. Metoda 162% od poprzedniego dołka jest racjonalną taktyką spekulacyjną (zob. wyjaśnienie tu). Racjonalność wynika z faktu, że nie łamie teorii efektywnego rynku, ponieważ po pierwsze liczba ta jest średnią, od której zawsze będą odchylenia, a zatem gracze o mniejszej awersji do ryzyka mogą sprzedawać trochę wyżej, a o większej awersji - niżej. Po drugie 162% stanowi jedynie punkt minimalny sprzedaży, wynikający z minimalizowania czasu trzymania waloru. Nie chodzi więc o to, że od tego punktu zaczną się spadki, ale o to, że że jest to dobry moment, by sprzedać akcje, jeśli nie mamy żadnego innego punktu zaczepienia (np. gdy nie potrafimy dobrze wycenić). Jeśli pojawią się spadki od tego punktu, to bardziej wynikają z samospełniającej się przepowiedni.

Spójrzmy więc na sytuację WIG-NIERUCHOMOŚCI:

Indeks właśnie wchodzi w rejony 161,8%, dlatego zapewne pojawi się spora grupa traderów chętnych do sprzedaży. Warto tu sobie jeszcze zadać pytanie, czy Fibo indeksu ma znaczenie, jeśli nie jest on handlowany? Indeks nazywany jest benchmarkiem nie bez powodu - jeżeli traderzy wiedzą, że inni też wiedzą, że oni patrzą na to Fibo indeksu, to przewidują, że paru graczy zareaguje tylko w oparciu o indeks, a więc oni też zareagują i tak sprawdzi się samospełniająca się przepowiednia. Oczywiście taka reakcja może być chwilowa. 

Z drugiej strony znajdzie się inny argument przemawiający za porzuceniem tej strategii na samych indeksach. Mianowicie twierdzenie, na którym się opieramy zakłada cykl życia firmy - tzn. że kurs osiągnie maksimum globalne, od którego będzie spadał do zera. Indeksy ważone kapitalizacją nie spełniają oczywiście tego warunku. Stąd powodzenie ekspansji Fibo na indeksie będzie wynikało albo z tego, że akcje o dużym udziale w indeksie znajdują się w obszarze zwrotnym, albo znowu tylko z wiary samych traderów, czyli psychologii.

Dwie prognozy przemysłu na kolejny rok (gretl)

Od początku, gdy analizowałem dynamikę przemysłu z Eurostatu, wiedziałem, że są one mniej aktualne od tych podawanych przez GUS. Trzeba przyznać, że to opóźnienie jest spore, bo na dziś w Eurostacie nadal ostatnie wyniki są ze stycznia, a w GUSie od paru tygodni jest już luty. Można by więc zapytać, po co biorę dane z Eurostatu, skoro w GUSie są zawsze bardziej aktualne. Problem z GUS-owskim przemysłem jest m.in. taki, że nie jest on wyrównany sezonowo, a więc nie mogą służyć do prognozy giełdowej. Drugi problem to niekonsekwencja w podawaniu wartości za poszczególne lata: dane, które można na szybko wyciągnąć z GUSu w formie arkusza kalkulacyjnego zaczynają się (w przypadku przemysłu) od 2005 r. Tymczasem roczniki statystyczne zaczynają się od roku 2000 (Rocznik 2001 za rok 2000), a co więcej, możemy tam znaleźć roczne zmiany w produkcji przemysłowej nawet od 1994 r. Fakt, że taki brak spójności w ogóle występuje, świadczy o dalekiej niedoskonałości instytucjonalnej naszego państwa. Mówię o tym także z tego względu, że artykuł z bankiera, który omawiałem poprzednio, przedstawia wykres z danymi sprzed 2000 r.

Jeżeli GUS nie podaje danych wyrównanych sezonowo, możemy dość łatwo sami je wyrównać za pomocą X-13ARIMA-SEATS lub TRAMO-SEATS. Dodatkowy atut jaki zyskamy dzięki temu zabiegowi, będzie automatyczna prognoza na kolejny rok. Oba te narzędzia stanowią standardowe dodatki do gretla. O tym jak je zainstalować i znaleźć model napisałem w tym art. Wtedy pominąłem możliwość użycia prognozy wprost z tych dodatków. Tym razem wykorzystam wbudowane prognozy. Ale po kolei. 

X-13ARIMA-SEATS:

Algorytm X-11

Algorytm ten powstał już w 1965 r.  w amerykańskim GUS, tj. The Census Bureau. Ponieważ powstawały różne eksperymentalne warianty, nadawano im kolejne numery - od 1 do 11. W USA z jakiegoś powodu mają hopla na punkcie litery "X". W tym wypadku "X" ma wywodzić się od ekstremalności ("eXtremal") sytuacji, z jakimi program miał sobie poradzić, a więc z nietypowymi obserwacjami lub brakującymi obserwacjami. Kolejne innowacje doprowadziły do powstania X-12, który potrafił np. adaptować model do zmian parametrów w czasie [1]. W końcu połączono ten nowy wariant z zupełnie innym (stworzonym przez innych naukowców) programem - SEATS (ang. Signal Extraction in ARIMA Time Series), który filtruje i wygładza dane. W sumie to SEATS dekomponuje dane na trend i resztę [3, 5]. Tak powstał X-13.

Mamy więc roczne zmiany przemysłu (jako miesięczny szereg czasowy) z GUS od 2006: 

W opcjach 'Zmienna' szukamy X-13ARIMA:

Wybieramy X11, a więc nie X13. Wydaje się, że nie jest to błąd, ponieważ X13 dotyczy całego modelu (czyli z SEATS, a nie tylko ARIMA). Zaznaczamy co trzeba i klikamy OK. Dostajemy wykresy z wyrównaniem sezonowym i trendem:

Okazuje się, że wyrównanie sezonowe niewiele zmienia oryginalne dane. Zwróćmy uwagę jednak na szczegóły i ostatnie wartości w zbliżeniu:

Ostatnie dane wyrównane za luty wskazują na spadek, nie na wzrost, jak sugerują te oryginalne. Trend jednak jest pozytywny:

Kolejne zadanie to utworzenie prognozy. Mnóstwo statystyk, które otrzymujemy na starcie może do tego zniechęcać, ponieważ nie sa dobrze opisane. Po ich analizie doszedłem do wniosku, że interesuje nas część B 1.A  Forecasts of (prior adjusted) original series:

Zauważmy, że program wykonał prognozę na cały kolejny rok. Niewątpliwie w Excelu łatwiej się robi ładne wykresy, dlatego przeniosłem potrzebne dane do Excela i otrzymałem wykres:

Algorytm SEATS

Nieco inny model dostaniemy, zaznaczając w opcjach SEATS:

I dostajemy analogiczne wykresy

Okazuje się, że obecny trend nie wygląda już na rosnący:

Ponownie otrzymujemy długą tablicę ze statystykami. Prognozę danych skorygowanych sezonowo znajdziemy w części FORECAST OF FINAL COMPONENT 

Fajne jest tu, że mamy jasny podział na prognozę danych surowych (oryginalnych), wyrównanych sezonowo oraz samego trendu. Mimo iż naocznie trend wydaje się nierosnący, to prognoza wskazuje co innego:

Prognozowane dane wyrównane sezonowo są niemal identyczne jak te w X11, dlatego ich tu nie rysowałem. Oznacza to także, że trend można traktować po prostu jako trend w całej analizie X13-ARIMA-SEATS. Innymi słowy X-11 (czy X-12) używamy do prognozy np. dynamiki wyrównanej sezonowo, natomiast SEATS do wyłuskania samego trendu.

TRAMO/SEATS:

TRAMO (Time Series Regression with ARlMA Noise, Missing Observations, and Outliers) jest programem, który, jak sama nazwa wskazuje, estymuje parametry regresji z ARIMA, wykrywa brakujące obserwacje oraz obserwacje odstające i odpowiednio dostosowuje dla nich model. SEATS wykorzystuje TRAMO do dalszej reestymacji modelu. Jest to więc także użyteczne narzędzie w klasycznej ekonometrii. Formalnie SEATS powstał jeszcze przed TRAMO, choć oba w 1994 r. [3, 4], a potem je połączono w jedno [5]. 

Wybieramy Zmienna -> TRAMO/SEATS. Pojawia się okno.

 

Wybieramy w opcjach 'Wyniki' i zaznaczamy jak poniżej

Klikamy OK. Pojawią się analogiczne wykresy: 

Statystyki, jakie dostajemy używając procedury TRAMO/SEATS, są bardzo niejasne. Z mojej analizy wynika, że prognozy TRAMO znajdziemy pod opisem DECOMPOSITION OF THE SERIES: FORECAST

Natomiast prognozę SEATS znajdziemy na samym końcu pod opisem FORECAST OF FINAL COMPONENT:

Ten podział na dwie prognozy wynika z faktu, że mamy tu dwie procedury: najpierw TRAMO, a następnie SEATS, które dokonuje reestymacji modelu TRAMO. Wydaje się więc, że ta ostatnia prognoza jest poszukiwaną. Wydruk jest bardzo źle opisany - poszczególne prognozy nie przypisano nawet odpowiednich okresów. W każdym razie poniżej efekt, tego co wyłuskałem z TRAMO-SEATS:

Wg dodatkowej analizy, jaką przeprowadziłem, TRAMO-SEATS najlepiej sprawdza się w krótkim terminie, dlatego nie ma sensu brać na poważnie prognozy na cały rok. X13-ARIMA sprawdza się trochę lepiej w dłuższym terminie. Patrząc na oba wykresy danych wyrównanych sezonowo z prognozami, dochodzimy do wniosku, że kolejne miesiące powinny być w miarę pozytywne dla przemysłu. X-13ARIMA sugeruje jakiś szczyt w sierpniu. Niemniej trend pozostaje rosnący - ale oczywiście to wszystko może się zmieniać z miesiąca na miesiąc.

Literatura:

[1] Findley, D. F., B. C. Monsell, W. R. Bell, M. C. Otto, and B. C. Chen (1998). New capabilities of the X-12-ARIMA seasonal adjustment program (with discussion). Journal of Business and Economic Statistics 16, 127–177.

[2] Findley, D. F., D. P. Lytras, and A. Maravall (2016). Illuminating ARIMA Model-Based Seasonal Adjustment with Three Fundamental Seasonal Models. pp. 1–42.

[3] Gomez, V., Maravall, A. (1994a), Program SEATS Signal Extraction in ARIMA Time Series: Instructions for the User, Working Paper ECO 94/28, European University Institute, Florence.

[4] Gomez, V., Maravall, A. (1994b), Program TRAMO 'Time Series Observations With ARIMA Noise, Missing Observations and Outliers': Instructions for the User, Working Paper ECO 94/31, European University Institute, Florence.

[5] Gomez, V., Maravall, A. (1996). Programs TRAMO (Time Series Regression with Arima noise, Missing observations, and Outliers) and SEATS (Signal Extraction in Arima Time Series). Instruction for the User. Working Paper 9628.

Czy strategia Bankiera działa?

Od paru lat portal bankier.pl promuje strategię inwestycyjną opisaną w artykule Prosta instrukcja obsługi warszawskiej giełdy. Została ona sprowadzona do dwóch haseł: 

1) kupuj na dołku produkcji przemysłowej

2) kupuj, gdy leje się krew na rynku złotego.

Zanim zajmę się ich testowaniem, zwrócę uwagę na niezbyt poprawny tytuł. O ile "inwestowanie" rzeczywiście oznacza kupowanie i trzymanie jakiś czas aktywów, to słowo "instrukcja" jest pojęciem szerszym i powinna zawierać oprócz strategii kupowania także sposób sprzedawania. Oczywiście można by dwie przytoczone strategie odwrócić dla momentu sprzedaży, ale intuicja nie wystarczy i autor powinien przynajmniej wspomnieć coś o zbywaniu. W sumie, to wystarczyłoby zmienić tytuł np. na "Prosta strategia inwestowania na GPW" i nie byłoby tematu. Ale zajmijmy się już istotą.

Ad 1) Kupuj na dołku produkcji przemysłowej

Dokładniej autor proponuje trzymiesięczną średnią rocznej dynamiki produkcji przemysłowej, która spada poniżej zera. Autor nie wskazuje źródła pobranych danych, nie pisze też wprost czy chodzi o polski przemysł - tego trzeba się domyślić z kontekstu. Dane o polskim przemyśle pobrałem ze strony Eurostatu (odnośnik - należy zmienić ustawienia na "Poland", "Seasonally and calendar adjusted data" oraz "Percentage change on previous period"). Aby uzyskać wymagane dane, musiałem surowe dane obrobić:

1. pobrałem miesięczne zmiany procentowe produkcji przemysłowej dla Polski,

2. zamieniłem na logarytmiczne stopy zmian,

3. aby dostać roczną dynamikę, dodałem do siebie kolejnych 12 stóp,

4. uśredniłem 3 kolejne stopy.

Dwie uwagi. Pierwsza dotyczy polskich odpowiedników nazw angielskich. Otóż industry oraz manufacturing są uznawane przez słowniki za synonimy i tłumaczone jako produkcja przemysłowa czy po prostu przemysł, ale jest to ich potoczne rozumienie. W rzeczywistości ten drugi termin oznacza przemysł przetwórczy, inaczej przetwórstwo przemysłowe, które ogólnie mówiąc stanowi produkcję pośrednią między dostarczaniem surowców a wyrobami gotowymi. Tak więc industry jest szerszym pojęciem, a manufacturing jest tylko jego składnikiem. 

Po drugie powinniśmy zdać sobie sprawę, że wielkość produkcji nie jest podawana na określony moment, ale w przedziale czasowym. Wynika to zarówno z faktu, że jest ona procesem, a nie stanem oraz z tego, że stanowi ona przychód przedsiębiorstwa produkcyjnego w danym okresie. Natomiast kurs giełdowy można porównać z wartością bilansową, czyli powstającą na dany moment, czyli formalnie jest nieporównywalny z wielkością produkcji. Aby sprowadzić obie wielkości do porównywalności, badamy ich zmiany, np. miesiąc do miesiąca lub rok do roku. Warto więc sobie uświadomić, że analizując produkcję, badamy jakby "zmiany ze zmiany" kapitału gospodarki. Zmiana między jednym momentem a drugim powoduje, że możemy przyrównać okresy do siebie dla dwóch różnych serii i w ten sposób je porównać. Choć to jedynie kwestie techniczne, to rzutują na poprawność wniosków. Z jednej strony informacje ze sfery makro przychodzą z opóźnieniem, z drugiej strony należy uważać, na który okres obliczamy zmiany. Jeśli popełnimy błąd, to błędnie wnioskować, że produkcja wyprzedza giełdę. W kontekście tego badania taki błąd nie jest jednak wielką przeszkodą, bo nie chodzi o idealne przewidzenie momentu dołka, ale okolicy. Zawsze jednak lepiej zachować precyzję.

Eurostat podaje dynamikę (procentową zmianę) produkcji na dany miesiąc i rok, a więc okres dynamiki to (początek miesiąca - koniec miesiąca). Gdybyśmy chcieli te wartości porównać ze stopami zwrotu z WIG, musielibyśmy odjąć zamknięcie w ostatnim dniu miesiąca od analogicznego zamknięcia z poprzedniego miesiąca. Jednak wiemy, że nie chodzi w naszym badaniu o stopy zwrotu, tylko o szacowanie okresu zmian trendu. Wynika z tego, że możemy wykorzystać kod R, jaki napisałem w tym artykule, gdzie porównywałem lokalne ekstrema WIG20 z sWIG80.

To po kolei. Nałożenie WIG i dynamiki przemysłu od marca 2001 do stycznia 2024:

 

Następnie filtruję oba wykresy (filtr lowess) i oznaczam lokalne ekstrema, przy czym dla przemysłu minimum ma być poniżej zera, a maksimum powyżej. 

Lokalne minima:

Pionowe linie zaznaczają dołki - czerwone przemysłu, niebieskie WIGu. Niektóre się nakładają, dlatego nie widać dwóch czerwonych, ale łatwo je usytuować, bo to najbardziej skrajne dna. W sumie mamy tu 4 czerwone dołki i 7 niebieskich. Przeanalizujmy je:

- pierwszy dołek WIGu jest w zasadzie opóźniony w stosunku do przemysłu o ponad pół roku (przemysł dobrze go prognozuje),

- drugi dołek WIGu spotyka się z dołkiem przemysłu,

- trzeci dołek WIGu wyprzedza dołek przemysłu o 1 rok,

- czwarty dołek giełdowy pokrywa się z przemysłem,

- piąty dołek WIGu tak naprawdę nie jest dołkiem, a tempo przemysłu - chociaż spada - nie wskazuje na żadne zwroty (czyli przemysł prawidłowo mówi, aby trzymać),

- szóste minimum WIGu zrównuje się z minimum przemysłu

- siódmy dołek nie koresponduje z przemysłem, który jeszcze trochę spada. 

Sumarycznie okazuje się, że dołki produkcji przemysłowej nie wyprzedzają WIGu, ponieważ jeden raz go wyprzedziły, drugi raz na odwrót, w pozostałych występują w tym samym czasie lub w ogóle ich nie ma. Niemniej rzeczywiście dołek przemysłu ma walor informacyjny i pozwala przewidzieć dalszy kierunek WIGu. Dodatkowa informacja wynika z faktu, że filtr WIGu nie może zastąpić tutaj filtru przemysłu, ponieważ w jednym przypadku dołek WIGu dał fałszywy sygnał - pod koniec roku 2018, kiedy niby było dno, ale wzrost prawie żaden - a tymczasem tempo przemysłu było ciągle dodatnie. Wcześniej, w roku 2016 także filtr WIGu osiąga lokalne minimum, a przemysł nie (tzn. spada, ale jest powyżej zera), co można byłoby uznać za argument przeciwko przemysłowi, gdybyśmy nastawiali się na krótki lub średni termin - była to mocniejsza korekta trwająca niecały rok. W dłuższym terminie nie ma tu błędu i rzeczywiście przemysł prognozuje hossę. De facto dynamika przemysłu także spadała w tym czasie, więc również można było tu zmniejszać pozycję, chociaż opieranie się tylko na tym nie miałoby sensu, choćby dlatego, że mamy inne przykłady, kiedy dynamika ta także spadała, a jednak WIG rósł. 

Popatrzmy teraz na ostatnie dołki. Minimum giełdy w 2022 r. wyprzedza lokalne minimum przemysłu o niecały rok. Tzn. dynamika produkcyjna znajduje się nadal pod kreską (choć minimalnie) i nie wiadomo czy w ogóle możemy nazwać dołkiem ten okres, ale istnieje na to szansa. Zobaczmy w zbliżeniu ten obszar:

Lokalne maksima:

Oznaczenia analogicznie jak przy dołkach. Spójrzmy:

- pierwszą górkę WIGu zapowiadała wyraźnie zmiana trendu przemysłu (jednak pierwsza górka przemysłu została zignorowana przez WIG), 

- drugą górkę WIGu podobnie wskazywała wcześniejsza górka przemysłu, chociaż bessa giełdowa była wyjątkowo krótka, a kolejny dołek przemysłu powstał trochę późno,

- trzecia górka WIGu była sygnalizowana przez sferę realną już ok. rok wcześniej, chociaż sfera ta szybko zawinęła w górę, co z pewnością dawało schizofreniczny obraz; to ten przypadek, gdy bardziej liczyła się umiejętność spekulacji (krótki trend spadkowy),

- czwarte maksimum giełdowe wyprzedziło minimalnie maksimum przemysłowe,

- przy piątym to przemysł nieco szybciej zareagował od WIGu.

Przy szczytach sytuacja się o tyle komplikuje, że raz giełda zignorowała duże spadki dynamiki produkcji, a przy wielu innych czas spadków indeksu był krótki, tak że w praktyce można było nawet nie zdążyć dobrze wyjść, kiedy trzeba było za chwilę wejść.

Ad 2) Kupuj, gdy leje się krew na rynku złotego

Kiedy się wczytamy w treść propozycji, możemy się zawieść. Autor nie wskazuje już warunku zakupu (ani tym bardziej sprzedaży), a jedynie luźno sugeruje, że jeżeli euro w ciągu roku stosunkowo mocno się umocniło do złotego, to możemy rozważać moment zakupu akcji. Podaje pewne liczby, ale niewiele one wnoszą. Jedynie co możemy zrobić, to połączyć poprzedni warunek z warunkiem, aby euro się wzmocniło w ciągu roku. Autor sugeruje min. 10%, ale jakoś nie przekonuje mnie to - nie ma powodu, żeby nie mogłoby być to 5%. Popatrzmy.

I sfiltrujmy też, ale tym razem dołki WIG będą korespondowały z górkami euro.

Lokalne minima:

Niebieskie proste wskazują dołki WIGu, czerwone - górki EUR/PLN. Prawdą jest, że jeżeli maksima euro znajdują się części dodatniej, WIG zaczyna hossę. Tak więc na obecną chwilę warunek ten nie zostaje spełniony, wręcz przeciwnie.

Lokalne maksima:

Odwrotna sytuacja, tj. górki na indeksie (niebieska prosta) vs. dołki na walucie (czerwona):

Maksima okazują się dużo bardziej problematyczne. Nawet, jeśli symetrycznie weźmiemy po uwagę tylko dołki walutowe z części ujemnej, to nie dostajemy już sygnału sprzedaży przed rokiem 2013. Gdyby przyjąć, że dołki euro zwiastują bessę, to należałoby się przygotowywać na jej nadejście. Obiektywnie jednak nie ma na to argumentu. 

Podsumowanie

Teraz jasne staje się, dlaczego autor omawianego artykułu pominął kwestię sprzedawania. Obie jego strategie mają potencjał tylko przy znajdowaniu dołków na rynku akcji (tym bardziej należałoby zmienić tytuł). Jako całość nie jest to nawet strategia sensu stricte. Można powiedzieć, że połowicznie poprawna, ale niedokończona. Jeśli chodzi o prognozę na przyszłość, to nie ma jednoznacznego sygnału. Ostatni lokalny dołek przemysłu sugeruje dalsze wzrosty giełdowe, ale nie jest spełniony warunek słabości złotego. W dodatku jeśli przemysł się teraz zagnie i wróci do dołu, to czeka nas co najmniej korekta, a prędzej czy później rynek niedźwiedzia. Wydaje mi się, że lepiej zamiast na polski przemysł spoglądać na zagraniczną produkcję. Na gospodarkę w Polsce wpływa produkcja UE (zob. tu).  

Dodatek. Kod w R:

# wykorzystanie funkcji Evan'a Friedlanda (https://stackoverflow.com/questions/6836409/finding-local-maxima-and-minima)

inflect <- function(x, czy_przemysl, threshold){

  up   <- sapply(1:threshold, function(n) c(x[-(seq(n))], rep(NA, n)))

  down <-  sapply(-1:-threshold, function(n) c(rep(NA,abs(n)), x[-seq(length(x), length(x) - abs(n) + 1)]))

  a    <- cbind(x,up,down)

  

  if (czy_przemysl==TRUE) {

    list(minima = which(apply(a, 1, min) == a[,1] & a[,1] < 0), maxima = which(apply(a, 1, max) == a[,1] & a[,1] > 0))

  } else {

    list(minima = which(apply(a, 1, min) == a[,1]), maxima = which(apply(a, 1, max) == a[,1]))

  }

}

linieEkstremum <- function(filtr, maks) {

  

  n <- 2

  if (deparse(substitute(filtr)) == "filtrPrzemysl") {

    czy_przemysl <- TRUE

    kolor <- "red"

  } else {

    czy_przemysl <- FALSE

    kolor <- "blue"

  }

  

  bottoms <- lapply(1:n, function(th) inflect(filtr, czy_przemysl, threshold = 10)$minima)

  tops <- lapply(1:n, function(th) inflect(filtr, czy_przemysl, threshold = 20)$maxima)

  

  for(i in n:n){

    if (maks==1) {

      #abline(v=tops[[i]], col=kolor)

      return(tops[[i]])

    } else {

      #abline(v=bottoms[[i]], col=kolor)

      return(bottoms[[i]])

    }

  }

}

# potrzebne pakiety

if (require("zoo")==FALSE) {

  

  install.packages("zoo")

}

library(zoo)

if (require("data.table")==FALSE) {

  

  install.packages("data.table")

}

library(data.table)

if (require("lubridate")==FALSE) {

  

  install.packages("lubridate")

}

library(lubridate)

#zamieniam "\" na "/"

sciezka <- r"(C:\Documents\R\testy)"

sciezka <- gsub("\\", "/", sciezka, fixed=T)

setwd(sciezka)

# Przemysł

# link: https://ec.europa.eu/eurostat/databrowser/view/sts_inpr_m__custom_10635116/default/table?lang=en

nazwaPliku1 <- "sts_inpr_m_page_linear.csv"

tabela_przemysl <- fread(nazwaPliku1, tz="")

tabela_przemysl$TIME_PERIOD <- paste0(tabela_przemysl$TIME_PERIOD, "-01")

tabela_przemysl$TIME_PERIOD <- ceiling_date(as.IDate(tabela_przemysl$TIME_PERIOD), "month") - days(1)

tabela_przemysl$TIME_PERIOD <- format(tabela_przemysl$TIME_PERIOD, "%Y-%m")

tabela_przemysl <- tabela_przemysl[, c("TIME_PERIOD", "OBS_VALUE")]

przemysl <- tabela_przemysl[, "OBS_VALUE"]

przemysl_rocznie <- c(rep(NA, 11), rollsum(log(1+przemysl/100)*100, k=12, align="right"))

tabela_przemysl$Przemysl_rocznie <- przemysl_rocznie

tabela_przemysl$SMA3 <- c(rep(NA, 2), round(rollmean(przemysl_rocznie, k=3, align="right"), 4))

# WIG, dane ze stooq.pl

nazwaPliku <- "wig_d.csv"

#tabela_wig <- fread(nazwaPliku, tz="")

# jeżeli jest ostrzeżenie, to prawdopodobnie sa jakieś braki, stosujemy ich wypełnienie

tabela_wig <- fread(nazwaPliku, fill=TRUE, tz="")

tabela_wig$Data <- as.Date(tabela_wig$Data)

# Zakładając, że twoja ramka danych nazywa się tabela_wig, a kolumna z datą to Data

tabela_wig$Data <- as.Date(tabela_wig$Data) 

tabela_wig$RokMies <- format(tabela_wig$Data, "%Y-%m")

# Oblicz ostatni dzień każdego miesiąca

ost_dzien_mies_tab <- aggregate(tabela_wig$Data ~ tabela_wig$RokMies, FUN=max)

colnames(ost_dzien_mies_tab) <- c("RokMiesiac", "OstatniaData")

# Wyfiltruj wiersze, gdzie Data to ostatni dzień miesiąca

tabela_wig <- tabela_wig[tabela_wig$Data %in% ost_dzien_mies_tab$OstatniaData, c(1, 5)]

tabela_wig$Data <- format(tabela_wig$Data, "%Y-%m")

# Połączenie danych

dane <- merge(x=tabela_przemysl, y=tabela_wig, by.x="TIME_PERIOD", by.y="Data")

dane <- as.data.frame(na.omit(dane[, c("TIME_PERIOD", "SMA3", "Zamkniecie")]))

colnames(dane) <- c("Okres", "Przemysl_rr", "WIG")

daty <- as.Date(paste0(dane$Okres, "-01"))

rok = year(daty)

mc = month(daty)

przemysl <- ts(dane$Przemysl_rr, start=c(rok[1], mc[1]), frequency=12)

wig <- ts(dane$WIG, start=c(rok[1], mc[1]), frequency=12)

# Wykres

par(mar=c(5, 5, 3, 5))

# Daty, indeks dat

datyWykres <- seq(from=daty[1], to=daty[length(daty)], length.out=year(daty[length(daty)])-year(daty[1]))

indeksWykres <- seq(from=1, to=length(daty), length.out=year(daty[length(daty)])-year(daty[1]))

# Wykres przemysłu

plot(x=index(index(przemysl)), y=przemysl, xlab="Okres", ylab="Dynamika przemysłu (Polska)", lwd=2, xaxt="n", type="l", col="darkred")

axis(side=1, at=indeksWykres, labels=format(datyWykres, "%Y-%m"), cex.axis=0.8)

abline(h=0, col="darkgray")

# Wykres WIG

par(new=TRUE)

plot(x=index(index(przemysl)), y=wig, log="y", axes=FALSE, xlab="", ylab="", col="darkblue", lwd=2, type="l")

axis(side=4)

mtext("WIG", side=4, line=3)

legend(x="topleft", legend=c("Przemysł (Polska)\n log-zmiany % r/r", "WIG"), lty=c(1,1),

       lwd=c(2,2), col=c("darkred", "darkblue"), cex=0.7)

# filtr lowess

filtrPrzemysl <- lowess(przemysl, f=0.03)$y

filtrPrzemysl <- ts(filtrPrzemysl, start=c(rok[1], mc[1]), frequency=12)

filtrWig <- lowess(wig, f=0.05)$y

filtrWig <- ts(filtrWig, start=c(rok[1], mc[1]), frequency=12)

# Wykres filtrów

## Filtr przemysłu

plot(x=index(index(filtrPrzemysl)), y=filtrPrzemysl, xlab="Okres", ylab="Filtr przemysłu (Polska)", lwd=2, xaxt="n", type="l", col="darkred")

axis(side=1, at=indeksWykres, labels=format(datyWykres, "%Y-%m"), cex.axis=0.8)

abline(h=0, col="darkgray")

## Filtr WIG

par(new=TRUE)

plot(x=index(index(filtrPrzemysl)), y=filtrWig, log="y", axes=FALSE, xlab="", ylab="", col="darkblue", lwd=2, type="l")

axis(side=4)

mtext("Filtr WIG", side=4, line=3)

## Dołki przemysł

abline(v=linieEkstremum(filtrPrzemysl, maks=0), col="red")

## Dołki WIG

abline(v=linieEkstremum(filtrWig, maks=0), col="blue")

legend(x="topleft", legend=c("Przemysł (Polska)\n log-zmiany % r/r", "WIG"), lty=c(1,1),

       lwd=c(2,2), col=c("darkred", "darkblue"), cex=0.5)

## Górki przemysł

abline(v=linieEkstremum(filtrPrzemysl, maks=1), col="red")

## Górki WIG

abline(v=linieEkstremum(filtrWig, maks=1), col="blue")

legend(x="topleft", legend=c("Przemysł (Polska)\n log-zmiany % r/r", "WIG"), lty=c(1,1),

       lwd=c(2,2), col=c("darkred", "darkblue"), cex=0.5)

# Zbliżenie

datyWykres <- seq(from=daty[1], to=daty[length(daty)], length.out=length(daty))

indeksWykres <- seq(from=1, to=length(daty), length.out=length(daty))

pocz <- index(index(przemysl))[length(przemysl)-24]

kon <- index(index(przemysl))[length(przemysl)]

plot(x=index(index(przemysl)), y=filtrPrzemysl, xlab="Okres", ylab="Filtr przemysłu (Polska)", lwd=2, xaxt="n", type="l", col="darkred", xlim=c(pocz, kon))

axis(side=1, at=indeksWykres, labels=format(datyWykres, "%Y-%m"), cex.axis=0.8)

abline(h=0, col="grey")

abline(v=indeksWykres, col="grey")

par(new=TRUE)

plot(x=index(index(filtrPrzemysl)), y=filtrWig, log="y", axes=FALSE, xlab="", ylab="", col="darkblue", lwd=2, type="l", xlim=c(pocz, kon))

axis(side=4)

mtext("Filtr WIG", side=4, line=3)

legend(x="bottomleft", legend=c("Przemysł (Polska)\n log-zmiany % r/r", "WIG"), lty=c(1,1),

       lwd=c(2,2), col=c("darkred", "darkblue"), cex=0.8)

########

# EUR/PLN, dane ze stooq.pl

nazwaPliku <- "eurpln_d.csv"

#tabela_eur <- fread(nazwaPliku, tz="")

# jeżeli jest ostrzeżenie, to prawdopodobnie sa jakieś braki, stosujemy ich wypełnienie

tabela_eur <- fread(nazwaPliku, fill=TRUE, tz="")

tabela_eur$Data <- as.Date(tabela_eur$Data)

# Zakładając, że twoja ramka danych nazywa się tabela_eur, a kolumna z datą to Data

tabela_eur$Data <- as.Date(tabela_eur$Data) 

tabela_eur$RokMies <- format(tabela_eur$Data, "%Y-%m")

# Oblicz ostatni dzień każdego miesiąca

ost_dzien_mies_tab <- aggregate(tabela_eur$Data ~ tabela_eur$RokMies, FUN=max)

colnames(ost_dzien_mies_tab) <- c("RokMiesiac", "OstatniaData")

# Wyfiltruj wiersze, gdzie Data to ostatni dzień miesiąca

tabela_eur <- tabela_eur[tabela_eur$Data %in% ost_dzien_mies_tab$OstatniaData, c(1, 5)]

tabela_eur$Data <- format(tabela_eur$Data, "%Y-%m")

# 

eur <- tabela_eur$Zamkniecie

eur_rocznie <- c(rep(NA, 12), 100*diff(eur, lag=12) / head(eur, -12)) 

tabela_eur$Eur_rocznie <- eur_rocznie

# Połączenie danych

dane <- merge(x=tabela_eur, y=tabela_wig, by="Data")

dane <- as.data.frame(na.omit(dane[, c("Data", "Eur_rocznie", "Zamkniecie.y")]))

colnames(dane) <- c("Okres", "EURPLN_rr", "WIG")

daty <- as.Date(paste0(dane$Okres, "-01"))

rok = year(daty)

mc = month(daty)

eur <- ts(dane$EURPLN_rr, start=c(rok[1], mc[1]), frequency=12)

wig <- ts(dane$WIG, start=c(rok[1], mc[1]), frequency=12)

# Wykres

par(mar=c(5, 5, 3, 5))

# Daty, indeks dat

datyWykres <- seq(from=daty[1], to=daty[length(daty)], length.out=year(daty[length(daty)])-year(daty[1]))

indeksWykres <- seq(from=1, to=length(daty), length.out=year(daty[length(daty)])-year(daty[1]))

# Wykres eurpln

plot(x=index(index(eur)), y=eur, xlab="Okres", ylab="Stopa zmian EUR/PLN", lwd=2, xaxt="n", type="l", col="darkred")

axis(side=1, at=indeksWykres, labels=format(datyWykres, "%Y-%m"), cex.axis=0.8)

abline(h=c(-10, -5, 0, 5, 10), col="darkgray")

par(new=TRUE)

# Wykres WIG

plot(x=index(index(eur)), y=wig, log="y", axes=FALSE, xlab="", ylab="", col="darkblue", lwd=2, type="l")

axis(side=4)

mtext("WIG", side=4, line=3)

legend(x="topleft", legend=c("EUR/PLN zmiany % r/r", "WIG"), lty=c(1,1),

       lwd=c(2,2), col=c("darkred", "darkblue"), cex=0.7)

# filtr lowess

filtrEur <- lowess(eur, f=0.03)$y

filtrEur <- ts(filtrEur, start=c(rok[1], mc[1]), frequency=12)

filtrWig <- lowess(wig, f=0.05)$y

filtrWig <- ts(filtrWig, start=c(rok[1], mc[1]), frequency=12)

# Wykres filtrów

## Filtr eur

plot(x=index(index(filtrEur)), y=filtrEur, xlab="Okres", ylab="Filtr zmian EUR/PLN", lwd=2, xaxt="n", type="l", col="darkred")

axis(side=1, at=indeksWykres, labels=format(datyWykres, "%Y-%m"), cex.axis=0.8)

abline(h=0, col="darkgray")

## Filtr WIG

par(new=TRUE)

plot(x=index(index(filtrEur)), y=filtrWig, log="y", axes=FALSE, xlab="", ylab="", col="darkblue", lwd=2, type="l")

axis(side=4)

mtext("Filtr WIG", side=4, line=3)

## Górki eur

abline(v=linieEkstremum(filtrEur, maks=1), col="red")

## Dołki WIG

abline(v=linieEkstremum(filtrWig, maks=0), col="blue")

legend(x="topleft", legend=c("EURPLN zmiany % r/r", "WIG"), lty=c(1,1),

       lwd=c(2,2), col=c("darkred", "darkblue"), cex=0.7)

## Dołki eur

abline(v=linieEkstremum(filtrEur, maks=0), col="red")

## Górki WIG

abline(v=linieEkstremum(filtrWig, maks=1), col="blue")

legend(x="topleft", legend=c("EURPLN zmiany % r/r", "WIG"), lty=c(1,1),

       lwd=c(2,2), col=c("darkred", "darkblue"), cex=0.7)

Tuskowy efekt flagi to kiepska taktyka

Ostatnie 2 tygodnie charakteryzowały się po raz kolejny swego rodzaju dywergencją między polską giełdą a S&P 500, DAX-em i CAC40:

WIG spada, reszta rośnie. Analitycy raczej nie zauważyli przyczyny takiego stanu rzeczy. Z dużym prawdopodobieństwem przyczyną były wypowiedzi rządu, w szczególności D. Tuska na temat zagrożenia ze strony Rosji.

Jest zadziwiające, że media w ogóle nie biorą pod uwagę możliwości, że straszenie Tuska wojną z Rosją w ostatnich tygodniach może mieć źródło w nadchodzących wyborach. Dziennikarze i publicyści bezkrytycznie poddają się tej narracji jak małe dzieci. Tusk zapewne opracował sobie taktykę, że będzie straszył tylko do pewnego stopnia, aby nie wywołać paniki, ale to paradoksalnie działa na media niczym scena w horrorze, w której niby nic się nie dzieje, ale każdy wie, że zaraz coś się wydarzy, co wywołuje automatycznie napięcie i strach. Niby wywoływanie strachu ma komercyjny sens, ale widząc co mówią i jak się dziennikarze zachowują, to można odnieść wrażenie, że naprawdę w to wierzą, że za chwilę Rosja zaatakuje Polskę. Wspólna wyprawa Dudy i Tuska do Białego Domu podgrzewa całą atmosferę. 

Może najpierw parę faktów. Wiele miesięcy przed wojną USA ostrzegało świat przed atakiem Rosji na Ukrainę. Obecnie żadnych tego typu informacji nie ma, a wręcz przeciwnie - właśnie powstał raport amerykańskiego wywiadu na temat zagrożenia ze strony Rosji. Jak czytamy "Rosja prawie na pewno nie chce bezpośredniego konfliktu zbrojnego z USA i NATO" (link). Pełen raport znajdziemy tutaj i tam na str. 14 pada to zdanie. Zauważmy, że jak tylko informacja o tym raporcie się pojawiła, WIG gwałtownie odbił.

Z tym efektem flagi jest jednak kłopot. Bo po pierwsze jest to reakcja krótkoterminowa, a po drugie po silnej reakcji w górę, często następuje silny spadek w dół. PO nie będzie w stanie utrzymać społeczeństwa w takim pół-napięciu, więc będzie musiała albo zwiększyć "dawkę", albo zacząć akcję uspokajania.

Ale czy to w ogóle ma sens? Seo i Horiuchi zbadali wpływ militarnych sporów na poparcie liderów w 27 krajach między 2008 a 2014. Okazuje się, że przywódcy państw, które bronią się przed obcą agresją, nie odczuwają żadnej zmiany w poparciu społecznym, podczas gdy ci, którzy inicjują konflikty ofensywne, spotykają się z negatywną reakcją opinii publicznej. Poniższy rysunek z przytoczonej pracy przedstawia dokładniej te wyniki:

Kraj broniący się przed agresją jest po lewej. Poparcie na plus praktycznie się nie zmienia, ale dezaprobata wzrasta o 1.09 pkt proc. Po prawej jest kraj agresora - poparcie dla lidera spada z każdej strony prawie o 3 pkty proc. Czerwony kolor oznacza, że wyniki są istotne statystycznie.

W bardzo krótkim terminie może nastąpić jednak mała poprawa, na co wskazuje to samo badanie (chociaż to zaledwie parę tygodni).

Polska byłaby więc po lewej. Chociaż są to wyniki nieistotne statystycznie, to mimo wszystko dezaprobata raczej wzrasta. 

Wynika z tego, że "efekt flagi" jest albo iluzoryczny, albo zbyt krótkoterminowy, żeby opierać na nim predykcje. Oryginalnie był on zauważony tylko dla USA, dokładniej dotyczył poparcia dla bieżącego prezydenta. W polskim systemie poparcie to dotyczyłoby więc premiera lub prezydenta. 

Gdyby chcieć zażartować, powiedziałbym, że obecna władza po prostu chce dokopać inwestorom: co wypowiedź, to jakieś tąpnięcie. To też pokazuje, jaką nieprzewidywalnością cechuje się giełda wbrew technikom. Pozytywnym sygnałem jest pojawienie się tego amerykańskiego raportu, który powinien doprowadzić do zaniku wspomnianej dywergencji.

 

Literatura:

Seo, T,, Horiuchi, Y., Natural Experiments of the Rally ’Round the Flag Effects Using Worldwide Surveys, Jan 2022.

Dlaczego oddzielne badanie autokorelacji jest błędem?

 W wielu naukach, szczególnie w finansach, dużą wagę przywiązuje się do autokorelacji. Jeżeli stopy zmian jakoś autokorelują, oznacza to możliwość lepszej prognozy dalszego kierunku. Niewielu jednak wie, że proste badanie autokorelacji jest błędem. Nawet jeżeli testy sugerują istotność statystyczną autokorelacji, to w rzeczywistości może być ona pomylona z innym efektem, jeśli jest sprawdzana osobno.

W poniższych przykładach używam języka R, w tym dwa pakiety: forecast oraz rugarch:

if (require("forecast")==FALSE) {

  install.packages("forecast")

  library("forecast")

}

if (require("rugarch")==FALSE) {

  install.packages("rugarch")

  library("rugarch")

}

Przykład 1. 

Najbardziej oczywisty - może się zmienić średnia. Budujemy dwie próbki o rozkładzie normalnym z różnymi średnimi, i łączymy w jedną próbkę. Aby mieć właściwe porównanie wszędzie będę stosował to samo ziarno nr 44.

set.seed(44)

rnorm_test1 <- rnorm(n=500, mean=1)

rnorm_test2 <- rnorm(n=500, mean=4)

rnorm_test <- c(rnorm_test1, rnorm_test2)

Pacf(rnorm_test)

plot(x=rnorm_test[-length(rnorm_test)], y=rnorm_test[-1])

abline(lsfit(x=rnorm_test[-length(rnorm_test)], y=rnorm_test[-1]), lwd=2, col="blue")

 

Pozorna autokorelacja cząstkowa sięga 8 rzędów. Poniższy wykres wskazuje czarno na białym, z czego to wynika.

Przykład 2. 

A) Mniej oczywisty przykład dotyczy efektu ARCH. Problem ten analizował m.in. A. Bera [1]. Zrobię symulację procesu w pakiecie rugarch - najpierw wskazuję o jaki model mi chodzi oraz ustalam mu parametry.

Specyfikację modelu robimy za pomocą funkcji ugarchspec(). Po pierwsze wybieram rodzaj modelu. Będzie to standardowy GARCH, oznaczony tu "sGARCH". Po drugie ustawiam liczbę rzędów procesu. Ponieważ interesuje mnie tylko ARCH(1) ustawiam garchOrder = c(1, 0), co oznacza tworzenie efektu ARCH pierwszego rzędu - warunkowa wariancja zależy tylko od składnika losowego modelu. Części warunkowej średniej, tj. ARMA, nie trzeba już ustawiać, ale można - czyli armaOrder = c(0, 0), aby nie mieć żadnych autokorelacji w szeregu czasowym. Można też dodać distribution.model = "norm", aby zaznaczyć rozkład normalny. Nie jest to konieczne, ale dodam też include.mean = FALSE, czyli nie zawieraj wyrazu wolnego dla części średniej warunkowej, bo nie jest potrzebny. Ostatnią sprawą do wpisania są parametry modelu za pomocą fixed.pars. Ponieważ mamy tylko ARCH(1) w rozkładzie normalnym wystarczy ustawić wyraz wolny w równaniu wariancji (oznaczony jako omega, np. na poziomie 0.1) oraz współczynnik kierunkowy ARCH(1) (tj. efekt MA(1) w równaniu wariancji, oznaczony jako alpha1, np. równy 0.9). 

Po tym etapie wykonujemy symulację przy pomocy funkcji ugarchpath(). Musi zawierać co najmniej 2 argumenty - rodzaj modelu, który utworzyliśmy za pomocą ugarchspec() oraz liczbę obserwacji (oznaczona n.sim). Tak zbudowana symulacja jest obiektem, który zostanie nazwany simArchModel. Żeby wyciągnąć sam szereg czasowy możemy wtedy użyć komendy simArchModel@path$seriesSim.

set.seed(44)

myGarch1 <- ugarchspec(variance.model = list(model = "sGARCH", 

                                             garchOrder = c(1, 0)), 

                       distribution.model = "norm", mean.model = list(armaOrder = c(0, 0), include.mean = FALSE), 

                       fixed.pars=list(omega=0.1, alpha1=0.9))

simArchModel = ugarchpath(myGarch1, n.sim = 1000)

simArch1 <- simArchModel@path$seriesSim

Pacf(simArch1, plot=F)

Pacf(simArch)

simArchLm <- lm(tail(simArch1, -2) ~ head(simArch1, -2))

plot(x=head(simArch1, -2), y=tail(simArch1, -2))

abline(simArchLm)

PACF sugeruje 2 pierwsze rzędy ujemnej korelacji. Statystyki:  

summary(simArchLm)

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)       -0.00932    0.02310   -0.40  0.68672    
head(simArch1, -2) -0.10622    0.03153   -3.37  0.00078 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.73 on 996 degrees of freedom
Multiple R-squared:  0.0113,	Adjusted R-squared:  0.0103 
F-statistic: 11.3 on 1 and 996 DF,  p-value: 0.000784

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

Jedna uwaga techniczna co do użycia head(x, -k) i tail(x, -k). Ta funkcja działa nieintuicyjnie. Zwyczajne head(x, k) bierze pierwsze k wartości, a resztę wyrzuca. Wydawałoby się, że ujemne k będzie oznaczać wyrzucenie pierwszych k wartości i pozostawienie reszty. Jednak funkcja zachowuje się inaczej: ucina ostatnie k wartości. Wyjaśnienie może jest trudne, ale wydaje mi się sensowne: mamy tu wektor skierowany w danym kierunku, powiedzmy, że w lewo. Niech k = 5. Funkcja head(x, k=5) oznacza, że pierwszych pięć wartości wektora zostaje na miejscu, a resztę ucinamy. Można też powiedzieć, że te pięć wartości wysuwamy do przodu, a resztę zostawiamy. Te pięć pierwszych wartości tworzy nam nowy wektor. Teraz niech k = -5. W tej sytuacji kierunek wektora staje się przeciwny - czyli zostaje skierowany w prawo. Bierzemy więc znowu pięć pierwszych wartości, ale z punktu widzenia kierunku wektora będzie to ostatnich 5 wartości, które znowu wysuwamy, ale tym razem w prawo. Jednak head() zachowa się tak samo jak wcześniej - pierwsze wartości zostają zachowane. Funkcja więc bierze znowu od początku zachowane wartości, ale w pierwotnej kolejności. Ponieważ "wysunięte daleko w prawo" zostało pięć ostatnich, to są one niewidoczne dla head() i dlatego ucięte zostaje 5 ostatnich.   

Analogicznie dzieje się w sytuacji tail(x, -k), tylko w przeciwnym kierunku. Ponieważ tail oznacza czytanie od końca, to k oznacza wysunięcie k wartości z prawej strony, a więc -k wysunięcie k wartości od lewej strony (tj. od początku) i czytamy wektor tak jakby tych wysuniętych wartości nie było.

W sumie w head / tail są zawarte dwa rodzaje instrukcji: (1) wysuwanie k wartości w kierunku zależnym od znaku k, (2) czytanie wartości pozostałych po "wysunięciu". Oczywiście "wysunięcie" jest tylko moim sposobem na zrozumienie tej mało intuicyjnej funkcji, gdy znak k jest ujemny.

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

Właściwie jedynie bliski zera R^2 mówi, że realna korelacja nie istnieje. Niemniej standardowe testy na autokorelacje jak Box-Pierce czy Ljung-Box wprowadzają w błąd:

Box.test(x=simArch, lag=2, type = "Box-Pierce")

Box.test(x=simArch, lag=2, type = "Ljung-Box")

	Box-Pierce test

data:  simArch1
X-squared = 20, df = 2, p-value = 0.00004


	Box-Ljung test

data:  simArch1
X-squared = 20, df = 2, p-value = 0.00004

Warto dodać, że gdybyśmy chcieli użyć funkcji auto.arima z pakietu forecast, dostaniemy także błędny model mimo, że optymalny. Najlepiej w tej funkcji używać kryterium BIC, ponieważ najmocniej "karze" za złe lub nadmierne dopasowanie. Ponadto warto dodać parametry: stepwise = FALSE oraz approximation = FALSE, które domyślnie są ustawione na TRUE. Te dwa parametry zostały wprowadzone po to, by efektywniej znaleźć optymalne modele, ale kosztem jest mniejsza dokładność. Opcje te mają sens tylko dla bardzo długich szeregów czasowych lub w których występują sezonowości.

arma_test_bic <- auto.arima(simArch1, stepwise = FALSE, approximation = FALSE, ic="bic")

arma_test_bic 

ARIMA(0,0,2) with zero mean 

Coefficients:
         ma1     ma2
      -0.111  -0.110
s.e.   0.031   0.031

sigma^2 = 0.525:  log likelihood = -1096
AIC=2198   AICc=2198   BIC=2212

Jak widać algorytm wskazał błędny model MA(2).

B) Jeszcze silniejszy efekt możemy dostać, jeśli zwiększymy rząd ARCH do 2. Ustawmy alpha1 na 0.49 i alpha2 na 0.5:

set.seed(44)

myGarch2 <- ugarchspec(variance.model = list(model = "sGARCH", 

                                             garchOrder = c(2, 0)), 

                       distribution.model = "norm", mean.model = list(armaOrder = c(0, 0), 

                                                                     include.mean = FALSE), 

                       fixed.pars=list(omega=0.1, alpha1=0.49, alpha2=0.5))

simArchModel = ugarchpath(myGarch2, n.sim = 1000)

simArch2 <- simArchModel@path$seriesSim

Pacf(simArch2, plot=F)

Pacf(simArch2)

simArchLm <- lm(tail(simArch2, -2) ~ head(simArch2, -2))

plot(x=head(simArch2, -2), y=tail(simArch2, -2))

abline(simArchLm)

summary(simArchLm)

 Autokorelacja drugiego rzędu sięga niemal 30%

Coefficients:
                  Estimate Std. Error t value            Pr(>|t|)    
(Intercept)        -0.0373     0.0571   -0.65                0.51    
head(simArch2, -2)  -0.2853     0.0304   -9.38 <0.0000000000000002 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.8 on 996 degrees of freedom
Multiple R-squared:  0.0812,	Adjusted R-squared:  0.0803 
F-statistic: 88.1 on 1 and 996 DF,  p-value: <0.0000000000000002

C) Co by się stało, gdybyśmy usunęli efekt ARCH i zostawili tylko stałą wariancję? Ustawmy garchOrder na c(1, 0), ale alpha1 niech równa się prawie 0:

set.seed(44)

myGarch3 <- ugarchspec(variance.model = list(model = "sGARCH", 

                                             garchOrder = c(1, 0)), 

                       distribution.model = "norm", mean.model = list(armaOrder = c(0, 0), 

                                                                     include.mean = FALSE), 

                       fixed.pars=list(omega=0.1, alpha1=0.0000000001))

simArchModel = ugarchpath(myGarch3, n.sim = 1000)

simArch3 <- simArchModel@path$seriesSim

Pacf(simArch3, plot=F)

Pacf(simArch3)

simArchLm <- lm(tail(simArch3, -2) ~ head(simArch3, -2))

plot(x=head(simArch3, -2), y=tail(simArch3, -2))

abline(simArchLm)

summary(simArchLm)

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)
(Intercept)       -0.00818    0.01002   -0.82     0.41
head(simArch3, -2)  0.01541    0.03179    0.48     0.63

Residual standard error: 0.316 on 996 degrees of freedom
Multiple R-squared:  0.000236,	Adjusted R-squared:  -0.000768 
F-statistic: 0.235 on 1 and 996 DF,  p-value: 0.628

Dostajemy po prostu błądzenie losowe, dlatego testy odrzucają możliwość autokorelacji. Czyli dopiero po usunięciu efektu ARCH dostajemy prawidłowe wyniki.

Sprawdźmy jeszcze, jaki model auto.arima wskaże tym razem:

arma_test_bic <- auto.arima(simArch3, stepwise = FALSE, approximation = FALSE, ic="bic")

arma_test_bic 

ARIMA(0,0,0) with zero mean 

sigma^2 = 0.0999:  log likelihood = -267.1
AIC=536.3   AICc=536.3   BIC=541.2

Prawidłowo: brak autoregresji.

Z BIC jest ten problem, że czasami karze za mocno niedopasowanie, co prowadzi do błędu odwrotnego - mimo, że efekt ARMA istnieje, to odrzuca go (zob. ten wpis). Ponieważ AIC z kolei prowadzi do "przeuczenia", rozsądnym rozwiązaniem wydaje się posługiwanie HQC, jeżeli próba jest stosunkowo duża (por. zob. tu). Można to zrobić w samym rugarch, który zawiera funkcję autoarfima. Jest podobna do tej z forecast, ale bardziej rozbudowana pod niektórymi względami, ale jednocześnie czasochłonna. Aby znaleźć model dla HQC stosujemy np. kod:

autoarfima(data=simArch3, criterion = "HQIC", ar.max=5, ma.max=5, method = "partial")

Określenie "method" jest tu niefortunne, bo nie chodzi o wybór np. LSE czy MLE, ale po prostu o to czy szukać wszystkich możliwych kombinacji rzędów parametrów od zera do ustawionego max, a więc jeśli mamy max 3, to szuka tylko pierwszego i pozostałych nie, potem tylko drugiego, a pozostałych nie, potem tylko trzeciego, a pozostałych nie, potem pierwszego i drugiego, a nie trzeciego, potem pierwszego i trzeciego, a nie drugiego, potem drugiego i trzeciego, w końcu wszystkie trzy. Jest to method = "full". Jest to wyjątkowo czasochłonna metoda przeszukiwania, która wydaje mi się w większości przypadków zbędna. Jeżeli pierwszy rząd nie zawiera cząstkowej autokorelacji, a tylko drugi, to ustawiając method = "partial", parametr pierwszy po prostu (teoretycznie) wyniesie zero, a więc nie trzeba usuwać pierwszego rzędu z obliczeń. Jeżeli występuje tylko na trzecim rzędzie, to wg "partial"" powinniśmy dostać teoretyczne zera dla dwóch pierwszych rzędów. Tak więc tutaj algorytm oblicza parametr dla każdego pośredniego rzędu. 

Otrzymane statystyki z autoarfima:

*----------------------------------*
*          ARFIMA Model Fit        *
*----------------------------------*
Mean Model	: ARFIMA(1,0,0)
Distribution	: norm 

Optimal Parameters
------------------------------------
       Estimate  Std. Error  t value Pr(>|t|)
ar1   -0.053005    0.031698  -1.6722 0.094481
sigma  0.315631    0.007058  44.7214 0.000000

Robust Standard Errors:
       Estimate  Std. Error  t value Pr(>|t|)
ar1   -0.053005    0.031006  -1.7095  0.08736
sigma  0.315631    0.007138  44.2202  0.00000

LogLikelihood : -265.8 

Information Criteria
------------------------------------
                    
Akaike       0.53551
Bayes        0.54533
Shibata      0.53551
Hannan-Quinn 0.53924

Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
                          statistic p-value
Lag[1]                  0.000001256  0.9991
Lag[2*(p+q)+(p+q)-1][2] 0.082771917  1.0000
Lag[4*(p+q)+(p+q)-1][5] 0.882737349  0.9628

H0 : No serial correlation

Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
                        statistic p-value
Lag[1]                     0.3014  0.5830
Lag[2*(p+q)+(p+q)-1][2]    0.4997  0.6937
Lag[4*(p+q)+(p+q)-1][5]    1.5118  0.7370


ARCH LM Tests
------------------------------------
             Statistic DoF P-Value
ARCH Lag[2]     0.7098   2  0.7012
ARCH Lag[5]     2.3098   5  0.8048
ARCH Lag[10]    6.9816  10  0.7272

Nyblom stability test
------------------------------------
Joint Statistic:  0.3984
Individual Statistics:            
ar1   0.1721
sigma 0.2294

Asymptotic Critical Values (10% 5% 1%)
Joint Statistic:     	 0.61 0.749 1.07
Individual Statistic:	 0.35 0.47 0.75


$rank.matrix
   AR MA Mean ARFIMA     HQIC converged
1   1  0    0      0   0.5392         1

 

Najlepszym modelem okazał się AR(1). Jednocześnie jednak model jest nieistotny stat. (p-value = 0.095), co sumarycznie przynosi poprawny wynik. HQIC lepiej działa, jeśli wiemy, że jakaś zależność istnieje (podobnie AIC), dlatego jego użycie na losowych danych może przynosić niewłaściwie rezultaty. Najlepiej używać razem z BIC.

 

Powyższe statystyki z autoarfima zawierają oprócz samych parametrów również informacje o stabilności modelu, autokorelacjach reszt i wariancji reszt. Wszystkie te testy posiadają p-value > 0.05, co oznacza, że model wydaje się stabilny i nie zawiera autokorelacji reszt. Oczywiście wiemy, że tak właśnie jest, więc działają tutaj prawidłowo.  

Z tych testów warto poświęcić chwilę testowi Nybloma na stabilność. Joint Statistic: 0.4 dotyczy stabilności całego modelu. porównujemy z teoretycznymi wartościami dla kolejnych p-value, szczególnie patrzymy na wartość korespondującą z 5% - tutaj jest to 0.75. Ponieważ ta liczba jest większa od 0.35 , nie ma podstaw do odrzucenia hipotezy zerowej o stabilności modelu. Dodatkowa analiza dotyczy każdego współczynnika osobno. Patrzymy na Individual Statistic. Mamy tu dwa współczynniki: dla AR(1) 0.17 oraz dla wariancji 0.23. Wartości te porównujemy z wartością teoretyczną odpowiadającą 5% - wynosi 0.47. Ponieważ 0.47 > 0.17 i 0.23, to znaczy, że ani parametr AR(1) ani wariancja całego modelu nie zmienia się.

Widać więc, że test pozwala sprawdzić nie tylko stabilność średniej, ale także wariancji. 

To może częściowo zastąpić badanie struktury modelu za pomocą pakietu strucchange. Ograniczeniem strucchange - przynajmniej na ten moment - jest brak wyróżnienia i oddzielenia efektów GARCH. W końcu wykryta heteroskedastyczność może być spowodowana zarówno niestabilnością parametrów, jak i efektem GARCH.

Rozwiązanie problemu

Zestaw statystyk, który otrzymujemy w autoarfima, nasuwa trop do rozwiązania problemu. Rzeczywiście, jest to szybki sposób na ocenę autokorelacji. 

Ad 1) Była zmienna rnorm_test złożona z dwóch próbek losowych z rozkładu normalnego o innych średnich. Użyjmy dla niej autoarfima:

autoarfima(data=rnorm_test, criterion = "HQIC", ar.max=5, ma.max=5, method = "partial")

*----------------------------------*
*          ARFIMA Model Fit        *
*----------------------------------*
Mean Model	: ARFIMA(4,0,3)
Distribution	: norm 

Optimal Parameters
------------------------------------
       Estimate  Std. Error    t value Pr(>|t|)
mu     1.664246    0.000848    1961.68        0
ar1   -0.895300    0.000021  -43261.14        0
ar2    0.748268    0.000017   43452.27        0
ar3    1.032321    0.000024   42584.16        0
ar4    0.074309    0.000022    3395.11        0
ma1    0.956703    0.000017   56730.59        0
ma2   -0.612226    0.000006 -100923.73        0
ma3   -0.896438    0.000023  -39376.32        0
sigma  1.043292    0.004795     217.58        0

Robust Standard Errors:
       Estimate  Std. Error     t value Pr(>|t|)
mu     1.664246    0.001353    1229.637        0
ar1   -0.895300    0.000033  -26992.452        0
ar2    0.748268    0.000038   19613.716        0
ar3    1.032321    0.000059   17386.314        0
ar4    0.074309    0.000041    1829.927        0
ma1    0.956703    0.000027   34910.145        0
ma2   -0.612226    0.000006 -103777.853        0
ma3   -0.896438    0.000054  -16747.143        0
sigma  1.043292    0.027597      37.805        0

LogLikelihood : -1462 

Information Criteria
------------------------------------
                   
Akaike       2.9427
Bayes        2.9869
Shibata      2.9425
Hannan-Quinn 2.9595

Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
                         statistic p-value
Lag[1]                      0.6171  0.4321
Lag[2*(p+q)+(p+q)-1][20]    7.2453  1.0000
Lag[4*(p+q)+(p+q)-1][34]   15.9922  0.6699

H0 : No serial correlation

Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
                        statistic p-value
Lag[1]                      2.513 0.11294
Lag[2*(p+q)+(p+q)-1][2]     3.476 0.10432
Lag[4*(p+q)+(p+q)-1][5]     6.116 0.08444


ARCH LM Tests
------------------------------------
             Statistic DoF P-Value
ARCH Lag[2]      4.208   2 0.12197
ARCH Lag[5]      8.450   5 0.13313
ARCH Lag[10]    16.301  10 0.09133

Nyblom stability test
------------------------------------
Joint Statistic:  9.36
Individual Statistics:             
mu    0.07491
ar1   0.01170
ar2   0.01451
ar3   0.04257
ar4   0.06617
ma1   0.07728
ma2   0.07368
ma3   0.05045
sigma 0.24675

Asymptotic Critical Values (10% 5% 1%)
Joint Statistic:     	 2.1 2.32 2.82
Individual Statistic:	 0.35 0.47 0.75


Elapsed time : 1.476 


$rank.matrix
   AR MA Mean ARFIMA    HQIC converged
1   4  3    1      0   2.959         1

Model wskazał błędnie ARMA(4, 3), w dodatku wszystkie parametry istotne. Popatrzmy jednak na test stabilności. 9.36 > 2.32, czyli testowa przekracza wielokrotnie wartość krytyczną. Wniosek: model jest na pewno niestabilny. Test dla osobnych współczynników nie przekracza teoretycznych, co dowodzi, że ten dodatkowy test nie jest bardziej szczegółową wersją tego pierwszego. Dlaczego dla osobnych nie ma niestabilności? Zauważmy, że ta niestabilność dotyczy w rzeczywistości jedynie stałej, czyli mu. Nie wiem dlaczego test nie wskazuje tego, ale na pewno ma to związek z tym, że model został wybrany błędnie z dużą liczbą parametrów. Aby się o tym przekonać, wystarczy, że ustawimy ar.max=1 i ma.max=0:

 autoarfima(data=rnorm_test, criterion = "HQIC", ar.max=1, ma.max=0, method = "partial")

*----------------------------------*
*          ARFIMA Model Fit        *
*----------------------------------*
Mean Model	: ARFIMA(1,0,0)
Distribution	: norm 

Optimal Parameters
------------------------------------
       Estimate  Std. Error  t value Pr(>|t|)
mu      2.46629    0.131594   18.742        0
ar1     0.68144    0.023154   29.430        0
sigma   1.33144    0.029772   44.721        0

Robust Standard Errors:
       Estimate  Std. Error  t value Pr(>|t|)
mu      2.46629    0.170801   14.440        0
ar1     0.68144    0.015803   43.121        0
sigma   1.33144    0.029462   45.191        0

LogLikelihood : -1705 

Information Criteria
------------------------------------
                   
Akaike       3.4164
Bayes        3.4311
Shibata      3.4164
Hannan-Quinn 3.4220

Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
                        statistic p-value
Lag[1]                       90.5       0
Lag[2*(p+q)+(p+q)-1][2]     104.1       0
Lag[4*(p+q)+(p+q)-1][5]     132.9       0

H0 : No serial correlation

Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
                        statistic p-value
Lag[1]                      4.450 0.03489
Lag[2*(p+q)+(p+q)-1][2]     4.473 0.05628
Lag[4*(p+q)+(p+q)-1][5]     5.004 0.15278


ARCH LM Tests
------------------------------------
             Statistic DoF P-Value
ARCH Lag[2]      4.540   2  0.1033
ARCH Lag[5]      5.816   5  0.3246
ARCH Lag[10]    13.804  10  0.1821

Nyblom stability test
------------------------------------
Joint Statistic:  11.13
Individual Statistics:              
mu    10.89280
ar1    0.09111
sigma  0.14545

Asymptotic Critical Values (10% 5% 1%)
Joint Statistic:     	 0.846 1.01 1.35
Individual Statistic:	 0.35 0.47 0.75


Elapsed time : 0.02069 


$rank.matrix
  AR MA Mean ARFIMA  HQIC converged
1  1  0    1      0 3.422         1

Tym razem zarówno cały model, jak i stała mu zostały prawidłowo uznane za niestabilne. A skoro tak, to nie można traktować poważnie takich autokorelacji. Na marginesie - efekt ARCH nie występuje, na co wskazuje test ARCH LM (p-value>0.05).

Ad 2A) W drugim przykładzie mieliśmy wprowadzony ARCH do wariancji. Ponownie możemy sprawdzić go za pomocą autoarfima

autoarfima(data=simArch, criterion = "HQIC", ar.max=5, ma.max=5, method = "partial", distribution.model="norm")

*----------------------------------*
*          ARFIMA Model Fit        *
*----------------------------------*
Mean Model	: ARFIMA(0,0,2)
Distribution	: norm 

Optimal Parameters
------------------------------------
       Estimate  Std. Error  t value Pr(>|t|)
ma1    -0.11131    0.031434  -3.5412 0.000398
ma2    -0.10991    0.031448  -3.4950 0.000474
sigma   0.72384    0.016184  44.7245 0.000000

Robust Standard Errors:
       Estimate  Std. Error  t value Pr(>|t|)
ma1    -0.11131    0.080460  -1.3835  0.16653
ma2    -0.10991    0.107543  -1.0220  0.30677
sigma   0.72384    0.083658   8.6524  0.00000

LogLikelihood : -1096 

Information Criteria
------------------------------------
                   
Akaike       2.1976
Bayes        2.2123
Shibata      2.1976
Hannan-Quinn 2.2032

Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
                        statistic p-value
Lag[1]                   0.001732  0.9668
Lag[2*(p+q)+(p+q)-1][5]  0.273995  1.0000
Lag[4*(p+q)+(p+q)-1][9]  2.276878  0.9695

H0 : No serial correlation

Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
                        statistic p-value
Lag[1]                      311.0       0
Lag[2*(p+q)+(p+q)-1][2]     423.5       0
Lag[4*(p+q)+(p+q)-1][5]     585.5       0


ARCH LM Tests
------------------------------------
             Statistic DoF P-Value
ARCH Lag[2]      348.1   2       0
ARCH Lag[5]      373.2   5       0
ARCH Lag[10]     380.7  10       0

Nyblom stability test
------------------------------------
Joint Statistic:  1.714
Individual Statistics:            
ma1   0.1458
ma2   0.9028
sigma 1.1342

Asymptotic Critical Values (10% 5% 1%)
Joint Statistic:     	 0.846 1.01 1.35
Individual Statistic:	 0.35 0.47 0.75



$rank.matrix
   AR MA Mean ARFIMA  HQIC converged
1   0  2    0      0 2.203         1

 

Algorytm wybrał MA(2), co koresponduje z ujemną autokorelacją cząstkową dwóch rzędów. Ale znowu jak spojrzymy na test Nybloma, okaże się, że jest to niestabilny model. W dodatku test ARCH LM informuje o występowaniu efektu ARCH w składnikach losowych, tak samo test Weighted Ljung-Boxa dla wariancji reszt ma p-value = 0. Natomiast test Weighted Ljung-Boxa dla samych reszt daje p-value dużo powyżej 0.05, co wynika z faktu, że wszystkie autokorelacje zostały uwzględnione w modelu MA(2).

Mamy więc podstawy do stworzenia modelu ARCH. Robimy to dwóch krokach.

Krok 1. Budujemy najprostszy model teoretyczny, ale tym razem bez fixed.pars:

myGarch1 <- ugarchspec(variance.model = list(model = "sGARCH", 

                                             garchOrder = c(1, 0)), 

                       distribution.model = "norm", mean.model = list(armaOrder = c(0, 0)))

Krok 2. Znajdujemy model empiryczny: bierzemy model (myGarch1) z pierwszego kroku i wstawiamy symulację procesu ARCH z p. 2A (simArch1)

ugarchfit(spec = myGarch1, data = simArch1)

*---------------------------------*
*          GARCH Model Fit        *
*---------------------------------*

Conditional Variance Dynamics 	
-----------------------------------
GARCH Model	: sGARCH(1,0)
Mean Model	: ARFIMA(0,0,0)
Distribution	: norm 

Optimal Parameters
------------------------------------
        Estimate  Std. Error  t value Pr(>|t|)
mu     -0.006221    0.011653 -0.53384  0.59345
omega   0.109187    0.008094 13.49034  0.00000
alpha1  0.803593    0.071384 11.25730  0.00000

Robust Standard Errors:
        Estimate  Std. Error  t value Pr(>|t|)
mu     -0.006221    0.011105 -0.56021  0.57534
omega   0.109187    0.008383 13.02461  0.00000
alpha1  0.803593    0.073067 10.99798  0.00000

LogLikelihood : -735.2 

Information Criteria
------------------------------------
                   
Akaike       1.4765
Bayes        1.4912
Shibata      1.4764
Hannan-Quinn 1.4820

Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
                        statistic p-value
Lag[1]                      3.057 0.08039
Lag[2*(p+q)+(p+q)-1][2]     3.169 0.12628
Lag[4*(p+q)+(p+q)-1][5]     3.974 0.25734
d.o.f=0
H0 : No serial correlation

Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
                        statistic p-value
Lag[1]                     0.2462  0.6197
Lag[2*(p+q)+(p+q)-1][2]    0.9056  0.5299
Lag[4*(p+q)+(p+q)-1][5]    1.9559  0.6289
d.o.f=1

Weighted ARCH LM Tests
------------------------------------
            Statistic Shape Scale P-Value
ARCH Lag[2]     1.313 0.500 2.000  0.2518
ARCH Lag[4]     2.120 1.397 1.611  0.4142
ARCH Lag[6]     2.382 2.222 1.500  0.5961

Nyblom stability test
------------------------------------
Joint Statistic:  0.4109
Individual Statistics:              
mu     0.06427
omega  0.14734
alpha1 0.27249

Asymptotic Critical Values (10% 5% 1%)
Joint Statistic:     	 0.846 1.01 1.35
Individual Statistic:	 0.35 0.47 0.75

Sign Bias Test
------------------------------------
                   t-value   prob sig
Sign Bias           1.1690 0.2427    
Negative Sign Bias  0.2163 0.8288    
Positive Sign Bias  0.6087 0.5429    
Joint Effect        1.4816 0.6865    


Adjusted Pearson Goodness-of-Fit Test:
------------------------------------
  group statistic p-value(g-1)
1    20     11.28       0.9141
2    30     15.32       0.9824
3    40     23.12       0.9796
4    50     32.20       0.9695

Przypomnę, że simArch1 miał ustawioną średnią wariancję (omega) na 0.1 oraz parametr ARCH (alpha1) na 0.9. Funkcja poprawnie dopasowała omegę, a nieco gorzej alpha1, ale oba są istotne. Teraz dokonujemy diagnostyki:

- autokorelacje w resztach modelu: nie występują (Ljung-Box Test on Standardized Residuals, p-value>0.05),

- autokorelacje w wariancji reszt modelu (Ljung-Box Test on Standardized Squared Residuals): nie występują, co oznacza, że efekt GARCH nie występuje w szeregu czasowym,

- Efekt ARCH w resztach modelu (ARCH LM Test): nie występuje, co oznacza, że uwzględniliśmy cały efekt ARCH w modelu

- stabilność (Nyblom stability test): model jest stabilny, ponieważ żadna statystyka nie przekracza wartości krytycznych,

- czy występują szoki nieuwzględnione przez model: nie występują (Sign Bias Test, p-value > 0.05)

- czy przyjęty rozkład reszt jest odpowiedni: rozkład normalny jest odpowiedni (Goodness-of-Fit Test, p>0.05 dla wszystkich kwantyli).

Powyższa diagnostyka dowodzi, że model jest poprawny.

Została ostatnia sprawa do sprawdzenia. Na samym początku robiliśmy test, z którego wynikało, że parametr AR(1) jest istotny statystycznie, choć wiedzieliśmy od początku, że żadnej autokorelacji nie ma. Ale teraz, kiedy mamy testy uwzględniające efekt ARCH, możemy sprawdzić jaka będzie ocena modelu, który zakładałby, że istnieje nie tylko ARCH(1), ale i AR(1). Wpisujemy więc armaOrder = c(1, 0):

myGarch1 <- ugarchspec(variance.model = list(model = "sGARCH", 

                                             garchOrder = c(1, 0), 

                                             submodel = NULL), 

                       distribution.model = "norm", mean.model = list(armaOrder = c(1, 0)))

ugarchfit(spec = myGarch1, data = simArch1)

*---------------------------------*
*          GARCH Model Fit        *
*---------------------------------*

Conditional Variance Dynamics 	
-----------------------------------
GARCH Model	: sGARCH(1,0)
Mean Model	: ARFIMA(1,0,0)
Distribution	: norm 

Optimal Parameters
------------------------------------
        Estimate  Std. Error  t value Pr(>|t|)
mu     -0.006138    0.011482 -0.53458  0.59294
ar1    -0.024017    0.028336 -0.84758  0.39667
omega   0.109324    0.008095 13.50565  0.00000
alpha1  0.804331    0.071556 11.24056  0.00000

Robust Standard Errors:
        Estimate  Std. Error  t value Pr(>|t|)
mu     -0.006138    0.011090 -0.55346  0.57995
ar1    -0.024017    0.031426 -0.76424  0.44472
omega   0.109324    0.008389 13.03140  0.00000
alpha1  0.804331    0.072694 11.06468  0.00000

LogLikelihood : -734.9 

Information Criteria
------------------------------------
                   
Akaike       1.4777
Bayes        1.4974
Shibata      1.4777
Hannan-Quinn 1.4852

Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
                        statistic p-value
Lag[1]                      1.457  0.2273
Lag[2*(p+q)+(p+q)-1][2]     1.552  0.4090
Lag[4*(p+q)+(p+q)-1][5]     2.389  0.5982
d.o.f=1
H0 : No serial correlation

Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
                        statistic p-value
Lag[1]                     0.1987  0.6558
Lag[2*(p+q)+(p+q)-1][2]    1.0633  0.4781
Lag[4*(p+q)+(p+q)-1][5]    2.1068  0.5934
d.o.f=1

Weighted ARCH LM Tests
------------------------------------
            Statistic Shape Scale P-Value
ARCH Lag[2]     1.722 0.500 2.000  0.1894
ARCH Lag[4]     2.373 1.397 1.611  0.3640
ARCH Lag[6]     2.579 2.222 1.500  0.5545

Nyblom stability test
------------------------------------
Joint Statistic:  0.7708
Individual Statistics:              
mu     0.06371
ar1    0.40036
omega  0.15814
alpha1 0.27546

Asymptotic Critical Values (10% 5% 1%)
Joint Statistic:     	 1.07 1.24 1.6
Individual Statistic:	 0.35 0.47 0.75

Sign Bias Test
------------------------------------
                   t-value   prob sig
Sign Bias           1.1068 0.2687    
Negative Sign Bias  0.1919 0.8478    
Positive Sign Bias  0.5713 0.5680    
Joint Effect        1.3374 0.7203    


Adjusted Pearson Goodness-of-Fit Test:
------------------------------------
  group statistic p-value(g-1)
1    20      9.28       0.9685
2    30     22.70       0.7901
3    40     21.44       0.9899
4    50     26.40       0.9966

Skąd wiadomo, że ten drugi model jest gorszy? Po pierwsze widzimy, że ar1 jest nieistotny stat., ale po drugie wartości wszystkich kryteriów (AIC, BIC, HQC, SIC) zwiększyły się, co oznacza, że model się pogorszył. W sumie więc dostajemy informację, że dodanie autokorelacji obserwacji pogarsza jedynie model, a to znaczy, że autokorelacje faktycznie nie występują. 

Ad 2B) W przykładzie 2B robiłem symulację ARCH dwóch rzędów. Co by się stało, gdybym użył jednak do dopasowania modelu dla jednego rzędu, czyli z przykładu 2A? Do funkcji ugarchfit wstawię więc myGarch1, ale symulację simArch2:

 ugarchfit(spec = myGarch1, data = simArch2)

*---------------------------------*
*          GARCH Model Fit        *
*---------------------------------*

Conditional Variance Dynamics 	
-----------------------------------
GARCH Model	: sGARCH(1,0)
Mean Model	: ARFIMA(0,0,0)
Distribution	: norm 

Optimal Parameters
------------------------------------
        Estimate  Std. Error  t value Pr(>|t|)
omega    0.35066    0.023680  14.8085        0
alpha1   0.80988    0.091434   8.8575        0

Robust Standard Errors:
        Estimate  Std. Error  t value Pr(>|t|)
omega    0.35066     0.11587   3.0264 0.002475
alpha1   0.80988     0.19997   4.0500 0.000051

LogLikelihood : -1204 

Information Criteria
------------------------------------
                   
Akaike       2.4114
Bayes        2.4212
Shibata      2.4114
Hannan-Quinn 2.4151

Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
                        statistic   p-value
Lag[1]                    0.06781 0.7945595
Lag[2*(p+q)+(p+q)-1][2]   6.96722 0.0121615
Lag[4*(p+q)+(p+q)-1][5]  15.92422 0.0002506
d.o.f=0
H0 : No serial correlation

Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
                        statistic            p-value
Lag[1]                    0.01798 0.8933177180472207
Lag[2*(p+q)+(p+q)-1][2]  48.79191 0.0000000000001286
Lag[4*(p+q)+(p+q)-1][5]  79.48513 0.0000000000000000
d.o.f=1

Weighted ARCH LM Tests
------------------------------------
            Statistic Shape Scale P-Value
ARCH Lag[2]     97.16 0.500 2.000       0
ARCH Lag[4]     98.80 1.397 1.611       0
ARCH Lag[6]     99.68 2.222 1.500       0

Nyblom stability test
------------------------------------
Joint Statistic:  0.4387
Individual Statistics:             
omega  0.2553
alpha1 0.2381

Asymptotic Critical Values (10% 5% 1%)
Joint Statistic:     	 0.61 0.749 1.07
Individual Statistic:	 0.35 0.47 0.75

Sign Bias Test
------------------------------------
                   t-value   prob sig
Sign Bias           0.3514 0.7254    
Negative Sign Bias  0.1107 0.9118    
Positive Sign Bias  0.3619 0.7175    
Joint Effect        0.3531 0.9497    


Adjusted Pearson Goodness-of-Fit Test:
------------------------------------
  group statistic      p-value(g-1)
1    20     97.32 0.000000000001632
2    30     98.24 0.000000001867482
3    40    114.00 0.000000002827673
4    50    105.30 0.000005407768343

Ocena modelu rzeczywiście jest zła: występują autokorelacje w resztach i wariancji, efekt ARCH nadal się pojawia, a rozkład reszt nie pasuje. Musismy więc rzeczywiście poprawić model do ARCH dwóch rzędów:

myGarch2 <- ugarchspec(variance.model = list(model = "sGARCH", 

                                             garchOrder = c(2, 0)), 

                       distribution.model = "norm", mean.model = list(armaOrder = c(0, 0)))

ugarchfit(spec = myGarch2, data = simArch2)

*---------------------------------*
*          GARCH Model Fit        *
*---------------------------------*

Conditional Variance Dynamics 	
-----------------------------------
GARCH Model	: sGARCH(2,0)
Mean Model	: ARFIMA(0,0,0)
Distribution	: norm 

Optimal Parameters
------------------------------------
        Estimate  Std. Error  t value Pr(>|t|)
mu     -0.001293    0.013566  -0.0953  0.92408
omega   0.110512    0.011278   9.7993  0.00000
alpha1  0.384594    0.055023   6.9897  0.00000
alpha2  0.541667    0.062540   8.6611  0.00000

Robust Standard Errors:
        Estimate  Std. Error   t value Pr(>|t|)
mu     -0.001293    0.012938 -0.099924  0.92041
omega   0.110512    0.010805 10.227909  0.00000
alpha1  0.384594    0.044567  8.629512  0.00000
alpha2  0.541667    0.074556  7.265275  0.00000

LogLikelihood : -963 

Information Criteria
------------------------------------
                   
Akaike       1.9340
Bayes        1.9536
Shibata      1.9339
Hannan-Quinn 1.9414

Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
                        statistic p-value
Lag[1]                      3.253  0.0713
Lag[2*(p+q)+(p+q)-1][2]     3.341  0.1134
Lag[4*(p+q)+(p+q)-1][5]     4.082  0.2440
d.o.f=0
H0 : No serial correlation

Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
                        statistic p-value
Lag[1]                      1.329  0.2491
Lag[2*(p+q)+(p+q)-1][5]     2.025  0.6126
Lag[4*(p+q)+(p+q)-1][9]     3.470  0.6795
d.o.f=2

Weighted ARCH LM Tests
------------------------------------
            Statistic Shape Scale P-Value
ARCH Lag[3]    0.4403 0.500 2.000  0.5070
ARCH Lag[5]    0.9204 1.440 1.667  0.7567
ARCH Lag[7]    2.2697 2.315 1.543  0.6603

Nyblom stability test
------------------------------------
Joint Statistic:  0.4937
Individual Statistics:             
mu     0.1045
omega  0.2665
alpha1 0.1167
alpha2 0.1563

Asymptotic Critical Values (10% 5% 1%)
Joint Statistic:     	 1.07 1.24 1.6
Individual Statistic:	 0.35 0.47 0.75

Sign Bias Test
------------------------------------
                   t-value   prob sig
Sign Bias           0.6479 0.5172    
Negative Sign Bias  0.5606 0.5752    
Positive Sign Bias  0.1738 0.8621    
Joint Effect        1.0146 0.7977    


Adjusted Pearson Goodness-of-Fit Test:
------------------------------------
  group statistic p-value(g-1)
1    20     11.72       0.8972
2    30     22.46       0.8007
3    40     25.68       0.9502
4    50     36.60       0.9046

Nastąpiła pełna poprawa. Podobnie jak poprzednio można porównać model dopisując ar1. Zobaczymy wtedy, że wszystkie kryteria wskazują na pogorszenie, a więc że pierwotna autokorelacja okazała się fałszem (już tego nie pokazuję).

Podsumowanie

Podsumowując, zarówno samo sprawdzanie autokorelacji, jak i konstrukcja modelu autoregresji bez sprawdzenia stabilności i efektów ARCH, jest błędne. Autokorelacja może się okazać złudzeniem.  Należy wykonać wszystkie testy i na podstawie kryteriów, takich jak BIC czy HQC, poszukać najlepszego modelu.    

Literatura:

[1] Bera, A.K., M.L. Higgins and S. Lee (1993) "Interaction Between Autocorrelation and Conditional Heteroskedasticity: A Random Coefficients Approach", Journal of Business and Economic Statistics, 10, 133-142.