Krótkookresowa prognoza WIG - na marzec 2018

WIG zaliczył w lutym solidny spadek wartości - niecałe 7%. Jest to powyżej odchylenia standardowego miesięcznej stopy zwrotu 5,8% (w okresie 2006-2017). Może to oznaczać, że indeks powróci do średniej w marcu i kwietniu - są to zresztą statystycznie dobre miesiące dla giełdy. Przeprowadziłem krótką analizę w gretlu w oparciu X-13-ARIMA (zob. Usuwanie sezowości i prognozowanie przy pomocy X-13-ARIMA i TRAMO/SEATS) oraz ARIMAX z inflacją w przemyśle jako zmienną egzogeniczną.

1) X-13-ARIMA dla miesięcznej stopy zwrotu WIG: automatycznie uzyskany trend wygląda następująco:

Wykres sugeruje, że stopa WIG znajduje się poniżej trendu, a sam trend wyprostował się. Dokładna prognoza wg X-13-ARIMA na marzec 2018 to +0,13%.

2) ARMAX dla miesięcznej stopy zwrotu WIG z inflacją w przemyśle m/m jako zmienną egzogeniczną: aby wykonać taką prognozę, należy wykonać trochę więcej obliczeń. Najpierw bierzemy inflację w przemyśle z GUSu (dla okresu 01.2006-01.2018). Dlaczego akurat inflacja w przemyśle? Okazuje się ona ujemnie skorelowana ze stopą zwrotu WIG. Od 2006 r. korelacja między miesięcznym tempem WIG a inflacją przemysłu m/m wyniosła -0,26 (por. z art. Czy inflacja jest dobra dla akcji?). Jest to korelacja istotna stat. na p = 0,01. Oznacza to, że WIG zachowuje się średnio biorąc odwrotnie do inflacji w przemyśle.

Następnie także za pomocą dodatku armax szukam optymalnego ARMA dla samej inflacji przemysłu. AIC wskazał ARMA(2,3), a HQC ARMA(0,1). Który bierzemy pod uwagę? Wg Venus Liew HQC przewyższa wszystkie kryteria, gdy próba jest stosunkowo duża i wynosi co najmniej 120 (dla modelu ARMA) [1]. Przy mniejszej próbie lepszy jest AIC. W naszym przypadku jest ponad 140 danych, czyli korzystamy z HQC. Stąd można wyznaczyć prognozę tej inflacji na luty, marzec:

Teraz możemy wstawić prognozę z lutego i marca do danych, a dodatkowo ustawić zakres próby na 01.2006 - 02.2018, po to, aby gretl uznał marzec za prognozę.

Przedostatnim krokiem będzie optymalizacja modelu stopy zwrotu WIG za pomocą dodatku armax. Poszukujemy modelu ARMAX(p, q) z inflacją przemysłu w roli zmiennej egzogenicznej. Dodatek wskazał optimum dla ARMAX(6, 6) - zarówno dla AIC jak i HQC. Uzyskujemy parametry:

 Ostatni krok to prognoza na marzec:

Stopa w marcu może być więc ciągle ujemna. Ale gdybyśmy wybrali model tylko minimalnie gorszy ARMAX(8, 2):

gdzie HQC jest praktycznie identyczny, to dostalibyśmy prognozę w marcu na plus:

Literatura:
[1] Liew, Venus Khim−Sen,  Which Lag Length Selection Criteria Should We Employ?, May 2004.

Usuwanie sezowości i prognozowanie przy pomocy X-13-ARIMA i TRAMO/SEATS

Różne branże, np. budownictwo, charakteryzuje sezonowość, która zniekształca obraz ich sytuacji finansowej. Jeśli dokonujemy analizy fundamentalnej albo sporządzamy prognozę ekonometryczną, to powinniśmy usunąć zmienność kalendarzową. Nazywa się to wyrównaniem sezonowym. Weźmy GUS-owski Wskaźnik ogólnego klimatu koniunktury (który jest średnią arytmetyczną z odpowiedzi na pytania dyrektorów przedsiębiorstw o aktualną oraz przewidywaną ogólną sytuację gospodarczą danego przedsiębiorstwa), który jest aktualizowany co miesiąc. Dla budownictwa w okresie 1.2001-7.2017 wygląda on tak:

Wskaźnik ten wydaje się być podany jako niewyrównany sezonowo, co nawet GUS zaznacza w tabeli. Ale zarówno w "Uwagach ogólnych" jak i "Wyjaśnieniach metodycznych" GUS stwierdza, że dane dotyczące wskaźnika koniunktury zostały wyrównane sezonowo metodą TRAMO/SEATS. Skąd taka niekonsekwencja? Nie wiem. Powinny być to dane niewyrównane, bo występuje ewidentna cykliczność. Np. gdy powiększymy ten wykres dla lat 2014-2017 to zobaczymy idealnie ten sam wzór:

Obecnie mamy do czynienia ze szczytem cyklu, tzn. że trzeba spodziewać się pogorszenia nastrojów koniunkturalnych. Czy to oznacza, że giełda powinna reagować spadkiem? Oczywiście nie, bo efektywny rynek powinien dyskontować takie informacje o sezonowości, a więc racjonalny inwestor użyje narzędzia do wyrównania sezonowego.

Są dwie znane metody do tego: pierwsza to X-12-ARIMA, która niedawno została zastąpiona przez X-13-ARIMA-SEATS, a druga to TRAMO/SEATS. Obydwie bardzo łatwo użyć w Gretlu. Tutaj wyjaśniam skąd takie nazwy i co te programy robią, a także pokazuję jak wyciągnąć z nich prognozy.

Nie chcę się skupiać na teoretycznych aspektach tych metod, przejdę więc od razu do rzeczy, czyli do praktycznego zastosowania obydwu metod. Pakiety X-13-ARIMA i TRAMO/SEATS nie są automatycznie instalowane razem z Gretlem, dlatego musimy na stronie Gretla je najpierw ściągnąć i zainstalować, co nie jest żadnym problemem, bo po instalacji powinniśmy mieć te pakiety w Gretlu. Czasami mogą wystąpić problemy ze znalezieniem ścieżki. Żeby sprawdzić czy Gretl poprawnie łączy się z programem, wchodzimy w gretlu w: Narzędzia -> Ustawienia -> Ogólne -> Programy. Wyszukujemy tutaj "Ścieżka do programu X-12-ARIMA" (przynajmniej u mnie tak jest, że nie zaktualizowano nazwy na X-13-ARIMA) i "Ścieżka do programu TRAMO". Jeśli wyświetlona obok ścieżka jest błędna, to musimy ją poprawić.

Aby użyć metody, musimy mieć zaznaczoną naszą zmienną w oknie gretla,a sama zmienna musi być szeregiem czasowym o określonej wcześniej częstości, tutaj 12 miesięcy (jednak nie każda częstość jest dozwolona). Same programy działają po wejściu w Zmienna -> Analiza ... i w zależności od typu wybieramy program.

X-13-ARIMA-SEATS:
Przy standardowych opcjach, klikamy OK; wówczas GUS-owski Wskaźnik koniunktury budownictwa zostanie podzielony na 3 części: wyrównany sezonowo, trend/cykl oraz odchylenia przypadkowe:

Wyrównanie sezonowe nie jest filtrowaniem z odchyleń przypadkowych, a jedynie ze zmian sezonowych. Mówiąc prosto: zmienna 'wyrównany sezonowo' to zmienna złożona z dwóch części: 'trendu' oraz zmian przypadkowych.
Wynika z tego, że filtrowanie z danych sezonowych oraz ze zmian przypadkowych zawiera trend. Dlatego inwestor będzie najbardziej zainteresowany drugim wykresem, który wskazuje w jakiej fazie cyklu koniunkturalnego rzeczywiście się znajdujemy.

Oprócz wykresu dostajemy także bardzo duży zestaw statystyk, których ilość może zniechęcać do ich analizy. Ogólnie to interesuje nas tylko wybór modelu, jego parametry i istotność statystyczna. W zestawie otrzymanych danych szukamy terminu "Final automatic model choice". W tym przypadku dostajemy coś takiego:

Final automatic model choice : (0 1 1)(0 1 1)

Co oznacza, po kolei w 1 nawiasie:
AR(0) - czyli nie występuje autoregresja liniowa w oryginalnych danych (z miesiąca na miesiąc)
d = 1 - czyli oryginalne dane są niestacjonarne i trzeba je zamienić na zmiany pierwszego rzędu
MA(1) - występuje średnia ruchoma pierwszego rzędu (z miesiąca na miesiąc) w oryginalnych danych.

W 2 nawiasie:
SAR(0) - nie występuje sezonowa autoregresja
Sd = 1 - oryginalne dane są niestacjonarne i trzeba je zamienić na zmiany 12-miesięczne
SMA(1) - występuje sezonowa średnia ruchoma pierwszego rzędu w oryginalnych danych (nie mylić z Simple Moving Average z AT)

Trzeba pamiętać, że AR i MA odnoszą się do zmian z okresu na okres (tutaj z miesiąca na miesiąc). Natomiast SAR i SMA odnoszą się do zmian danego miesiąca w stosunku do zmian tego samego miesiąca z poprzedniego roku.

Istotność statystyczną sezonowości oceniamy szukając wyników do 3 testów:
- test na stabilność sezonowości: Test for the presence of seasonality assuming stability
Wynik:
Seasonality present at the 0.1 per cent level.
tj. sezonowość występuje poziomie istotności 0,1%.

- nieparametryczny test na stabilność sezonowości:Nonparametric Test for the Presence of Seasonality Assuming Stability.
Wynik:
Seasonality present at the one percent level.

- test na zmienność amplitudy sezonowości: Test na Moving Seasonality Test
Wynik:
Moving seasonality present at the one percent level
tj. wykryta sezonowość zmienia swoją wielkość (amplitudę) z okresu na okres przy 1% istotności.
Ta własność nie występowałaby gdyby np. w danym miesiącu zmiana byłaby zawsze taka sama.

Jak powiedziałem inwestor jest zainteresowany głównie trendem, a więc także jego prognozą. Jeżeli na początku zaznaczymy opcję 'Zapisywanie danych' - składnik trendowy, to możemy wykorzystać sam trend jako zmienną do prognozy. Jakie jednak przyjąć dla niej parametry? Skoro X-13-ARIMA poszukuje optymalnego modelu, to moglibyśmy ponownie wykonać całą procedurę już dla trendu - spodziewając się, że wskaże model bez sezonowości, np. (1, 1, 1)(0, 0, 0). Jak najbardziej jest to możliwe, jednak okazuje się, że gdy to zrobiłem, to dostałem:
Final automatic model choice : (3 2 1)(1 0 1) ,
co sugeruje, że sezonowość występuje także w trendzie!

Aby sprawdzić prognozę trendu, możemy wykorzystać zwykłą procedurę ARIMA (Model -> Modele szeregów czasowych -> model ARIMA) i zaznaczyć AR=3, d=2, MA=1, SAR=1, Sd=0, MA=1.

Wyniki:

Następnie sprawdzamy prognozę (Analiza -> Prognoza). Przy wyborze 12 miesięcy do przodu, można spodziewać się ożywienia w budownictwie:

Metoda ta ma ograniczenia w postaci maksymalnych rzędów części AR i MA. Jest też inna metoda. Instalujemy dodatek o nazwie armax (Narzędzia -> pakiety funkcji -> na serwerze gretla, tam wybieramy 'armax'). Dodatek ten jest fajny, bo szuka najlepszego możliwego modelu ARMAX (literka X oznacza tylko, że możemy dodać zmienną zewnętrzną do modelu) przy wybranym maksymalnym rzędzie AR i MA (nie ma tutaj ograniczeń). Ze względu na miesięczną częstotliwość, wybrałem max rząd 12 dla obydwu. I tu trzeba się zastanowić czy tak wysoki pułap chcemy sprawdzić, bo będziemy długo czekać zanim program przeliczy parametry AR i MA od 0 do 12. Jeśli zależy nam na tej liczbie, a chcemy przyspieszyć działanie, możemy zmienić metodę estymacji. Program używa dwóch metod: standardowa to dokładna metoda największej wiarogodności (exact ML), druga to metoda warunkowej wiarogodności (conditional ML).

Która metoda jest lepsza? Obydwa estymatory są zgodne, ale exact ML jest efektywniejsza od conditional ML w tym sensie, że jej wariancja jest najmniejsza [1], tzn. także lepsza od uogólnionej metody najmniejszych kwadratów - w sytuacji gdy występuje proces ARMA [2]. Ta druga jest jednak obliczeniowo bardzo szybka przy dużej liczbie parametrów.

Pamiętać należy o dwóch rzeczach: po pierwsze dodatek armax nie sprawdza sezonowości, co w naszym przypadku może być zaletą, bo już ją usunęliśmy; po drugie jeśli wiemy już, że nasz szereg jest niestacjonarną zmienną, to musimy przekształcić je w różnice (Dodawanie zmiennych -> pierwsze różnice), bo armax bada tylko model ARMAX, a nie ARIMAX.

Optimum zostaje osiągnięte, gdy kryteria AIC, BIC lub HQC osiągną najmniejszą wartość (tak jak entropia - im mniejsza, tym mniejsza niepewność co do modelu). Tutaj znów powstaje pytanie, które kryterium wybrać? Najczęściej wybierane jest AIC, jednak nie jest to wybór absolutny, co więcej nie jest optymalny w przypadku procesów ARMA. Hannan i Quinn [3] pokazali empirycznie, że ich kryterium, które potem zostało nazwane od ich nazwisk HQC, jest prawdopodobnie lepsze gdy proces jest autoregresyjny, a próba jest stosunkowo duża (np. powyżej 100). W przeciwnym wypadku AIC lub BIC może być lepsze. AIC minimalizuje błąd prognozy, a BIC razem z HQC dążą do wartości oczekiwanej prawdziwego modelu. Jeżeli wszystkie 3 kryteria wskazują najniższą wartość dla tego samego modelu, to nie mamy wątpliwości, że ten właśnie jest najlepszy z danego zestawu modeli. Dodatek armax umożliwia ocenę tych 3 kryteriów jednocześnie.

Najpierw użyłem metody exact ML, w której AIC, BIC i HQC wskazały jednoznacznie ten sam model: AR(4)-MA(2). Aby sprawdzić prognozę, robię to samo co poprzednio - używam standardowego modelu ARIMA, w której wpisuję odpowiednio rzędy parametrów (oczywiście bez sezonowości). Tym razem dostałem następujący wykres prognozy:

 

A więc zupełnie inna prognoza niż ta wzięta z X-13-ARIMA. Która jest bardziej prawdopodobna? Tutaj idealnie sprawdzą się kryteria AIC, BIC i HQC. Trzeba tylko wiedzieć, że:
AIC = Kryt. inform. Akaike'a
BIC = Kryt. bayes. Schwarza
HQC = Kryt. Hannana-Quinna

Zauważmy, że każdy z tych kryteriów jest niższy dla ostatniego modelu, co oznacza, że ożywienie w budownictwie nie jest wcale takie pewne.

Analogicznie sprawdźmy conditional ML. Wyniki tym razem nie są jednoznaczne, bo każde kryterium przynosi inne rezultaty. Biorąc pod uwagę, że mamy 199 obserwacji, wykorzystam wnioski Hannana i Quinna i skorzystam z HQC, które wskazało model AR(8)-MA(0). Użycie prognozy z tymi parametrami przyniosło bardzo podobną prognozę, co dla exact ML, tylko nieco bardziej stromą.


TRAMO/SEATS:
Drugą procedurą jest TRAMO/SEATS ("Time Series Regression with ARIMA Noise Missing Observations, and Outliers" / "Signal Extraction in ARIMA Time Series"). Przy standardowych opcjach, klikamy OK; GUS-owski Wskaźnik koniunktury budownictwa został podzielony na 3 części: wyrównany sezonowo, trend/cykl oraz odchylenia przypadkowe:

Znów dostaniemy skomplikowane statystyki, ale najbardziej interesuje nas to:
MODEL FITTED

      NONSEASONAL     P= 0     D= 1     Q= 1
         SEASONAL    BP= 0    BD= 1    BQ= 1
      PERIODICITY    MQ= 12

Przekładając na X13-ARIMA zapisalibyśmy to tak: (0, 1, 1)(0, 1, 1)

Czyli dostaliśmy ten sam schemat co w X-13-ARIMA. Nie zawsze jednak tak będzie.

Tutaj także możemy zapisać składnik trendowy (Wynik -> Zapisz...), który chcielibyśmy prognozować. Aby to zrobić możemy wybrać - tak jak poprzednio - jedną z dwóch metod. A więc pierwsza to zaznaczenie zmiennej trendu i użycie dla niej TRAMO-SEATS. Otrzymamy statystyki, wśród których szukamy:

 MODEL FITTED

      NONSEASONAL     P= 3     D= 1     Q= 1
      PERIODICITY    MQ= 12

Zniknęła sezonowość (SEASONAL), co świadczy, że jesteśmy na dobrej drodze. Uzyskane parametry wykorzystujemy do prognozy, a więc znów wchodzimy w Modele szeregów czasowych -> ARIMA. Tam wpisujemy te parametry AR=3, I=1, MA=1. Wyniki:

I prognoza:

W przeciwieństwie do X-13-ARIMA, Tramo sygnalizuje szczyt.

Druga metoda to identycznie jak przy X-13-ARIMA użycie dodatku armax. Tworzymy więc pierwsze przyrosty zmiennej trendu, zaznaczamy tę nową zmienną i możemy porównać dwie pod-metody: exact ML i conditional ML. Wyniki dla exact ML nie są jednoznaczne. AIC wskazuje AR(6)-MA(10), BIC i HQC - AR(1)-MA(4). W dodatku conditional ML prowadzi do jeszcze innych wniosków: dla AIC AR(12)-MA(1), BIC AR(3)-MA(1) i HQC AR(6)-MA(1). Idąc tropem HQC, exact ML i conditional ML przynoszą podobną prognozę. Poniżej wykres dla pierwszej z nich.

Ta prognoza współgra z tą, przeprowadzoną na Tramo. Sugeruje stabilność w branży budowniczej, ale nie ożywienie. W sumie więc zaprzecza początkowej intuicji, że trend jest rosnący.

Uwaga na fałszywe trendy
Niedawno opisywałem filtr Hodrika-Prescotta. Przestrzegałem, że może być mylący, bo zawsze wskaże "trend" nawet, jeśli jest to tylko przypadek. Tutaj jest podobnie. Samą zmienną trendu nie można się sugerować. Może wspomagać, ale musi być powiązana z istotnością statystyczną modelowania pierwotnej zmiennej. Np. wygenerujmy 199 obserwacji procesu Browna i zastosujmy X-13-ARIMA. Dostaniemy:

Jaki otrzymałem model optymalny?

Final automatic model choice : (0 1 1)

Wprawdzie nie pokazał żadnej sezonowości, ale przecież powinniśmy mieć (0 1 0). Przekonamy się o tym dopiero gdy sprawdzimy istotność statystyczną takiego modelu; wtedy otrzymamy p = 0,63.
Tylko że ta wartość p dotyczy zmiennej ruchu Browna, a nie zmiennej "trendu". Ponieważ procedura filtruje dane, to zaczyna tworzyć sztuczny trend.

Dlatego nie można się sugerować tutaj nazwą "trend", bo prawidłowo powinno brzmieć po prostu "filtr". Dlatego zarówno X-13-ARIMA jak i TRAMO/SEATS wymagają dużej ostrożności przy wnioskowaniu co do przyszłości.

Literatura:

[1] Andersen, E. B., Asymptotic Properties of Conditional Maximum-Likelihood Estimators, 1970,
[2] Beach, C. M., MacKinnon, J. G., A Maximum Likelihood Procedure for Regression with Autocorrelated Errors, Jan. 1978,
[3] Hannan, E. J., Quinn, B. G., The Determination of the Order of an Autoregression, 1979.

P. S. Więcej informacji na temat tych programów (m.in. ich historii) napisałem tutaj.

Historyczne wyceny WIG20 okazują się potwornie żmudne - nowe wnioski z analizy

 Im dłużej poświęcam czas na historyczne wyceny WIG20, tym więcej zauważam przeszkód. Ogólnie można je podzielić na 5 grup. Pierwsza to wpływ takich operacji jak np. odcięcie dywidend i splity. Problem jest następujący. Powiedzmy, że spółka ogłasza, że wypłaci w danym roku bardzo dużą dywidendę. Kurs oczywiście rośnie. Ale wtedy wycena musiałaby dodać tę dywidendę bieżącą, dopóki nie ma odcięcia. W sukurs przychodzi stooq.pl, który ma domyślną opcję korygowania ceny przed zdarzeniem odcięcia - historycznie wszystkie ceny przed dniem odcięcia są pomniejszane o tę odcinaną dywidendę, dzięki czemu nie ma nagłych przeskoków cen, tak że wycena w oparciu o te skorygowane ceny przed odcięciem będzie uwzględniać już tylko przyszłe dywidendy. Analogicznie z innymi operacjami. 

Skorygowane ceny na pierwszy rzut oka wydają się pełnym rozwiązaniem tego problemu. Ale pomyślmy głębiej. Cena w przeszłości jest korygowana o wszystkie dywidendy z przyszłości. Gdyby odcięcie dotyczyło tylko danego roku, nie byłoby problemu, bo wartość księgowa jeszcze w tym roku zostanie także skorygowana o tyle samo. Ale gdy cena jest korygowana o dywidendę za wszystkie kolejne lata, to powstaje kłopot, bo wartość księgowa będzie zawyżona w stosunku do ceny. Zarówno C/WK jak i C/Z będą za niskie. Skoro EPS uzyskuję przez odwrócenie C/Z i pomnożenie przez ceny WIG20 na dany dzień, to znaczy, że EPS będzie zawyżony, tak że wycena może okazać się zawyżona.

Z tym problemem wiąże się też inny - dostępność danych - jeżeli spółka została wycofana z giełdy, na Stooq już jej nie znajdziemy. Ceny starych spółek składu indeksu można znaleźć na gpwbenchmark.pl, ale tylko te faktyczne. Wszelkie operacje, jakie zaszły międzyczasie, są trudne do znalezienia.

Drugi (mniejszy) problem: trzeba pamiętać, że niektóre spółki nadal istnieją, ale zmieniły nazwę. Przykładem jest BZWBK, który obecnie nazywa się Santander Bank Polska. Na Stooq i w Bankierze więc będzie pod nową nazwą, ale w historii GPW będzie widniała stara nazwa. Tu powiem od razu, że irytuje mnie, że te serwisy nie konsolidują tych starych nazw z nowymi. Wystarczy przekierowanie ze starych na nowe. To chyba nie takie trudne?

Trzeci problem: wycena na bieżącym ROE okazuje się zgubna, gdy jest ono znacznie poniżej kosztu kapitału własnego. Okazuje się konieczne wyłuskanie średniej rentowności z przeszłości i na jej podstawie szacować przyszłość. Jeżeli jednak operujemy średnią ROE, to prowadzi to do kolejnego spostrzeżenia: aby wycena była spójna, koszt kapitału własnego musi też być średnią z przeszłości. Ja dotychczas operowałem bieżącymi szacunkami Damodarana. 

Czwarty problem, który sobie niedawno uświadomiłem to fakt, że szacunki Damodarana dotyczą wszystkich branż, czyli powiedzmy WIG, a nie WIG20. Czyli dotychczas zakładana stopa zwrotu była za wysoka, co zaniżało wyceny. 

Piąty problem: gdy największe spółki mają straty w danym roku, bieżący EPS zostaje tak zaniżony, że wycena okazuje się mało wiarygodna. Jeżeli spółka ma jednorazową stratę, to inwestorzy nie mogą zakładać, że we wszystkich kolejnych okresach będzie strata. Muszą zakładać, że będzie zysk, który będzie rósł. Dzieje się tak z dwóch powodów. Po pierwsze jednorazowa strata nie zwiastuje kolejnych. Oczywiście wiele zależy tu od określonych przyczyn i warunków sytuacji. Po drugie nawet jeśli spółka będzie miała wiele kolejnych strat z rzędu, to zostanie ona wykreślona z głównego indeksu: jej kapitalizacja spadnie, a zastąpi ją inna spółka. To jest zresztą istotny powód, dla którego lepiej inwestować w ETFy niż oddzielnie w każdą spółkę, jeśli nie mamy chęci / czasu na aktywne zarządzanie portfelem.

Teraz krótko omówię "rozwiązanie" każdego problemu.

Pierwszy problem rozwiązuję w ten sposób, że liczę wycenę dwukrotnie - raz z oryginalnymi cenami (bez korekt cenowych), a raz z korektą cen - jeżeli w ogóle jeszcze istnieją takie spółki na giełdzie (ponieważ niestety Stooq usuwa wycofane spółki ze swojej strony).  

Z tymi cenami bez korekt występuje też problem techniczny ze strony serwisów, jak bankier.pl czy biznesradar.pl. W przypadku splitów lub scaleń zmieniają wstecznie historyczne liczby akcji w wynikach finansowych, co okazuje się w zasadzie błędem z ich strony, a przynajmniej dużą przeszkodą dla ścisłych (księgowych) analiz. Przykładowo PZU miał split w 2015 r., a więc liczba akcji powinna się zmienić, ale tego nie widać w tabelach wyników w tych serwisach. Dlatego  tam gdzie trzeba wyjątkowo ręcznie zmieniłem liczbę akcji przed splitem, żeby dostosować do ceny sprzed splitu. Zaznaczyłem takie operacje purpurowym kolorem z komentarzem.

Jeśli chodzi o drugi problem, to nie jest on trudny dla mnie, bo akurat pamiętam, które spółki zmieniały nazwę.

Trzeci problem wymagał ode mnie dużej pracy - najpierw wyłuskałem wszystkie ROE WIG20 od 2009 do 2024 roku. Oto efekt końcowy:

Średnie ROE wyniosło 10,1% i wydaje się mniej więcej utrzymywać na tym poziomie. 

Następnie wyciągnąłem historyczne premie za ryzyko dla Polski od 2006 do 2024 wg szacunków Damodarana. Na koniec została najprostsza sprawa - premia za czas, tj. historyczne rentowności 10-letnich obligacji skarbowych. Suma obydwu premii daje koszt kapitału własnego.

No i tu dochodzimy do czwartego problemu. Szacunki Damodarana dotyczą kosztu kapitału dla całego kraju - wszystkie branże są brane pod uwagę. Czyli bardziej prawidłowo to jest dla WIG, a nie WIG20. Wydawałoby się, że skoro WIG20 zawiera więcej ryzyka politycznego (które raczej nie jest dywersyfikowalne), to suma sumarum oba indeksy mogą mieć zbliżoną oczekiwaną stopę zwrotu. Okazuje się jednak, że ryzyko polityczne nie kompensuje w całości ryzyka małych spółek. WIG daje większy zwrot. Dokładniej rzecz biorąc, porównałem WIG20TR (czyli WIG20 + dywidendy) z WIG i oto co dostałem:

Czyli WIG rośnie średniorocznie szybciej od WIG20TR o prawie 2 pkt proc. To rzeczywiście sporo, trzeba przyznać. Po prostu małe i średnie spółki są bardziej ryzykowne i inwestorzy wymagają większego zwrotu kapitału. 

Tak więc okazuje się, że dotychczas zawyżałem koszt kapitału WIG20 o te 1,9 pkt proc. 

Gdyby odjąć aż 1,9 pkt proc, to powinienem przyjąć zaledwie 9% oczekiwanej stopy zwrotu To nie za mało? Na wszelki wypadek sprawdziłem za pomocą testu CUSUM zmianę średniej różnicy, ale ewidentnie mamy tu losowe odchylenia:

Średnia różnica jest stabilna, czyli muszę po prostu przyjąć, że WIG20 ma niższy koszt kapitału (r) z powodu blue chipów. Dlatego ustalam, że r jest równy kosztowi kapitału dla Polski (lub WIG) minus 1,9 pkt proc. Wtedy wykres dla r kształtuje się następująco:

Jeżeli WIG20 płaci tylko 9%, to znaczy, że ROE > r, a to wyjaśnia dlaczego C/WK > 1 (średnia z lat 2006-2024 to 1,49). 

Trzeba zauważyć, że średnia ROE i r będą się zmieniać z roku na rok, bo wycena jest na dany dzień historyczny, a więc zawiera dane do ostatniego dostępnego roku. Chodzi o to, aby uzyskać faktyczne wyceny z przeszłości, a nie sztucznie dopasowywać do nowych danych.

Zostaje piąty problem. Proponuję robić 3 warianty. Pierwszy jest najważniejszy i polega na uśrednieniu dwóch kolejnych wariantów. Nie chodzi tu o sztuczną średnią z wycen, ale z założeń. I tak, drugi wariant uwzględnia wszystkie bieżące zyski i straty, a także wszystkie bieżące parametry (ROE, r itd.). Dlatego ten wariant nazywam "Wg danych bieżących". Trzeci wariant nazywa się "Wg danych uśrednionych": zeruje straty dla obliczenia EPS oraz uśrednia parametry z poprzednich lat. Usuwanie strat nie polega na usuwaniu wartości księgowej ani kapitalizacji spółek ze stratami, a jedynie wynik netto zamiast być stratą, stanie się zerem, co technicznie wdrażam przez sumę warunkową (wartość większa od zera). 

Przy tych założeniach mam następujące wnioski:

1. Im bliżej teraźniejszości, tym znacznie lepiej używać cen bez korekt, a im dalej w przeszłość, tym lepiej sprawdza się wycena po cenach z korektami. Zrobiłem wyceny na 30.10.2015, 18.11.2016, 17.11.2017, 19.11.2021 i 18.11.2022. Prześledzę krótko te wyceny:

30.10.2015: Cena zamknięcia 2060,0 pkt. Powiedzmy środek bessy. Wycena bez korekt (wariant 1): 1667 pkt, czyli całkowite dno, które zostanie wydrążone po prawie roku. Wycena z korektami: wariant 1 daje 2286 pkt, czyli bieżący kurs byłby poniżej tej wartości. Najbliżej kursu był wariant 2, bo wskazał   1954 (odchylenie ledwo -5%).

18.11.2016: Cena zamknięcia 1725,74 pkt. Było to prawie dno bessy. Wycena bez korekt (wariant 1): 941,78 pkt, czyli mocna przesada. Wycena z korektami (wariant 1):   1179 pkt, czyli bieżący kurs byłby powyżej tej wartości o 32%. Innymi słowy można było się spodziewać dalszych spadków, a mimo to rozpoczęła się hossa. Jedynie wariant 3 wskazał prawidłową wycenę (1788 pkt).

17.11.2017: Cena zamknięcia 2443,83 pkt. Był to szczyt hossy. Wycena bez korekt (wariant 1):   2354 pkt, czyli niemal w punkt. Wycena z korektami (wariant 1): 2891, tzn. zawyżona co najmniej o 18%. 

19.11.2021: Cena zamknięcia 2248,18 pkt. Był to szczyt hossy. Wycena bez korekt (wariant 1): 3851 pkt, czyli mocne przeszacowanie. Wycena z korektami (wariant 1): nie ma nawet sensu wspominać, bo grubo ponad 4000. Można najwyżej wspomnieć o wariancie 3: 3281 pkt, czyli co najmniej 46% za dużo. Powinniśmy się więc byli spodziewać kontynuacji hossy, tymczasem dostaliśmy za chwilę gwałtowny zjazd. Trzeba tu jednak przypomnieć, że przyczyna była bardziej przypadkowa i niezależna od czynników ekonomiczno-finansowych: wojna na Ukrainie. Gdyby nie to, zapewne rynek byka by się utrzymał. To też nauczka, że "czarny łabędź" wcale nie musi być taki rzadki (półtora roku wcześniej była pandemia).

18.11.2022: Cena zamknięcia 1707,23 pkt. Był to początek hossy. Wycena bez korekt (wariant 1):   2601 pkt, czyli 52% wyżej. Kurs dotarł tam, ale dopiero po 2,5 latach. Trochę bardziej wiarygodna była wycena wariantu 3: 2334 pkt, którą indeks osiągnął po ok. roku. Wyceny z korektami są zdecydowanie przesadzone.

2. Przy tych nowych założeniach dużą rolę zaczyna odgrywać współczynnik zatrzymania. Za bieżący współczynnik przyjmowałem średnią z ostatnich 2 lat, a więc można nim trochę kalibrować i dostosowywać, tak by wycena się bardziej "zgadzała". Jednak przypomnę, że, jak pokazałem w tym wpisie, EPS WIG20 stoi w miejscu od 2007 r.:

Na marginesie nie ma znaczenia, że Stooq błędnie wylicza C/Z dla indeksu (zob. tu), z którego wyliczałem EPS, bo chodzi tu wyłącznie o zmianę poziomów, a nie sam poziom. Ponieważ wzrost zależy od tego, ile zatrzymamy zysku, to widzimy, że nie możemy dowolnie ustawiać współczynnika zatrzymania. 

Oczywiście ktoś obeznany w temacie może zauważyć, że zatrzymany zysk może być kompensowany większym nabyciem akcji własnych, ale nie ma znaczenia co jak nazwiemy czy jaką mechanikę zastosujemy, bo liczy się odpowiednio rosnąca dywidenda, a nie precyzyjna proporcja zatrzymania zysku. A dywidenda z kolei rośnie mniej więcej tak jak WIG20, co można sprawdzić analizując stopę dywidendy WIG20 (dane historyczne dostępne w Stooq). Spójrzmy na ten rysunek:

Skoro średnia - niecałe 3% - jest dość stabilna (przy okazji w następnych latach stopa będzie pewnie spadać w ramach regresji do średniej, przynajmniej stopniowo), to znaczy, że dywidenda zmienia się właśnie w tempie WIG20. Problem polega na tym, że główny indeks nie pokonał nawet szczytu z 2007 r. Jednak od 2000 r. średnio rośnie on w tempie 3,4% (w sumie to tempo inflacji). Dlatego, skoro wiemy, że średnie ROE = 10,1%, to ok. 34% tej rentowności można traktować jako substytut współczynnika zatrzymania (10%*0.34 = 3,4%). Tak więc moje szacunki pokrywają się z pośrednio wyliczoną wartością tego współczynnika. Dla niektórych lat przyjmowałem 47-50%, co może być trochę za dużo (bo daje to 5% wzrostu), ale akurat taką średnią 2-letnią dostałem i myślę, że można do takiego poziomu podnieść w czasach hossy.

3. Chociaż trudno mi się pogodzić z r = 9%, to mam z tyłu głowy, że towarzyszy mu właśnie niski wzrost dywidendy. W modelu Gordona mamy r = d + g, gdzie d to stopa dywidendy i powinna być równa 2,9%, a g to wzrost dywidendy, który wychodzi 3,4%. Razem to 6,3%, czyli faktycznie mamy nadwyżkę 2,7%. Ta nadwyżka wynika z faktu, że parametry, których używam w modelu nie są zupełnie stałe, ale są zmiennymi losowymi, tak że do wzoru r = d + g należy dodać korektę związaną z wariancją. Na ten temat chciałbym bardziej szczegółowo napisać oddzielny artykuł.     

Link do pełnej analizy znajduje się tutaj. Niedługo dodam tam bieżącą wycenę WIG20. Nie chcę mieszać jej z historycznymi wycenami, dlatego napiszę o tym osobny wpis.

Prognozowanie ROE S&P 500

Jest oczywiste, że nawet jeśli dałoby się dokładnie przewidzieć księgowe zyski czy rentowności spółek, to nie oznacza jeszcze, że możemy przewidzieć choćby kierunek ich ceny. Jak się zastanowić, to paradoksalnie właśnie ich losowość pomaga w prognozie. Bo to jak usuwanie sezonowości z danych - jeśli wiadomo, że w danym sezonie zawsze są lepsze/gorsze wyniki, to nie będzie on miał wpływu na ceny, bo rynek uwzględnia w cenie takie zależności. Warunkiem jest pewność tej sezonowości, bo jeżeli sezonowość zaczyna się wahać, okazuje się nieregularnym cyklem, to nie usuwamy jej i wchodzi w grę możliwość wykorzystania jej do przewidywania kierunku. Powiedzmy, że mamy model wskazujący, że ROE wzrośnie w następnym roku. W teorii oznacza to, że rynek powinien to uwzględnić poprzez wyższą wycenę akcji, ale w praktyce występują różne korekty, w tym przypadku spadkowe, które można wykorzystać do kupienia, aby sprzedać później po tej wyższej cenie.  

Warto na początku pokazać relację między rocznymi stopami zwrotu a zmianami ROE S&P 500 od 1981 r. do 2024:

Dane ROE od 1980 r. wyłuskałem z dwóch źródeł: 

1) raport "Economic and Capital Markets Analysis: 2023 Year in Review U.S. Chartbook" opracowany przez batesgroup.com, s. 19;

2) https://www.multpl.com/sitemap

Ten pierwszy pokazuje na wykresie historyczne P/BV indeksu od 1980. Aby uzyskać ROE, podzieliłem P/BV przez P/E, które z kolei dostałem z tego drugiego źródła. Z niego także dostaniemy P/BV, ale od 1999 r. Dlatego do 1998 mam szacunkowe ROE, a od 1999 już dokładne. To opracowanie historyczne ROE wraz z analizą udostępniam tutaj. Poniżej końcowy wykres:

To co się rzuca w oczy, to niezwykle niskie ROE w latach 80-tych (5-6%). Jednak jak sprawdzimy, P/E wynosił wtedy ok. 10-11, czyli 3 razy mniej niż obecnie. Ale ROE obecne jest też ponad 3x większe. Czyli wszystko się zgadza.

Tak jak pisałem poprzednio, ROE bywa rożnie definiowane. Z punktu widzenia wyceny interesuje nas postać EPS / BV * (1+w), gdzie w to tempo wzrostu EPS w następnym roku (lub oczekiwane w danym roku). Mógłbym oczywiście rozbić te wszystkie czynniki i sprawdzać oddzielnie, może nawet do jakichś interesujących wniosków bym doszedł, ale nie chcę komplikować analizy, dlatego będę prognozować w całości taką zmienną i będę ją nazywał po prostu ROE. Oprócz wspomnianej poprawności teoretycznej zaletą tego podejścia jest uwzględnienie jednocześnie rentowności i wzrostu zysku.

Chociaż rok 2025 jeszcze nie minął, to są prognozy EPS S&P 500 na ten rok. Np. na tej stronie podane są kwartalne prognozy i sumarycznie przewiduje się 239 dol. Byłby to wzrost o 14,5%. Zakładając obecne wskaźniki P/E i P/BV, w = 14,5%, ROE(2025) = 20,3%. W sumie dostalibyśmy następujący wykres:

Wszystkie poniższe kody są zapisane w języku R. Ładowane są pakiety: forecast, tseries, lmtest oraz strucchange.

Luka PKB jako zmienna niezależna

ROE można próbować prognozować za pomocą zmiennych różnych ekonomicznych, jak wzrost PKB, inflacja lub stopa procentowa. Wszystkie te zmienne są ze sobą powiązane, zatem nie możemy ich wrzucać po prostu do modelu, jakby były niezależne. Zamiast tworzyć skomplikowane modele, możemy zauważyć, że pewnym pośrednikiem tych zmiennych jest luka PKB. Ostatnio pisałem o tym, że FED powinien podnosić stopę procentową, aby zmniejszyć lukę, która jest wyraźnie dodatnia. Wiemy, że FED zachował się politycznie pod dyktando Trumpa - w związku z tym inflacja wzrośnie, choćby z powodu nakładanych ceł, które nie zostaną skompensowane niższą stopą proc.

Nakładając wykres luki PKB na ROE, dostajemy zaskakującą korelację: 

W sumie korelacja wynosi 55%. Jednak obie zmienne mogą być niestacjonarne (test ADF tak sugeruje), a wtedy taka korelacja jest wątpliwa. Trzeba je sprowadzić do stacjonarności. Okazuje się, że po tym przekształceniu korelacja zostaje zachowana (54,4%). W dodatku ujawnia się wtedy ujemna korelacja (-0,37) z opóźnionym ROE:

ccf(zmianyLuka, zmianyRoe)

Oznaczałoby to, że jeżeli w danym roku luka rośnie, to ROE rośnie i w następnym roku ROE spada albo rośnie wolniej. Analogicznie w przypadku zmniejszenia się luki. Dodatnia korelacja jest dość oczywista - wyższe lub niższe ROE odzwierciedla stan gospodarki (wzrost PKB oraz w jakimś stopniu inflację). Trudniejsza do wyjaśnienia jest ta ujemna opóźniona reakcja. Najpierw wykonamy test Grangera:

 > grangertest(x=zmianyLuka, y=zmianyRoe, order=1)

Granger causality test

Model 1: zmianyRoe ~ Lags(zmianyRoe, 1:1) + Lags(zmianyLuka, 1:1)
Model 2: zmianyRoe ~ Lags(zmianyRoe, 1:1)
  Res.Df Df     F Pr(>F)  
1     40                  
2     41 -1 4.439 0.0415 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

, który potwierdza, że to nie przypadkowa korelacja, ale przyczynowość (luka jest wcześniej). Następnie zbudujemy liniowy model z opóźnieniem zmian ROE i zbadamy i jego stabilność

> zmianyRoe_ts <- ts(zmianyRoe, start = 1982)  
> zmianyLuka_ts <- ts(zmianyLuka, start = 1982)
> zmianyRoe_ts1 <- window(zmianyRoe_ts, start = 1983)  
> zmianyLuka_ts1 <- window(zmianyLuka_ts, start = 1982, end = 2024)
> modelRoe <- zmianyRoe_ts1 ~ zmianyLuka_ts1
> bai = breakpoints(modelRoe)
> bai

	 Optimal 2-segment partition: 

Call:
breakpoints.formula(formula = modelRoe)

Breakpoints at observation number:
24 

Corresponding to breakdates:
2006 

Test Bai-Perrona wskazuje na zmianę struktury w 2006 r. Sprawdźmy jak zmieniają się współczynniki:

> coef(bai)
            (Intercept) zmianyLuka_ts1
1983 - 2006    0.562118      0.0796585 

2007 - 2025 0.241999 -2.1666214 

Wg testu do 2006 r. ujemna zależność w ogóle nie istniała. Żeby się upewnić czy przyczyną nie jest jakiś inny efekt jak autokorelacja reszt, kwadratów reszt czy heteroskedastyczność wykonujemy kilka testów:

> reszty <- resid(lm(modelRoe))
> 
> # autokorelacje reszt:
> bgtest(modelRoe, order = 2)

	Breusch-Godfrey test for serial correlation of order up to 2

data:  modelRoe
LM test = 3.643, df = 2, p-value = 0.162

> Box.test(reszty, type = "Ljung-Box")

	Box-Ljung test

data:  reszty
X-squared = 0.1114, df = 1, p-value = 0.739

> 
> # autokorelacje kwadratów reszt:
> Box.test(reszty^2, type = "Ljung-Box")

	Box-Ljung test

data:  reszty^2
X-squared = 0.4259, df = 1, p-value = 0.514

> 
> # heteroskedastyczność
> bptest(modelRoe)

	studentized Breusch-Pagan test

data:  modelRoe
BP = 2.145, df = 1, p-value = 0.143

> gqtest(modelRoe) 

	Goldfeld-Quandt test

data:  modelRoe
GQ = 0.9562, df1 = 20, df2 = 19, p-value = 0.54
alternative hypothesis: variance increases from segment 1 to 2

> hmctest(modelRoe)

	Harrison-McCabe test

data:  modelRoe
HMC = 0.4871, p-value = 0.504

Żaden z tych testów nie wykrył problemów z resztami, a zatem zmiana struktury wydaje się faktem. Ostatni test, jaki warto przeprowadzić, aby się upewnić, że zmiana struktury dotyczy współczynnika kierunkowego, a nie wyrazu wolnego, to Score-CUSUM, oddzielający graficznie możliwe załamanie w każdym współczynniku:

Wynik pokrywa się z testem Bai-Perrona, który wskazał rok 2006 jako rok zmiany struktury. Teraz dodatkowo mamy dowód, że chodzi wyłącznie o zmiany luki PKB, a nie wyrazu wolnego czy wariancji. 

Nadal nie odpowiedziałem skąd może wynikać ta zmiana. Biorąc pod uwagę, że rok 2007 był przełomowy w historii finansów, możliwe, że tu coś "pękło". Danych jest trochę mało, ale taki model nadal jest istotny:

> zmianyRoe_ts2 <- window(zmianyRoe_ts, start = 2007)  
> zmianyLuka_ts2 <- window(zmianyLuka_ts, start = 2006, end = 2024)
> modelRoe <- zmianyRoe_ts2 ~ zmianyLuka_ts2
> summary(lm(modelRoe))

Call:
lm(formula = modelRoe)

Residuals:
    Min      1Q  Median      3Q     Max 
-11.044  -0.403   0.184   1.762   4.152 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)       0.242      0.758    0.32  0.75338    
zmianyLuka_ts2   -2.167      0.463   -4.68  0.00022 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.3 on 17 degrees of freedom
Multiple R-squared:  0.563,	Adjusted R-squared:  0.537 
F-statistic: 21.9 on 1 and 17 DF,  p-value: 0.000215

Rysunek dobrze naświetla co tu się dzieje:

Rok 2007-2008 był wstrząsem który mógł zmienić strukturę zachowania samego PKB (być może jest to związane z działaniami państwa w celu zminimalizowania kryzysu). Rzeczywiście, od 2007-2008 r. pojawiła się ujemna autokorelacja 1 rzędu dla zmian wzrostu PKB (do 2005 była ujemna, ale nieistotna stat., po tym roku wynosi ponad 0,6). Logika podpowiada, że to właśnie ta autokorelacja jest przyczyną ujemnej opóźnionej korelacji z ROE. Wyżej widzieliśmy, że ROE koreluje dodatnio ze zmianami PKB dla tego samego roku. Zatem skoro zmiany PKB ujemnie autokorelują, to jeśli w jednym roku są wyższe, to w następnym roku będą niższe lub ujemne, a z tego wynika, że w tym drugim roku zmiana ROE też będzie średnio biorąc niższa lub ujemna.  

Zmiana luki PKB w 2025 najprawdopodobniej będzie lekko ujemna, co sugeruje wzrost ROE w 2026. Usuwając stałą z modelu, która w tym przypadku raczej przeszkadza, dostajemy:

> modelRoe <- zmianyRoe_ts2 ~ zmianyLuka_ts2 - 1
> summary(lm(modelRoe))

Call:
lm(formula = modelRoe)

Residuals:
    Min      1Q  Median      3Q     Max 
-10.803  -0.155   0.425   1.994   4.410 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
zmianyLuka_ts2   -2.161      0.451   -4.79  0.00015 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.22 on 18 degrees of freedom
Multiple R-squared:  0.561,	Adjusted R-squared:  0.536 
F-statistic:   23 on 1 and 18 DF,  p-value: 0.000146

I teraz prognoza:

> predict.lm(lm(modelRoe), newdata = data.frame(zmianyLuka_ts2 = zmianyLuka[length(zmianyLuka)]))
       1 
0.359036

Na 2026 prognoza ROE = 20,3 + 0,36 = 20,7%. 

Uzyskaliśmy więc prognozę ROE na poziomie ok. 21%. Skoro stopa zwrotu koreluje dodatnio z ROE, to możemy się spodziewać w sumie wzrostu S&P 500 jeszcze w 2026... 

Nazwałbym osiągnięciem to co uzyskałem, bo nie jest to mechaniczny model oparty na ARMA. Żaden system AI czegoś podobnego nie zbuduje. Wymaga namysłu, wiedzy i prawdziwego rozumowania ekonomicznego. Ale nie znaczy to, że z ARMA można zupełnie zrezygnować. Zawsze warto podeprzeć się paroma modelami, bo prognoza z kilku perspektyw jest bardziej obiektywna. 

Autoregresja i średnia krocząca

Najpierw jednak trzeba sprawdzić rozkład zmian ROE. Właściwie powinienem to sprawdzić na samym początku, ale normalność nie ma aż takiego znaczenia w przypadku zwykłej regresji liniowej. Jednak tam gdzie stosuje się metodę największej wiarygodności, lepiej dostosować odpowiednio dane do rozkładu albo rozkład do danych. Robimy test Shapiro-Wilka:

 > shapiro.test(zmianyRoe)

	Shapiro-Wilk normality test

data:  zmianyRoe
W = 0.8992, p-value = 0.00102

Czyli rozkład jest daleki od normalnego. Popatrzmy w ogóle na te zmiany w końcu:

Wygląda to tak, jakby wariancja się zmieniła. Wykonajmy gqtest i hmctest na samych zmianach:

> modelRoe <- zmianyRoe_ts ~ 1
> gqtest(modelRoe)

	Goldfeld-Quandt test

data:  modelRoe
GQ = 2.002, df1 = 21, df2 = 21, p-value = 0.0598
alternative hypothesis: variance increases from segment 1 to 2

> hmctest(modelRoe)

	Harrison-McCabe test

data:  modelRoe
HMC = 0.3332, p-value = 0.056

Rzeczywiście jest blisko heteroskedastyczności, co stanowi dodatkowy argument, żeby dokonać korekty obserwacji. Idziemy więc do auto.arima i tam dostosujemy dane.

Aby znormalizować niegaussowskie szeregi w auto.arima stosujemy transformację Boxa-Coxa. Są tu dwa argumenty do użycia: lambda oraz biasadj. Pierwszy to parametr transformacji. Domyślnie ustawiony na NULL, a teraz ustawiamy go na "auto". Drugi to parametr odwrotnej transformacji. Jeżeli ustawimy go na FALSE, to estymatory dopasowania oraz prognozy - jeżeli lambda nie jest NULL - są medianą (przekształcenie monotoniczne funkcji nie zmienia jej wartości środkowej, ale średnią owszem). Żeby sprowadzić medianę do średniej, potrzebna jest korekta - to robi biasadj. Możemy ustawić biasadj = TRUE, jeśli chcemy mieć średnią, czyli wartość oczekiwaną. Istotne, żeby w takim przypadku używać samego ROE jako zmiennej prognozowanej, dlatego że transformacja Boxa-Coxa działa tylko na dodatnich wartościach. 

Użyję kryterium AIC. BIC okazało się za surowe - wskazało błądzenie losowe. 

> roe_ts <- ts(roe, start = 1981)
> arima_model <- auto.arima(roe,
+                           ic = "aic", max.p = 10, max.q = 10,
+                           stepwise = FALSE, approximation = FALSE, lambda = "auto", biasadj = TRUE)
> arima_model
Series: roe 
ARIMA(0,1,2) with drift 
Box Cox transformation: lambda= 1.18272 

Coefficients:
         ma1     ma2  drift
      -0.365  -0.268  0.516
s.e.   0.145   0.147  0.322

sigma^2 = 32:  log likelihood = -137.31
AIC=282.62   AICc=283.65   BIC=289.76
> prognoza <- forecast(arima_model, h=1)$mean
> prognoza
Time Series:
Start = 46 
End = 46 
Frequency = 1 
    out 
19.9176 
> roe_ts <- ts(roe, start = 1981)
> arima_model <- auto.arima(roe_ts,
+                           ic = "aic", max.p = 10, max.q = 10,
+                           stepwise = FALSE, approximation = FALSE, lambda = "auto", biasadj = TRUE)
> arima_model
Series: roe_ts 
ARIMA(0,1,2) with drift 
Box Cox transformation: lambda= 1.18272 

Coefficients:
         ma1     ma2  drift
      -0.365  -0.268  0.516
s.e.   0.145   0.147  0.322

sigma^2 = 32:  log likelihood = -137.31
AIC=282.62   AICc=283.65   BIC=289.76
> prognoza <- forecast(arima_model, h=1)$mean
> prognoza
Time Series:
Start = 2026 
End = 2026 
Frequency = 1 
    out 
19.9176

    

Optymalny model to MA(2) i prognozuje on lekki spadek ROE w 2026 do poziomu ok. 20%. Gdybyśmy pominęli transformację, wartość niewiele by się zmieniła, bo dostalibyśmy równo 20%.

W zasadzie nie musimy nawet zastanawiać się nad tym czy transformować czy nie, tylko od razu użyć z tego samego pakietu modelu BATS (akronim od: Box-Cox transform, ARMA errors, Trend, and Seasonal components), który automatycznie sprawdza czy transformować dane i jeśli potrzeba to robi; poza tym wyławia trend, sezonowość i ARMA w resztach. 

> roe_model_bats <- bats(roe_ts, biasadj = TRUE)
> prognoza <- forecast(roe_model_bats, 1)$mean
> prognoza
Time Series:
Start = 2026 
End = 2026 
Frequency = 1 
[1] 19.6489

Tak samo jak poprzednio, jeśli biasadj = FALSE (domyślnie), dostaniemy medianę. Właściwie to w ogóle można pominąć ten argument, bo mediana często jest lepsza od średniej  dla niegaussowskich rozkładów (zob. ten wpis). 

BATS nie jest to bardziej zautomatyzowana wersja auto.arima, tylko bardziej złożony model w przestrzeni stanów (state space model). Stanem może być średni poziom ROE, nachylenie trendu albo pewna cykliczność. Jak sama nazwa wskazuje, to stan, który obowiązuje w danym czasie, a więc w następnych okresach może się zmienić. Mimo to współczynniki ARMA pozostają stałe. Intuicyjnie: jeżeli mamy tłumioną sinusoidę, to funkcja jest stała, w tym zawiera parametr tłumienia; stanem jest wtedy amplituda (wielkość ruchu), która się zmniejsza wraz z każdym okresem. 

ARMA + luka PKB

W końcu można połączyć model czysto ekonomiczny z ARMA, czyli użyć ARMAX. Do auto.arima wpisujemy argument xreg, który będzie zmianami luki PKB. Powrócę więc też do zmian ROE - i to od 2006 r., bo, jak pamiętamy, regresja z luką do tego roku była niestabilna.  Oznacza to, że lepiej wtedy nie transformować danych, bo występują ujemne wartości. Muszą wystarczyć surowe:

> arima_model <- auto.arima(zmianyRoe_ts2,
+                           ic = "bic", max.p = 10, max.q = 10,
+                           stepwise = FALSE, approximation = FALSE, xreg = zmianyLuka_ts2)
> arima_model
Series: zmianyRoe_ts2 
Regression with ARIMA(0,1,1) errors 

Coefficients:
         ma1    xreg
      -0.619  -2.478
s.e.   0.197   0.381

sigma^2 = 10.7:  log likelihood = -46.07
AIC=98.14   AICc=99.86   BIC=100.81

> prognoza <- forecast(arima_model, h=1, xreg = zmianyLuka_ts2[length(zmianyLuka_ts2)])$mean
> prognoza
Time Series:
Start = 2026 
End = 2026 
Frequency = 1 
[1] 1.29567

Nie ma znaczenia którego kryterium użyję - daje to samo. W sumie prognoza jest mocno dodatnia. Na 2026 mamy 20,3 + 1,3 = 21,6%. Pamiętać oczywiście trzeba, że model powstał w oparciu o bardzo małą próbę (ledwo 19 obserwacji), więc nie jest to szczególnie wiarygodna wartość. Raczej traktować ją należy jak drogowskaz.

Podsumowując, mechaniczne metody prognozują delikatny spadek ROE, ale efekt ekonomiczny zmienia ten kierunek na dodatni. I wstępnie dodaje giełdzie amerykańskiej pretekstu do dalszych zwyżek. 

Link do danych z tej analizy