Jak szybko, ale prawidłowo, porównać zarobki w Polsce z innymi krajami?

Kiedy porównujemy zarobki w różnych krajach, nie możemy jedynie opierać się na danych przekształconych o kurs walutowy, bo ceny w krajach różnią się niezależnie od waluty. Kraje mają inne warunki geograficzne, a stąd jedne surowce są bardziej dostępne, inne mniej. Kraje mają różną kulturę i w jednym może być większy popyt na jakiś produkt, w drugim mniejszy. Ponadto nie wszystkie dobra są handlowalne, a to znaczy, że kurs walutowy nie może ich uwzględnić. Stąd kraje bogatsze będą mieć droższe np. usługi publiczne, a to z kolei przekłada się na wyższe ceny także dóbr handlowalnych.

W poprzednim wpisie pokazałem, jak szybko porównać ceny między krajami OECD. Teraz pokażę, jak porównać dochody między nimi. Najbardziej aktualne dane dotyczą roku 2021, więc pamiętajmy, że są opóźnione w stosunku do porównań cenowych.

1. wchodzimy na stronę OECD z tabelą PKB per capita PPP (w dolarach) . Na tablicy powinna widnieć etykieta:

(jeśli jest inna klikamy w Measure i wybieramy właściwą).

2. W pierwszej kolumnie szukamy kraju, czyli Polski (Poland). Na 2021 r. było to 37 711 USD:

3. Dla porównania szukamy Niemiec (Germany) - jest to 58 784 USD. W USA 70 181 USD, a średnio w OECD 48 958 USD.

4. Czyli średnio realnie zarabiają więcej: 

- Niemcy o 56% więcej (58,8 / 37,7 = 1,56),

-  USA o 86%,

- OECD o 30%.

Porównajmy teraz dochody krajów bez korekty o PPP:

1. Klikamy w Measure:

2. Wybieramy Per head, current prices, current exchange rates

3. Szukamy odpowiednich krajów do porównań na 2021. Polska miała wtedy 17 804 USD. Zauważmy, że to odpowiadało przeciętnemu wynagrodzeniu: średni kurs USD wynosił ok. 3,8. Stąd miesięcznie dostajemy ponad 5,6 tys. zł (17804/12*3,8 = 5634). ZUS podaje, że w 2021 było to 5663 zł.

4. Niemcy mieli 51 204, USA 70 181, OECD 42 581 USD. Zatem nominalnie Niemcy zarabiają prawie 3 razy więcej, USA prawie 4 razy więcej, OECD o 40% więcej.

Jednak w 2022 przedstawione statystyki nieco się zmienią, bo mamy sporo wyższą inflację od Niemiec i USA. Nasz sąsiad prawdopodobnie zarabiał realnie w tym roku o 60% więcej, a USA o 90%. Natomiast w ujęciu nominalnym nieco ich nadgonimy. To też pokazuje, że nominalne porównania nie mają dużego sensu.

Dodatek

Chociaż obliczenia PKB PPP są skomplikowane (zob. dokumentację OECD), to sama idea jest prosta:

PKB [w PLN] = wartość wydatków w gospodarce = średnia cena [w PLN] * ilość dóbr i usług.

Średnia cena w PLN najpierw jest zamieniona w USD:

PKB [USD] = średnia cena [PLN] * ilość dóbr i usług * kurs złotego w USD = średnia cena [USD] * ilość dóbr i usług

Ten PKB to właśnie current prices, current exchange rates.

Żeby przekształcić ten PKB w PKB PPP musimy pomnożyć pierwszy z nich przez odpowiedni współczynnik korygujący:

PKB PPP [USD] = średnia cena [USD] * ilość dóbr i usług * WK,

WK = średnia faktyczna cena [USD] / średnia cena [USD]

Nowa, faktyczna cena w USD jest to cena produktu w USA, a więc w dolarach, a nie jak dotychczas jedynie przekształcona z PLN na USD. 

Przykład:

cena chleba w Polsce = 4 zł

cena chleba w USA = 1,6 USD.

Bochenek chleba starczy na 3 dni, tj. w ciągu roku kupujemy 122 bochenków. Wtedy:

wydatki w Polsce na chleb [PLN] na osobę w ciągu roku = 4 zł*1*122 = 488 zł.

Analogicznie w USA:

wydatki w USA na chleb [USD] na osobę w ciągu roku = 1,6 USD*1*122 = 195 USD.

W 2021 r. kurs złotego wyniósł ok. 0,26 USD. Stąd:

wydatki w Polsce na chleb [USD] na osobę w ciągu roku = 4 PLN*0,26 USD/PLN*122 = 127 USD.

WK = cena chleba w USA w USD / cena chleba w Polsce w USD

WK = 1,6 USD / (4*0,26 USD) = 1,538

127 USD*1,538 = 195 USD.

Dostaliśmy tyle samo wydatków na chleb w Polsce i USA (195 USD rocznie), ponieważ założyliśmy, że ludzie w obu krajach zjadają tyle samo chleba. Różnice będą występować dopiero przy mniej podstawowych produktach. Powiedzmy jednak, że Amerykanie kupują co dwa dni chleb, tak że w ciągu roku kupują 183 bochenków chleba:

wydatki w USA na chleb [USD] na osobę w ciągu roku = 1,6 USD*1*183 = 293 USD.

Dane dla Polski się nie zmieniają, czyli nadal wydajemy 195 dol. na chleb. Porównamy oba kraje:

293 / 195 = 1,5.

Czyli USA wydają o 50% więcej na chleb, ale to również znaczy, że ich dochody przeznaczane na chleb są o 50% większe. Gdybyśmy jednak uwzględnili wszystkie produkty i usługi, dostalibyśmy 1,86, tzn. 86% więcej - to właśnie PKB PPP wykazane przez statystyki OECD.

Jak szybko porównać ceny w Polsce z innymi krajami?

Często chcielibyśmy porównać obecne ceny w Polsce z innym krajem. Wyszukiwarka na pierwszych stronach odniesie nas najprawdopodobniej do artykułów, które nie dostarczają najbardziej aktualnych danych, a poza tym są informacją już przetworzoną, z drugiej ręki. Dlatego najlepszy sposób jest następujący: 

1. wchodzimy na stronę OECD z tabelą porównawczych indeksów cenowych .

2. Mamy tabelę, w której: 

* etykietami wierszy są nazwy krajów, które chcemy porównać z krajem bazowym,

* etykietami kolumn są nazwy krajów wraz z ich walutą - chodzi o kraj bazowy. Zaczynamy od kolumn - szukamy etykiety naszego kraju z walutą, czyli Poland PLN. 

3. Po znalezieniu kolumny szukamy wiersza z etykietą Poland.

4. Zauważmy, że przecięcie wiersza i kolumny daje 100. To jest indeks porównawczy a 100 poziom odniesienia.

5. W kolumnie Poland PLN znajdują się indeksy porównawcze krajów OECD w stosunku do Polski. Np. chcemy sprawdzić o ile różnią się ceny w Niemczech. Szukamy wiersza Germany i sprawdzamy liczbę na przecięciu z kolumną Poland. Obecnie (10.2022) wynosi 177:

 To znaczy, że w październiku 2022 ceny towarów i usług ogólnie były wyższe w Niemczech o 77% niż w Polsce. Jeżeli kupujemy w Polsce coś za 1 zł, to w Niemczech kupimy to samo - średnio biorąc - za 1,77 zł. Oczywiście oznacza to, że najpierw musimy zamienić 1 zł na pewną ilość euro - tutaj byłoby to ok. 0,21 EUR. Czyli gdyby w obu krajach ceny były identyczne, to zapłacilibyśmy 0,21 EUR. W rzeczywistości jednak zapłacimy 0,37 EUR (0,21*1,77 = 0,372), co wynika m.in. z wyższych zarobków i większego bogactwa Niemców.

Możemy to nawet sprawdzić zaczynając porównanie od strony Niemiec. Na przecięciu kolumny Germany EUR i wiersza Poland dostaniemy w tym samym miesiącu 57. To znaczy w Polsce jest o 43% taniej (57/100-1 = -0,43). Dlaczego nie -77%? Bo procentowy wzrost nie odpowiada procentowi spadku. Wiemy, że cena wzrosła z 0,21 EUR do 0,37 EUR (+77%). Powiedzmy teraz, że ją obniżamy z powrotem do 0,21 EUR. To znaczy, że cena spadnie o 0,21/0,37 - 1 = -43%. Czyli dokładnie tyle, ile przewidywaliśmy.  

6. Łatwiej porównamy ceny po ściągnięciu pliku w formacie .CSV lub Excel (wybieramy opcję Export w lewym górnym rogu i wybieramy format).

7. W arkuszu kalkulacyjnym (nie musi być to oczywiście Excel) szukamy pola Country currency i filtrujemy po Poland.

8. Szukamy pola Country i znowu Germany. 

9. Szukamy pola Value. Widzimy, że 177, czyli to samo, co w wersji online.

10. Inny przykład: poszukajmy USA (United States) ceny są wyższe o 130% (Value = 230).

11. Posortujmy pole Value rosnąco. Okazuje się, że w tej liście tylko dwa kraje mają niższe ceny od Polski: Turcja (-36%) i Kolumbia (-24%).

Dodatek

Zainteresowanych dokładną metodologią obliczania PPP odsyłam do dokumentacji OECD. Pokazany w tym wpisie porównawczy indeks poziomu cen (ang. Comparative Price Index)  nazywany jest tam Indeksem poziomu cen (Price Level Index - PLI). 

PLI = PPP / kurs waluty y wyrażony w walucie x,

PPP = cena dóbr w walucie x / cena dóbr w walucie y,

Przykład:

PPP = cena chleba w euro w Niemczech / cena chleba w złotych w Polsce = 1,5 EUR / 4 PLN = 0.375 EUR / PLN.

kurs złotego w euro = 0,21 [EUR/PLN].

PLI = 0,375 [EUR/PLN] / 0,21 [EUR/PLN] = 1,79.

PKN - aktualizacja wyceny

Ponieważ Orlen znowu pokazał większe zyski, tak że EPS wynosi obecnie ok. 40 zł, to wycena się podniesie. Jednak dodatkowo zrobiłem 3 poprawki. Po pierwsze podniosłem współczynnik wypłaty do 15%, tak że w 2023 dywidenda może wynieść 6 zł brutto. Po drugie podniosłem koszt kapitału dla zysków rezydualnych z 9,6% do 10,5%. Wg moich obliczeń klasyczny CAPM faktycznie daje 9,6% dla PKN, ale po uwzględnieniu kapitału ludzkiego (wzrostu wynagrodzeń - zob. Czy trój-betowy CAPM wyjaśnia stopy zwrotu?) oraz ryzyka międzyokresowego (spreadu między rentownością długoterminową a krótkoterminową - zob. Spread między długo- a krótkoterminową stopą procentową) koszt kapitału własnego rośnie o ok. 1 pkt proc. Po trzecie obniżyłem rezydualny wzrost zysków w Scenariuszu 2 z 4 do 3%. Ktoś może powiedzieć, że sobie sztucznie obniżyłem, żeby wycena mi pasowała. Ale de facto chodzi o coś innego, mianowicie cel inflacyjny NBP powinien wynosić 2,5% (podkreślam - powinien wg założeń polityki tej instytucji). Dodatkowe 0,5 pktu proc. wynika z tego, że cel inflacyjny może być rozumiany jako mediana, a średnia może być większa. Odpowiedź na pytanie, dlaczego zakładam max wzrost na poziomie średniej inflacji zawarłem w poprzednim wpisie.   

Jeszcze kwestia dlaczego zwiększyłem współczynnik wypłaty do 15%. Jest to subiektywne odczucie, że spółka będzie naciskana by wypłacać większą dywidendę i 6 zł w kolejnym roku powinno być w zasięgu możliwości. Jednocześnie zapowiadano, że 4 zł dywidendy to minimum. 

W sumie dostałem średnią wycenę 65 zł. W Scenariuszu 1 dla założenia rezydualnego wzrostu g = 1%, wartość wyniosła 60 zł, a dla g = 3% 71 zł. W Scenariuszu 2, tj. z prowizją 0,2%, jest to odpowiednio 58 i 68 zł. Obecnie kurs znajduje się w tych rejonach (64-65 zł), więc racjonalnie. 

Pamiętajmy, że jest to wycena tylko na dziś, a więc uwzględnia bieżącą sytuację. Jutro wycena może być już zupełnie inna. Przykładowo, gdyby współczynnik wypłaty zwiększyć do 25%, to wartość rośnie do ponad 100 zł. Jest to ciągle "spółka polityczna", a to oznacza, że każda wypowiedź polityka będzie pociągać kurs w jakimś kierunku. Niestety może być wykorzystywane przez nich samych, więc ryzyko jest tu podwójne.

Link do obliczeń: Wycena DDF - PKN

Czy PKN Orlen powinien być wycofany z giełdy?

Tak się zastanawiam, po co właściwie PKN jest na giełdzie? Wydaje się to coraz mniej sensowne w kontekście tego czym staje się koncern. Spółki akcyjne, a tym bardziej giełdowe, z definicji mają jak najwięcej zarabiać dla właścicieli. Tymczasem PKN jest spółką strategiczną dla państwa. Istnieje ogromny nacisk polityczny, aby jej zasoby wykorzystać do minimalizacji cen dla odbiorców. Już w sierpniu skrytykowałem PO za populistyczny przekaz, że PKN nie powinien w ogóle zarabiać, że powinien służyć społeczeństwu.  Ale ta narracja jest kontynuowana przez rząd. Kilka dni temu premier M. Morawiecki oświadczył, że "zyski spółek Skarbu Państwa są dziś wysokie i na tym polega zmiana sposobu kształtowania cen, aby zysk netto był niewielki, bo ten zysk pokrywa wszystkie koszty" i dodał że ten sektor powinien się zamknąć w okolicy zera (źródło). W dodatku sam prezes Obajtek w każdym wywiadzie na pierwszym miejscu wymienia zawsze bezpieczeństwo energetyczne, unika tematu zysków, a jedynie na pytania o ceny paliw, opowiada, że koncern dąży do jak najniższych cen dla odbiorców końcowych. No i jeszcze ostatnio z różnych stron słyszymy o wprowadzeniu dodatkowego podatku od "nadwyżkowych zysków" firm energetycznych.

To wszystko razem wpływa oczywiście na ceny tych akcji, ale nie o to chodzi. Samo istnienie takich firm na giełdzie wydaje się w obecnych warunkach kompletnie bez sensu. 

Wcześniej wielu wydawało się, że kryzys energetyczny i wyższe ceny surowców będą sprzyjać właśnie wycenom takich firm, jak PKN. Dziś jednak widać, ta szansa została całkowicie przykryta i to z nawiązką przez ryzyko polityczne.

Właściwie nie ma znaczenia, czy PKN będzie znacjonalizowany czy nie, jeśli jest kontrolowany przez rząd. Ale wycofanie z giełdy, a następnie nacjonalizacja spółki byłoby zwyczajnie bardziej uczciwym posunięciem niż udawanie, że działa się dla dobra akcjonariuszy.

Zróbmy teraz małą analizę. W art. Czy opłaca się kupować PKN Orlen? pokazałem, że wycena PKN mieści się od 65 do 75 zł. Jednak za ta wyceną stało założenie, że PKN będzie rozwijał się średnio w tempie gospodarki, ok. 6%. Wiemy już jednak, że celem zarówno rządu jak i całej opozycji (może za wyjątkiem Konfederacji) jest, aby zyski były zbliżone do zera. Możemy więc przyjąć, że będą rosły średnio w tempie 1-4%, tak aby dać szansę na długoterminowy wzrost o inflację. Jednak dodatkowo założymy, że przez kolejne 4 lata spółka nie będzie w ogóle zarabiać. Po tym okresie będzie rosła w tempie 1 lub 4%.

Zacznę od najprostszej wyceny, z zerową prowizją maklerską. Ostatnio XTB wprowadził zerowe prowizje dla transakcji do pewnego progu, więc również taka wycena może dotyczyć inwestora krótkoterminowego. EPS wg BiznesRadar wynosi prawie 32 zł. Wiemy również, że PKN chce wypłacać co najmniej 3,5 zł na dywidendę. 3,5 zł dotyczyło 2021 r. i aby uzyskać taki sam stosunek w kolejnych latach, dostaniemy 4,3 zł (32*0,135 = 4,3). Przy inflacji 3-4% koszt kapitału własnego (r) wynosi 10-11%. Powiedzmy, że inflacja w 2023 wyniesie 10%. Utrzymując ten sam stosunek r do inflacji dostaniemy r = 27,5% (11/4*10% = 27,5%). W kolejnym roku niech będzie inflacja 5% i wtedy r =13,8%, Potem 4% i r = 11%, a rezydualnie odpowiednio 3,5 i 9,6%. Koszt kapitału 9,6% jest niższy niż poprzednio zastosowany 11%. Jednak sam koncern zdywersyfikował działalność, łącząc się z energetyką, co nieco go przybliży do portfela rynkowego. Poza tym w okresie funkcjonowania PKN na giełdzie 2000-2021 średnia stopa zwrotu z WIG wyniosła tylko 9,2%. W tym samym czasie średnia PKN równa się 11,5%, ale tak jak mówię spadnie z powodu tej dywersyfikacji (również agencja Fitch podniosła ocenę kredytową dla spółki).

Przy tych założeniach dostałem przedział 47-65 zł  (tzn. 47 dla g = 1% i 65 dla g = 4%), średnio 56 zł. Jeśli ktoś uważa, że nawet rząd przesadza i należy przyjąć co najmniej 2% wzrostu, to i tak minimalna cena wyniosłaby 51,5 i średnia 58. Ta ostatnia rzeczywiście niedawno obowiązywała. Ale czy to koniec? Wydaje się bardziej prawdopodobne, że spadnie do 50 zł albo nawet niżej. To jest Scenariusz A.

W Scenariuszu B zakładam, że inwestor płaci koszty transakcyjne. Chociaż wyżej zasugerowałem, że wycena bez tych kosztów będzie prawidłowa, to jednak przy większych wolumenach z pewnością się pojawi. Na XTB - z tego co wiem - wynosi ona 0,2% od aktywów. Przy tej prowizji otrzymałem wycenę 45 zł dla g = 1% i 62 zł dla g = 4%. Średnio 54 zł.

Różnice nie są wielkie między scenariuszami A i B z powodu małego g. Warto zauważyć, że akcje PKN niedawno znajdowały się na tym poziomie:

Obecny poziom 63 zł mieści się jeszcze w przedziale wyceny. Nie zdziwiłbym się, gdyby kurs doszedł do 65 zł zgodnie ze scenariuszem 2A. Mimo wszystko kupowanie powyżej 60 zł na dziś to już czysta spekulacja.

W każdym razie fundamentalni wg mnie powinni całkowicie wyeliminować tego typu spółki z portfela. Mówię o typie, bo nie chodzi tylko o ten przypadek, ale też np. energetyczne. 

Wielokrotna zmiana struktury modelu - funkcja breakpoints (pakiet strucchange dla R)

Dynamika środowiska giełdowo-gospodarczego powoduje, że parametry modeli ekonometrycznych są niestabilne i nie można ich traktować jako coś pewnego. Wiadomo choćby, że występują cykle nieokresowe, tzn. powtarzające się wzrosty i spadki o zmiennej długości i amplitudzie (i stąd niełatwe do prognozowania). Nie ma więc mowy o sezonowości, ale można badać zmienne korelacje, średnie oraz lokalne trendy. W tym art. analizowałem testy na zmianę struktury modelu, ale ich ograniczeniem była pojedyncza zmiana parametru. Jeżeli w próbie występują dwie albo więcej zmian, testy tego typu mogą okazać się niewystarczające, a może nawet wprowadzające w błąd. Poniższe przykłady i cała analiza będzie wykonywana w języku R. Najlepiej używać nakładki na suche środowisko R, RStudio, które znacznie ułatwia wprowadzanie kodów. Rzeczy, które tu omawiam są kontynuacją serii artykułów: 1, 2, 3, 4. Jeżeli coś jest niejasne, to prawdopodobnie w  którymś z nich jest wyjaśnione. Jeżeli nie, należy sięgnąć do pomocy technicznej języka R. 

Zanim przejdę do przykładów nasz skrypt musi mieć wstępny kod. A właściwie przedwstępny, dlatego że w późniejszych przykładach będą kolejne wstępne kody. Dla jasności zapisów kody będę oznaczał jako: Przedwstęp, Wstęp, Pod-wstęp (o ile jest), Część główna. 

Potrzebne są nam cztery pakiety: tseries, strucchange, forecast oraz lmtest:

# Przedwstęp

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

  install.packages("tseries")

  library("tseries")

}

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

  install.packages("strucchange")

  library("strucchange")

}

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

  install.packages("forecast")

  library("forecast")

}

if (require(lmtest)==FALSE) {

  install.packages("lmtest")

  library("lmtest")

}

Jedna uwaga. Raz na jakiś czas pakiety są aktualizowane, a wtedy i tak musimy je na nowo zainstalować, stąd jeśli nie zależy nam na efektywności kodu, to możemy zawsze używać po prostu:

  install.packages("tseries")

  library("tseries")

  install.packages("strucchange")

  library("strucchange")

  install.packages("lmtest")

  library("lmtest")

Przykład 1. Błądzenie losowe

Symulujemy ruch Browna z podziałem na 3 części: 50 danych ze średnią 1, 50 danych ze średnią 1.5 oraz 50 danych ze średnią 1. Połączoną próbę nazwiemy sim123:

# Wstęp

  set.seed(550)

  sim1 = arima.sim(list(order=c(0,0,0)), n=50, mean=1)

  sim2 = arima.sim(list(order=c(0,0,0)), n=50, mean=1.5) 

  sim3 = arima.sim(list(order=c(0,0,0)), n=50, mean=1)

  sim123 = c(sim1, sim2, sim3)

# Część główna

Testujemy stacjonarność (KPSS):

  # test stacjonarności

  if (length(stopa)>200) {

    kpss.test(sim123, lshort=FALSE)

  }  else {

    kpss.test(sim123)

  }

Efekt:

KPSS Test for Level Stationarity

data:  sim123

KPSS Level = 0.41947, Truncation lag parameter = 4, p-value = 0.06877

KPSS wykrył prawidłowo niestacjonarność (p < 0.1). Teraz jednak chcemy znaleźć punkty przełamań. To zastosujmy jak dotychczas QLR. Aby dobrze pokazać zmienną sim123 razem ze statystykami F, stosujemy funkcję par:

par(mfrow=c(2,1), mar=c(2,5,2,2))

plot(sim123, type="l") 

qlr = Fstats(sim123 ~ 1)

plot(qlr, alpha=0.05)

Efekt:

Graficzny test QLR pokazał tylko jeden punkt zmiany, czyli wprowadza w błąd. A wynika on właśnie z tego, że test dzieli próbę na dwie części i poszukuje istotnych zmian. Bai i Perron [1] zaproponowali procedurę uogólniającą QLR na testowanie wielu punktów zmian struktury regresji liniowej. 

Od teraz zamiast QLR będziemy stosować funkcję breakpoints(). De facto już ją stosowałem do znajdowania momentu zmiany, ale wtedy ograniczałem się do jednego punktu. Np. żeby dostać konkretny punkt zmiany w sim123, który odpowiada statystyce supF wpisuję tylko:

breakpoints(qlr)

Efekt:

Breakpoints at observation number:

107 

Corresponding to breakdates:

0.7066667 

Widzimy, że to odpowiada mniej więcej maksimum na wykresie statystyk F.

I uwaga, pokażę, jaki błąd początkowo popełniałem. Funkcja breakpoints() ma taki argument jak breaks, który oznacza maksymalną liczbę wybranych punktów "przełamań". Wydawało się, że wystarczy wpisać:  breakpoints(qlr, breaks = 2) , aby znaleźć 2 takie momenty. Ale przecież argument qlr dotyczy tylko jednego takiego punktu, dlatego to nie zmieni efektu. Co więc zrobić? Okazuje się, że te funkcje są dość elastyczne i można je zastosować na różnych funkcjach. I tak zamiast argumentu qlr, czyli zamiast statystyk F, można wpisać sam model sim123 ~ 1. Mówiąc inaczej breakpoints częściowo zastępuje funkcję Fstats. Dlaczego częściowo? Bo funkcja plot(breakpoints()) działa wtedy inaczej. 

Po kolei. Pierwsza część nowego kodu będzie wyglądać następująco. Dla ułatwienia przekształcamy szybko dane w szereg czasowy:

sim123 = ts(sim123)

Następnie stosujemy na tej zmiennej breakpoints: 

bai = breakpoints(sim123 ~ 1, breaks = 2)

summary(bai)

Efekt:

Breakpoints at observation number:

m = 1      107

m = 2   49 102

Corresponding to breakdates:                               

m = 1                     0.713333333333333

m = 2   0.326666666666667 0.68             

Fit:  

m   0     1     2    

RSS 148.9 136.5 124.3

BIC 434.6 431.6 427.6

Ponieważ sami wybraliśmy dwa punkty, to dostaliśmy analizę dla dwóch punktów. Najpierw wskazano jakie to punkty. Warto zwrócić uwagę, że algorytm inaczej oblicza je zależnie od ustawienia breaks, co wyraża tu indeks m. Dla m = 1 algorytm oszacował 107 (na datach - ponieważ nie ma ustawień dat - w częściach ułamkowych - 0,71), a dla m = 2 49 (0,33) i 102 (0,68). 

Najważniejsza statystyka jest na dole. Patrzymy jak zachowuje się BIC. Dla m = 2 BIC jest najmniejszy w porównaniu z m = 0 i m = 1. Zatem wnioskujemy, że algorytm prawidłowo wskazał istotne 2 przełamania. 

Co by się stało, gdybyśmy zwiększyli m do 4? Popatrzmy:

bai = breakpoints(sim123 ~ 1, breaks = 4)

summary(bai)

Efekt:

Breakpoints at observation number:

m = 1            107

m = 2      49    102

m = 3      49 75 102

m = 4   24 49 75 102

Corresponding to breakdates:  

m = 1            107

m = 2      49    102

m = 3      49 75 102

m = 4   24 49 75 102

Fit:                   

m   0     1     2     3     4    

RSS 148.9 136.5 124.3 122.7 122.5

BIC 434.6 431.6 427.6 435.7 445.4

 

Widać, że dla m >= 3 BIC wzrósł, a więc wynik się pogorszył. Wniosek: nie ma więcej niż 2 punkty zmiany modelu - prawidłowo.

Jeżeli chcemy zobaczyć to graficznie, to używamy plot:

bai = breakpoints(sim123 ~ 1, breaks = 4)

plot(bai)

Efekt:

Plot możemy jeszcze użyć do dwóch celów. 

W połączeniu z funkcją lines możemy na wykresie czasowym pokazać momenty zmian w czasie: 

par(mfrow=c(2,1), mar=c(2,5,2,2))

sim123 = ts(sim123)

plot(sim123) 

bai = breakpoints(sim123 ~ 1)

lines(bai, breaks = 2)

Efekt:

Grafika wyżej zawiera pewien błąd związany z ucięciem 15% skrajnych danych na Fstat. Wcześniej to poprawiałem dla prawdziwych danych, ale w tym wypadku to byłaby bardziej sztuka dla sztuki. 

Jeszcze ciekawszą grafikę można pokazać wykorzystując funkcję breakfactor:

sim123.model = lm(sim123 ~ breakfactor(bai)-1)

dopas = fitted(sim123.model)

plot(sim123)

lines(dopas, col=4, lwd=3)

Efekt:

Parametr lwd kalibruje grubość linii, col - kolor. 

Zwróćmy uwagę, że zapisując funkcję breakfactor usnąłem wyraz wolny poprzez odjęcie 1 na końcu. W języku R jest pewnego rodzaju umowa, że jeśli nie usuniemy sami wyrazu wolnego, to on zostanie sam dodany. Moim zdaniem to jest zbędne i trochę mylące, bo nie trudno przecież dodać 1, jeśli się chce mieć wyraz wolny, a jeśli mamy tylko wyraz wolny bez zmiennej, to tym bardziej . W każdym razie odjęcie 1 pozwoli nam tutaj wyłuskać prawidłowe parametry w postaci trzech segmentów: bez zmiany parametru, po zmianie i po drugiej zmianie. Aby wyłuskać zmienne współczynniki, używamy funkcji coef():

coef(sim123.model)

Efekt:

breakfactor(bai, breaks = 2)segment1 breakfactor(bai, breaks = 2)segment2 breakfactor(bai, breaks = 2)segment3 

                           0.9340913                            1.6430498                            0.7015123 

W sumie cały nasz nowy kod może wyglądać następująco:

#Wstęp

#symulacja procesu losowego  - tu ruchu Browna - zmienne parametry

set.seed(199)

sim1 = arima.sim(list(order=c(0,0,0)), n=50, mean=1)

sim2 = arima.sim(list(order=c(0,0,0)), n=50, mean=1.5) 

sim3 = arima.sim(list(order=c(0,0,0)), n=50, mean=1)

sim123 = c(sim1, sim2, sim3)

# Część główna

# test stacjonarności

if (length(sim123)>200) {

  kpss.test(sim123, lshort=FALSE)

}  else {

  kpss.test(sim123)

}

#ustawienie parametru graf

par(mfrow=c(3,1), mar=c(5,5,2,2))

sim123 = ts(sim123)

plot(sim123) 

#Procedura Bai-Perrona:

bai = breakpoints(sim123 ~ 1, breaks = 4)

lines(bai, breaks = 2)

summary(bai)

plot(bai)

sim123.model = lm(sim123 ~ breakfactor(bai)-1)

dopas = fitted(sim123.model)

plot(sim123)

lines(dopas, col=4, lwd=3)

coef(sim123.model)

#przywracanie parametru graf

par(mfrow=c(1,1), mar=c(5, 4, 4, 2) + 0.1)

Przykład 2.

Wykonajmy ten sam kod, ale zmieńmy jedną linię. Zamiast

sim123.model = lm(sim123 ~ breakfactor(bai)-1)

wpiszmy

sim123.model = lm(sim123 ~ breakfactor(bai, breaks = 1)-1)

Efekt:

breakfactor(bai, breaks = 1)segment1 breakfactor(bai, breaks = 1)segment2 

                            1.292695                             0.655961 

Czyli widzimy, że jeśli nie wstawimy parametru breaks do breakfactor, to algorytm automatycznie wybiera optymalną liczbę punktów. Dopisanie tego parametru wymusza taką liczbę zmian, jaką chcemy. 

Przykład 3. ARMA(1, 1)

Znów zmieniamy parametry procesu w połowie próby dwukrotnie, ale tym razem średnią zachowuję taką samą, a zmienia się część AR i MA:

#Wstęp

set.seed(1980)

sim1 = arima.sim(list(order=c(1,0,1), ar=0.6, ma=0.3), n=50, mean=1)

sim2 = arima.sim(list(order=c(1,0,1), ar=0, ma=0), n=50, mean=1) 

sim3 = arima.sim(list(order=c(1,0,1), ar=0.6, ma=0.3), n=50, mean=1)

sim123 = c(sim1, sim2, sim3)

# Część główna

# test stacjonarności

if (length(sim123)>200) {

  kpss.test(sim123, lshort=FALSE)

}  else {

  kpss.test(sim123)

}

#ustawienie parametru graf

par(mfrow=c(3,1), mar=c(5,5,2,2))

sim123 = ts(sim123)

#Procedura Bai-Perrona:

bai = breakpoints(sim123 ~ 1)

plot(bai)

Efekt:

KPSS Test for Level Stationarity

data:  sim123

KPSS Level = 0.31029, Truncation lag parameter = 4, p-value = 0.1

Tym razem dostaliśmy niespójność: dane są stacjonarne mimo wskazania punktów zmian. Ponadto test wskazał 3 punkty, choć powinien dwa.

Błąd, o którym już pisałem w art. o trzech technikach testu F wynika z tego, że sprawdziliśmy sobie samą średnią (wyraz wolny), a nie cały model, który tutaj dodatkowo (oprócz wyrazu wolnego) zawiera zmienne zależne. Poprawiamy:

# Wstęp

set.seed(1980)

sim1 = arima.sim(list(order=c(1,0,1), ar=0.6, ma=0.3), n=50, mean=1)

sim2 = arima.sim(list(order=c(1,0,1), ar=0, ma=0), n=50, mean=1) 

sim3 = arima.sim(list(order=c(1,0,1), ar=0.6, ma=0.3), n=50, mean=1)

sim123 = c(sim1, sim2, sim3)

# Część główna

# test stacjonarności

if (length(sim123)>200) {

  kpss.test(sim123, lshort=FALSE)

}  else {

  kpss.test(sim123)

}

sim123 = ts(sim123)

#tworzenie modelu arma(1,1)

sim = ts(sim123[-1])

sim_ar1 = ts(sim123[-length(sim123)])

sim_arma11 = arima(sim123, order=c(1,0,1))

sim_ma = ts(residuals(sim_arma11))

sim_ma1 = ts(sim_ma[-length(sim_ma)])

#Procedura Bai-Perrona:

bai = breakpoints(sim ~ sim_ar1 + sim_ma1)

plot(bai)

Efekt:

Tym razem jest OK.

To dokończmy kod:

plot(sim123) 

lines(bai, breaks = 2)

summary(bai)

sim.model = lm(sim ~ breakfactor(bai)-1)

dopas = fitted(sim.model)

lines(dopas, col=4, lwd=3)

coef(sim123.model)

#przywracanie parametru graf

par(mfrow=c(1,1), mar=c(5, 4, 4, 2) + 0.1)

Efekt:

Kolejne przykłady będą robione na danych rzeczywistych. Dane w formie plików najpierw ściągam z internetu. I tu pokażę jedną nową rzecz. Żeby otworzyć plik w R, musimy podać ścieżkę i tu napotykamy niekompatybilność R z Windowsem: w tym pierwszym ścieżki używają slasha, a w tym drugim backslasha. Coś, co np. w VBA jest łatwe do podmiany, w R okazuje się zadziwiająco trudne. W każdym razie jeśli mamy w miarę nową wersję R, to stosujemy następujący kod. Powiedzmy, że nasza ścieżka z plikami (folder roboczy) to C:\R\testy

# Wstęp

#zamieniam "\" na "/"

sciezka = r"(C:\R\testy)"

sciezka = gsub("\\\\", "/", sciezka)

Gdyby zdefiniować sciezka = "C:\R\testy" dostalibyśmy błąd, dlatego takie przekształcenie jest potrzebne. Funkcja gsub() jest ogólną funkcją do zamiany znaków tekstowych. Dlaczego trzeba aż 4 razy użyć backslasha? Ogólnie backslash ma znaczenie specjalne - wyłącza znaczenie znaków specjalnych, także takich jak on sam. A więc żeby program go rozumiał jako zwykły znak tekstowy, trzeba jego samego użyć - dopisać po pierwszym "\", czyli dostać "\\". Ale to znaczenie występuje na dwóch poziomach. Pierwszy poziom to sam program - interpreter R. Drugi poziom to funkcja typu regex (tutaj gsub). Oba są jakby niezależne. Stąd trzeba aż 4 razy użyć backslasha - dla samego programu oraz dla regexa. Bardziej detaliczny opis poniżej.

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

 O ile np. w VBA interpreter najpierw przetwarza to co wpiszemy na dosłowne znaki, o tyle R od razu interpretuje znaki specjalne i wtedy nie wie co w konsoli wpisać (nie da się wpisać "semantycznego" znaczenia znaku specjalnego). Co by jednak się stało, gdyby w kodzie czy konsoli wpisać sam znak "\"? Błędu nie będzie, bo interpreter widzi, że po backslashu nic już nie ma, więc nie ma też specjalnego znaku. Co ciekawe, błędu nie będzie też, gdy wpiszemy po "\" cyfrę powyżej 0. Np. "\1", "\2" itp. nie wywoła błędu, ale "\0" już tak.

W sumie, aby zamienić "\" na "/" trzeba wyłączyć znak specjalny "\" poprzez zapisanie "\\". Zapiszmy taki kod:

sciezka = r"(C:\R\testy)"

sciezka

W konsoli dostaniemy napis "C:\\R\\testy". Widzimy, że program sam doda jeden backslash. Teraz funkcja gsub() ma znaleźć te znaki i zamienić na slash. Kłopot polega na tym, że mamy do czynienia z podwójną interpretacją znaku "\": najpierw w samym programie (pierwszy poziom), a następnie w funkcji gsub (drugi poziom). Pierwszy argument gsub() to poszukiwane znaki, w tym przypadku szukamy "\". Ale wiemy już, że sam interpreter R zakazuje wpisywać do kodu pojedynczy "\" i trzeba go zastąpić "\\". A więc funkcja regex  (czyli gsub) widzi tylko jeden "\". I tu zaczyna się drugi poziom problemu, bo regex nie wie jak interpretować samotny backslash: jako dosłowny tekst czy jako wyłącznik znaku specjalnego. Aby uzyskać pewność trzeba dodać kolejny "\" przed "\". W ten sposób mielibyśmy trzy backslashe. Jednakże trzy backslashe są znowu niedopuszczalne dla samego interpretera R (czyli dla poziomu pierwszego), ponieważ widzi on nowy (drugi) backslash, który jak pamiętamy był znakiem niejednoznacznym i potrzebował dodatkowego "\".  Dopiero dodanie czwartego backslasha stabilizuje sytuację. Skomplikowana sprawa, ale myślę, że w miarę logicznie wyjaśniona.

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

Tak naprawdę problem ten możemy ominąć wstawiając do gsub dodatkowo opcję fixed = T:

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

Podsumujmy część wstępną:

# Wstęp

#zamieniam "\" na "/"

sciezka = r"(C:\R\testy)"

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

# albo

# sciezka = gsub("\\\\", "/", sciezka)

# ustawiamy folder roboczy

setwd(sciezka)

I teraz idziemy z żywymi przykładami.

Przykład 4. Ceny węgla (ICE Rotterdam Coal)

Sprawdzamy notowania węgla od 31.08.2006 do 31.10.2022 - dane miesięczne, źródło stooq.pl: kontrakty na węgiel (oznaczenie LU.F). Mam plik o nazwie coal.csv. Początek pliku wygląda tak:

 Data Otwarcie Najwyższy Najniższy Zamknięcie

31.08.2006 70 70.5 67.9 68.45

30.09.2006 67.5 67.5 63.75 66.75

31.10.2006 65 68.55 63.5 67.75

30.11.2006 68.25 68.75 66.2 68.75

31.12.2006 67 69.25 67 68

I kod:

# Pod-wstęp

nazwa = "Coal.csv"

  

  plik = read.csv(nazwa, sep=";")

  

  if (ncol(plik)==1) {

    

    plik = read.csv(nazwa, sep=",")

    

  }

  

Żeby zawsze mieć prawidłowo odczytaną datę, używamy opcji TryFormats funkcji as.Date():

  daty = as.Date(plik[,1], tryFormats = c("%d.%m.%Y", "%Y.%m.%d", "%d-%m-%Y", "%Y-%m-%d", "%d/%m/%Y", "%Y/%m/%d"))

  

Wadą danych stooq są ostatnie dane, których daty są błędnie oznaczane. Aby nie usuwać ręcznie zbędnych danych, stosujemy prosty warunek

 

  #jeżeli ostatnia data przekracza wczorajszą datę, to usuwamy ostatnią obserwację

  

  if (daty[length(daty)]>=Sys.Date()) {

    

    plik = plik[-nrow(plik),]

    

    daty = daty[-length(daty)]

  }

  

Dopasowujemy format do szeregu czasowego i tworzymy stopy zwrotu (zmienna o nazwie stopa)

  

  

  rok = as.numeric(format(daty, "%Y"))

  mc = as.numeric(format(daty, "%m"))

  dz = as.numeric(format(daty, "%d"))

  

  cena = as.numeric(plik[,5])

  cena = ts(cena, start=c(rok[1], mc[1]), frequency=12)

  

  stopa = diff(cena)/cena[-length(cena)]

Dane gotowe. Kurs prezentuje poniższa grafika (kod: plot(cena) , plot(stopa))

Badamy stacjonarność zmiennej:

# Część główna

  # test stacjonarności

  if (length(stopa)>200) {

    kpss.test(stopa, lshort=FALSE)

  }  else {

    kpss.test(stopa)

  }

  Efekt:

KPSS Test for Level Stationarity

data:  stopa

KPSS Level = 0.22501, Truncation lag parameter = 4, p-value = 0.1

Zmiany procentowe są stacjonarne. Oznacza to, że raczej nie mamy tutaj do czynienia ze zmienną średnią. Sprawdzamy autokorelację cząstkową. Możemy ją sprawdzić za pomocą funkcji pacf(). Niestety ona sama nie zawiera informacji o istotności stat. Dlatego do pacf() dodamy kod na skumulowany rozkład odwrotny. W sumie:

  pacf(stopa, plot=FALSE, lag.max=10)

  ist_pacf = qnorm((1 + 0.95)/2)/sqrt(sum(length(stopa)))

  noquote(c("Istotna autokorelacja będzie od poziomu: ", ist_pacf))

Efekt:

Partial autocorrelations of series ‘stopa’, by lag

0.0833 0.1667 0.2500 0.3333 0.4167 0.5000 0.5833 0.6667 0.7500 0.8333 

 0.242  0.089 -0.081  0.086 -0.005 -0.010  0.100 -0.004 -0.047  0.055 

Istotna autokorelacja będzie od poziomu:  0.140717213329746 

Ułamkowe rzędy mogą być mylące, bo są to po prostu kolejne części roku: 1/12 = 0.0833 = 1 miesiąc, 1/6 = 0.1667 = 2 miesiące itd.  Dla kwartałów byłoby to 0.25 = kwartał  ; 0.5 = 2 kwartał itd. Istotna korelacja występuje od poziomu 0.14, a to znaczy, że mamy do czynienia z korelacją 1 rzędu. Prawdopodobnie zastosujemy AR(1). Upewniamy się za pomocą funkcji auto.arima:

  arma_stopa = auto.arima(stopa,ic="bic", max.P=7, max.Q=7)

  arma_stopa

Efekt:

ARIMA(1,0,0) with zero mean 

Coefficients:

         ar1

      0.2612

s.e.  0.0705

sigma^2 = 0.01162:  log likelihood = 157.29

AIC=-310.59   AICc=-310.53   BIC=-304.05

Istotnie występuje AR(1), więc tworzymy ten model. Zwróćmy uwagę, że stopa zwrotu już zaczyna się od drugiej obserwacji, więc zmienna opóźniona 1 rzędu będzie zaczynać się od 3-ciej. Choć moglibyśmy użyć podobnej techniki jak przy symulacji, to użyjemy innej, w sumie prostszej techniki, używając funkcji lag. Mimo prostoty funkcja nie jest zbyt intuicyjna, więc opiszę ją nieco dokładniej. 

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

Zacznijmy od pytania czym jest zmienna opóźniona o 1 okres? Powiedzmy, że zmienna to nasz szereg czasowy. Wtedy opóźniona o 1 będzie to:

zmienna1 = zmienna[-length(zmienna)]

i aby porównać ją z bieżącym okresem, musieliśmy utworzyć nową bieżącą zmienną , którą definiowaliśmy jako 

zmienna0 = zmienna[-1]

Zatem zmienna1 to zmienna kończąca się na przedostatniej obserwacji, a zmienna0 to zmienna zaczynająca się od drugiej obserwacji. Jeżeli więc tworzymy model typu:

zmienna0 ~ zmienna1

to patrzymy jak pierwsza obserwacja (zmiena1) wpływa na drugą (zmiena0), druga na trzecią itd. 

I teraz najważniejsze: opisana technika wcale nie tworzy zmiennej opóźnionej sensu stricte. Czas płynie do przodu, a to znaczy, że opóźniona obserwacja to taka, która występuje później niż bieżąca. A więc w rzeczywistości tworzyliśmy tutaj nie zmienną opóźnioną, ale przyspieszoną (pospieszoną)! Zmienna0 to zmienna obserwowana w czasie teraźniejszym i patrzymy jak na nią wpływa zmienna1 z przeszłości. Skoro więc porównujemy dzisiejsze obserwacje z tym co było, to w pewnym sensie traktujemy obserwacje z przeszłości jako teraźniejsze, bo skoro mają teraz wpływ, to znaczy, że ciągle istnieją tu i teraz. Z technicznego punktu widzenia bierzemy dane z przeszłości i przesuwamy do przodu. To przesunięcie do przodu oznacza właśnie przyspieszenie, a nie opóźnienie. 

Dlaczego o tym mówię? Bo funkcja lag tworzy zmienną opóźnioną sensu stricte. Wielu tego nie rozumie i w poradnikach możemy przeczytać, że lag działa na odwrót i trzeba używać ujemnego kroku, aby dostać opóźnienie. W rzeczywistości ten ujemny lag da nam właśnie pospieszenie okresu, ale to właśnie o to nam chodzi. Dlatego do intuicyjnie rozumianych "opóźnień" używamy ujemnego kroku, a do "pospieszeń" dodatniego. W każdym razie nie ma tu żadnego błędu, a jedynie musimy pamiętać, że jeżeli wyciągamy zmienną z przeszłości do teraźniejszości, to de facto ją pospieszamy, choć w modelu nazywamy ją zmienną opóźnioną.

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

Tworzymy zmienną opóźnioną:

stopa1 = lag(stopa, k = -1)

To niestety nie wystarcza do stworzenia AR(1), bo przesuwa całą zmienną razem z okresem, podczas gdy my chcemy zachować stałość okresu, a jedynie przesuwać obserwacje w czasie. Żeby uzyskać żądany efekt, musimy użyć funkcji intersect (albo ts.intersect aby działać od razu na ts), która działa jak funkcja INNER JOIN w SQL. Nazwa intersect jest myląca, powinna działać jak przecięcie zakresów w VBA albo mieć nazwę np. "innerjoin". Ale nieważne, tworzymy nową zmienną:

dane = ts.intersect(stopa, stopa1)

którą wykorzystujemy do przekształcenia naszych zmiennych:

stopa0 = dane[, "stopa"]

stopa1 = dane[, "stopa1"]

I dalej już to co wcześniej:

  #test Bai-Perrona

  bai = breakpoints(stopa0 ~ stopa1, breaks = 3)

  summary(bai)

  plot(bai)

Efekt:

Breakpoints at observation number:

m = 1         164

m = 2   29    164

m = 3   29 58 164

Corresponding to breakdates:           

m = 1                   2020(5)

m = 2   2009(2)         2020(5)

m = 3   2009(2) 2011(7) 2020(5)

Fit:             

m   0        1        2        3       

RSS    2.228    2.149    2.082    2.041

BIC -297.590 -288.798 -279.130 -267.116

A więc brak zmian struktury. Z jednej strony to pozytywne, bo znaczy, że można lepiej prognozować kolejne zmiany dzięki autoregresji, która pozostaje stabilna. Z drugiej strony ostatnie lata przynoszą o wiele większe wahania, więc możliwe, że to wariancja stopy zwrotu się zmieniła, a nie średnia czy autokorelacja.

Zmianę wariancji, jak pamiętamy, wykonujemy np. za pomocą takich testów jak Goldfelda-Quandta, Breuscha-Pagana czy Harrisona-McCabe. Kod:

  czas = c(1:length(stopa))

  bptest(stopa ~ czas)

  gqtest(stopa ~ 1)

  hmctest(stopa ~ 1)

Normalnie w każdym z tych testów użyłbym wyrazu wolnego za zmienną niezależną, ale dla pierwszego nie da się - pokazuje się błąd, dlatego ustawiłem czas. Wszystkie z tych testów p-value < 0,01, zatem nie ma wątpliwości, że wariancja stopy zwrotu jest zmienna w czasie. Skoro test Bai-Perrona był niewrażliwy na jej zmianę, to znaczy, że można go spokojnie używać na danych heteroskedastycznych.

Cały kod:

# Pod-wstęp

nazwa = "Coal.csv"

plik = read.csv(nazwa, sep=";")

if (ncol(plik)==1) {

  plik = read.csv(nazwa, sep=",")

}

daty = as.Date(plik[,1], tryFormats = c("%d.%m.%Y", "%Y.%m.%d", "%d-%m-%Y", "%Y-%m-%d", "%d/%m/%Y", "%Y/%m/%d"))

#jeżeli ostatnia data przekracza wczorajszą datę, to usuwamy ostatnią obserwację

if (daty[length(daty)]>=Sys.Date()) {

  plik = plik[-nrow(plik),]

  daty = daty[-length(daty)]

}

rok = as.numeric(format(daty, "%Y"))

mc = as.numeric(format(daty, "%m"))

dz = as.numeric(format(daty, "%d"))

cena = as.numeric(plik[,5])

cena = ts(cena, start=c(rok[1], mc[1]), frequency=12)

stopa = diff(cena)/cena[-length(cena)]

# Część główna

#test stacjonarności

if (length(stopa)>200) {

  kpss.test(stopa, lshort=FALSE)

}  else {

  kpss.test(stopa)

}

pacf(stopa, plot=FALSE, lag.max=10)

ist_pacf = qnorm((1 + 0.95)/2)/sqrt(sum(length(stopa)))

noquote(c("Istotna autokorelacja będzie od poziomu: ", ist_pacf))

arma_stopa = auto.arima(stopa,ic="bic", max.P=7, max.Q=7)

arma_stopa

stopa1 = lag(stopa, k = -1)

dane = ts.intersect(stopa, stopa1)

stopa0 = dane[, "stopa"]

stopa1 = dane[, "stopa1"]

#test Bai-Perrona

bai = breakpoints(stopa0 ~ stopa1, breaks = 3)

summary(bai)

plot(bai)

czas = c(1:length(stopa))

bptest(stopa ~ czas)

gqtest(stopa ~ 1)

hmctest(stopa ~ 1)

Przykład 5. Inflacja USA

Inflacja roczna USA w latach 1775-2021 (247 obs), średnioroczna. Źródło: https://www.measuringworth.com/

Zanim przejdę do części głównej, zmodyfikuję nieco część pod-wstępną, tak aby uwzględniała zarówno roczne, jak i miesięczne, a także aby można było łatwo przełączać się między danymi pierwotnymi a przekształconymi. W tym ostatnim chodzi o to, czy nasza zmienna już jest gotowa (stopa inflacji jest podana), czy może musimy ją uzyskać (przekształcić ceny w stopy zwrotu). Oprócz tego do read.csv dodam opcję dec , aby prawidłowo czytać ułamki dziesiętne (w moim pliku używam przecinków).  Dzięki temu nie muszę tworzyć dwóch różnych skryptów:

# Pod-wstęp

nazwa = "inflation US.csv"

plik = read.csv(nazwa, sep=";", dec=",")

if (ncol(plik)==1) {

  plik = read.csv(nazwa, sep=",", dec=".")

}

if (is.numeric(plik[1,1])) {

  daty = paste("01.01.", plik[,1], sep="")

} else { 

  daty = plik[,1] 

}

  daty = as.Date(daty, tryFormats = c("%d.%m.%Y", "%Y.%m.%d", "%d-%m-%Y", "%Y-%m-%d", "%d/%m/%Y", "%Y/%m/%d"))

#jeżeli ostatnia data przekracza wczorajszą datę, to usuwamy ostatnią obserwację

  if (daty[length(daty)]>=Sys.Date()) {

    plik = plik[-nrow(plik),]

    daty = daty[-length(daty)]

  }

fr = 1

kol = 2

st = FALSE

rok = as.numeric(format(daty, "%Y"))

mc = as.numeric(format(daty, "%m"))

dz = as.numeric(format(daty, "%d"))

data_startowa = c(rok[1], mc[1])

cena = as.numeric(plik[,kol])

cena = ts(cena, start = data_startowa, frequency = fr)

if (st == TRUE) {

  stopa = diff(cena)/cena[-length(cena)]

} else {

  stopa = cena

}

Jak widać, mimo, że okres pierwotnie podany jest w latach, to przekształcam je w pełne daty od pierwszego stycznia i z tego już wyciągam lata. Nowością są też 3 nowe zmienne: fr, kol oraz st, gdzie odpowiednio oznaczają częstość (tu roczna, czyli 1), nr kolumny szeregu czasowego (tu 2, wcześniej na danych stooq.pl było to 5) i zmienna logiczna typu TRUE / FALSE , gdzie TRUE oznacza, że zamieniamy pierwotny szereg w stopy oraz FALSE przeciwnie (tu FALSE).

W części głównej kodu jest bardzo podobnie jak w poprzednim przykładzie:

# Część główna

# test stacjonarności

if (length(stopa)>200) {

  kpss.test(stopa, lshort=FALSE)

}  else {

  kpss.test(stopa)

}

# sprawdzanie autokorelacji

pacf(stopa, plot=FALSE, lag.max=10)

ist_pacf = qnorm((1 + 0.95)/2)/sqrt(sum(length(stopa)))

noquote(c("Istotna autokorelacja będzie od poziomu: ", ist_pacf))

# optymalne arima

arma_stopa = auto.arima(stopa,ic="bic", max.P=7, max.Q=7)

arma_stopa

Efekt:

KPSS Test for Level Stationarity

data:  stopa

KPSS Level = 0.43043, Truncation lag parameter = 15, p-value = 0.06404

Partial autocorrelations of series ‘stopa’, by lag

     1      2      3      4      5      6      7      8      9     10 

 0.390  0.158 -0.153  0.048  0.063 -0.063 -0.040  0.116  0.009 -0.049 

[1] Istotna autokorelacja będzie od poziomu:  0.124709521934413  

                      

ARIMA(1,1,1) 

Coefficients:

         ar1      ma1

      0.3758  -0.9817

s.e.  0.0614   0.0143

sigma^2 = 34.33:  log likelihood = -784.25

AIC=1574.5   AICc=1574.6   BIC=1585.02

Inflacja USA okazuje się procesem niestacjonarnym, co podpowiada nam już, że dostaniemy model niestabilny. Najlepszym ARIMA okazał się jednak model bez wyrazu wolnego, ale jest to proces zintegrowany stopnia 1 z AR(1) i MA(1) - w sumie ARIMA(1,1,1). KPSS i ARIMA świetnie się tutaj uzupełniają.

Nie tworzymy w tym momencie modelu ARIMA, tylko zwyczajny ARMA(1,1), co pozwoli nam zobaczyć niejako, jak ta średnia się zmienia w czasie. Ale żeby to było możliwe, musimy dodać wyraz wolny. Następnie testujemy stabilność tego modelu:

# tworzenie optymalnego arma

stopa_ar1 = lag(stopa, k = -1)

stopa_ma = arma_stopa$residuals

stopa_ma1 = lag(stopa_ma, k = -1)

dane = ts.intersect(stopa, stopa_ar1, stopa_ma1)

stopa0 = dane[, "stopa"]

stopa_ar1 = dane[, "stopa_ar1"]

stopa_ma1 = dane[, "stopa_ma1"]

#test Bai-Perrona

bai = breakpoints(stopa0 ~ stopa_ar1 + stopa_ma1 + 1, breaks = 3)

summary(bai)

if (st ==TRUE) {

  par(mfrow=c(3,1), mar=c(2,5,2,2))

  plot(cena)

} else {

  par(mfrow=c(2,1), mar=c(2,5,2,2))

}

plot(bai)

nowy_model = lm(stopa0 ~ breakfactor(bai)-1)

rok0 = as.numeric(format(as.yearmon(time(stopa0)[1]), "%Y"))

mc0 = as.numeric(format(as.yearmon(time(stopa0)[1]), "%m"))

dopas = ts(fitted(nowy_model), start=c(rok0, mc0), frequency = fr)

coef(nowy_model)

plot(stopa0, ylab=colnames(plik)[kol])

lines(dopas, col=4, lwd=3)

par(mfrow=c(1,1), mar=c(5, 4, 4, 2) + 0.1)

Efekt:

Breakpoints at observation number:

m = 1      85    

m = 2   39 85    

m = 3   39 85 145

Corresponding to breakdates:          

m = 1        1860     

m = 2   1814 1860     

m = 3   1814 1860 1920

Fit:            

m   0    1    2    3   

RSS 8250 7141 6729 6577

BIC 1584 1571 1578 1595

breakfactor(bai)segment1 breakfactor(bai)segment2 

               0.3842353                2.3413665 

Zwracam uwagę, że w linii 

bai = breakpoints(stopa0 ~ stopa_ar1 + stopa_ma1 + 1, breaks = 3)

jest dodany wyraz wolny, a w linii

nowy_model = lm(stopa0 ~ breakfactor(bai)-1)

odjęty. W tym ostatnim tworzymy nowy model, więc tworzony byłby nowy wyraz wolny - a tego już nie chcemy.

Istotna zmiana modelu nastąpiła w 1860 r. Problem z jego interpretacją polega na tym, że nie widzimy w ogóle, czy niestabilność dotyczy wyrazu wolnego, części AR, MA, a może wszystkich na raz. To jest ograniczenie tego testu. Są ścisłe sposoby na jego rozwiązanie, ale na razie oprę się jedynie na KPSS. To znaczy wiem, że KPSS wskazał na niestacjonarność, czyli niestabilność właśnie wyrazu wolnego. Dlatego przetestujmy sam wyraz wolny. Kod będzie prawie ten sam, jedynie zmieni się od testu Bai-Perrona:

#test Bai-Perrona

stopa0 = stopa

bai = breakpoints(stopa0 ~ 1, breaks = 3)

summary(bai)

if (st ==TRUE) {

  par(mfrow=c(3,1), mar=c(2,5,2,2))

  plot(cena)

} else {

  par(mfrow=c(2,1), mar=c(2,5,2,2))

}

plot(bai)

nowy_model = lm(stopa0 ~ breakfactor(bai) - 1)

rok0 = as.numeric(format(as.yearmon(time(stopa0)[1]), "%Y"))

mc0 = as.numeric(format(as.yearmon(time(stopa0)[1]), "%m"))

dopas = ts(fitted(nowy_model), start=c(rok0, mc0), frequency = fr)

coef(nowy_model)

plot(stopa0, ylab=colnames(plik)[kol])

lines(dopas, col=4, lwd=3)

par(mfrow=c(1,1), mar=c(5, 4, 4, 2) + 0.1)

Efekt:

Breakpoints at observation number:

m = 1         166

m = 2   40    135

m = 3   40 77 166

Corresponding to breakdates:               

m = 1             1940

m = 2   1814      1909

m = 3   1814 1851 1940

Fit:              

m   0    1    2    3   

RSS 9864 9317 9110 8906

BIC 1623 1620 1625 1631

breakfactor(bai)segment1 breakfactor(bai)segment2 

               0.5975301                3.7677778 

Teraz już dostaliśmy jasną interpretację, choć zmieniła się data breakpoint: do 1939 r. inflacja w USA wynosiła średnio ok. 0.6%, a po tym roku 3.77%. Data jest nieprzypadkowa z trzech powodów:

1) między 1933 a 1939 wprowadzono szereg programów i reform, zwanych New Deal. Program był tak głęboki, że prawdopodobnie wpłynął na trwały wzrost cen,

2) druga wojna św. wprowadziła poważne zmiany w polityce fiskalno-pieniężnej w wielu krajach. Np. na Wikipedii ang. czytamy, że po 2 wojnie św. idee Keynesa zostały powszechnie przyjęte w krajach zachodnich (aż do lat 70-tych, gdy pojawiła się stagflacja),

3) Bank centralny USA - FED - powstał formalnie w 1913, a zaczął działać w 1914 r. Chociaż jego celem była m.in. kontrola inflacji, to de facto stał się maszyną do jej tworzenia. Jak czytamy w "Tajnikach bankowości" system bankowy zmieniono w ten sposób, że tylko Banki Rezerwy Federalnej mogły drukować banknoty [2]. Pierwszej poważnej ekspansji monetarnej dokonano podczas I wojny św., a drugiej podczas Wielkiego Kryzysu od 1929 do 1933 r. 

Na koniec możemy wykonać test homoskedastyczności:

# test homoskedastyczności składnika losowego 

czas = c(1:length(stopa))

bptest(stopa ~ czas)

gqtest(stopa ~ 1)

hmctest(stopa ~ 1)

Efekt:

studentized Breusch-Pagan test

data:  stopa ~ czas

BP = 18.209, df = 1, p-value = 1.98e-05

Goldfeld-Quandt test

data:  stopa ~ 1

GQ = 0.37915, df1 = 123, df2 = 122, p-value = 1

alternative hypothesis: variance increases from segment 1 to 2

Harrison-McCabe test

data:  stopa ~ 1

HMC = 0.71341, p-value = 1

Jedynie BP wskazał na zmienność wariancji - zakładam, że wynika to z występowania warunkowej heteroskedastyczności, być może ARCH. Na razie jeszcze tym efektem się nie zajmuję. 

Zupełnie z boku: analizę dla tego przykładu warto porównać z artykułem FED do likwidacji?! Wówczas porównałem okres inflacji przed powstaniem FEDu z okresem po jego powstaniu. Chociaż dostajemy różnicę ponad 25 lat, to należy pamiętać, że test Bai-Perrona znajduje punkt zmiany struktury. W przypadku inflacji taki punkt naprawdę nie istnieje, bo przecież ten wpływ FEDu był stopniowy. Test uśrednił więc ten proces do punktu. 

Przykład 6. NewConnect

Miesięczne stopy zwrotu od 01.2012 do 10.2022 (źródło: stooq.pl)

Wbrew pozorom nie jest wcale takie proste znaleźć dane finansowe, które mają więcej niż jedną zmianę struktury. Jednak takim przykładem okazuje się osławiony indeks NewConnect. Kod generalnie jest ten sam co dla inflacji, choć zmieniają się założenia (fr, kol, st) i pojawia się nazwa_zmiennej (znaczenie tylko dla wykresu):

# Pod-wstęp

nazwa = "ncindex_m.csv"

plik = read.csv(nazwa, sep=";", dec=",")

if (ncol(plik)==1) {

  plik = read.csv(nazwa, sep=",", dec=".")

}

if (is.numeric(plik[1,1])) {

  daty = paste("01.01.", plik[,1], sep="")

} else { 

  daty = plik[,1] 

  daty = as.Date(daty, tryFormats = c("%d.%m.%Y", "%Y.%m.%d", "%d-%m-%Y", "%Y-%m-%d", "%d/%m/%Y", "%Y/%m/%d"))

#jeżeli ostatnia data przekracza wczorajszą datę, to usuwamy ostatnią obserwację

  if (daty[length(daty)]>=Sys.Date()) {

    plik = plik[-nrow(plik),]

    daty = daty[-length(daty)]

  }

#początkowe założenia i ustawienia

fr = 12

st = TRUE

kol = 5

if (st ==TRUE) {

  par(mfrow=c(3,1), mar=c(2,5,2,2))

  plot(cena, ylab= paste("cena", substr(nazwa, 1, nchar(nazwa)-4)))

  nazwa_zmiennej = paste("stopa zwrotu", substr(nazwa, 1, nchar(nazwa)-4))

} else {

  par(mfrow=c(2,1), mar=c(2,5,2,2))

  nazwa_zmiennej = paste(dimnames(plik[kol])[2], substr(nazwa, 1, nchar(nazwa)-4))

}

rok = as.numeric(format(daty, "%Y"))

mc = as.numeric(format(daty, "%m"))

dz = as.numeric(format(daty, "%d"))

data_startowa = c(rok[1], mc[1])

cena = as.numeric(plik[,kol])

cena = ts(cena, start = data_startowa, frequency = fr)

if (st == TRUE) {

  stopa = diff(cena)/cena[-length(cena)]

} else {

  stopa = cena

}

# Część główna

# test stacjonarności

if (length(stopa)>200) {

  kpss.test(stopa, lshort=FALSE)

}  else {

  kpss.test(stopa)

}

# sprawdzanie autokorelacji

pacf(stopa, plot=FALSE, lag.max=10)

ist_pacf = qnorm((1 + 0.95)/2)/sqrt(sum(length(stopa)))

noquote(c("Istotna autokorelacja będzie od poziomu: ", ist_pacf))

# optymalne arima

arma_stopa = auto.arima(stopa,ic="bic", max.P=3, max.Q=3)

arma_stopa

# ODTĄD ZBĘDNE

# tworzenie optymalnego arma

stopa_ar1 = lag(stopa, k = -1)

stopa_ar2 = lag(stopa, k = -2)

stopa_ma = arma_stopa$residuals

stopa_ma1 = lag(stopa_ma, k = -1)

stopa_ma2 = lag(stopa_ma, k = -2)

dane = ts.intersect(stopa, stopa_ar1, stopa_ar2, stopa_ma1, stopa_ma2)

stopa0 = dane[, "stopa"]

stopa_ar1 = dane[, "stopa_ar1"]

#stopa_ar2 = dane[, "stopa_ar2"]

stopa_ma1 = dane[, "stopa_ma1"]

#stopa_ma2 = dane[, "stopa_ma2"]

# DOTĄD ZBĘDNE

stopa0 = stopa

#test Bai-Perrona

#bai = breakpoints(stopa0 ~ stopa_ar1 + stopa_ar2 + stopa_ma1 + stopa_ma2, breaks = 4)

bai = breakpoints(stopa0 ~ 1, breaks = 4)

summary(bai)

plot(bai)

nowy_model = lm(stopa0 ~ breakfactor(bai) - 1)

rok0 = as.numeric(format(as.yearmon(time(stopa0)[1]), "%Y"))

mc0 = as.numeric(format(as.yearmon(time(stopa0)[1]), "%m"))

dopas = ts(fitted(nowy_model), start=c(rok0, mc0), frequency = fr)

coef(nowy_model)

plot(stopa0, ylab= paste(nazwa_zmiennej, substr(nazwa, 1, nchar(nazwa)-4)))

lines(dopas, col=4, lwd=3)

par(mfrow=c(1,1), mar=c(5, 4, 4, 2) + 0.1)

#test homoskedastyczności składnika losowego 

czas = c(1:length(stopa))

bptest(stopa ~ czas)

gqtest(stopa ~ 1)

hmctest(stopa ~ 1)

Efekt:

KPSS Test for Level Stationarity

data:  stopa

KPSS Level = 0.12962, Truncation lag parameter = 4, p-value = 0.1

Warning message:

In kpss.test(stopa) : p-value greater than printed p-value

Partial autocorrelations of series ‘stopa’, by lag

0.0833 0.1667 0.2500 0.3333 0.4167 0.5000 0.5833 0.6667 0.7500 0.8333 

 0.189  0.028  0.043 -0.024  0.104  0.189  0.088 -0.240 -0.074  0.117 

[1] Istotna autokorelacja będzie od poziomu:  0.172565206655387     

                   

Series: stopa 

ARIMA(0,0,0) with zero mean 

sigma^2 = 0.004339:  log likelihood = 167.85

AIC=-333.7   AICc=-333.67   BIC=-330.84

Breakpoints at observation number:

m = 1            102

m = 2         83 102

m = 3      62 83 102

m = 4   41 62 83 102

Corresponding to breakdates:

m = 1                            2020(7)

m = 2                   2018(12) 2020(7)

m = 3           2017(3) 2018(12) 2020(7)

m = 4   2015(6) 2017(3) 2018(12) 2020(7)

Fit:               

m   0         1         2         3         4        

RSS    0.5597    0.5440    0.4640    0.4562    0.4514

BIC -325.9824 -319.9147 -330.7315 -323.1976 -314.8319

breakfactor(bai)segment1 breakfactor(bai)segment2 breakfactor(bai)segment3 

            -0.008154642              0.063809295             -0.021800546 

studentized Breusch-Pagan test

data:  stopa ~ czas

BP = 7.701, df = 1, p-value = 0.005519

Goldfeld-Quandt test

data:  stopa ~ 1

GQ = 6.0217, df1 = 64, df2 = 63, p-value = 1.005e-11

alternative hypothesis: variance increases from segment 1 to 2

Harrison-McCabe test

data:  stopa ~ 1

HMC = 0.14108, p-value < 2.2e-16

Kwestia techniczna: zaznaczyłem pewną część kodu jako zbędną, którą można usunąć. Nie usunąłem, aby pokazać w miarę ogólny kod, który może służyć do modelowania ARMA(p, q). Ponieważ w tym przypadku zmienna nie jest ograniczona przez zmienne opóźnione, to stopa stanowi zmienną bezpośrednio użytą do modelu. Jednak, żeby zachować ogólność zapisu, zachowuję zmienną stopa0 i to na niej działam. Stąd tutaj wyjątkowo przypisuję stopa0 = stopa. Jeśli model zawierałby części AR lub MA, to po prostu usuwam to przypisanie. 

Przykład NewConnect jest o tyle wyjątkowy, że KPSS wskazuje na stacjonarność stóp zwrotu, optymalny jest ARMA(0,0,0) z zerową średnią, a jednocześnie wyraz wolny jest niestabilny: 

* do końca 2018 r. wynosił -0,8%, 

* do lipca 2020 +6,4% 

* do dziś -2,2%. 

Tak się zmieniała się średnia miesięczna stopa zwrotu. Dlaczego więc KPSS nie potwierdza tego? Odpowiedź już była podana przy badaniu inflacji - nie wiemy czy przypadkiem między 2018 a 2020 nie powstał proces ARMA, który następnie znowu zanikł tak że sumarycznie BIC wskazuje optimum przy zerowych autoregresjach. Można by jeszcze pomyśleć o wpływie zmian wariancji na stabilność - wszystkie testy na wariancję dowodzą jej zmienności w czasie. Jednakże zarówno KPSS jak i test Bai-Perrona są odporne na heteroskedastyczność. Stąd KPSS nie wykrył zmian w rozkładzie stóp, natomiast Bai-Perron wykrył, bo jest wrażliwy na budowę modelu, który w dodatku musi być liniowy.

 

Analiza, którą przedstawiłem ma szczególne znaczenie przy prognozowaniu, gdyż widzimy, że niektóre procesy są niestabilne na poziomie zarówno średniej, jak i wariancji. W przypadku inflacji USA powinniśmy po pierwsze ustawić "uczenie się" parametrów przez program od roku 1941, a po drugie rozważyć model klasy (G)ARCH zamiast ARMA. Te pierwszy obejmuje drugi i czasem niepotrzebnie nazywany jest ARMA-GARCH. Z kolei poświęcanie czasu na NewConnect wydaje się złym pomysłem dopóki średnia jest ujemna - i to nie do końca z powodu bessy, bo to giełda, która od początku swojej historii spadła o 70% (od 2012 o ok. 30%).

Literatura:

[1] Bai J., Perron P. (2003), Computation and Analysis of Multiple Structural Change Models, Journal of Applied Econometrics, 18, 1-22 ;

[2] Rothbard, M. N.., Tajniki bankowości, Warszawa 2007, s. 254.