środa, 1 stycznia 2025

Jak giełda powinna reagować na praworządność w Polsce?

Od wyborów parlamentarnych w 2015 do kolejnych w 2023 r. WIG urósł zaledwie 31%. W tym samym czasie nominalny PKB wzrósł o 88% (źródło). Realny PKB wzrósł o 33% (źródło). Wynika z tego, że ceny wzrosły w sumie o 55%. Inwestowanie w nasze spółki było zatem przez 8 lat rządów PISu nieopłacalne. Pomyślmy czysto logicznie: skoro PKB w sumie wzrósł prawie 3 razy więcej niż główny indeks, który w jakimś stopniu powinien odzwierciedlać stan gospodarki, to znaczy, że cały koszt, że tak powiem, musiał pochodzić od innych czynników niż gospodarka. Takim czynnikiem może być strach przed wojną z Rosją, który eksplodował w 2022 r. Byłby to dobry argument w rękach obrońców PISu, gdyby nie to, że można łatwo udowodnić, że wojna była jedynie pretekstem do przeceny. Jak bowiem wyjaśnić, że Litwa, która teoretycznie jest w gorszej pozycji niż Polska - chociażby dlatego, że stoi na przeszkodzie Rosji do połączenia z obwodem królewieckim (oczywiście upraszczam, bo jeszcze jest Białoruś, ale trudno nazwać ją przeszkodą) - praktycznie nie spadła w 2022? Łotwa niewiele gorzej od Litwy, za to w całym okresie PISu urosła mocniej. Popatrzmy na indeksy w skali logarytmicznej:



Litwa w całym okresie poszła w górę o 94%, a Łotwa 134% (w skali liniowej). 

Skoro wojna nie może być odpowiedzią na słabe zachowanie polskiej giełdy, to trzeba szukać w innych obszarach. Takimi obszarami są praworządność, stopień demokratyczności oraz poziom populizmu. Wszystkie one prawdopodobnie miały pewien wpływ, ponieważ ostatnie rządy PISu powszechnie uważa się za populistyczne i mało demokratyczne (tzn. dążące do centralizacji władzy), a przez większość prawników także za łamiące prawo (źródło). Kwestię populizmu w kontekście wpływu na gospodarkę omawiałem tutaj, praworządności / demokratyczności na rynek akcji tutajtutaj, i na gospodarkę tutaj.

Kiedy władzę przejęła znowu PO, od początku byłem wstrzemięźliwy co do dalszej przyszłości naszego rynku, bo widziałem, że politycy z tego obozu przejęli populistyczny język PISu. Szczególnie narażone stały się spółki skarbu państwa, a te jak wiadomo dużo ważą w indeksie. 

Okazało się, że to dopiero początek. Rok temu pisałem o bezprawnym, siłowym przejęciu TVP i PR, nawet zrobiłem analizę prawną (zob. tu). Był to pewien szok, ale giełda przyjęła tę akcję bardzo spokojnie, można by pomyśleć, że nie będzie miała na nią wpływu. Nie takie rzeczy giełda przełykała - pamiętam Smoleńsk, miałem obawy, że będzie tąpnięcie w poniedziałek, ale reakcja była słaba i chwilowa. Jest jednak różnica między tymi zdarzeniami. Smoleńsk był zdarzeniem przypadkowym, a przejęcie mediów publicznych było zaplanowane z założeniem, że cel uświęca środki.  

Ktoś mógłby powiedzieć, że przecież już wcześniej złamali Konstytucję, powołując 2 sędziów TK na zapas w 2015 r. Tylko że wtedy szybko się wycofali pod wpływem bezpiecznika, jakim był TK. Dziś tego bezpiecznika praktycznie nie ma. Wtedy PIS sprytnie wykorzystał tę zagrywkę jako pretekst do wprowadzenia swoich dodatkowych sędziów, nazywanych potem "dublerami", rozpoczynając tym samym kryzys TK. 

I teraz dochodzimy do meritum. PIS wykorzystał nielegalność przejęcia mediów publicznych do uznania wszystkich innych kontrowersyjnych działań rządu za nielegalne. Co więcej, udało im się narzucić narrację, że w Polsce narasta chaos prawny, tak że obywatel nie może mieć pewności czy wydany wyrok w jego sprawie jest wiążący. Dla przedsiębiorców czy inwestorów to duży problem, bo 3 razy się zastanowią, czy podjąć jakieś decyzje. Jest to zachęta do nieuczciwych praktyk czy nawet oszustw, bo oszust wiedząc, że sprawa przez lata nie zostanie rozwiązana, chętniej podejmie ryzyko.

Z drugiej strony PIS już raz ośmieszył się ze sprawą Kamińskiego i Wąsika. Dla każdego normalnie myślącego człowieka jest jasne, że nie można ułaskawić kogoś, kto nie został uznany prawomocnie za winnego. Ośmieszyli się także dlatego, że Sąd Najwyższy rozstrzygnął ten spór ostatecznie (zob. tu). Pytanie czy dalej się ośmieszają? W zasadzie odpowiedziałem - udało im się narzucić narrację, więc nie ośmieszyli się w tym wypadku. Ale nie dlatego, że się poprawili, tylko dlatego, że sprawa jest bardziej skomplikowana prawnie, a przez to łatwiej nią manipulować.

Weźmy tylko dwa przykłady z brzegu. Mamy postanowienie SN (zob. II CSKP 556/22) potwierdzające, że orzeczenia tzw. neosędziów należy uznać za nieistniejące, ponieważ sędziowie je wydający byli wadliwie obsadzeni. Żeby nie było - art. 187 p. 1 Konstytucji zawiera lukę, którą PIS wykorzystał. Chodzi o zdanie:  Krajowa Rada Sądownictwa składa się z: (...) piętnastu członków wybranych spośród sędziów Sądu Najwyższego, sądów powszechnych, sądów administracyjnych i sądów wojskowych. Nie ma w tym zdaniu nic o tym, kto wybiera tych sędziów. A ponieważ nie ma, to - wydawałoby się słusznie - powinna to precyzować jakaś ustawa (a tę PIS bez żadnych problemów uchwalił i wprowadził w życie). Dopóki tej interpretacji nie podważył SN, krzyki o nielegalność KRS były jedynie opiniami prawnymi mieszającymi się z opiniami politycznymi. SN wyszedł od spostrzeżenia, że skoro zgodnie z art. 186 Konstytucji KRS stoi na straży niezależności sądów i niezawisłości sędziów, to sama musi być niezależna od władzy ustawodawczej i wykonawczej. A zatem jej członkowie nie mogą być wybierani przez którąś z tych władz. Wprawdzie są wyjątki od tej reguły (w art. 187 też jest o nich mowa), ale wyjątki są właśnie wyjątkami, gdy nie są rozszerzane na regułę. Podobnie jest z sędziami TK, którzy nawet w całości są wybierani przez sejm, bo taki właśnie wyjątek wprowadzono. Ustawa nie może sama narzucać tego wyjątku. Niemniej, można spotkać głosy, podważające stanowisko SN i zapytać, dlaczego z góry należy zakładać, że sędzia wybierany przez sejm czy senat będzie nieobiektywny, skoro sędzia zasadniczo jest nieusuwalny? Stąd, powtarzam, mógłbym przyjąć interpretację PISu za prawidłową i wiążącą, gdyby nie wyrok SN, a także TSUE i ETPC (trzeba więc też brać pod uwagę prawo międzynarodowe), które uchwała SN zresztą przywołuje. Po to właśnie są sądy, aby ograniczać możliwości władzy ustawodawczej i wykonawczej. 

Drugi przykład jest bardziej zagmatwany. Chodzi o sprawę prokuratora krajowego, który wg PISu zasiada obecnie bezprawnie. Mamy tu dwie warstwy. Pierwsza to pozbycie się prokuratora krajowego od Ziobry, a druga to powołanie nowego. PIS próbuje narzucić narrację, że ich prokurator krajowy był nielegalnie odwołany, a nowy nielegalnie powołany. Strona rządowa natomiast twierdzi, że ten pierwszy został wcześniej (za PISu) wadliwie powołany, a więc nie mógł być też później odwołany. Dokładniej, politycy wygłaszają taką oto formułkę, żeby przeciętny Kowalski mógł zrozumieć: prokurator krajowy nie został w ogóle przywrócony ze stanu spoczynku, bo zastosowano błędne przepisy. Niestety, diabeł tkwi w szczegółach. Bo to zdanie jest tak proste, że wydawałoby się, że PIS czegoś nie zauważył i popełnił oczywisty błąd. W rzeczywistości, żeby do takiej konkluzji dojść, trzeba przeprowadzić skomplikowaną analizę prawną. Ministerstwo Sprawiedliwości zamówiło 3 opinie prawne (zob. tu) i warto przeczytać choćby opinię prof. G. Kucy, żeby zrozumieć z czym mamy do czynienia. Z ustawy o prokuraturze nie wynika wcale, że PIS powołał nieprawidłowo prokuratora krajowego. Wynika to ze szczegółowej analizy prawno-logicznej. Ale ta analiza została poddana krytyce przez Izbę Karną Sądu Najwyższego (zob. tu), która przyjęła uchwałę, zgodnie z którą PISowski prokurator krajowy został przywrócony do służby skutecznie. Normalnie rząd znalazłby się w ogromnych kłopotach po tym orzeczeniu, ale w pewnym sensie uratował go sam PIS. Art. 179 Konstytucji stwierdza, że to KRS wybiera sędziów do SN. I tu wracamy do przykładu pierwszego - to neo-sędziowie z KRS wskazali sędziów do tejże Izby Karnej SN. Tym samym - logicznie biorąc - sami sędziowie wybrani do tej Izby byli wadliwie wybrani i odziedziczyli z automatu status neo-sędziów. Tak więc PO ma mocny pretekst do podważania ich orzeczeń i wyroków, niezależnie od tego czy mają oni rację czy nie.

No ale jest jeszcze druga warstwa. PIS stawia sprawę jasno - POwski prokurator krajowy nie jest żadnym prokuratorem krajowym, bo jego powołania wymaga opinia prezydenta. Rzeczywiście, art. 14 par. 1 ustawy o prokuraturze brzmi:

Prokuratora Krajowego (...) powołuje spośród prokuratorów Prokuratury Krajowej i odwołuje z pełnienia tych funkcji Prezes Rady Ministrów na wniosek Prokuratora Generalnego. Prokuratora Krajowego oraz pozostałych zastępców Prokuratora Generalnego powołuje się po uzyskaniu opinii Prezydenta Rzeczypospolitej Polskiej, a odwołuje za jego pisemną zgodą.

Zwróćmy jednak uwagę, że artykuł nie mówi o pisemnej opinii. Oznacza to, że dowolna opinia prezydenta  na ten temat może posłużyć za podkładkę do powołania prokuratora. W dodatku nawet gdybyśmy przyjęli, że opinia musi być formalna, to również w tej kwestii SN się wypowiedział (zob. tu) i zauważył, że po zwróceniu się przez premiera o opinię ws. D. Korneluka na stanowisko prokuratora krajowego, prezydent odpowiedział 14 dni później, że Korneluk nie został skutecznie powołany na to stanowisko z uwagi na brak jego opinii. Można więc to czytać jako opinię negatywną, która jednak nie ma żadnego znaczenia oprócz formalnego. Jest to biurokratyczna bzdura, żeby uzyskanie opinii było obligatoryjne, ale nie miało żadnego znaczenia dla decyzji premiera. A to prowadzi do prostego jej ominięcia - cokolwiek A. Duda odpisze, będzie traktowane jako opinia. Co więcej, SN doszedł do wniosku, że nawet gdyby Duda nic nie odpisał, to też należy to traktować jako opinię. Ważne, że została wysłana prośba o nią od premiera. I trzeba przyznać, że ma to sens - skoro mamy do czynienia z biurokratyczną bzdurą, to należy ją traktować jak na to zasługuje.   

Sporów jest więcej, ale te dwa przykłady są chyba najważniejsze i dowodzą, że to strona rządowa ma przewagę w argumentacji. Cały problem leży w komunikacji. Ma on dwojaką naturę. Z jednej strony wydaje się, że politycy koalicji rządzącej nie potrafią narzucić swojej narracji w sposób przekonywający dla ogółu. Częściowo wynika to z faktycznie bezprawnego przejęcia mediów publicznych, które opozycja wykorzystuje do rozszerzenia na pozostałe działania. Polityk PO nie może przecież powiedzieć: przyznaję, że to akurat było nielegalne, ale cała reszta już zgodna z prawem. Z drugiej strony mediom komercyjnym zależy na podsycaniu konfliktu, bo to zawsze będzie oglądane. To prowadzi zazwyczaj do politycznej jatki zamiast merytorycznej dyskusji. Ponadto z jakiegoś powodu brakuje (chyba, bo może nie mam o tym wiedzy) dyskusji między ekspertami prawnymi o odmiennych poglądach. Jeśli coś takiego jest, to raczej są to niepogłębione, krótkie wywiady, w których nie daje się szansy na pełną odpowiedź na krytykę. Takie dyskusje powinny być organizowane przez media publiczne właśnie, ale wydaje się to niemożliwe, gdy są one pod pełną kontrolą jednej strony politycznej.

W sumie, subiektywnie praworządność w Polsce może się wydawać nadal łamana, ale obiektywnie biorąc na dziś mamy do czynienia z jej poprawą w stosunku do tego co było rok temu i wcześniej. Potwierdzają to statystyki opracowane przez World Justice Project:


Między 2021 a 2023 byliśmy na 36 miejscu, w 2024 skoczyliśmy na 33.

Miejsce 33 to ciągle daleko za Litwą (18) i Łotwą (21). Co chwilę też wybuchają nowe spory prawne. Pewne jest jedno - dopóki sprawa neo-sędziów i TK nie zostanie w pełni uregulowana (rozwiązana), dopóty nie można mówić o przywróceniu praworządności.  

Jak to wpłynie na giełdę? Może być tak, że negatywne czynniki instytucjonalne zostaną zastąpione przez negatywne czynniki gospodarcze. Zakładając, że działa tu psychologia, czyli subiektywizm, możemy mieć do czynienia z przesadnymi spadkami w kolejnych miesiącach, a gdy inwestorzy zauważą, że chaos prawny jest sztucznie wykreowany w mediach, nastąpi ożywienie. 

niedziela, 8 grudnia 2024

Prognoza krocząca z ruchomym oknem w rugarch (język R)

 Tytuł brzmi koślawo, bo słowo "krocząca" jest synonimem do "ruchomej" i np. średnia krocząca i średnia ruchoma oznaczają to samo. Z tym że ta średnia oznacza średnią z ostatnich x okresów, a prognoza jest zawsze do przodu. W ogóle określenie prognoza krocząca nie jest najlepsze i powinno się raczej mówić o prognozie nawijanej (rolującej) lub rekurencyjnej. W języku angielskim mówi się o recursive forecast. Natomiast rolling forecast zwykło się nazywać prognozę kroczącą z ruchomym oknem będącą analogią do średniej ruchomej. Tym przypadkiem zajmę się obecnie.

Krótko na temat samej funkcji w rugarch. Interesuje nas postać:

ugarchroll(spec, data, forecast.length = 500, n.start = NULL, refit.every = 25, refit.window = c("recursive", "moving"), window.size = NULL, solver = "hybrid", cluster = NULL)

gdzie:

spec - teoretyczny model GARCH,

data - dane najlepiej w xts,

forecast.length - długość prognozy dla próby testowej (tj. out-of-sample), która liczona jest do końca całej próby,

n.start - można użyć zamiast forecast.length, zamiast liczyć od końca można po prostu określić punkt startowy dla prognozy,

refit.every - określa, co którą obserwację w próbie testowej ponownie estymować; jeśli więc chcemy aby dla każdej następnej obserwacji reestymować, wpisujemy 1,

refit.window - wybór czy reestymować od początku próby do ostatniej dostępnej obserwacji ("recursive"), tzn. wybieramy tu prognozę nawijaną na kolejny okres, czy też reestymować w oknie o stałej liczebności, czyli tak jak np. średnia 30-okresowa ("moving"),

window.size - jeżeli wybraliśmy refit.window = "moving", wtedy ustawiamy tu długość okna, w którym estymowane są parametry; w przeciwnym wypadku dajemy NULL,

solver - zazwyczaj wybieram "hybrid", bo jest najbardziej wszechstronny,

cluster - również tu stosujemy przetwarzanie równoległe, aby przyspieszyć obliczenia (na podstawie pakietu parallel).


W tym badaniu ustawiam forecast.length = 30, refit.every = 3, refit.window = "moving", window.size = 300.

Temat ogólnie stanowi kontynuację prognozy za pomocą pakietu rugarch. Stanęło na tym, że modelowanie na dużej próbie jest utrudnione z uwagi na niestabilność parametrów wariancji warunkowej - przynajmniej dla tygodniowych stóp zwrotu. Wychwycił to specjalny test na proces TV-GARCH. Od teraz zawsze będę stosował ten test w pierwszej kolejności.

W końcu udało mi się znaleźć w miarę optymalną próbę - nie za małą i nie za dużą dla GARCH - jest to próba między 300 a 350 tygodniowych stóp zwrotu. Wybrałem 330, z czego ostatnie 30 będą out-of-sample. 

Tym razem będę modelował WIG zamiast WIG20, a potem S&P 500 (źródło stooq.pl). Wszystkie potrzebne pakiety:

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

  install.packages("tidyverse")

}

library("dplyr")

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

  install.packages("xts")

}

library("xts")

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

  install.packages("rugarch")

}

library("rugarch")

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

  install.packages("foreach")

}

library("foreach")

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

  install.packages("doParallel")

}

library("doParallel")

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

  install.packages("tvgarch")

}

library("tvgarch")


WIG

#Przygotowanie danych:
nazwaPliku = "wig_d.csv"
Obs <- 331
plik <- read.csv(nazwaPliku)
plik$Data <- as.Date(plik$Data)
tabela <- tail(to.weekly(as.xts(plik), indexAt="endof", OHLC=FALSE), Obs)
cena <- tabela$Zamkniecie
stopaZwrotu <- diff.xts(log(cena))
tabela$stopaZwrotu <- stopaZwrotu
tabela <- na.omit(tabela)
stopaZwrotu <- tabela$stopaZwrotu


Tak jak powiedziałem, zaczynam od testowania TV-GARCH. Testuję ostatnie 330 stóp zwrotu WIG (zmienna stopaZwrotu):


# Sprawdzamy TV-GARCH
tvGarchTest <- tvgarchTest(y = as.numeric(stopaZwrotu))
tvGarchTest

Testing GARCH(1,1), i.e., the model under H0, against 
  TV-GARCH(1,1), i.e., the model under H1: 

Estimation results for model under H0:
Model: GARCH(1,1) 

            intercept.h  arch1  garch1
Estimate:     0.0002362 0.2062 0.51256
Std. Error:   0.0003463 0.1080 0.03947
                     
Log-likelihood: 723.2


Transition variable in TV-GARCH(1,1): time 

Results from the Non-Robust TR^2 Test: 

                 NonRobTR2 p-value
H0:B3=B2=B1=0       2.4250  0.4890
H03:B3=0            1.4251  0.2326
H02:B2=0|B3=0       0.7273  0.3938
H01:B1=0|B3=B2=0    0.2775  0.5984

Results from the Robust TR^2 Test: 

                 RobTR2 p-value
H0:B3=B2=B1=0     3.749  0.2898
H03:B3=0          1.122  0.2895
H02:B2=0|B3=0     2.245  0.1340
H01:B1=0|B3=B2=0  1.370  0.2419

                                Single
No. of locations (alpha = 0.05)      0


Żaden test nie sugeruje nawet, że może wystąpić niestabilność. 


Jeśli chodzi o ogólną postać modelu wariancji, to przyjmuję EGARCH(1, 1), która okazywała się dotychczas najlepszą pod względem dopasowania. Kłopot zaczyna się dla części ARMA(p, q). Gdybym wiedział od początku jaką postać powinienem przyjąć, wtedy stosuję dwie funkcje: ugarchspec() oraz ugarchroll(). Nie jest już potrzebna ugarchfit(), którą częściowo zastępuje ugarchroll(). A nie jest potrzebna, bo normalnie służy nie tylko do utworzenia modelu empirycznego, ale i też do diagnostyki modelu. A skoro uznajemy, że prawidłową postać modelu znam od początku, to nie jest konieczna diagnostyka. Ale nawet jeżeli nie znam modelu, a interesuje mnie jedynie prognoza out-of-sample, to diagnostyka także jest zbędna. Może inaczej - nie jest niezbędna. Można jej użyć, ale tylko żeby się dowiedzieć, czy model jest dobrze dopasowany do próby in-sample.

Funkcja ugarchfit() przyda się jednak do prognozowania naprawdę przyszłych, czyli nieistniejących jeszcze danych. To znaczy do samej prognozy użyję ugarchforecast(), której argumentem jest ugarchfit(). 

W sumie, tak jak poprzednio, ugarchroll() posłuży do testowania prognozy out-of-sample, a ugarchfit() + ugarchforecast() do prognozy na kolejne okresy.

Kryterium wyboru najlepszego modelu będzie się tym razem opierać wyłącznie na out-of-sample: kierunkowej dokładności (ang. directional accuracy criterion - DAC) oraz na sumie odchyleń (SO) od średniej prognozy. DAC musi być jak największe, SO jak najmniejsze. Zwrócę uwagę, że nie będę korygował tych miar o liczbę współczynników modelu, jak to robią AIC, BIC i HQIC, dlatego należy traktować to badanie ad hoc.

Ostatnia sprawa to rozkład reszt. Tym razem dla przyspieszenia obliczeń przyjmę rozkład normalny, chociaż zdaję sobie sprawę, że mogę dostać przez to niedopasowaną wariancję warunkową. Jednak będzie to punkt wyjścia do dalszych analiz. 

rozklad <- "norm"
hqc <- list()
model_Garch <- "eGARCH"
submodel_Garch <- NULL
wyrzuc <- 30
q <- 1
p <- 1

# Ustawienie równoległego backendu do użycia wielu procesorów
ileRdzeni <- detectCores() - 1  # Użyj jednego mniej niż całkowita liczba rdzeni
klaster <- makeCluster(ileRdzeni)
clusterExport(klaster, varlist = c("stopaZwrotu", "rozklad", "model_Garch", "submodel_Garch", "q", "p", "wyrzuc"))
clusterEvalQ(klaster, library(rugarch))

# Utworzenie siatki parametrów dla pętli
siatka_parametrow <- expand.grid(k = 0:1, i = 0:20, j = 0:20)

# Iteracja po siatce parametrów bez foreach
wyniki <- lapply(1:nrow(siatka_parametrow), function(param) {
  k <- siatka_parametrow$k[param]
  i <- siatka_parametrow$i[param]
  j <- siatka_parametrow$j[param]
  
  rzad_modelu <- sprintf("ARMA(%d, %d)%s", i, j, ifelse(k == 1, "-ARCH_In_Mean", ""))
  
  tryCatch({
    prognoza <- NULL
    
    # Define the GARCH specification
    mojGarch <- ugarchspec(variance.model = list(model = model_Garch, submodel = submodel_Garch,
                                                 garchOrder = c(q, p)), 
                           mean.model = list(armaOrder = c(i, j), archm = as.logical(k)), 
                           distribution.model = rozklad)
    
    # Try ugarchroll and handle potential errors by attempting to use resume()
    prognoza <- try(
      ugarchroll(spec = mojGarch, data = stopaZwrotu, forecast.length = wyrzuc,
                 n.start = NULL, refit.every = 3, refit.window = "moving",
                 window.size = 300, solver = "hybrid", calculate.VaR = TRUE, 
                 VaR.alpha = c(0.01, 0.05),
                 cluster = klaster),
      silent = TRUE
    )
    
    # Check if ugarchroll returned an error
    if (inherits(prognoza, "try-error")) {
      message("Błąd w ugarchroll, przechodzę do resume(): ", attr(prognoza, "condition")$message)
      return(list(error = attr(prognoza, "condition")$message, rzad_modelu = rzad_modelu))
    }
    
    # Use resume() if no error occurred
    prognoza <- resume(prognoza)
    
    # Extract forecasts and calculate metrics
    prognozaSrednia <- prognoza@forecast$density$Mu # średnia prognoza
    prognozaSigma <- prognoza@forecast$density$Sigma # warunkowe odchylenie st
    stopaZwrotuZgodna <- prognoza@forecast$density$Realized # faktyczne stopy zwrotu
    odchylenie <- sum((prognozaSrednia - stopaZwrotuZgodna)^2)
    
    # Initialize DAC value
    dac_wartosc <- NA
    
    # Skip DACTest for ARMA(0, 0) with k = 0
    if (i + j + k > 0) {
      dac_wartosc <- DACTest(forecast = as.numeric(prognozaSrednia), 
                                     actual = as.numeric(stopaZwrotuZgodna), 
                                     test = "PT")$DirAcc  # Liczenie DAC dla prognozy
    }
    
    return(list(
      rzad_modelu = rzad_modelu,
      dac = dac_wartosc,
      odchylenie = odchylenie
    ))
    
  }, error = function(e) {
    # Catch block for unexpected errors
    message("Błąd w ARMA(", i, ", ", j, ") z k = ", k, ": ", conditionMessage(e))
    return(list(error = conditionMessage(e), rzad_modelu = rzad_modelu))
  })
})


# Zatrzymaj klaster po przetwarzaniu
stopCluster(klaster)


dac <- list()
odchylenie <- list()

# Przetwórz wyniki i przechowaj je w oryginalnych listach
for (wynik in wyniki) {
  if (!is.null(wynik)) {
    for (nazwa_modelu in names(wynik)) {
      dac[[wynik$rzad_modelu]] <- wynik$dac
      odchylenie[[wynik$rzad_modelu]] <- wynik$odchylenie
      #var5[[wynik$rzad_modelu]] <- var5
    }
  }
}


# DAC w próbie testowej
dacDf <- as.data.frame(do.call(cbind, dac))
rownames(dacDf) <- "DAC (out of sample)" 

# odchylenie w próbie testowej
odchylenieDf <- as.data.frame(do.call(cbind, odchylenie))
rownames(odchylenieDf) <- "odchylenie (out of sample)" 


porownanie <- bind_rows(dacDf, odchylenieDf)

# Znajdź wszystkie kolumny z maksymalną wartością w trzecim wierszu (DAC w próbie testowej)
max_dac <- max(round(porownanie[1, ], 4), na.rm = TRUE)
kryterium_dac <- porownanie[, which(round(porownanie[1, ], 4) == max_dac), drop = FALSE]

# Znajdź wszystkie kolumny z minimalną wartością w pierwszym wierszu (HQC)
min_odchylenie <- min(round(porownanie[2, ], 4), na.rm = TRUE)
kryterium_odchylenie <- porownanie[, which(round(porownanie[2, ], 4) == min_odchylenie), drop = FALSE]

kryterium_dac
kryterium_odchylenie



Z punktu widzenia DAC najlepszy okazał się ARMA(1, 3)-ARCH_In_Mean (DAC = 0,667 , SO = 0,196). Z punktu widzenia sumy odchyleń najlepszy model to ARMA(7, 8)-ARCH_In_Mean (DAC = 0,633, SO = 0,178).

Te dwa modele sprawdzimy.

Kryterium DAC

Wyciągamy parametry i ustalamy je w formie zmiennych:

# parametry - DAC

wyciagnij_liczby <- function(wyr, x) {
  as.numeric(gsub("[^0-9]", "", substr(wyr, x, x+2), ))
}

rzadAR_DAC <- wyciagnij_liczby(colnames(kryterium_dac), 6)[1]
rzadMA_DAC <- wyciagnij_liczby(colnames(kryterium_dac), 9)[1]
czy_archm_DAC <- grepl("ARCH_In_Mean", colnames(kryterium_dac))[1]


Specyfikacja otrzymanego modelu:

mojGarch_DAC <- ugarchspec(variance.model = list(model = model_Garch, submodel = submodel_Garch,
                                                 garchOrder = c(p, q)), 
                          mean.model = list(armaOrder = c(rzadAR_DAC, rzadMA_DAC), archm = czy_archm_DAC), 
                           distribution.model = rozklad)
 

Następnie stosujemy właśnie prognozę kroczącą - rolującą (robię to drugi raz, pierwszy raz w pętli tylko w celu wyszukania optimum, teraz pełna analiza - chociaż można by to wszystko zapisać w pętli). 

# Prognoza z reestymacją - prognoza krocząca
klaster <- makeCluster(detectCores() - 1)
prognozaGarch_DAC <- ugarchroll(spec = mojGarch_DAC, data = stopaZwrotu, forecast.length = wyrzuc,
                                     n.start = NULL, refit.every = 3, refit.window = "moving",
                                     window.size = 300, solver = "hybrid", calculate.VaR = TRUE, 
                                     VaR.alpha = c(0.01, 0.05),
                                     cluster = klaster)
prognozaGarch_DAC <- resume(prognozaGarch_DAC)
stopCluster(klaster)

prognozaMu_DAC <- prognozaGarch_DAC@forecast$density$Mu
prognozaSigma_DAC <- prognozaGarch_DAC@forecast$density$Sigma


Aby sprawdzić czy rzeczywiście dostałem ten model co chciałem, robię test DAC:
DACTest(forecast=prognozaReest_DAC, 
        actual=tail(as.numeric(stopaZwrotu), wyrzuc),
        test="PT")

Stat
[1] 1.975

$p.value
[1] 0.02411

$H0
[1] "Independently Distributed"

$Decision
[1] "Reject  H0"

$DirAcc
[1] 0.6667


Czyli rzeczywiście to co powinno być. 


Ostatnia część analityczna to utworzenie prognozy na przyszłe okresy, których jeszcze nie ma. Wykorzystamy całą próbę i funkcje ugarchfit() i ugarchforecast():

# Prognoza na podstawie całej próby na 1 przyszły okres
dopasGarchCalaProba_DAC <- tryCatch(
  {
    ugarchfit(spec = mojGarch_DAC, data = stopaZwrotu, solver = "hybrid")
  }, warning = function(w) {
    message("Ostrzeżenie: ", conditionMessage(w))
    ugarchfit(spec = mojGarch_DAC, data = stopaZwrotu, solver = "hybrid")
  }
)

Normalnie w tym miejscu przeprowadziłbym diagnostykę czyli wpisał dopasGarchCalaProba_DAC, ale tego nie robię, bo nie jest to teraz najważniejsze.

Tworzę prognozę na 1 okres do przodu zarówno dla średniej jak i wariancji warunkowej (ryzyka):

prognozaGarchNowa_DAC <- ugarchforecast(dopasGarchCalaProba_DAC, n.ahead = 1)
# średnia
prognozaMuNowa_DAC <- as.numeric(fitted(prognozaGarchNowa_DAC))
# ryzyko
prognozaSigmaNowa_DAC <- as.numeric(sigma(prognozaGarchNowa_DAC))



Dalsza część to już dostosowanie do wykresu:


# Dla wykresu - tworzenie daty, której jeszcze nie ma - 1 okres na przyszłość
datyTest <- tail(index(stopaZwrotu), wyrzuc)
czest <- "week"
nowaData <- seq(from = datyTest[length(datyTest)], by = czest, length.out = 2)[-1]

# Dla wykresu - cały zakres dat
datyWykres <- c(datyTest, nowaData)

# Dla wykresu - średnia prognoza i ryzyko, cały zakres (z reestymacją + na przyszłość)
prognozaMu_xts <- xts(c(prognozaMu_DAC, prognozaMuNowa_DAC), order.by = datyWykres)
prognozaSigma_xts <- xts(c(prognozaSigma_DAC, prognozaSigmaNowa_DAC), order.by = datyWykres)

# wykres
plot(prognozaMu_xts, col = "blue", lwd = 2)
lines(prognozaSigma_xts, col = "red", lwd = 2)
lines(-prognozaSigma_xts, col = "red", lwd = 2)
lines(tail(stopaZwrotu, wyrzuc), col="gray", lwd=2)
addLegend(legend.loc = "topleft", legend.names = c("Średnia prognoza krocząca", 
                                                   "Średnie ryzyko prognozy kroczącej", 
                                                   "Tygodniowa stopa zwrotu WIG"),
          lwd=3, bty="o", cex=0.85, col=c("blue", "red", "gray"))



Sporo błędów, model DAC = 67%, ale nie widać tego specjalnie, szczególnie w ostatnich tygodniach.   

Kryterium sumy odchyleń

Kod jest analogiczny, jedynie zamiast kryterium_DAC wstawiamy kryterium_odchylenie i następnie tak samo z pozostałymi zmiennymi, nie będę więc już powtarzał kodu. 




Z pozoru dziwne się może wydać, że tak dużo odchyleń wychodzi poza obszar oczekiwanego ryzyka. Niestety może być to cena przyjęcia rozkładu normalnego dla przyspieszenia obliczeń.

Uśrednione prognozy

Oba kryteria można uśrednić. Wtedy dostaniemy:



Widać, że to uśrednienie ma sens - tam gdzie kryterium DAC było gorsze, tam lepiej wypadało kryterium odchyleń. Można powiedzieć, że prognoza na następny tydzień jest delikatnie wzrostowa, może wokół zera. 

Natomiast ryzyka nie udało się tu dobrze zamodelować, przypuszczalnie z powodu tego nieszczęsnego rozkładu normalnego. 

S&P 500

Jedyną zmianą w kodzie jest nazwaPliku = "^spx_d.csv". 

Kryterium DAC

Otrzymałem 2 modele:  ARMA(12, 20) (DAC = 0,8 i SO = 0,008) oraz ARMA(9, 14)-ARCH_In_Mean (DAC = 0,8 i SO = 0,0081).

ARMA(12, 20)


ARMA(9, 14)-ARCH_In_Mean


Kryterium sumy odchyleń (SO)

W tym wypadku najlepsze okazały się ARMA(11, 7)-ARCH_In_Mean (DAC = 0,67 i SO = 0,0079) oraz ARMA(7, 9)-ARCH_In_Mean (DAC = 0,67 i SO = 0,00791)

ARMA(11, 7)-ARCH_In_Mean


ARMA(7, 9)-ARCH_In_Mean


Uśrednione prognozy

Wszystkie 4 prognozy możemy uśrednić i dostajemy:


W sumie oczekiwać można spadków amerykańskiego indeksu w następnym tygodniu przy średnim ryzyku. Niemniej trzeba mieć na względzie, że badanie zawiera wiele uproszczeń, jak nieuwzględnienie poprawnego rozkładu reszt, brak korekty o liczbę współczynników modelu oraz wnioskowanie na podstawie zaledwie 30 ostatnich obserwacji, tzn. out-of-sample, wyłączając dopasowanie in-sample.