Giełdowy Racjonalista

Jak bardzo AI przyspieszy wzrost zysków w kolejnych latach?

2026-02-18T21:46:00.038+01:00

To co się obecnie dzieje w świecie AI, szczególnie LLM (Large Language Model), przyprawia o zawrót głowy. Tempo rozwoju tej technologii jest oszałamiające, zarówno wertykalnie jak i horyzontalnie. To znaczy serwisy branżowe co chwilę donoszą o nowych funkcjach lub poprawionej efektywności jakiegoś modelu, a jednocześnie można korzystać z coraz większej ilości narzędzi do różnych celów. Dla gigantów z S&P 500 rozwój ich własnych AI to priorytet. Dla inwestorów oczywiście najważniejsze to pytanie o ich zyskowność, a więc czy zwrot będzie odpowiedni w stosunku do kosztów. Technologia ta pozwala zarówno zmniejszać koszty (zastępować ludzi) jak i zwiększać przychody (firmy kupują AI zamiast zatrudniać), jednak jej utrzymanie kosztuje bardzo dużo. Sprzęt to jedno, ale istnieje też ryzyko niepowodzenia. Nie chodzi tu nawet o realne możliwości danego modelu, ale o to, że konkurencja może przejąć klientów. Np. Microsoft i Google (czy Alphabet - choć ta nazwa jest taka jakaś...) wydadzą miliardy dolarów na AI (w 2025 jest to odpowiednio 65 i 60 mld dol - źródło), a w tym samym czasie ktoś trzeci "niszczy system", przejmuje klientów, a giganci przepalają fortunę.

Pytanie więc brzmi: jak może wyglądać wzrost EPS S&P 500 w kolejnych 5-10 latach? I teraz ciekawa sprawa: wydawałoby się, że będzie to powyżej średniej; ale okazuje się, że średni wzrost od 1981 do 2022 r. wynosi 15%. Dopiero mediana daje bardziej umiarkowany, choć nadal bardzo wysoki, bo 10%. A takie są mniej więcej oczekiwania co do kolejnych lat. J.P. Morgan szacuje wzrost 13-15% na lata 2026-2027 (źródło), a Goldman Sachs 12% w 2026 i 10% w kolejnym (źródło).

Sam Standard and Poor's szacuje jednak ogromny wzrost ok. 20% na 2026 (źródło). Biorąc pod uwagę, że mamy do czynienia ciągle z hossą, przyjmę ten optymistyczny scenariusz, natomiast na 2027 średnią z wyżej wymienionych instytucji 14 i 10%, czyli 12%.

Skłania to do pytania o lokalny szczyt samego impetu zysków. Coraz wolniejsze zmiany, a w końcu spadek zysku oznaczałoby występowanie cykliczności. Aby to sprawdzić najlepiej użyć periodogramu lub analizy spektralnej.

W R najłatwiejszą funkcją (w sensie najmniejszej liczby ustawień) jest spectrum():

spectrum(x, method = c("pgram", "ar"))

Z metod mamy do wyboru periodogram lub spektrum dla modelu autoregresyjnego. Zaletą tego pierwszego jest to, że stanowi metodę nieparametryczną - nie nakłada specjalnych warunków na sygnał (chociaż musi być stacjonarny). Jego wadą jest jednak duża wrażliwość na wielkość próby, wykres jest poszarpany (duża wariancja) i piki mogą się układać przypadkowo niezależnie od wielkości próby (estymator jest niezgodny). Zaletą tego drugiego jest gładkość i stabilność (estymator zgodny), ale posiada dużą wadę: sygnał musi być liniowym procesem autoregresyjnym.

Pomimo prostoty to spectrum() jest tylko nakładką (wrapper) na odpowiednio spec.pgram() i spec.ar(). Lepiej ich używać, bo od razu widać jakie są dodatkowe argumenty i można nimi manipulować. Tak więc dla periodogramu najlepiej użyć:

spec.pgram(x, spans = NULL, kernel, taper = 0.1,
pad = 0, fast = TRUE, demean = FALSE, detrend = TRUE,
plot = TRUE, na.action = na.fail, ...)

Dzięki kalibracji spans, kernel i taper możemy zminimalizować przytoczone wyżej wady surowego periodogramu.

Na ten moment zwrócę tylko uwagę na demean i detrend. Już samo detrend = TRUE oznacza usunięcie zarówno średniej jak i (liniowego) trendu, dlatego demean (usunięcie średniej) jest zbędne. Jeśli używamy czystego EPS, to jednak detrend nie wystarczy, bo usuwa tylko zwykły trend liniowy, a EPS rośnie przecież nieliniowo, powiedzmy wykładniczo (link do danych):

Po zlogarytmowaniu:

W tym miejscu kłania się problem skrajnego minimalizmu dokumentacji R, bo okazuje się, że nie wystarczą argumenty z funkcji spec.pgram, ale trzeba uwzględnić argumenty funkcji plot.spec, którą tamta funkcja wykorzystuje. Niestety domyślnie funkcja zamienia dane na logarytmy i dlatego trzeba dodać specjalny argument log = "no".

W sumie tworzymy surowy periodogram z dwoma niezbędnymi argumentami (zmienną oznaczamy jako pg):

pg <- spec.pgram(x = log_eps, log = "no")

Powyższy wykres periodogramu w zasadzie udowadnia dużą przydatność tego narzędzia. Maksimum wyróżnia się od całej reszty, a jego położenie względem osi X wskazuje długość cyklu. Aby ją wyciągnąć używamy dwóch zmiennych obiektu pg:

pg$freq oraz pg$spec.

freq_values <- pg$freq

spec_values <- pg$spec

Oznaczamy dalej:

freq_values <- pg$freq
spec_values <- pg$spec
idx_max <- which.max(spec_values)

Szukaną częstotliwością odpowiadającą maksimum oznaczymy fr_max i będzie to:

fr_max <-freq_values[which.max(spec_values)]

Otrzymałem fr_max = 0.125.

Długość cyklu jest odwrotnością częstotliwości. Najprostsze wyjaśnienie może być na takim przykładzie. Powiedzmy, że mamy tylko 3 punkty w czasie, a każdy punkt to koniec roku. Mamy więc tylko 2 pełne lata. Pierwszy rok wzrost, a drugi rok spadek. To oznacza, że okres lub cykl trwa 2 lata i rok to jednostka czasu. Czyli w jednym roku mieści się pół cyklu, tj. 1/2 - odwrotność okresu. Dostajemy więc długość cyklu:

> 1/fr_max
[1] 8

A więc dostaliśmy "magiczną" liczbę 8 lat - tyle wynosi pełen cykl EPS: wzrost + spadek.

Kolejnym krokiem jest użycie wykrytej cykliczności do budowy filtra Butterwortha (FB). Dotychczas używałem go w gretlu, ale teraz zrobię to w R, w bardziej zaawansowany sposób. Użyję pakietu gsignal, w którym potrzebujemy dwóch funkcji: butter() oraz filtfilt(). Pierwsza oblicza współczynniki filtra, a druga podstawia je bezpośrednio do filtra.

Najważniejsze argumenty w funkcji butter() to rząd filtra oraz graniczna częstotliwość (znormalizowana częstotliwość odcięcia). Pozostałe argumenty dotyczą typu FB (w większości przypadków interesuje nas dolnoprzepustowy, to jest domyślne ustawienie), rodzaju płaszczyzny (wybór między filtrem cyfrowym a analogowym, interesuje nas zawsze cyfrowy - domyślne ustawienie) oraz sposobu konstrukcji (tu wybieramy "Sos", skrót od Second-order sections - mechanizm ten zwiększa stabilność numeryczną filtrowania poprzez podział jednej wielkiej operacji na wiele mniejszych i ich połączenie w kaskadową całość).

Chociaż pakiet ma ciekawe funkcje optymalizujące zarówno rząd filtra jak i częstość odcięcia, to dla naszych potrzeb wystarczy dotychczasowa wiedza, że ustawiamy rząd filtra = 2, natomiast znormalizowana częstotliwość odcięcia jest to f_cut / (fs/2), gdzie f_cut - częstotliwość odcięcia odpowiadająca maksimum periodogramu plus pewna nadwyżka (która bywa potrzebna, żeby uwzględnić niestabilność cyklu), fs to częstotliwość próbkowania (por. z opisem funkcji SciPy, który jest pierwowzorem). W naszym przypadku dane są roczne, więc pobieramy 1 na rok, tj. fs = 1, tak że zostaje 2*f_cut. To mnożenie przez 2 wydaje się być związane z twierdzeniem Nyquista o próbkowaniu. Zgodnie z tym twierdzeniem, aby sygnał został odtworzony bez zniekształceń, częstotliwość próbkowania musi być co najmniej dwa razy większa od częstotliwości odcięcia. Filtrowanie jest podobną operacją do próbkowania, dlatego można traktować taką częstość jako odpowiednią do zachowania sygnału. Od razu jednak zaznaczę, że trochę zgaduję, bo może być tak, że to tylko przyjęta konwencja konstrukcji tej funkcji. Mimo wszystko nie ma wątpliwości, że konwencja ta wywodzi się wprost z pojęcia unormowanej częstotliwości [zob. Loy, G., Musimathics. The Mathematical Foundations of Music, Volume 2, 2007, s. 19].

Żeby to zilustrować, pokażę jak wygląda FB z częstością odcięcia bez normalizacji, a jak wygląda po normalizacji. Przykład dla sinus z szumem:

# --- Parametry sygnału ---

fs <- 100 # częstotliwość próbkowania w Hz

t <- seq(0, 5, by=1/fs)

f_sin <- 5 # częstotliwość sinusa w Hz

# --- Sygnał sinusoidalny + biały szum ---

set.seed(123)

x <- sin(2*pi*f_sin*t) + rnorm(length(t), sd=0.5)

# --- Widmo sygnału ---

pg <- spec.pgram(x, log="no", plot=FALSE)

freq_values <- pg$freq # jednostki Nyquista [0,0.5]

spec_values <- pg$spec

idx_max <- which.max(spec_values)

fr_max <- freq_values[idx_max] # główny pik

# --- Dwie wersje częstotliwości odcięcia ---

fr_cut1 <- fr_max # zwykłe

fr_cut2 <- fr_cut1*2 # "podwójne" dla efektu wizualnego

# --- Filtry Butterwortha z gsignal ---

order_filtr <- 2

coef1 <- butter(order_filtr, fr_cut1, type="low", output="Sos")

coef2 <- butter(order_filtr, fr_cut2, type="low", output="Sos")

x_filt1 <- filtfilt(coef1, x)

x_filt2 <- filtfilt(coef2, x)

# --- Porównanie na wykresie ---

plot(t, x, type='l', col='gray', main='Sygnał sinus + szum')

lines(t, x_filt1, col='blue', lwd=2)

lines(t, x_filt2, col='red', lwd=2)

legend("topright", legend=c("Oryginalny", "fr_cut", "fr_cut*2"),

col=c("gray","blue","red"), lwd=2)

Widzimy, że czerwony filtr (z normalizacją) dużo lepiej odzwierciedla przebieg faktycznych cykli. Nawet jednak ten filtr nie dostosowuje się idealnie do głównego sygnału - tak jak pisałem tutaj ważny jest kierunek, a nie wartość. Chociaż niebieska linia też wskazuje ten sam kierunek, to spójrzmy na końcówkę po prawej stronie - czerwona linia prawidłowa już zakręca, a niebieska dalej spada, wprowadzając w błąd. Przy większej zmienności także może zmylić badacza/inwestora/trejdera. Dlatego wybór częstości ma duże znaczenie.

Wspomniałem wyżej, że samo podwojenie częstości może nie wystarczyć i przy niestabilności cykli trzeba dodać nadwyżkę. Pytanie ile dodać? Intuicyjnym sposobem jest dodanie takiej nadwyżki, która koresponduje z rozpoczęciem się piku spektralnego lub jego otoczeniem. Przy dużej zmienności piki nie będą bowiem linią pionową, ale zawsze będą mieć kształt trójkąta albo paraboli. Aby wykonać to zadanie wykorzystamy funkcję findpeaks() w tym samym pakiecie, która służy do szukania lokalnych ekstremów. Kod:

# 1. Wyszukiwanie lokalnego ekstremum i dopasowanie paraboli do struktury piku
peak <- findpeaks(pg$spec, MinPeakHeight = max(pg$spec)*0.99)

# 2. Pobranie precyzyjnego, ułamkowego indeksu prawego punktu przecięcia z podstawą (a1)
idx_a1 <- peak$roots$a1

# 3. Zaokrąglenie wyniku do najbliższego całkowitego indeksu istniejącego binu w periodogramie
idx_s <- round(peak$roots$a1)

# 4. Odczytanie dyskretnej częstotliwości (cykle/rok) odpowiadającej wybranemu binowi szumu
fr_s <- pg$freq[idx_s]

# 5. Obliczenie statystycznej odległości (pasma przejścia) między sygnałem a poziomem tła
d_peak <- abs(fr_max - fr_s)

# 6. Wyznaczenie finalnej częstotliwości granicznej pasma zaporowego dla projektu filtra
fr_cut <- fr_max + d_peak

Użycie progu 99% wysokości maksymalnej jako MinPeakHeight jest celowym zabiegiem. Choć mogłoby się wydawać, że szukamy po prostu max(pg$spec), to funkcja findpeaks robi coś znacznie ważniejszego: wykonuje interpolację (dopasowanie paraboli) do struktury piku (stąd potrzebne przybliżenie). Dzięki temu nie otrzymujemy tylko surowego, najwyższego binu z widma, ale cały obiekt (peak) zawierający precyzyjne parametry, takie jak ułamkowe indeksy podstawy piku (roots). To właśnie te dodatkowe dane, a nie tylko samo maksimum, są kluczowe do poprawnego wyznaczenia marginesu filtra.

Otrzymane fr_cut wyniosło 0.208. Zobaczmy jak to wygląda na periodogramie:

W końcu dochodzimy do funkcji butter(). tworzymy dolnoprzepustowy filtr Butterwortha:

order_filtr <- 2
fr_norm <- fr_cut*2
coef_butter <- butter(order_filtr, fr_norm, type = "low", output = "Sos")
sygnal_filtr <- filtfilt(coef_butter, log_eps)
rysuj_wykres(x=lata, y1=log_eps, y2=sygnal_filtr, y1_nazwa="EPS S&P 500")

Trzeba sobie teraz zadać pytanie czy przewidywania przedstawione na początku nie są przesadnie optymistyczne biorąc pod uwagę, że w ciągu ostatnich10 lat filtr rósł w linii prostej. Oczywiście można twierdzić, że era AI wszystko zmienia, że game changer" itd., ale era Internetu na początku 2000 r. też miała wszystko zmienić.

Teraz zostaje nam prognoza naszego filtra. Standardem jest rozpoczęcie od ARIMA. Znajdźmy optymalny model - możemy użyć pakietu forecast albo rugarch. Ponieważ jednak bardziej nas interesuje długoterminowa prognoza, to użyjemy tego pierwszego, bo ma on opcję użycia warunkowej metody największej wiarygodności (w tym przypadku jest ona tożsama z warunkową metodą najmniejszych kwadratów, bo dla reszt zakładany jest rozkład Gaussa), która jest lepsza dla takich celów (zob. np. tu). Używając funkcji auto.arima() wpiszemy więc argument method = "CSS" (ang. Conditional Sum of Squares).

Warto dodać, że dzięki zastosowaniu filtracji Butterwortha w dziedzinie logarytmów, uzyskujemy sygnał o ustabilizowanej wariancji i wyeliminowanych ekstremach, jak to w 2008-2009. W takim układzie dodatkowa transformacja Box-Coxa w modelu ARIMA staje się zupełnie zbędna.

Stosujemy kod:

library("forecast")

arima_model <- auto.arima(sygnal_filtr, ic = "aic", max.p = 10, max.q = 10, method = "CSS", stepwise = FALSE, approximation = FALSE)
summary(arima_model)
dopas <- arima_model$fitted
prognoza <- forecast(arima_model, h=6)$mean
rysuj_wykres(x=lata, y1=log_eps, y2=dopas, y_prognoza=prognoza, y1_nazwa="EPS S&P 500")

Otrzymałem model ARIMA(3,1,2)


Coefficients:
        ar1     ar2    ar3    ma1    ma2
      2.106  -2.001  0.770  1.857  0.926
s.e.  0.111   0.175  0.107  0.101  0.082

Oraz prognozy na kolejne 6 lat:


[1] 5.62094 5.65171 5.58799 5.51704 5.51587 5.57493

To samo dostałem dla kryterium "bic", więc to raczej nie przypadek.

Wg modelu prognoza na 2026 r. wynosi 5.62 i na 2027 5.65. Żeby porównać ze wskazanymi na początku 20% i 10%, musielibyśmy retrasformować uzyskane liczby exp(5.62) = 275.9 i exp(5.65) = 284.3. W 2025 mamy 244.5, stąd prognoza wzrostu na 2026 wynosi 275.9 / 244.5 - 1 = 13%. analogicznie na 2027 będzie to 284.3 / 275.9 - 1 = 3%. Czyli dostajemy tutaj mniejsze prognozy niż pierwotnie przyjęte (20% i 12%), ale należy uwzględnić fakt, że retrasformacja z postaci liniowej do nieliniowej powoduje pewien błąd w sytuacji, gdy wartość oczekiwana parametrów jest nieznana (zob. tu). Ta dodatkowa zmienność może jeszcze podnieść wartość. Niemniej model sugeruje szczyt w 2027 r., w efekcie spadki od 2028 r. W sumie możemy spodziewać się spowolnienia już od przyszłego roku.

Poziom WIG20 bardziej racjonalny niż sądziłem

2025-12-03T22:25:00.016+01:00

Dopóki nie miałem dość motywacji, by wnikliwie przeanalizować parametry wyceny WIG20 byłem absolutnie pewien, że mamy do czynienia z przewartościowaniem. Wychodzi jednak na to, że sprawa jest dużo bardziej złożona. Od razu mówię, że obecna wycena może nie być ostateczna, bo możliwe, że pewne (mniejsze lub większe) drobiazgi będę jeszcze poprawiał. Wycena jest na dzień 01.12.2025, kurs = 3007.41 pkt. Wartość wewnętrzna wychodzi w przedziale 2440-3010, tzn. indeks osiągnął górny pułap racjonalnej wyceny.

Przyczyny poprawy wyceny

Przyczyną tej zmiany jest, najkrócej mówiąc, zrozumienie, że należy przyjąć niższą oczekiwaną stopę zwrotu. Sprawę tę wyjaśniłem w poprzednim artykule, więc nie będę się tu o tym rozwodził. Jednocześnie myślę, że później trzeba będzie głębiej przyjrzeć się związkowi między oczekiwaną stopą zwrotu a oczekiwanym wzrostem dywidendy, bo dotychczas traktowałem je raczej jako niezależne od siebie. W każdym razie WIG20 ma niższy koszt kapitału od WIG o ok. 2 pkt proc. I te 2 procenciki wszystko zmieniają...

No dobra, może nie wszystko. Pisałem, że sprawa jest bardziej skomplikowana. Stopa dyskontowa to jedno, ale założenia ROE i EPS to drugie. To po kolei. Szczegółowa analiza wszystkich spółek WIG20 pozwala nie tylko sprawdzić, na ile obliczenia wskaźników przez GPW pokrywają się z moimi, ale też zobaczyć, które spółki są problematyczne (mają straty) i co się stanie, jeśli usuniemy te straty. Wyjaśniałem już, że w długim terminie spółki WIG20 nie mogą ponosić strat, bo zostałyby wyrzucone z indeksu i zastąpione zyskownymi. A nawet jeśli nie zostaną wyrzucone, to ich udział spadnie z powodu niskiej kapitalizacji i nie będą mieć większego wpływu na indeks. Z drugiej strony całkowite wyzerowanie strat też się mija z celem. Ktoś może argumentować, że jeśli straty są jednorazowe, to nie powinniśmy ich uwzględniać, bo wycenę tworzą jedynie przyszłe zyski. Ale przecież te przyszłe zyski są generowane przez kapitał własny, a ten został właśnie obniżony o te straty. Jeśli teraz wyzerujemy straty, to zawyżymy ROE, bo spadek KW będzie większy niż spadek sumarycznego EPS. Mój pomysł polegał na utworzeniu tych oddzielnych wariantów wyceny, a następnie uśrednieniu założeń obu, uzyskując w ten sposób unikalny wariant. Nazywam go pierwszym, bo jest ustawiony w pierwszej kolejności i jest najważniejszy. Pozostałe dwa służą głównie do porównań.

Porównując wskaźniki z mojej analizy do tych oficjalnych, zauważymy, że C/WK są bardzo zbliżone (u mnie 1.55, oficjalne - 1.61), natomiast C/Z już trochę oddalone (14.93 vs. 15,67). Różnica w ROE jest zaniedbywalna (10.4 vs. 10.3). Do tego momentu wszystko się zgadza. Ale co się stanie, jeśli stworzymy czwarty wariant polegający na użyciu bieżącego EPS z tych oficjalnych danych, ale użyciu pozostałych parametrów? Otóż dostaniemy wartość wewnętrzną 2538 pkt. Natomiast pierwszy wariant wskazał 3013.

Różnica między pierwszym a czwartym wariantem jest taka, że pierwszy mówi coś w rodzaju: "biorę EPS ze wszystkich spółek, następnie biorę EPS tylko spółek bez strat, i uśredniam oba", a czwarty mówi: "biorę EPS taki, jaki wychodzi z oficjalnego wskaźnika C/Z". Potem oba warianty biorą już te same uśrednione parametry.

2 techniki czy 2 metody?

Fakt, że C/Z jest mocniej odchylone w mojej analizie od GPW oznacza, że poszczególne spółki mają wpisane inne zyski. Aby się upewnić, że poprawnie konstruuję indeks oraz jego wskaźniki, analizę (a właściwie syntezę) wykonałem dwiema różnymi technikami. W pierwszej wyłuskałem z portalu bankier.pl liczbę akcji, księgową kapitalizację (kapitał własny) i zysk z ostatnich 4-ch kwartałów. W drugiej technice wykorzystałem wskaźniki podane przez stooq.pl. Najpierw spójrzmy sumarycznie: C/WK jest identyczne u obu, ale C/Z "bankierowe" jest znacznie niższe od "stooqowego" (14.9 vs. 17).

Głównym winowajacą tej różnicy jest Orlen. Otóż wg Stooq oraz GPW C/Z Orlenu wynosi 17, a u mnie 7.7. Ewidentnie coś tu jest nie tak. Zacząłem więc drążyć. Przejrzałem rzeczywiste sprawozdania. Okazało się, że zarówno Bankier, jak i Stooq mają rację. Bo tu wchodzimy w temat kreatywnej księgowości. Czwarty kwartał 2024 w sprawozdaniu faktycznie pokazuje duży zysk netto 4 664 000 tys. zł. Ale gdybyśmy dodali do niego pozostałe kwartały, dostalibyśmy więcej niż zysk za cały rok 2024! Ten paradoks powstaje w wyniku ogromnych odpisów ujawnionych dopiero na etapie rocznego sprawozdania (6,2 mld zł).

Jak porównamy sobie np. kwartalne wyniki w Biznesradar, to zobaczymy, że zyski za 1-3 kw 2024 są takie same jak w Bankierze, ale za 4-ty kw jest strata 1,55 mld zł, a w Bankierze wspomniane 4,7 mld zysku. Bankier pokazuje oryginalne dane z raportu kwartalnego, a Biznesradar wrzuca odpis do ostatniego kwartału, aby wyrównać z danymi rocznymi.

Z księgowego punktu widzenia to Bankier prezentuje poprawne wyniki. Ale to Biznesradar oddaje pełny obraz sytuacji. Czyli fajnie tu widać różnicę między księgowym a inwestorskim punktem widzenia. Stooq oraz GPW podają więc C/Z dla Orlenu taki, jaki wynika inwestorskiej perspektywy.

I teraz pozornie wydaje się, że skoro podejście Stooq jest lepsze dla inwestora, to powinniśmy porzucić bankierową technikę. Ale znowu - zastanówmy się czy rzeczywiście te odpisy są tak ważne, by je umieszczać w EPS. To są jednorazowe zdarzenia, które nie mogą być traktowane na równi ze zwykłymi kosztami. Zresztą kurs Orlenu dowodzi, że inwestorzy już dawno "przetrawili" te straty. Jeśli mamy prognozować przyszłe C/Z, to prędzej to bankierowe (7,7) będzie przybliżać jego wartość.

Ale, tak jak pisałem, nie można zupełnie pomijać jednorazowych strat, które mają negatywny wpływ na wartość poprzez niższy KW - niezależnie od tego czy ta strata jest czysto księgowa (papierowa / memoriałowa) czy kasowa. Zanim te odpisy nastąpiły, kupiono aktywa za pieniądze. Te pieniądze ktoś dał w dobrej wierze. Jeśli aktywa nie ma, to i zwrotu pieniędzy (z inwestycji) nie ma. Przecież cały kapitał własny opiera się na dobrej wierze inwestorów, że spółka zwróci w przyszłości im pieniądze i jeszcze wypłaci nadwyżkę - w postaci dywidend. Skoro tej przyszłości nie ma, to znaczy, że te pieniądze przepadły. Stąd papierowe zyski mają znaczenie dla wyceny.

Podsumowując tę część, zrobiłem dwa arkusze z wyceną na ten sam dzień liczone inną techniką... czy może metodą? Żeby było jasne - technika to nie metoda. Przez technikę rozumiem sposób zdobywania danych, a przez metodę sposób obliczania wyników i wyceny, którego rdzeń stanowi jakaś filozofia podejścia. Metoda będzie się różnić w sensie założeń, a technika w sensie źródła danych. Przykładowo, bankierową techniką obliczałem C/Z mnożąc cenę zamknięcia przez liczbę akcji, i dzieląc to przez sumę zysków z 4-ch ostatnich kwartałów (dane z bankiera). Stooqowa technika pozwala od razu wpisać C/Z wprost z portalu. Natomiast metoda jest tu inna tylko w tym sensie, że zysk TTM (Trailing Twelve Months, tj. z 12-tu ostatnich miesięcy) uwzględnia (Stooq) lub nie uwzględnia (Bankier) różnic między raportem rocznym a kwartalnym.

Uśrednienie prognoz

Po głębszym przeanalizowaniu tej sprawy myślę, że jest sensowne wyciągnięcie średniej z tych dwóch arkuszy. Tego już tam nie pokazuję, więc tutaj o tym napiszę. W każdym arkuszu za najważniejszy uważam wariant 1 i 4. W technice bankierowej są to odpowiednio wartości 3013 oraz 2538. W technice stooqowej 2585 oraz 2443. Najpierw uśredniłbym 3013 z 2585, czyli dostajemy:

(3013 + 2585)/2 = 2800.

Następnie wyciągnę średnią 2538 i 2443:

(2538 + 2443)/2 = 2490.

Czy jest sens wyciąganie jeszcze średnich z tych liczb? Wydaje mi się to już trochę sztuczne, bo pierwsza wartość koncentruje się bardziej na przyszłości, natomiast druga (czyli oparta o oficjalne wskaźniki GPW) bardziej na teraźniejszości. Ta pierwsza jest jakby bardziej nastawiona na to co będzie za jakiś czas, a ta druga bardziej na wycenę bieżącą (chociaż to taka luźna interpretacja).

Przedział prognozy

Uzyskany przedział 2490-2800 jest bardzo umowny i należy go prędzej traktować jak średnią z przedziałów. Bo taki najbardziej konserwatywny wynik byłby na poziomie tego 2440, a najbardziej optymistyczny 3310.

Rzecz jasna nie oznacza to, że WIG20 nie pójdzie jeszcze w górę, a jedynie, że każdy ruch powyżej obecnego poziomu rodzi coraz wyższe ryzyko (dla inwestora przesadnie duże). Przypomnę może wycenę vs. indeks w listopadzie roku 2017 (przy okazji - jakieś deja vu - znowu rok po rządach Trumpa):

Jest to jeden z możliwych scenariuszy.

Link do wycen - ten sam co poprzednio. Po prostu dorzuciłem dwa nowe arkusze z wycenami.

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

2025-11-30T12:44:00.012+01:00

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

2025-10-10T05:25:00.015+02:00

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

FED powinien podnosić, nie obniżać stóp

2025-09-16T07:45:00.007+02:00

Normalnie uważa się, że nie należy się kłócić z rynkiem. Ale w dzisiejszych czasach, to raczej rynek nie chce się kłócić z Trumpem. Trump żąda obniżki stóp procentowych od FEDu, rynek mu wierzy i pewnie pod to tak ciągną indeksy. Wszyscy wierzą, że FED pokornie obniży stopy i wygląda na to, że tak się właśnie stanie. Będzie to dowód na upolitycznienie tej instytucji, albowiem nic nie wskazuje, aby stopy były niższe, wręcz przeciwnie - powinny być lekko podniesione.

Po pierwsze zgodnie z wytycznymi celem FEDu powinno być utrzymanie inflacji na poziomie 2% (zob. źródło). Tymczasem inflacja w USA sięga 3% (źródło):

Po drugie reguła Taylora, którą FED powinien się kierować, wprost wskazuje zbyt niski poziom (źródło):

Powiedzmy jednak, że krytyk tego podejścia uzna je za zbyt sztywne, bo nie bierze pod uwagę przyszłej koniunktury. Pretekstem do obniżek mają być gorsze dane z rynku pracy. Tyle że jakoś FED prognozuje normalny, nawet lepszy niż w ostatnich latach, wzrost (źródło):

Stąd reguła Taylora w wersji dynamicznej, którą prezentowałem kilka lat temu, a teraz zaktualizowałem, pokrywa się z klasyczną (link do danych):

I rzeczywiście, zgodnie z tą regułą, luka PKB jest już przesadnie dodatnia i powinna być jak najprędzej ograniczana (źródło):

Luka trochę się zmniejszyła, zatem powinni przynajmniej utrzymać stopy na obecnym poziomie. Więc to jest test dla FEDu, czy poddadzą się politycznym naciskom.

Ile powinien kosztować WIG20 trzy lata temu?

2025-09-03T07:35:00.030+02:00

Żeby ocenić, na ile model wyceny jest sensowny albo skuteczny, należy przetestować go na historii. Modelując wycenę WIG20 na drugą połowę 2022 r., spodziewałem się, że model Gordona wskaże niedowartościowanie na początku hossy. Jednak indeks okazał się być prawidłowo wyceniony. Możliwe, że model jest za prosty, może trzeba by założenia parametrów poprawić, a może po prostu rynek zgodził się z modelem.

Jeśli chodzi o założenia parametrów można je oczywiście zmieniać, ale nie znaczy to, że można wstawiać dowolne liczby. Stopa dyskontowa 11,6% to suma premii za ryzyko USA, krajowego (na podst. strony Damodarana) i stopy wolnej od ryzyka (przyjąłem 6%). Współczynnik wypłaty dywidendy / zatrzymania jest trudny do oszacowania. Dziś wg moich szacunków wsp. wypłaty wynosi ok. 60-70%, czyli zatrzymania 30-40%. Można nawet dla szybkiego sprawdzianu spojrzeć na PKO: EPS = 8, dywidenda = 5.5, czyli wsp. wypłaty = 5,5/8 = 69%, a więc wsp. zatrzymania 31%. Podobnie PEKAO: odpowiednio 25,75 i 18,35, wsp. zatrzymania = 1 - 18,35/25,75 = 0,29. Ostatnio jak robiłem, to wpisywałem wsp. wypłaty = 66,66%, czyli prawidłowo, ale zamieniłem teraz na równe 70% (praktycznie nie zmienia to wyceny). Taki współczynnik generuje wzrost zysku między 3,5 a 4,6%.

Współczynnik ten jednak się zmienia. W 2022 r. wyszło mi niecałe 50% dla wsp. zatrzymania, dokładniej 47%, i tyle przyjąłem.

Ostatni aspekt, jaki chciałem krótko poruszyć dotyczy ROE. Nigdy tego szczegółowo nie analizowałem, bo nie uważałem tego za zbyt ważne, ale teraz wyjaśnię dokładniej, że istnieją 3 różne metryki ROE:

1) podstawowa: Z(t) / KW(t) , gdzie Z(t) to zysk netto w okresie t i KW(t) to wartość księgowa kapitału własnego w okresie t,

2) stosowana w różnych analizach finansowych: Z(t) / [KW(t) + KW(t-1)]/2,

3) jedyna prawidłowa dla wyceny: Z(t) / KW(t-1), dlatego że ROE porównuje się z kosztem kapitału własnego, a ten jest stopą procentową / zwrotu.

ROE dla wyceny można też zapisać jako Z(t) * (1 + w)/ KW(t) = Podstawowe ROE * (1+w) , gdzie w to tempo wzrostu zysku. Chodzi nam o oszacowanie przyszłego ROE, czyli przyszłego zysku. Tej formuły używam, ale za w nie wstawiam wzrostu dywidendy, tylko tempo nominalnego wzrostu PKB. Zysk może rosnąć nie tylko dzięki zyskowi zatrzymanemu, ale także dzięki emisji akcji i obligacji. Ta emisja pochodzi z dochodów ludności, które rosną właśnie w tempie w. Tak więc ten wzrost nie pochodzi od strony podażowej, ale popytowej - innymi słowy spółka może wyemitować maksymalnie tyle ile ludzie dodatkowo zarabiają. Tak więc wzrost z tej emisji nowego kapitału jest ograniczony możliwościami gospodarki. Z Banku Światowego można uzyskać dane na temat statystyk wzrostu nominalnego PKB i właśnie stąd w = 7%.

Wycenę zrobiłem na dzień 18.11.2022, kiedy skład WIG20 został zaktualizowany. Dostałem wartość indeksu na poziomie 1783 zł. Tego dnia WIG20 kosztował na koniec dnia 1707,23 zł, zatem ściśle biorąc był niedowartościowany, ale zaledwie o 4%, więc trudno tu mówić o jakimś niedowartościowaniu.

Co by się stało, gdyby przyjąć wsp. zatrzymania na poziomie dzisiejszym, 30-40%? Przy 30% jest to 1461, a przy 40% 1612. Czyli powiedzmy 1500 - też nie jest to poziom odległy, bo jeszcze miesiąc wcześniej indeks tam się wspinał. A co powodowało dalszy wzrost? Zysk netto spółek mocno wzrósł w kolejnym roku - średnio gdzieś o 40%.

Link do wyceny znajdziemy tutaj.

Krótka uwaga na temat free floatu

2025-08-31T22:49:00.005+02:00

Ostatnia wycena ad hoc WIG20, czyli z modelu Gordona, okazuje się jeszcze niższa niż te 1950 zł. Problem polega na wstawieniu odpowiedniej kapitalizacji i liczby akcji. W poprzedniej wycenie wstawiłem nie tyle wyemitowaną liczbę akcji, co liczbę akcji w obrocie, tj. free float. Po uwzględnieniu wszystkich akcji (na podstawie danych z bankier.pl), wycena spada do 1882 zł (link).

Można by przekonywać, że wartość tworzą jedynie ci, którzy handlują akcjami i stąd wycena nie powinna dotyczyć części zamrożonej. Free float obejmuje jednak tylko liczbę akcji na daną chwilę, a pakiety kontrolne też mogą być sprzedawane.

Mniejszy free float rodzi zarówno przychody jak i koszty. Koszty to mniejsza płynność (trudniej kupić / sprzedać po własnej cenie) oraz związane z tym ryzyko manipulacji kursem (przez zmowę spekulantów). Ponadto duża koncentracja akcji w rękach paru podmiotów rodzi ryzyko wywłaszczenia. O problemie wywłaszczenia, czy w ogóle to jest legalne, pisałem w art. Wyciskanie akcjonariusza, czyli kiedy wycena nie ma sensu. W sumie wzrasta koszt kapitału własnego - a to obniża wartość.

Z drugiej strony większa koncentracja kapitału w rękach inwestorów instytucjonalnych sprawia, że zarząd spółki nie może robić co mu się żywnie podoba i może być łatwo odwołany. Duzi inwestorzy patrzą na ręce zarządowi i nie będą pozwalać np. na dodatkową emisję akcji lub obligacji. Oczywiście ta zaleta koncentracji kapitału powstaje jedynie przy warunku, że pakiety kontrolne nie należą do podmiotów powiązanych z zarządzającymi. Jeśli tak się dzieje, to nie można mieć pewności, że będą zatrudniani ludzie wg klucza kompetencji, a wszelkie emisje będą tylko ułatwiane.

W każdym razie kupując akcje powinniśmy też zwracać uwagę na poziom free floatu i kto ma kontrolę. Wydaje się, że najlepszą sytuacją jest, jeśli po pierwsze kontrola jest rozproszona między wiele funduszy niezależnych od siebie, ale sumarycznie mają powyżej 50%, a po drugie, żeby ta kontrola nie była powyżej powiedzmy 80-85%. Wzorową sytuacją jest Kruk (kontrola = 73%, najwięcej ma niecałe 13% fundusz), a przeciwieństwem mBank (kontrola 90,08%, a inny (niemiecki) bank ma aż 69% akcji).

Policzmy w końcu poprawnie wskaźniki indeksów

2025-08-24T15:52:00.009+02:00

Zawsze gdy analizujemy dane finansowe, staramy się je uzyskać z najbardziej rzetelnych i aktualnych źródeł. Gdy więc chcemy zobaczyć bieżące C/Z i C/WK indeksu lub spółki, pierwszym wyborem powinno być to oficjalne. W przypadku WIG20, będzie nim strona gpwbenchmark. Widzimy, że C/Z równa się 14,66, a C/WK 1,64 (stan na 21.08.2025). Gdy jednak chcemy znaleźć historyczne wskaźniki, to dużo łatwiej oprzeć się na źródłach pośrednich. Korzystamy więc np. ze stooq. I okazuje się, że portal pokazuje inne dane: C/Z = 16,19 i C/WK = 2,86. Z czego wynika ta różnica? Czy stooq.pl pokazuje błędne wartości?

W jaki sposób policzyć wskaźnik dla indeksu? Zajmijmy się C/WK. Na logikę, cały indeks jako portfel traktujemy jak jedną spółkę. Czyli:

C/WK = Kapitalizacja portfela / Wartość księgowa portfela

Całą kapitalizację oznaczymy Kap. Powiedzmy, że mamy 3 spółki w portfelu. Zapiszemy:

gdzie Kap_i oznacza kapitalizację spółki i, czyli

gdzie N_i - liczba akcji spółki i, C_i - cena akcji.

Wtedy C/WK indeksu ma postać:

WK_i - wartość księgowa na akcję spółki i.

Taką właśnie definicję wskaźników znajdziemy w dokumencie GPW "Podstawowe algorytmy wskaźników rynku kapitałowego". Jest to nic innego jak średnia ważona harmoniczna:

I dla C/Z analogicznie:

gdzie Z_i - EPS dla i-tej spółki.

Taki wzór powinno się stosować.

Dlaczego nazywa się to średnią ważoną, skoro nie widać tu żadnych wag? Wzór można tak przekształcić, żeby wagi się pojawiły. Wtedy będzie można porównać powyższy wzór z innymi średnimi ważonymi. Najpierw zdefiniujemy udział i-tej spółki w portfelu:

Czy za udział można by alternatywnie wstawić liczbę akcji podzieloną przez liczbę wszystkich akcji? Nie, bo wtedy split czy resplit zmieniałby udział. "Waga" jest więc ciut lepszym słowem niż udział, bo jak coś waży, to wiadomo o co chodzi. Definicja udziału wywodzi się jednak formalnie z teorii portfela, w której inwestor przeznacza określoną sumę pieniędzy na zakup danego waloru (zob. ten wpis, w którym jest to wyjaśnienie).

To teraz przekształcimy pierwotny wzór:

Czyli dostajemy postać:

Powiedzmy jednak, że jesteśmy otwarci na eksperymenty z innymi metodami. Porównajmy naszą formułę ze średnią ważoną arytmetyczną i geometryczną.

Średnia ważona arytmetyczna

To metoda intuicyjna w tym sensie, że wykorzystujemy naturalnie wagi (udziały) spółek z portfela. Każdą wagę mnożymy przez wskaźnik dla danej spółki:

Wskaźnik C/Z wygląda analogicznie. Teraz pomyślmy co się stanie, gdy zysk trzeciej spółki wynosi zero lub jest bardzo niski. C/Z skacze do nieskończoności albo do ogromnych wartości, a więc średnia ważona C/Z dla całego portfela będzie nienormalnie duża. Czyli z punktu widzenia C/Z widać, że nie należy tego sposobu używać. Ale wszystkie wskaźniki muszą być liczone tą samą metodą, aby zachować spójność i proporcjonalność (chociażby, żeby policzyć ROE ze stosunku C/WK do C/Z). Zatem dostajemy dowód, że C/WK też nie może być liczone w ten sposób.

Średnia ważona geometryczna

Trzecią możliwością jest średnia geometryczna z wagami w potęgach:

Tutaj można przeprowadzić identyczne rozumowanie co w przypadku średniej arytmetycznej. Jeśli zysk wyniesie zero, to średnia ważona geometryczna skacze do nieskończoności. Nie musi być to nawet dokładnie zero, wystarczy, że będzie bardzo mały. Gdyby była strata, to już w ogóle nie dałoby się nic policzyć. Jedna spółka całkowicie psuje wskaźnik. Trzeba więc porzucić tę metodę, ale żeby zachować spójność, także C/WK nie może być w ten sposób liczone.

Przykłady

Zobaczmy przykłady trzech wariantów-okresów. Dla uproszczenia liczba akcji wszędzie równa się 1.

Okres I – wszystkie spółki mają podobny zysk

Dane wejściowe:

Spółka 1: C1 = 100, Z1 = 10, Kap1 = 100, C1/Z1 = 10

Spółka 2: C2 = 150, Z2 = 10, Kap2 = 150, C2/Z2 = 15

Spółka 3: C3 = 120, Z3 = 10, Kap3 = 120, C3/Z3 = 12

Suma kapitalizacji i zysków:

Kap1 + Kap2 + Kap3 = 370

Z1 + Z2 + Z3 = 30

Wagi:

w1 = Kap1 / (Kap1 + Kap2 + Kap3) = 100 / 370 = 0,27

w2 = 150 / 370 = 0,41

w3 = 120 / 370 = 0,32

Średnie ważone:

Arytmetyczna = 12,7

Geometryczna = 12,5

Harmoniczna = 12,3

Okres II – spółka 3 ma spadek zysku o 95% i kapitalizacji o 40%
Dane wejściowe:

Spółka 1: C1 = 100, Z1 = 10, Kap1 = 100, C1/Z1 = 10

Spółka 2: C2 = 150, Z2 = 10, Kap2 = 150, C2/Z2 = 15

Spółka 3: C3 = 72, Z3 = 0,5, Kap3 = 72, C3/Z3 = 144

Suma kapitalizacji i zysków:

Kap1 + Kap2 + Kap3 = 322

Z1 + Z2 + Z3 = 20,5

Wagi:

w1 = Kap1 / (Kap1 + Kap2 + Kap3) = 100 / 322 ≈ 0,31

w2 = 150 / 322 ≈ 0,47

w3 = 72 / 322 ≈ 0,22

Średnie ważone:

Arytmetyczna = 42,3

Geometryczna = 21,9

Harmoniczna = 15,7

Okres III – spółka 3 ma zysk zero i spadek ceny o 50%

Dane wejściowe:

Spółka 1: C1 = 100, Z1 = 10, Kap1 = 100, C1/Z1 = 10

Spółka 2: C2 = 150, Z2 = 10, Kap2 = 150, C2/Z2 = 15

Spółka 3: C3 = 60, Z3 = 0, Kap3 = 60, C3/Z3 = nieskończoność

Suma kapitalizacji i zysków:

Kap1 + Kap2 + Kap3 = 310

Z1 + Z2 + Z3 = 20

Wagi:

w1 = Kap1 / (Kap1 + Kap2 + Kap3) = 100 / 310 = 0,32

w2 = 150 / 310 = 0,48

w3 = 60 / 310 = 0,19

Średnie ważone:

Arytmetyczna = nieskończoność

Geometryczna = nieskończoność

Harmoniczna = 15,5

Spójrzmy jak "oporna" jest ta średnia harmoniczna - nawet przy nieskończonej wartości C/Z dla jednego składnika, jest stabilna, a dla ogromnych wzrostów wskaźników niewiele się zmienia.

Co ze stratami?

Konstrukcja średniej harmonicznej pozwala też użyć jej dla spółek ze stratami, jeśli nie przeważają nad zyskami pozostałych. Ponieważ wskaźniki poszczególnych spółek sztucznie usuwa się z analizy, gdy występują straty, to lepiej używać wzoru ze stosunkiem sum, ponieważ ten z wagami wymaga użycia indywidualnych wskaźników.

I tu zaczyna się problem. Obecnie dwie z WIG20 są pod kreską (Pepco i PGE). Jeśli ich nie uwzględnimy, to sztucznie zaniżymy C/Z. Powoduje to nie tylko błędną analizę porównawczą, ale też dostaniemy nieprawidłowe ROE dzieląc C/WK przez C/Z.

Potrzebujemy więc odpowiednich danych, tj. strat dla tych spółek, których nie wyłuskam ze wskaźników (raczej nie podaje się ich, gdy są ujemne). Spojrzałem więc na dane w dwóch źródłach - bankier.pl i biznesradar.pl. Oba portale pokazują inne liczby.

Biznesradar: dla Pepco wskazuje ponad 3 mld zł strat, dla PGE ponad 1,5 mld strat.

Bankier: dla Pepco w ogóle nie wykazuje ani zysku ani strat, dla PGE to 1,64 mld strat.

Biorąc dane Bankiera dostałbym C/Z = 14,59, natomiast wg Bizneradar ok. 17. Zakładając, że bankier.pl ma bardziej poprawne dane, możemy uznać, że gpw uwzględnia straty - wtedy wszystko będzie prawidłowo (przypominam, że podali 14,66).

Straty trzeba koniecznie uwzględnić we wskaźnikach nie tylko dlatego, że błędnie omijamy gorsze spółki, ale także po to, aby utrzymać współmierność z innymi wskaźnikami. Jeżeli policzymy C/WK normalnie sumując wszystkie kapitalizacje i dzieląc je na sumę wartości księgowych (czyli średnią harmoniczną), dostaniemy dla WIG20 1,65 (stan na 21.08.2025), czyli praktycznie idealnie się zgadza z danymi gpw, a różnica wynika z zaokrągleń (przypominam, że na gpw podali 1,64).

Jak robią na stooq.pl?

Wg średniej ważonej arytmetycznej dostajemy C/Z 17,7 i C/WK 3,25. Wg średniej ważonej geometrycznej będzie to odpowiednio 13,97 i 2,26. W ogóle to nie odpowiada 16,19 i 2,86, które prezentował stooq na 21.08. Czyli stosują jeszcze jakąś inną metodę. Jedynie, co warto odnotować, to to, że (17,7 + 13,97) / 2 = 15,84, a (3,25 + 2,26) / 2 = daje 2,76, czyli naprawdę blisko. Sugeruje to, że używają jakiejś pośredniej metody.

Problem z ROE

Takie różnice w szacowaniu wskaźników prowadzą do dużych rozbieżności w określeniu rentowności kapitału własnego. Poprawna metoda, którą tu pokazałem, wskazuje, że dla WIG20 ROE = 1,66 / 15,11 = ok. 11%. Wg oficjalnych wskaźników GPW, byłoby to 1.64 / 14.66 = 11,2%. Ale wg stooq to jest 17,7%. Raczej nie jest przypadkiem, że wg śr. ważonej arytm. wychodzi 17,8%.

Wszystko wskazuje na to, że stooq.pl zawyża wszystkie te wskaźniki, co prowadzi do zaburzeń przy wycenie indeksu. Moje obliczenia we wpisie Jaka piękna katastrofa , w którym badałem na szybko historyczne EPS, ROE i WK okazują się błędne. Jednak sama wycena, którą wtedy pokazałem, może być nadal poprawna. Od tamtego czasu wzrósł EPS, zatem będzie wyższa. Na ten moment wyceniam WIG20 na ok. 1950, czyli nadal uważam, że indeks jest mocno przewartościowany.

Wszystkie obliczenia do tego artykułu zamieściłem w tym linku.

Przecięcie KSP z DSP to zapowiedź bessy?

2025-08-11T16:22:00.008+02:00

Ostatnio giełda mnie mocno zaskakuje, chyba nie tylko mnie (nie spodziewałem się tak szybkiego powrotu do hossy), dlatego nie będę się mądrzył, że zaraz to na pewno wszystko walnie. Jedynie zwrócę uwagę na ryzyko jakie warto wziąć pod uwagę w najbliższej przyszłości.

Parę lat temu napisałem ten wpis i omówiłem wyniki badań związku między spredem rentowności a zmianami PKB. Dla przypomnienia:

KSP - krótkoterminowa stopa procentowa

DSP - długoterminowa stopa procentowa

Badania Clintona [Clinton, K., The term structure of interest rates as a leading indicator of economic activity: A technical note, 1995] wskazują, że spread między DSP a KSP to "doskonały predyktor zmian aktywności gospodarczej w Kanadzie". Gdy DSP jest znacznie większa od KSP, silny wzrost PKB nastąpi mniej więcej w kolejnym roku. Gdy KSP > DSP, prawdopodobnie szykuje się recesja. Trzeba dodać, że to badanie jest dość stare, bo dotyczy lat 1962-93.

Przeanalizujemy sytuację w USA. Do niedawna KSP była rzeczywiście większa od DSP - i to po raz pierwszy w tak widoczny sposób. Co więcej, porównałem momenty lokalnych szczytów S&P 500, po których następował spadek nie mniej niż 20%:

Najczęściej szczyty korespondowały z przecięciami KSP (czerwona) z DSP (zielona). Oczywiście sam moment był losowy, ale widać, że zakres, gdy obie stopy były bardzo blisko siebie, wiązał się z największymi spadkami indeksu. W dwóch przypadkach mechanizm ten nie zadziałał (1987, 2022) i prawidłowo, bo przyczyny leżały poza gospodarką USA. Wprawdzie na początku 2023 nastąpiło przecięcie, ale było to już zakończenie krótkiej, wojennej, bessy. Chociaż raczej mówimy tu o pełnym przecięciu, to czasem wystarczyło stykanie się. W 2019 KSP i DSP dotykały się, ale nie przecięły i nie było bessy - a mimo wszystko indeks i tak zaliczył istotną korektę.

Z czego wynika ten mechanizm? W uproszczeniu powinniśmy kojarzyć tak: KSP - teraźniejszość i bliska przyszłość (do roku), DSP - niedaleka przyszłość (do roku) plus daleka przyszłość (po roku). Przez przyszłość mam na myśli sferę realną oraz inflację. Inflacja teraz jest nieistotna, zajmijmy się częścią realną. W uproszczeniu: większe KSP niż DSP to lepsza teraźniejszość niż przyszłość i vice versa - mniejsze KSP niż DSP to gorsza teraźniejszość niż przyszłość.

To jednak duże uproszczenie, a sam spred niewiele mówi, biorąc pod uwagę, że jest on ujemny od końcówki 2022. Przez to zresztą rok temu pomyliłem się, sądząc, że bessa tuż za rogiem. Trzeba podzielić przyszłość na dwie części:

- bieżące oczekiwania, bardziej pewne

- większą niepewność przyszłości.

Jeżeli KSP przegania DSP, to można sądzić, że następuje przesadny optymizm czasem dyskontowany od razu (np. rok 2000), a czasami dopiero gdy rynek orientuje się, że był to hurraoptymizm i następuje powrót KSP do DSP (rok 2007). W pierwszym przypadku inwestorzy wierzą, że inwestycje firm szybciej się zwrócą lub przyniosą większe zyski. Sprzedają więc krótkie terminy i kupują długie, podwyższając KSP. W drugim przypadku ich oczekiwania maleją albo niepewność rośnie i znów kupują krótkie terminy, obniżając KSP. Jednak niepewność przyszłości i bieżące oczekiwania powinniśmy oddzielić. Jeśli pierwsze rośnie, a drugie spada, to DSP może stać w miejscu (bo DSP zawiera obydwa czynniki), a KSP może ją wtedy przeciąć. Gdy niepewność zacznie rosnąć, a bieżące oczekiwania będą stać w miejscu, to KSP może stać w miejscu (bo KSP nie zawiera dalekiej niepewności), a DSP przeciąć KSP (z dołu w górę). I taka sytuacja byłaby pozytywnym sygnałem w dłuższym terminie mimo iż niepewność wzrasta. Spred staje się dodatni, a rosnącą niepewność inwestorzy uznają bardziej za szansę niż zagrożenie, wiedząc, że w równowadze stopa procentowa powinna odzwierciedlać produktywność gospodarki.

Poprawka do analizy premii za ryzyko S&P 500

2025-08-01T17:23:00.012+02:00

Wykonana ostatnio analiza premii za ryzyko S&P 500 wskazywała, że roczne zmiany można modelować za pomocą ARMA(7, 4). Dziś widzę, że zawierała 3 problemy. Po pierwsze metoda "full" w funkcji autoarfima (pakiet rugarch dla R), którą wówczas uważałem za bardziej precyzyjną, lepszą, wcale taka nie jest. Poprzez symulacje przetestowałem ją i okazało się, że "partial" dużo lepiej znajduje prawidłowe rzędy. Wydawało się, że skoro "full" jest tak potwornie wolna, bo sprawdza wszystkie możliwe kombinacje (a więc np. 1 rzędu nie ma, ale jest drugi i trzeci, albo drugiego nie ma, ale jest pierwszy i trzeci, albo pierwszego i drugiego nie ma, ale jest trzeci itd.), to będzie najdokładniejsza. Chociaż na obrazkach rzeczywiście dopasowanie było niezłe, to jednak wiadomo, że w próbie zawsze może być zbyt przypadkowe.

Drugi problem to założenie rozkładu t-studenta w resztach modelu. Właściwie dlaczego założyłem taki? Stopy zwrotu są zazwyczaj dalekie od normalności, a rozkład Studenta uwzględnia grubsze ogony. Proces MA można zapisać jako AR wielu rzędów z przeszłości (zob. ten wpis), a skoro AR to właśnie stopa zwrotu, tylko opóźniona ze stałymi współczynnikiami, to znaczy, że proces MA ma ten sam rozkład - zakładając stacjonarność stopy zwrotu. A skoro (stacjonarny) MA ma ten sam rozkład, to znaczy, że składnik losowy też ma ten sam rozkład (bo MA jedynie mnoży jego wartość z przeszłości przez stałą). W sumie więc, jeśli stopa zwrotu ma rozkład Studenta, to reszty w ARMA też muszą taki posiadać.

Tylko że okazało się, że nie można odrzucić hipotezy o gaussowskości ostatnich 80 rocznych stóp zwrotu lub premii:

> shapiro.test(as.numeric(stopa_zwrotu))

	Shapiro-Wilk normality test

data:  as.numeric(stopa_zwrotu)
W = 0.9804, p-value = 0.257

Nawet ostatnie 120 danych (czyli od 1904 r.) nie zmieniło tego wyniku:

	Shapiro-Wilk normality test

data:  as.numeric(stopa_zwrotu)
W = 0.9854, p-value = 0.221

Nienormaność pojawia się dopiero przy większych częstościach. Bardzo możliwe, że przy większej liczbie rocznych obserwacji należy odrzucić rozkład normalny, ale dla danych, które analizowałem, nie można tego zrobić.

Po trzecie porównałem autoarfima z auto.arima z pakietu forecast, zrobiłem różne symulacje ARMA i czasami auto.arima wskazywała prawidłowy model, czasami autoarfima. Stąd doszedłem do wniosku, że należy wykonywać testy obiema funkcjami. W dodatku najlepiej stosować zarówno kryterium nie tylko BIC, ale też AIC.

Premia za ryzyko (od samych akcji)

autoarfima:

klaster <- makeCluster(detectCores() - 1)

premia_model_arma <- autoarfima(data = premia, ar.max = 10, ma.max = 10, cluster = klaster, solver = "nlminb")

stopCluster(klaster)

Wszystkie kryteria wskazały ten sam model:

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

Optimal Parameters
------------------------------------
       Estimate  Std. Error  t value Pr(>|t|)
mu     0.077126    0.021203   3.6376 0.000275
sigma  0.189641    0.014992  12.6491 0.000000

Robust Standard Errors:
       Estimate  Std. Error  t value Pr(>|t|)
mu     0.077126    0.019321   3.9917 0.000066
sigma  0.189641    0.017623  10.7608 0.000000

LogLikelihood : 19.4946 

Information Criteria
------------------------------------
                     
Akaike       -0.43737
Bayes        -0.37781
Shibata      -0.43857
Hannan-Quinn -0.41349

Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
                        statistic p-value
Lag[1]                     0.7322  0.3922
Lag[2*(p+q)+(p+q)-1][2]    1.5233  0.3557
Lag[4*(p+q)+(p+q)-1][5]    4.0949  0.2425

H0 : No serial correlation

Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
                        statistic p-value
Lag[1]                      1.392  0.2381
Lag[2*(p+q)+(p+q)-1][2]     1.749  0.3081
Lag[4*(p+q)+(p+q)-1][5]     2.110  0.5928


ARCH LM Tests
------------------------------------
             Statistic DoF P-Value
ARCH Lag[2]      2.253   2  0.3241
ARCH Lag[5]      2.403   5  0.7911
ARCH Lag[10]     5.979  10  0.8171

Nyblom stability test
------------------------------------
Joint Statistic:  0.2461
Individual Statistics:            
mu    0.1522
sigma 0.0835

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

Model jest dobrze dopasowany: stabilny, nie występują autokorelacje w resztach ani ARCH. Czyli nie żadna ARMA, ale zwykła średnia. Średnioroczna premia za ryzyko wynosi 7,7% +/- 19% (nie mylmy tu błędu standardowego reszt, który wynosi 1,9% - on wskazuje jak średnia przypadkowo się odchyla, a nie sama stopa zwrotu).

To samo auto.arima:

premia_model_arma <- auto.arima(premia, max.p=10, max.q=10, stepwise = FALSE, approximation = FALSE)

Argumenty stepwise i approximation ustawiam na FALSE, bo zwiększają precyzję algorytmu. Niezależnie od kryterium dostaniemy model:

Series: premia 
ARIMA(0,0,0) with non-zero mean 

Coefficients:
       mean
      0.077
s.e.  0.021

sigma^2 = 0.0364:  log likelihood = 19.49
AIC=-34.99   AICc=-34.83   BIC=-30.23

Obydwie funkcje - na obydwu kryteriach BIC i AIC wskazały ten sam model - zwykłą średnią. To jest silny argument za odrzuceniem poprzedniego modelu ARMA.

Można pokazać to na wykresie:

Premia za ryzyko od akcji i obligacji skarbowych

Obligacje skarbowe są obarczone głównie ryzykiem stopy procentowej. Aby je uwzględnić metodą ad hoc, odejmujemy od stopy zwrotu z akcji rentowność bonów skarbowych. Nie jest to prawdziwa premia, bo nie uwzględnia zależności między akcjami a obligacjami, ale daje pewien obraz.

Test na normalność:


	Shapiro-Wilk normality test

data:  premia
W = 0.9845, p-value = 0.447

Czyli też zachowuje się normalnie.

Analogicznie jak poprzednio używamy autoarfima i auto.arima.

autoarfima:

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

Optimal Parameters
------------------------------------
       Estimate  Std. Error  t value Pr(>|t|)
mu     0.089739    0.019332   4.6421 0.000003
sigma  0.172907    0.013670  12.6491 0.000000

Robust Standard Errors:
       Estimate  Std. Error  t value Pr(>|t|)
mu     0.089739    0.017760   5.0529        0
sigma  0.172907    0.012238  14.1292        0

LogLikelihood : 26.885 

Information Criteria
------------------------------------
                     
Akaike       -0.62212
Bayes        -0.56257
Shibata      -0.62333
Hannan-Quinn -0.59825

Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
                        statistic p-value
Lag[1]                     0.4299  0.5121
Lag[2*(p+q)+(p+q)-1][2]    1.4963  0.3619
Lag[4*(p+q)+(p+q)-1][5]    4.3507  0.2134

H0 : No serial correlation

Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
                        statistic p-value
Lag[1]                     0.9464  0.3306
Lag[2*(p+q)+(p+q)-1][2]    3.1373  0.1288
Lag[4*(p+q)+(p+q)-1][5]    4.8670  0.1640


ARCH LM Tests
------------------------------------
             Statistic DoF P-Value
ARCH Lag[2]      5.497   2 0.06403
ARCH Lag[5]      5.442   5 0.36440
ARCH Lag[10]    10.107  10 0.43112

Nyblom stability test
------------------------------------
Joint Statistic:  0.1892
Individual Statistics:            
mu    0.1121
sigma 0.0580

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

auto.arima:

Series: premia 
ARIMA(0,0,0) with non-zero mean 

Coefficients:
       mean
      0.090
s.e.  0.019

sigma^2 = 0.0303:  log likelihood = 26.88
AIC=-49.77   AICc=-49.61   BIC=-45.01

Ponownie obie funkcje wskazały ten sam model zwykłej średniej: "pełna premia" wynosi 9% +/- 17,3%. Taki wynik uzyskałem niezależnie od użytego kryterium.

I wykres:

Chociaż brak zależności w rocznych stopach pogarsza sytuację prognostyczną, to z punktu widzenia inwestora może się ona nawet poprawić. Po pierwsze ułatwia mu wycenę - nie musi uwzględniać zmienności stopy dyskontowej. Po drugie zwiększa szansę, że wycena będzie szybko zbiegać do jego wyceny, a nie od niej uciekać. Na ten moment, nawet bez wyceny, możemy spekulować, że premia w 2025 spadnie, choć powinna pozostać dodania.

Krótka analiza i prognoza premii za ryzyko S&P 500 (język R)

2025-06-09T19:19:00.054+02:00

Powiedzmy, że chcemy oszacować premię za ryzyko amerykańskiej giełdy w roku 2025 i 2026. Określenie "premia za ryzyko" używam jako skrót na oznaczenie nadwyżki stopy zwrotu z indeksu (tu S&P 500) nad stopę zwrotu z 10-letnich obligacji skarbowych. Precyzyjnie przecież mówiąc o premii mamy na myśli oczekiwane lub wymagane wynagrodzenie za przyszłe ryzyko, czyli takie, którego okres jeszcze nie zaczął się. I może się zrealizować (wtedy albo stracimy, albo mniej zarobimy), albo nie zrealizować (zarobimy co najmniej tyle, ile wymagamy). Dodatnia premia za ryzyko wynika z tego, że inwestorzy są niechętni wobec ryzyka. To znaczy, że wyceniają niżej aktywa o wyższym ryzyku. Mechanizm ten pozwala łatwo zrozumieć, dlaczego oczekiwana (a więc i średnia) stopa zwrotu jest dodatnia, a nie zerowa, jak w hazardzie. Otóż cała tajemnica leży w tym, że to ludzie wyceniają aktywa, a nie natura. Jeżeli wielu inwestorów będzie wyceniać niżej akcje z powodu ich większego ryzyka, to przecież spółka nie przestaje zarabiać tyle co zarabiała, jej wartość księgowa się nie zmienia. Spadek ceny oznacza więc, że inwestorzy mogą więcej zarobić, sprzedając w przyszłości akcje. Ryzyko jest mimo wszystko miarą psychologiczną, a więc dość subiektywną - i to właśnie sprawia, że inwestowanie daje dużo większe możliwości niż hazard.

Historyczne roczne premie za ryzyko S&P 500 dla lat 1945-2024 biorę ze strony Damodarana. Średnia wyniosła 7,7%. To coś dużo biorąc pod uwagę, że średnia stopa wolna od ryzyka to ok. 5%. Porównałem więc te dane z S&P U.S. Equity Risk Premium Index. Obejmuje on zajęcie długiej pozycji w indeksie S&P 500 Futures Excess Return Index i krótkiej pozycji w indeksie S&P U.S. Treasury Bond Futures Excess Return Index [zob. S&P Factor Indices Methodology]. Od 2015 roku dostałem taki wykres:

A więc dane są poprawne, bo idealnie korelują. Oczywiście w tym okresie średnia była jeszcze wyższa i dla S&P U.S. Equity Risk Premium Index wyniosła 13,2%, a Damodarana 13,5%.

Stabilność średniej

Pytanie brzmi czy średnia się zmieniała w ciągu tych 80 lat. Przez pierwsze 40 lat wyniosła ok. 8,6%, a kolejne niecałe 6,8%. Sprawdźmy czy możemy uznać to za coś więcej niż przypadek.

Jeżeli dzielimy okres na dwie części, to najlepiej użyć testu Chowa. W tym artykule pokazałem, jak go wykonać w R.

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

install.packages("strucchange")

library("strucchange")

}

sctest(premia ~ 1, type="Chow", point=40)

	Chow test

data:  premia ~ 1
F = 0.18, p-value = 0.7

Jak widać nie ma podstaw do odrzucenia hipotezy, że średnia jest stała. Wykonałem dodatkowo najróżniejsze testy (ADF, KPSS, Bai-Perrona, CUSUM, MOSUM) i nie wykryłem niestacjonarności czy niestabilności.

Modelowanie

Jeżeli średnia się nie zmienia, to trzeba sprawdzić autokorelację. Najpierw ACF:

acf(premia, 15)

Do dziś nie mogę z tego, że R pokazuje bez sensu ten zerowy rząd. Głupota. W cząstkowej jest już normalnie od pierwszego rzędu:

pacf(premia, 15)

W Dodatku do tego artykułu zrobiłem analogiczną analizę dla nadwyżki nad bonami. Tam PACF 4-tego rzędu jest na granicy istotności. Tutaj natomiast nie można mówić o jakiejkolwiek istotności. Z tego z kolei wynika, że obligacje skarbowe mają swój udział w cyklu 4-letnim (bo je tam dodaliśmy).

Mimo wszystko sprawdzę model AR(4) i jego stabilność. W R kod będzie następujący:

# Pierwszy model - arima

premia_model_arma1 <- arima(

premia,

order = c(4, 0, 0), # ARMA(4, 0, 0)

include.mean = TRUE, # estymuj μ

fixed = c(0, 0, 0, NA, NA), # φ1=0, φ2=0, φ3=0, φ4=estymuj, μ=estymuj

transform.pars = FALSE # aby fixed=0 było literalne

)

summary(premia_model_arma1)

Coefficients:
      ar1  ar2  ar3    ar4  intercept
        0    0    0  0.205      0.079
s.e.    0    0    0  0.111      0.026

sigma^2 estimated as 0.0344:  log likelihood = 21.16,  aic = -36.33

Gdyby zastosować taki model, dostalibyśmy rzeczywiście średnią 7,9% (+/- 2,6%), ale dodatkowo współczynnik 4-tego rzędu 0,205 (+/- 0,11).

Aby sprawdzić stabilność tego modelu, trzeba niestety zamienić na formę typu:

premia4 <- tail(zmienna, -4)

premia_4 <- head(zmienna, -4)

premia.model <- premia4 ~ premia_4

i przykładowo, żeby użyć testu CUSUM-Score, wpisać:

par(mfrow=c(3,1), mar=c(2.5,5,1.5,3))

cus = efp(premia.model, type="Score-CUSUM")

plot(cus, functional=NULL)

Za jednym zamachem mamy tu sprawdzenie stabilności wszystkich współczynników, w tym wariancji. Brak jakichkolwiek wychyleń (byłyby zaznaczenia) świadczy o stabilności.

To jednak tylko punkt wyjścia, bo znacznie lepszym sposobem jest użycie autoarfima() w rugarch lub auto.arima() w forecast. Wybrałem ten pierwszy, bo dostarcza wszystkie potrzebne statystyki. Ten argument omówiłem w tym artykule. Sprawdziłem więc w autoarfima max ARMA(7, 7) w metodzie full. Wybrałem rozkład "std", bo dawał najlepsze rezultaty, jeśli chodzi o statystyki dopasowania. Wadą rugarch jest czasochłonność obliczeń, dlatego koniecznie trzeba używać argumentu cluster tam gdzie się da. Ponadto wybrałem solver "nlminb", ponieważ wydaje się dość efektywny. Użyłem funkcji dwukrotnie - wg kryterium BIC i AIC (HQIC nie używam, bo danych jest za mało, zob. tu). Okazało się, że oba dały ten sam model. Zwracam też uwagę, że nie stosuję tutaj modelu GARCH ze względu na małą liczbę danych. Efekt:

# Drugi model - autoarfima (rugarch)

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

install.packages("rugarch")

}

library("rugarch")

# Kryterium BIC

klaster <- makeCluster(detectCores() - 1)

premia_model_arma_bic <- autoarfima(data = premia, ar.max = 7, ma.max = 7, criterion = "BIC",

method = "full", arfima = FALSE, include.mean = NULL,

distribution.model = "std", cluster = klaster, external.regressors = NULL,

solver = "nlminb")

stopCluster(klaster)

premia_model_arma_bic

# Kryterium AIC

klaster <- makeCluster(detectCores() - 1)

premia_model_arma_aic <- autoarfima(data = premia, ar.max = 7, ma.max = 7, criterion = "AIC",

method = "full", arfima = FALSE, include.mean = NULL,

distribution.model = "std", cluster = klaster, external.regressors = NULL,

solver = "nlminb")

stopCluster(klaster)

premia_model_arma_aic

premia_model_arma2 <- premia_model_arma_aic

*----------------------------------*
*          ARFIMA Model Fit        *
*----------------------------------*
Mean Model	: ARFIMA(7,0,4)
Distribution	: std 

Optimal Parameters
------------------------------------
       Estimate  Std. Error    t value Pr(>|t|)
ar1     0.44901    0.000926    484.974        0
ar2     0.00000          NA         NA       NA
ar3     0.00000          NA         NA       NA
ar4     0.00000          NA         NA       NA
ar5    -0.22529    0.000465   -484.974        0
ar6     0.00000          NA         NA       NA
ar7     0.34287    0.000027  12785.413        0
ma1    -1.28980    0.000061 -21191.451        0
ma2     0.00000          NA         NA       NA
ma3     0.61606    0.000029  21191.451        0
ma4     0.44741    0.000021  21191.451        0
sigma   0.13607    0.005403     25.184        0
shape 100.00000    0.004585  21809.383        0

Robust Standard Errors:
       Estimate  Std. Error    t value Pr(>|t|)
ar1     0.44901    2.324983    0.19312  0.84686
ar2     0.00000          NA         NA       NA
ar3     0.00000          NA         NA       NA
ar4     0.00000          NA         NA       NA
ar5    -0.22529    1.166585   -0.19312  0.84686
ar6     0.00000          NA         NA       NA
ar7     0.34287    0.000385  891.10681  0.00000
ma1    -1.28980    0.005203 -247.88844  0.00000
ma2     0.00000          NA         NA       NA
ma3     0.61606    0.000766  804.54128  0.00000
ma4     0.44741    0.001619  276.38488  0.00000
sigma   0.13607    0.619235    0.21974  0.82608
shape 100.00000    1.325372   75.45050  0.00000

LogLikelihood : 46.12 

Information Criteria
------------------------------------
                     
Akaike       -0.95299
Bayes        -0.71479
Shibata      -0.97067
Hannan-Quinn -0.85749

Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
                         statistic p-value
Lag[1]                    0.000611  0.9803
Lag[2*(p+q)+(p+q)-1][32] 11.645908  1.0000
Lag[4*(p+q)+(p+q)-1][54] 17.816677  0.9976

H0 : No serial correlation

Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
                        statistic p-value
Lag[1]                     0.4785  0.4891
Lag[2*(p+q)+(p+q)-1][2]    0.7811  0.5750
Lag[4*(p+q)+(p+q)-1][5]    2.3214  0.5447


ARCH LM Tests
------------------------------------
             Statistic DoF P-Value
ARCH Lag[2]      1.519   2  0.4678
ARCH Lag[5]      2.647   5  0.7542
ARCH Lag[10]     9.429  10  0.4920

Nyblom stability test
------------------------------------
Joint Statistic:  231429
Individual Statistics:             
ar1   0.09969
ar5   0.05997
ar7   0.01355
ma1   0.03352
ma3   0.11951
ma4   0.09983
sigma 0.13154
shape 0.42456

Asymptotic Critical Values (10% 5% 1%)
Joint Statistic:     	 1.89 2.11 2.59
Individual Statistic:	 0.35 0.47 0.75

BIC i AIC wskazały ten sam model: AR o rzędach 1, 5 i 7; MA o rzędach 1, 3 i 4. W nieodpornych statystykach (tj. przy założeniu danego rozkładu, tu t-Studenta) wszystkie są istotne. Jednak statystyki odporne (na odchylenia od założonego rozkładu reszt) wskazują jedynie istotność rzędu 7 dla AR, a dla MA wszystkie są istotne. Pozostałe statystyki są okej, wszystkie współczynniki są stabilne, chociaż znowu cały model mocno niestabilny (pewnie niezbyt idealny rozkład).

Test MCS

Aby się upewnić, że na pewno drugi model jest lepszy od pierwszego, użyjemy testu MCS, który omówiłem w poprzednim artykule. Tak jak wtedy za stratę uznam sytuację, gdy znak prognozy in-sample nie zgadza się ze znakiem rzeczywistej wartości.

#MCS test

Funkcja_strat <- function(rzeczywiste, prognoza) {

roznica <- rzeczywiste - prognoza

strata <- ifelse(sign(rzeczywiste) != sign(prognoza), abs(roznica), 0)

return(strata)

}

# Obliczenie funkcji strat dla każdego modelu (in-sample)

straty1 <- Funkcja_strat(premia, as.numeric(fitted(premia_model_arma1)))

straty2 <- Funkcja_strat(premia, premia_model_arma2$fit@fit$fitted.values)

straty <- cbind(straty1, straty2)

mcs_stationary <- mcsTest(straty, alpha = 0.05, nboot = 500, nblock = 10, boot = "stationary")

print(mcs_stationary)

$includedR
[1] 1 2

$pvalsR
     [,1]
[1,] 0.23
[2,] 1.00

$excludedR
NULL

$includedSQ
[1] 1 2

$pvalsSQ
     [,1]
[1,] 0.23
[2,] 1.00

$excludedSQ
NULL

Test MCS nie pozwolił wykluczyć żadnego z dwóch modeli, ale wskazuje dużo bardziej na zachowanie drugiego, tj. z autoarfima, więc go wybierzemy.

Wykres modelu

Wobec tego sprawdźmy jak te dopasowania dla najlepszego modelu wyglądają:

# Wartości dopasowane

dopas <- premia_model_arma2$fit@fit$fitted.values

okres <- c(1945:2024)

# Wykres najlepszego modelu z prawdziwymi obserwacjami

plot(x=okres, y=premia, type="l", ylab = "")

lines(x=okres, y=dopas, col="red")

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

abline(h=mean(premia), col="blue", lwd=2)

legend("topright", legend=c("premia za ryzyko S&P 500", "prognoza in-sample", "średnia premia za ryzyko"),

col=c("black","red", "blue"), lwd=2, cex=0.6)

Prognoza na 2025-2026

Niestety nie mam wiedzy, czy można od razu w rugarch robić prognozy mając już najlepszy model z autoarfima. jestem niemal pewny, że jest sposób, ale szybciej jest mi tak zrobić. Dlatego robię to "na piechotę". Obliczam arfimaspec(), z tego arfimafit i dopiero z tego arfimaforecast(). W arfimaspec() wpisujemy normalnie maksymalne rzędy ARMA, ale za to dodajemy argument fixed.pars = list(ar2=0, ar3=0), który właśnie kontroluje ustawienie konkretnych współczynników.

rzadAR <- premia_model_arma2$fit@model$modelinc[2]

rzadMA <- premia_model_arma2$fit@model$modelinc[3]

premia_spec <- arfimaspec(mean.model = list(armaOrder = c(rzadAR, rzadMA), include.mean = TRUE,

arfima = FALSE, external.regressors = NULL), distribution.model = "std",

start.pars = list()

, fixed.pars = list(ar2=0,

ar3=0,

ar4=0,

ar6=0,

ma2=0)

)

premia_fit <- arfimafit(spec = premia_spec, data = premia, solver = "nlminb")

prognoza <- arfimaforecast(premia_fit, n.ahead = 2)

prognoza

*----------------------------------*
*        ARFIMA Model Forecast     *
*----------------------------------*

Horizon: 2
Roll Steps: 0
Out of Sample: 0

0-roll forecast: 
    T+1     T+2 
0.05436 0.07451

Na rok 2025 dostajemy premię niecałe 5,4% i na 2026 7,5%. Dodajemy nowe prognozy do wykresu:

nowa_prognoza <- as.numeric(prognoza@forecast$seriesFor)

# Wykres najlepszego modelu z prawdziwymi obserwacjami

plot(x=c(okres, 2025, 2026), y=c(premia, NA, NA), type="l", ylab = "")

lines(x=c(okres, 2025, 2026), y=c(dopas, NA, NA), col="red")

lines(x=c(okres, 2025, 2026), y=c(rep(NA, length(premia)-1), tail(dopas, 1), nowa_prognoza), col="green", lwd=2, type="o", cex=0.4)

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

abline(h=mean(premia), col="blue", lwd=2)

legend("topright", legend=c("premia za ryzyko S&P 500", "prognoza in-sample", "prognoza poza próbą", "średnia premia za ryzyko"),

col=c("black","red", "green", "blue"), lwd=2, cex = 0.6)

Prognoza krocząca

Czerwone linie są tak naprawdę jedynie dopasowaniem do danych, a nie prognozą. Aby uzyskać bardziej realny obraz dobroci modelu, wykonam prognozy kroczące za pomocą arfimaroll(), czyli prognozy out-of-sample. Jest to odpowiednik ugarchroll bez procesu GARCH (zob. ten art.). Zobaczmy 10 ostatnich kroczących-nawijanych prognoz, tj. od roku 2015:

# Prognoza krocząca-nawijana

klaster <- makeCluster(detectCores() - 1)

premia_prognoza_roll <- arfimaroll(spec = premia_spec, data = premia, forecast.length = 10,

refit.every = 1, refit.window = "recursive",

window.size = NULL, solver = "nlminb", fit.control = list(ar2=0, ar3=0),

cluster=klaster)

premia_prognoza_roll <- resume(premia_prognoza_roll)

stopCluster(klaster)

prognoza_roll <- premia_prognoza_roll@forecast$density$Mu

okres_prognozy_roll <- c(2015:2024)

names(prognoza_roll) <- okres_prognozy_roll

prognoza_roll

# Wykres z prognozą kroczącą

plot(x=c(okres, 2025, 2026), y=c(premia, NA, NA), type="l", ylab = "", xlim=c(2015, 2026))

lines(x=c(okres, 2025, 2026), y=c(dopas, NA, NA), col="red")

lines(x=c(okres, 2025, 2026), y=c(rep(NA, length(premia)-1), tail(prognoza_roll, 1), nowa_prognoza), col="green", lwd=2, type="o", cex=0.4)

lines(x=c(NA, okres_prognozy_roll), y=c(tail(dopas, 1), prognoza_roll), col="green", lwd=2, type="o", cex=0.4)

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

abline(h=mean(premia), col="blue", lwd=2)

legend("topright", legend=c("premia za ryzyko S&P 500", "prognoza in-sample", "prognoza poza próbą", "średnia premia za ryzyko"),

col=c("black","red", "green", "blue"), lwd=2, cex = 0.6)

Prognoza sugeruje spadek premii w 2025 i może wzrost 2026. Zauważmy tylko, że spadek nie poniżej zera, tylko poniżej średniej. No i oczywiście nie mówi to nic na temat zmian w obligacjach.

Rozkład częstości i gęstości

Skupiając się na prognozowaniu łatwo jednak zgubić szeroki obraz, tj. ekonomiczne znaczenie premii za ryzyko. Dlatego musimy mieć w pamięci, że prawdopodobnie tak będzie, ale faktycznie może się zdarzyć niespodzianka. Ostatnie lata przyzwyczajają nie tylko zwykłych ludzi, ale i ekspertów do tego, że co roku z indeksów USA jest dodatnia stopa zwrotu. Od 2009 r do 2024, czyli przez 16 lat zdarzyły się zaledwie 2 ujemne roczne (całkowite) stopy zwrotu! Tworzy to złudzenie, że ryzyka nie ma, bo wystarczy przeczekać parę lat. Jednak ekonomia robi swoje i musi to jakoś wyrównać. Dlatego muszą nie tylko przydarzać się szokowe, bardzo silne i szybkie spadki (prędzej czy później). Muszą występować dłuższe okresy, w których obligacje skarbowe przyniosą większe zwroty od akcji. Gdyby nie ich nie było, to rzeczywiście nie byłoby ryzyka.

Nawet historyczna ujemna skośność (średnia < mediana) i kurtoza w rozkładzie premii o tym świadczy, zobaczmy:

# Zapisujemy stare ustawienia graficzne

old_par <- par(no.readonly = TRUE)

# Dajemy więcej miejsca po prawej (dół, lewo, góra, prawo)

par(mar = c(3.5, 5, 4, 5) - 0.1, mgp = c(2, 1, 0))

# 1. Obliczamy histogram, ale nie rysujemy

h <- hist(premia,

breaks = 10,

plot = FALSE)

# 2. Estymujemy KDE (Kernel Density Estimation)

d <- density(premia,

kernel = "gaussian",

bw = "nrd0")

# 3. Rysujemy histogram na lewą oś Y (counts)

plot(h,

freq = TRUE,

col = "lightgray",

border = "white",

xlab = "Premia",

ylab = "Częstość",

xlim = range(h$breaks, d$x),

ylim = c(0, max(h$counts) * 1.1),

main = "Histogram i gęstość rozkładu – lewa oś: częstość, prawa oś: gęstość")

# 4. Nałóż KDE na ten sam wykres

par(new = TRUE)

plot(d$x, d$y,

type = "l",

lwd = 2,

axes = FALSE,

xlab = "",

ylab = "",

ylim = c(0, max(d$y) * 1.1))

# 5. Dodajemy prawą oś Y (density)

axis(side = 4,

at = pretty(c(0, max(d$y) * 1.1)))

mtext("Gęstość (KDE)", side = 4, line = 3)

# 7. Przywracamy poprzednie ustawienia

par(old_par)

Jeżeli rynek jest efektywny w długim okresie, to powinniśmy się spodziewać, że podobny rozkład wystąpi dla okresów wieloletnich. Zobaczmy dla 4-letnich. Zagregujmy stopy zwrotu:

grupy <- ceiling(1:length(premia)/4)

premia_4lata <- exp(tapply(X = log(1+premia), INDEX = grupy, FUN = sum)) - 1

I powtarzamy kod dla wykresu rozkładu częstości i gęstości, zastępując zmienną premia przez premia_4lata. Dostajemy:

Rozkład gęstości 4-letnich premii zawiera niewiele danych, ale też występuje na nim skośność.

Dodatek dla premii za ryzyko nad bonami skarbowymi 3-miesięcznymi

Powtórzmy całą procedurę dla nadwyżki ponad bony skarbowe. To znaczy premia zawiera również ryzyko długoterminowych obligacji skarbowych w postaci zmienności stopy procentowej.

Stabilność średniej

W całym okresie średnia wyniosła aż 9%. Odwrotnie niż dla obligacji pierwsze 40 lat dały mniej, bo 8%, a kolejne prawie 10%. Czy więc tym razem to istotna zmiana średniej? Test Chowa wskazuje na stabilność:

	Chow test

data:  premia ~ 1
F = 0.22, p-value = 0.6

ACF i PACF

acf(premia, 15)

pacf(premia, 15)

Okazuje się, że jedynie ta 4-tego rzędu jest na granicy istotności (przekracza 0,2). To potwierdzałoby istnienie 4-letnich cykli, chociaż wykorzystanie tego fenomenu na pewno nie jest proste.

Modelowanie

Sprawdzamy zatem AR(4). Efekt:

Coefficients:
      ar1  ar2  ar3    ar4  intercept
        0    0    0  0.236      0.090
s.e.    0    0    0  0.112      0.024

sigma^2 estimated as 0.0282:  log likelihood = 29.05,  aic = -52.11

Training set error measures:
                      ME  RMSE    MAE    MPE  MAPE   MASE     ACF1
Training set 0.000006594 0.168 0.1362 -266.2 433.2 0.6704 -0.06222

Gdyby zastosować taki model, dostalibyśmy średnią 9% (+/- 2,4%) - faktycznie, tyle co średnia arytmetyczna. Dodatkowo współczynnik 4-tego rzędu 0,24 (+/- 0,11).

To teraz autoarfima z tymi samymi ustawieniami.

*----------------------------------*
*          ARFIMA Model Fit        *
*----------------------------------*
Mean Model	: ARFIMA(4,0,7)
Distribution	: std 

Optimal Parameters
------------------------------------
       Estimate  Std. Error     t value Pr(>|t|)
mu     0.098811    0.000012   8073.8245  0.00000
ar1    0.000000          NA          NA       NA
ar2    0.000000          NA          NA       NA
ar3    0.945906    0.000100   9503.1496  0.00000
ar4    0.054088    0.000008   6695.8385  0.00000
ma1    0.279539    0.000033   8524.8266  0.00000
ma2    0.000000          NA          NA       NA
ma3   -1.685603    0.000062 -27299.6027  0.00000
ma4   -0.135260    0.000014  -9448.6750  0.00000
ma5   -0.118209    0.000014  -8582.6486  0.00000
ma6    0.000000          NA          NA       NA
ma7   -0.036416    0.000006  -5762.9626  0.00000
sigma  0.111497    0.009553     11.6718  0.00000
shape 16.523609   12.901423      1.2808  0.20028

Robust Standard Errors:
       Estimate  Std. Error      t value Pr(>|t|)
mu     0.098811    0.000566   174.510530  0.00000
ar1    0.000000          NA           NA       NA
ar2    0.000000          NA           NA       NA
ar3    0.945906    0.010638    88.917365  0.00000
ar4    0.054088    0.000072   750.823674  0.00000
ma1    0.279539    0.002076   134.669756  0.00000
ma2    0.000000          NA           NA       NA
ma3   -1.685603    0.002733  -616.737764  0.00000
ma4   -0.135260    0.000331  -408.478534  0.00000
ma5   -0.118209    0.000610  -193.739688  0.00000
ma6    0.000000          NA           NA       NA
ma7   -0.036416    0.000020 -1799.605712  0.00000
sigma  0.111497    0.093611     1.191060  0.23363
shape 16.523609 1287.342190     0.012835  0.98976

LogLikelihood : 62.2351 

Information Criteria
------------------------------------
                    
Akaike       -1.3059
Bayes        -1.0081
Shibata      -1.3327
Hannan-Quinn -1.1865

Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
                         statistic p-value
Lag[1]                       3.101 0.07826
Lag[2*(p+q)+(p+q)-1][32]    17.518 0.04329
Lag[4*(p+q)+(p+q)-1][54]    27.653 0.46556

H0 : No serial correlation

Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
                        statistic p-value
Lag[1]                     0.3235  0.5695
Lag[2*(p+q)+(p+q)-1][2]    0.6261  0.6373
Lag[4*(p+q)+(p+q)-1][5]    1.6151  0.7116


ARCH LM Tests
------------------------------------
             Statistic DoF P-Value
ARCH Lag[2]     0.6869   2  0.7093
ARCH Lag[5]     2.5554   5  0.7681
ARCH Lag[10]    5.7786  10  0.8335

Nyblom stability test
------------------------------------
Joint Statistic:  411.572
Individual Statistics:             
mu    0.03948
ar3   0.04534
ar4   0.02661
ma1   0.04062
ma3   0.15471
ma4   0.02460
ma5   0.03678
ma7   0.02549
sigma 0.13303
shape 0.08077

Asymptotic Critical Values (10% 5% 1%)
Joint Statistic:     	 2.29 2.54 3.05
Individual Statistic:	 0.35 0.47 0.75

Dostałem model ARMA(4, 7), ale przy usuniętych ar1, ar2, ma2 i ma6.

Jednak całkowity model okazuje się niestabilny, bo Joint Statistic = 412. Za to poszczególne współczynniki są stabilne shape, czyli kształt rozkładu reszt, jest niestabilny. Wynika to po prostu z tego, że rozkład t-Studenta nie jest idealnym rozkładem tutaj.

Test MCS

$includedR
[1] 2

$pvalsR
     [,1]
[1,]    0
[2,]    1

$excludedR
[1] 1

$includedSQ
[1] 2

$pvalsSQ
     [,1]
[1,]    0
[2,]    1

$excludedSQ
[1] 1

Test MCS potwierdza jednoznacznie, że drugi model (uzyskany z autoarfima) jest lepszy, albowiem każe odrzucić pierwszy model.

Prognozy kroczące i dwie nowe

*----------------------------------*
*        ARFIMA Model Forecast     *
*----------------------------------*

Horizon: 2
Roll Steps: 0
Out of Sample: 0

0-roll forecast: 
   T+1    T+2 
0.0309 0.1279

Podobnie jak poprzednio - rok 2025 spadkowa premia.

Rozkład częstości i gęstości

Roczne premie:

4-letnie premie:

Porównując ze sobą obydwa rodzaje premii (nad bonami i nad obligacjami), to co się rzuca w oczy, to wpływ obligacji skarbowych na 4-letnią cykliczność stopy zwrotu. Po odjęciu bonów premia wykazuje pewną autokorelację 4-tego rzędu. Po odjęciu obligacji skarbowych premia ta autokorelacja mieści się w granicach nieistotności, a więc przypadkowości (jednak zwracam uwagę, że to PACF to taka klasyczna metoda szukania cykli, raczej wstępna). Wydaje się, że taka dysproporcja wynika z tego, że długoterminowe obligacje odzwierciedlają oczekiwania koniunktury, a ta przecież porusza się właśnie w takich kilkuletnich cyklach. Dostajemy więc pewnego rodzaju "dowód", że indeks S&P 500 podąża w jakimś zakresie w kierunku zmian gospodarczych.

Modele w zbiorze ufności w rugarch (język R)

2025-05-03T13:46:00.006+02:00

Jeśli powiem, że każda prognoza gospodarcza czy giełdowa jest okraszona niepewnością, to będzie to banał. Ale jeśli powiem, że celem prognozy powinno być nie tyle odgadnięcie przyszłej wartości, co minimalizacja jej niepewności, to już banałem nie jest. Dla uproszczenia założymy, że niepewność to przedział ufności dla prognozy naszego modelu. Im więcej informacji zawiera model, tym przedział ufności staje się mniejszy, zbliżając się tym samym do prognozy punktowej. Widzimy więc, że idea nie jest wcale banalna: należy zebrać jak najwięcej istotnych informacji, aby zawężać przedział prognozy, zmniejszając tym samym ryzyko.

Analogicznie do przedziału ufności, Hansen et al.* zaproponowali zbiór ufności dla wielu modeli: z całej puli testowanych modeli, poszukiwane są te najlepsze, które wpadają do zbioru ufności. Im więcej posiadamy informacji, tym mniejszy będzie ten zbiór, zbliżając się tym samym do jednego najlepszego modelu.

Standardową metodą wyboru najlepszego modelu jest użycie kryterium informacyjnego, np. AIC, BIC czy HQIC. Niestety nie uwzględniają one niepewności czy przypadkowości parametrów. Z tym się ciągle zmagamy. Przykładowo, mam zawsze wątpliwości co do używania modeli parametrycznych, ponieważ wydają się one sztuczne, realizując uśrednianie na siłę. Dlatego tak ważne są testy stabilności parametrów. Ale stabilność to dopiero początek, bo:

(1) jeśli spośród wielu modeli różnica w kryterium informacyjnym jest niewielka, to może ona wynikać po prostu z przypadku,

(2) dodanie tylko jednej obserwacji bardzo często zmienia wybór modelu.

W rugarch jest funkcja mcsTest (Model Confidence Set - MCS), która testuje, czy dany model należy zachować w zbiorze ufności czy go odrzucić. Ma ona następujące argumenty:

- losses: macierz lub ramka danych z wartościami strat dla każdego modelu (każda kolumna odpowiada jednemu modelowi). To w jaki sposób definiujemy straty, zależy od nas. W przypadku rozkładów symetrycznych, najczęściej stosuje się kwadraty odchyleń od prognozy lub wartości absolutne tych odchyleń. Dla rozkładów niesymetrycznych stosować można semi-odchylenia albo miary oparte na wartości zagrożonej.

- alpha: poziom istotności (domyślnie 0.05).

- nboot: liczba replikacji bootstrapowych (domyślnie 100). Bootstrap to procedura wielokrotnego losowania próbki ze zwracaniem. nboot to właśnie liczba tych losowań.

- nblock: długość bloku dla bootstrappu blokowego (domyślnie 1). Jeżeli obserwacje mogą być zależne od siebie, to bootstrap musi losować blok występujących po sobie x obserwacji. Jeżeli uważamy, że każda obserwacja jest niezależna od innej, to x = 1. Im rząd autokorelacji większy, tym blok powinien być naturalnie też większy.

- boot: rodzaj bootstrapu ("stationary" lub "block"). Chociaż w obydwu przypadkach mamy do czynienia z blokami, to metoda blokowa, "block", oznacza, że całą próbę dzielimy na stałą liczbę przedziałów, tak że każdy przedział składa się bloku obserwacji o długości wyznaczonej dokładnie przez nblock. Potem dopiero losujemy blok, powtarzamy to losowanie dokładnie nboot razy. Natomiast metoda stacjonarna, "stationary", tworzy bloki o losowej długości ze średnią równą nblock. Jeśli występuje autokorelacja w szeregu czasowym, to blok nie może być za krótki. Najczęściej "stationary" jest lepszy od "block".

Co właściwie robi ten test? W skrócie, kiedy mamy macierz strat dla każdego modelu, liczymy różnice między każdym modelem (prognozą każdego modelu) i są one odpowiednio uśredniane. Tworzymy bloki bootstrapowe i dla każdego bloku obliczamy średnie różnice między średnią stratą w bloku, średnią stratą dla całego modelu i średnią stratą ze wszystkich modeli. Dla każdego modelu będą różne wartości tworzące rozkład różnic. Wartości zbyt odchylone na plus (tzn. prawostronne ogony rozkładu) zostają odrzucone i tak uzyskujemy zbiór ufności.

Oczywiście najpierw sami budujemy funkcję strat. Ja użyłem definicji najbardziej dla mnie intuicyjnej: jeżeli znak prognozy różni się od znaku wartości rzeczywistej, to bierzemy moduł z różnicy między nimi. Jeżeli znak jest ten sam, to strata wynosi zero.

Przetestujmy tę funkcję. Będziemy używać mcsTest na zbiorze uczącym. Reguły używania poniższych funkcji w rugarch omówiłem szczegółowo tutaj.

Stwierdziłem wcześniej, że próba między 300 a 400 tygodniowych danych giełdowych dla GARCH wydaje się już w miarę odpowiednia. Zrobię symulację EGARCH o próbie równej 330. Powiedzmy, że prawdziwy proces to EGARCH(1, 1)-ARMA(3, 1). Obok niego mamy też fałszywe - wszystkie EGARCH: ARMA(0, 0), ARMA(1, 0), ARMA(2, 0), ARMA(3, 0). Zobaczmy które z nich zostaną wykluczone przez mcsTest, a które przyjęte do zbioru ufności. Zrobimy test w dwóch wariantach: blokowy i stacjonarny. Ustawienia: nblock = 20 oraz nboot = 1000.

library("rugarch")

library("xts")

# 1. Definicje stałych parametrów

stale_3_1 <- list(

omega = 1,

alpha1 = 0.3,

beta1 = 0.2,

mu = 0.0001,

ar1 = -1,

ar2 = 0.5,

ar3 = 0.7,

ma1 = 0.4,

gamma1 = 0.05

)

stale_0_0 <- list(

omega = 1,

alpha1 = 0.3,

beta1 = 0.2,

mu = 0.0001,

gamma1 = 0.05

)

stale_1_0 <- list(

omega = 1,

alpha1 = 0.3,

beta1 = 0.2,

mu = 0.0001,

ar1 = -1,

gamma1 = 0.05

)

stale_2_0 <- list(

omega = 1,

alpha1 = 0.3,

beta1 = 0.2,

mu = 0.0001,

ar1 = -1,

ar2 = 0.5,

gamma1 = 0.05

)

stale_3_0 <- list(

omega = 1,

alpha1 = 0.3,

beta1 = 0.2,

mu = 0.0001,

ar1 = -1,

ar2 = 0.5,

ar3 = 0.7,

gamma1 = 0.05

)

stale_0_1 <- list(

omega = 1,

alpha1 = 0.3,

beta1 = 0.2,

mu = 0.0001,

ma1 = 0.4,

gamma1 = 0.05

)

# 2. Specyfikacje modeli EGARCH(1,1) z różnymi specyfikacjami średniej:

spec_3_1 <- ugarchspec(

variance.model = list(model = "eGARCH", garchOrder = c(1,1)),

mean.model = list(armaOrder = c(3,1), include.mean = TRUE),

distribution.model = "norm",

fixed.pars = stale_3_1

)

spec_0_0 <- ugarchspec(

variance.model = list(model = "eGARCH", garchOrder = c(1,1)),

mean.model = list(armaOrder = c(0,0), include.mean = TRUE),

distribution.model = "norm",

fixed.pars = stale_0_0

)

spec_1_0 <- ugarchspec(

variance.model = list(model = "eGARCH", garchOrder = c(1,1)),

mean.model = list(armaOrder = c(1,0), include.mean = TRUE),

distribution.model = "norm",

fixed.pars = stale_1_0

)

spec_2_0 <- ugarchspec(

variance.model = list(model = "eGARCH", garchOrder = c(1,1)),

mean.model = list(armaOrder = c(2,0), include.mean = TRUE),

distribution.model = "norm",

fixed.pars = stale_2_0

)

spec_3_0 <- ugarchspec(

variance.model = list(model = "eGARCH", garchOrder = c(1,1)),

mean.model = list(armaOrder = c(3,0), include.mean = TRUE),

distribution.model = "norm",

fixed.pars = stale_3_0

)

spec_0_1 <- ugarchspec(

variance.model = list(model = "eGARCH", garchOrder = c(1,1)),

mean.model = list(armaOrder = c(0,1), include.mean = TRUE),

distribution.model = "norm",

fixed.pars = stale_0_1

)

# 3. Symulacja 330 obserwacji według poprawnego modelu (ARMA(3,1))

set.seed(765)

symulacja <- ugarchpath(spec_3_1, n.sim = 330)

dane_symulowane <- as.numeric(fitted(symulacja))

# 4. Dopasowanie modeli do całego szeregu (in-sample)

dopasowanie_3_1 <- ugarchfit(spec = spec_3_1, data = dane_symulowane)

dopasowanie_0_0 <- ugarchfit(spec = spec_0_0, data = dane_symulowane)

dopasowanie_1_0 <- ugarchfit(spec = spec_1_0, data = dane_symulowane)

dopasowanie_2_0 <- ugarchfit(spec = spec_2_0, data = dane_symulowane)

dopasowanie_3_0 <- ugarchfit(spec = spec_3_0, data = dane_symulowane)

dopasowanie_0_1 <- ugarchfit(spec = spec_0_1, data = dane_symulowane)

# 5. Zapisanie wyników dopasowania w liście (dla przejrzystości)

dopasowania <- list(

model_3_1 = dopasowanie_3_1,

model_0_0 = dopasowanie_0_0,

model_1_0 = dopasowanie_1_0,

model_2_0 = dopasowanie_2_0,

model_3_0 = dopasowanie_3_0,

model_3_0 = dopasowanie_0_1

)

print("Dopasowane modele:")

print(names(dopasowania))

# 6. Definicja funkcji strat:

Funkcja_strat <- function(rzeczywiste, prognoza) {

roznica <- rzeczywiste - prognoza

strata <- ifelse(sign(rzeczywiste) != sign(prognoza), abs(roznica), 0)

return(strata)

}

# 7. Obliczenie funkcji strat dla każdego modelu (in-sample)

strata_3_1 <- Funkcja_strat(dane_symulowane, fitted(dopasowanie_3_1))

strata_0_0 <- Funkcja_strat(dane_symulowane, fitted(dopasowanie_0_0))

strata_1_0 <- Funkcja_strat(dane_symulowane, fitted(dopasowanie_1_0))

strata_2_0 <- Funkcja_strat(dane_symulowane, fitted(dopasowanie_2_0))

strata_3_0 <- Funkcja_strat(dane_symulowane, fitted(dopasowanie_3_0))

strata_0_1 <- Funkcja_strat(dane_symulowane, fitted(dopasowanie_0_1))

# Łączenie strat w macierz – wiersze odpowiadają obserwacjom, kolumny modelom

straty <- cbind(

model_3_1 = strata_3_1,

model_0_0 = strata_0_0,

model_1_0 = strata_1_0,

model_2_0 = strata_2_0,

model_3_0 = strata_3_0,

model_3_0 = strata_0_1

)

# 8. Usunięcie wierszy składających się tylko z zer

straty <- straty[rowSums(straty != 0) > 0, ]

straty <- coredata(straty)

# 9. Test MCS

# 9.1. Wariant 1: stacjonarny

mcs_stationary <- mcsTest(straty, alpha = 0.05, nboot = 1000, nblock = 20, boot = "stationary")

# 9.2. Wariant 1: blokowy

mcs_block <- mcsTest(straty, alpha = 0.05, nboot = 1000, nblock = 20, boot = "block")

# 10. Wyświetlenie wyników:

cat("\nWyniki mcsTest (ustawienie nboot = 1000, nblock = 20, boot = stationary):\n")

print(mcs_stationary)

Wyniki mcsTest (ustawienie nboot = 1000, nblock = 20, boot = stationary):

$includedR
[1] 1 5

$pvalsR
      [,1]
[1,] 0.000
[2,] 0.000
[3,] 0.000
[4,] 0.004
[5,] 0.697
[6,] 1.000

$excludedR
[1] 6 2 4 3

$includedSQ
[1] 1 5

$pvalsSQ
      [,1]
[1,] 0.000
[2,] 0.000
[3,] 0.000
[4,] 0.003
[5,] 0.697
[6,] 1.000

$excludedSQ
[1] 6 2 4 3

Należy zauważyć, że mamy tu 2 typy testów: statystyka R (Range statistic) oraz SQ (Semi-Quadratic statistic). Pierwszy z nich skupia się na różnicach w maksymalnych stratach między modelami, a drugi na różnicach w kwadratach strat. Najlepiej, jeśli obydwa wskazują ten sam zbiór, bo wtedy dostajemy potwierdzenie, że zbiór ufności zawiera odpowiednie modele.

Jak widać obydwa testy rzeczywiście wskazują to samo: included 1 oraz 5, czyli zbiór ufności zawiera tylko dwa modele - pierwszy oraz piąty. Pozostałe modele zostają odrzucone (excluded), zarówno przez R jak SQ. Możemy powiedzieć, że test zadziałał poprawnie, mimo przyjęcia błędnego modelu ARMA(3, 0). Model ten jest podobny do prawidłowego ARMA(3, 1), a co więcej poprawny model ma największe p-value = 1 (piąty 0.697). O tym może jeszcze powiem. Otóż wydrukowana lista może być trochę myląca, bo kolejność nie pokazuje p-value w kolejności od modelu 1 do 6, tylko w kolejności od najbardziej istotnego do najmniej istotnego wyniku dotyczącego odrzucenia modelu ze zbioru ufności. Im mniejsze p-value, tym większe prawdopodobieństwo, że dany model należy wyrzucić ze zbioru.

Żeby zawęzić zbiór ufności, należy zwiększyć próbę. Po zwiększeniu jej do 1600, dostałem:

Wyniki mcsTest (ustawienie nboot = 1000, nblock = 20, boot = stationary):

$includedR
[1] 1

$pvalsR
      [,1]
[1,] 0.000
[2,] 0.000
[3,] 0.000
[4,] 0.000
[5,] 0.047
[6,] 1.000

$excludedR
[1] 6 2 4 3 5

$includedSQ
[1] 1

$pvalsSQ
      [,1]
[1,] 0.000
[2,] 0.000
[3,] 0.000
[4,] 0.000
[5,] 0.047
[6,] 1.000

$excludedSQ
[1] 6 2 4 3 5

Przy tak dużej próbie test nie ma już wątpliwości, że jedyny prawidłowy model to ten pierwszy.

Bardzo podobne wyniki dostałem dla boot = "block" i nie ma sensu pokazywać ich. Podobnie jak dla "stationary" musiałem ustawić próbę liczną 1600, aby uzyskać wykluczenie wszystkich fałszywych modeli. Co jednak zauważyłem, to to, że w przypadku "block" wykluczenie modelu 5 było ledwo, ledwo, na styk, bo p-value wyniosło 0.049 (dla stationary też, ale niższe, bo 0.047). Stacjonarny test generalnei jest lepszą alternatywą.

Powyższe przykłady zakładały rozkład normalny reszt modelu. Dla rynków może być lepszy inny, np. skośny rozkład t-Studenta. W takiej sytuacji potrzebne są dwa dodatkowe parametry do stałych: skew oraz shape. Zamiast "norm" będzie "sstd". Przykładowo kod będzie taki:

# 1. Definicje stałych parametrów

stale_3_1 <- list(

omega = 1,

alpha1 = 0.3,

beta1 = 0.2,

mu = 0.0001,

ar1 = -1,

ar2 = 0.5,

ar3 = 0.7,

ma1 = 0.4,

gamma1 = 0.05,

skew = 1.5,

shape = 5

)

stale_0_0 <- list(

omega = 1,

alpha1 = 0.3,

beta1 = 0.2,

mu = 0.0001,

gamma1 = 0.05,

skew = 1.5,

shape = 5

)

stale_1_0 <- list(

omega = 1,

alpha1 = 0.3,

beta1 = 0.2,

mu = 0.0001,

ar1 = -1,

gamma1 = 0.05,

skew = 1.5,

shape = 5

)

stale_2_0 <- list(

omega = 1,

alpha1 = 0.3,

beta1 = 0.2,

mu = 0.0001,

ar1 = -1,

ar2 = 0.5,

gamma1 = 0.05,

skew = 1.5,

shape = 5

)

stale_3_0 <- list(

omega = 1,

alpha1 = 0.3,

beta1 = 0.2,

mu = 0.0001,

ar1 = -1,

ar2 = 0.5,

ar3 = 0.7,

gamma1 = 0.05,

skew = 1.5,

shape = 5

)

stale_0_1 <- list(

omega = 1,

alpha1 = 0.3,

beta1 = 0.2,

mu = 0.0001,

ma1 = 0.4,

gamma1 = 0.05,

skew = 1.5,

shape = 5

)

# 2. Specyfikacje modeli EGARCH(1,1) z różnymi specyfikacjami średniej:

spec_3_1 <- ugarchspec(

variance.model = list(model = "eGARCH", garchOrder = c(1,1)),

mean.model = list(armaOrder = c(3,1), include.mean = TRUE),

distribution.model = "sstd",

fixed.pars = stale_3_1

)

spec_0_0 <- ugarchspec(

variance.model = list(model = "eGARCH", garchOrder = c(1,1)),

mean.model = list(armaOrder = c(0,0), include.mean = TRUE),

distribution.model = "sstd",

fixed.pars = stale_0_0

)

spec_1_0 <- ugarchspec(

variance.model = list(model = "eGARCH", garchOrder = c(1,1)),

mean.model = list(armaOrder = c(1,0), include.mean = TRUE),

distribution.model = "sstd",

fixed.pars = stale_1_0

)

spec_2_0 <- ugarchspec(

variance.model = list(model = "eGARCH", garchOrder = c(1,1)),

mean.model = list(armaOrder = c(2,0), include.mean = TRUE),

distribution.model = "sstd",

fixed.pars = stale_2_0

)

spec_3_0 <- ugarchspec(

variance.model = list(model = "eGARCH", garchOrder = c(1,1)),

mean.model = list(armaOrder = c(3,0), include.mean = TRUE),

distribution.model = "sstd",

fixed.pars = stale_3_0

)

spec_0_1 <- ugarchspec(

variance.model = list(model = "eGARCH", garchOrder = c(1,1)),

mean.model = list(armaOrder = c(0,1), include.mean = TRUE),

distribution.model = "sstd",

fixed.pars = stale_0_1

)

# 3. Symulacja 330 obserwacji według poprawnego modelu (ARMA(3,1))

set.seed(55)