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

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)]
> fr_max
[1] 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.  lata i rok to jednostka czasu. Czyli w jednym roku mieści się pół cyklu, tj. 1/2 - odwrotność okresu.  (bo ta ostatnia mówi jak dużo cykli mieści się średnio w jednym roku, a więc odwrotność pokazuje ile lat mieści się w jednym średnim cyklu), oznaczmy cyc:

> cyc <- 1/fr_max
> cyc
[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 czysto techniczne ustawienie formatu wyniku (tu zawsze wybieramy "Sos", skrót od Second-order sections). 

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 2*f_cut / fs, 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 dodam, że trochę zgaduję, bo może być tak, że to tylko przyjęta konwencja konstrukcji tej funkcji.

Ż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ą miec 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 [8]

d_peak <- abs(fr_max - fr_s)

# 6. Wyznaczenie finalnej częstotliwości granicznej pasma zaporowego dla projektu filtra [9, 10]

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_max*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, która jest lepsza dla takich celów (zob. np. tu).

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

Stosujemy kod:

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

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

Wg modelu prognoza na 2026 r. wynosi  5.63 i na 2027 5.72. Żeby porównać ze wskazanymi na początku 20% i 10%, musielibyśmy retrasformować uzyskane liczby exp(5.63) = 279 i exp(5.72) = 305. W 2025 mamy 244.5, stąd 279 / 244.5 - 1 = 14%. W 2026 305 / 279 = 9,3%. Czyli dostajemy tutaj nieco mniejsze prognozy, 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. Ta dodatkowa zmienność może jeszcze podnieść wartość. Stąd uzyskany model potwierdza przewidywania J.P. Morgan i Goldman Sachs. Taki wynik zwiększa jego wiarygodność, a to z kolei uwiarygadnia prognozy na późniejsze lata - po 2027. Na dziś możemy więc spodziewać się spowolnienia od roku 2029. Byłby to okres trochę podobny do roku 2000.