czwartek, 13 lipca 2023

Mentzen ma rację co do euro

W debacie, jaka odbyła się niedawno pomiędzy R. Petru a S. Mentzenem, skupiono się głównie na wizerunkowych błędach tego drugiego, przez co niemal jednogłośnie uznano, że przegrał on tę potyczkę. Rzadko można znaleźć poważniejszą ocenę merytoryczną tego występu. Jako z jeden z nielicznych, portal oko.press podjął się tego w art. Debata Mentzen – Petru. Cztery bzdury o ZUS-ie, Orlenie, kodeksie pracy i euro . A więc już w tytule dowiadujemy się, że będą bzdury o euro. I analiza została przeprowadzona nie przez byle kogo, bo, cytuję, przez ekonomistę - Tomasza Makarewicza - juniora profesora na Uniwersytecie w Bielefeld i członka grupy eksperckiej Dobrobyt na Pokolenia.

Chciałbym się skupić na odpowiedzi Mentzena czy powinniśmy przystąpić do strefy euro. Mentzen: 

Nie, wspólna waluta może być tylko na optymalnym obszarze walutowym. Unia Europejska nigdy nie była optymalnym obszarem walutowym, co najwyżej Niemcy, Austria i kraje Beneluksu. Stąd właśnie problemy Grecji, Włoch czy Hiszpanii po wejściu do strefy euro - te kraje przestały się rozwijać.

Komentarz Makarewicza:

To jest temat, na który nie ma jednej odpowiedzi, ale zasadniczo wszyscy zgadzają się, że kwestia jest skomplikowana. Te kraje miały pewne strukturalne problemy także przed wejściem do strefy euro. W przypadku Grecji to był problem oligarchizacji gospodarki. Powiązanie kilku nielicznych klanów, do których należą ważniejsze zakłady gospodarcze, i świata polityki. Można się zastanawiać, czy euro spotęgowało te problemy, czy nie. Ekonomiści wciąż się o to spierają. Ale nie jest tak, że te kraje rozwijały się świetnie, po czym przystąpiły do strefy euro i padły.

A więc Makarewicz nie obala stwierdzenia czy opinii Mentzena. Uznaje jedynie sprawę za bardziej skomplikowaną. Ale uwaga - dopuszcza możliwość, że euro mogło spotęgować problemy. Tylko przecież tytuł stwierdza, że Menzten podał bzdurę. To gdzie ta bzdura? Jeżeli dziennikarz tak sobie dopisał, żeby wyszedł mocniejszy przekaz, to mamy do czynienia z manipulacją. Również źle to świadczy o Makarewiczu i jego organizacji, która na swoim profilu FB promuje ten artykuł. Żadnego sprostowania nie ma.

Dla przypomnienia jakiś czas temu napisałem art. Wejście do strefy euro nieopłacalne? , które wskazuje przesłanki, że euro rzeczywiście jest nieopłacalne dla Polski.

Kolejna sprawa. Makarewicz w ogóle nie odnosi się do kwestii optymalnego obszaru walutowego, o którym Mentzen wspomina dwukrotnie. Pytanie brzmi: dlaczego? Czy dlatego, że okazałoby się, że Mentzen powiedział prawdę, czy dlatego, że Makarewicz nie zna się na tym temacie? Jeżeli pomija, bo wie, że to niewygodna prawda, to traci wiarygodność obiektywnego eksperta. Jeżeli pomija, bo nie zna teorii optymalnego obszaru walutowego, to w ogóle traci miano eksperta. W każdym przypadku oko.press się kompromituje, atakując Mentzena, nie wiedząc o czym on mówi.

Zanim przejdę dalej, muszę od razu podkreślić, że ja też nie jestem ekspertem od tego tematu, ale powołam się na takich. Już w tym artykule Wejście do strefy euro nieopłacalne? napisałem, że P. Krugman - noblista w dziedzinie ekonomii - stwierdził w 2011 r., że Europa znajduje się w "głębokim kryzysie" właśnie z powodu używania jednej waluty. Czyli wprost potwierdza zdanie Mentzena. 

Ważne jest, aby zdać sobie sprawę, że zdanie Krugmana, a więc i Mentzena, nie jest jakąś tam opinią - parafrazując klasyka - ekonomistów przy kawie i ciastkach. Optymalny obszar walutowy to termin ukuty przez R. Mundella w 1961 r. [3].

Ogólnie są dwie skrajne sytuacje: cały świat posługuje się jedną walutą albo każdy człowiek ma swoją własną walutę. Ten drugi przypadek sprowadza się de facto do barteru. Teraz pomiędzy tymi dwiema skrajnościami możemy znaleźć optymalny obszar walutowy. Może to być np. rodzina, gmina, cały kraj, grupa krajów albo kontynent. Najbardziej naturalną formą obszaru walutowego jest oczywiście kraj, ale nie znaczy to, że jest ona optymalna. Barter przecież też jest naturalny, ale nikt powie, że jest lepszą formą wymiany od pieniądza w nowoczesnej gospodarce. Ale tak samo mało wygodne jak barter byłoby wprowadzenie osobnej waluty dla każdej gminy. Wynika to z ciągłego przemieszczania się towarów i ludzi między jednym regionem kraju a drugim, co sprawia, że oba regiony traktowane są jak całość. 

Mundell analizuje sprawę następująco. Powiedzmy, że popyt na jakieś dobra przesuwa się z gminy A do gminy B (w A spada, w B rośnie). Podaż pozostaje stała. Przedsiębiorstwa A mają dwa wyjścia: albo eksportować towary do B, albo obniżyć ceny towarów. Niestety pierwsze niesie za sobą dodatkowe koszty, a drugie obniża przychody, tak że warianty oba ograniczają zyski. Mniejsze zyski zmniejszą inwestycje lub dywidendy (a więc konsumpcję), a to z kolei pociągnie za sobą wzrost bezrobocia i PKB w gminie A. Jednocześnie w gminie B rośnie popyt, ale nie podaż, a więc jedynym efektem będzie wzrost cen w B. Jednak ekonomia nie znosi próżni: firmy będą starać się zasilić popyt, zwiększając produkcję i sprzedaż u siebie. Ale do tego potrzebni są pracownicy, których firma może jedynie ściągnąć z A, ponieważ tam są właśnie bezrobotni. Oczywiście zawsze jakieś bezrobocie w gminie B też będzie, ale skoro nie podejmują jeszcze pracy, to możliwe, że wymagają wyższej płacy, podczas gdy bezrobotni z A nie mają takich wymagań. Tak więc bezrobocie w A zostaje automatycznie wyeliminowane przez mechanizmy rynkowe, niezależnie od tego czy pracownicy będą dojeżdżać z A do B czy się przeprowadzą (mniej ludzi oznacza mniejszą stopę bezrobocia). Gorzej, gdyby gminy były od siebie sporo oddalone - wtedy nowo zatrudnieni musieliby się przeprowadzić, a nie byłoby pewności czy są tak mobilni. Bez tej mobilności w gminie A nastanie recesja, a w gminie B inflacja na skutek wymuszenia wyższych płac.

Co by się jednak stało, gdyby gminy miały inne waluty, w dodatku w systemie kursów płynnych? W gminie A jest waluta PLN A, w gminie B waluta PLN B. Powiedzmy, że zaraz po spadku popytu w A kurs PLN A spada. O ile na początku 1 PLN A = 1 PLN B, o tyle teraz 1 PLN A = 0,5 PLN B. Powoduje to, że w gminie B ten sam produkt można kupić o połowę taniej w gminie A. Gmina B będzie zmotywowana do kupowania tańszych towarów od A, tak że A nie musi ich już eksportować i ponosić dodatkowych kosztów. Oznacza to, że popyt w gminie A wróci do punktu wyjścia. Skąd wiadomo, że do punktu wyjścia? Z jednej strony wyższy popyt w gminie B zostaje zaspokojony dzięki nowej podaży powstałej w wyniku wzrostu importu z A do B (ceny w B pozostaną więc stałe mimo wzrostu popytu). Czyli nadwyżkowy popyt w B równa się nadwyżkowej podaży w B. Ale z drugiej strony popyt importowy równa się podaży na import z A, a to oznacza, że popyt importowy musi się równać nadwyżkowemu popytowi w B. W sumie logika nakazuje przyjąć, że cały popyt w A wróci do punktu wyjścia. Tym samym bezrobocie znika i PKB wraca do punktu wyjścia.

W powyższej analizie uwzględniliśmy jedynie dobra ruchome i płynne. Co jednak np. z nieruchomościami? W tej sytuacji eksport wykluczamy. Są jednak tu dwie możliwości. Pierwsza taka, że firmy z A przenoszą produkcję oraz dział sprzedaży do B. Druga taka, że firmy w B mając u siebie nadwyżkę popytu, zwiększą produkcję. Skupmy się na drugim wariancie. Problem w gminie B polega na znalezieniu rąk do pracy, ponieważ wszyscy gdzieś pracują, tzn. nie ma bezrobocia. Deweloper mógłby ściągnąć bezrobotnych z gminy A, ale przecież ciągle zakładamy, że nie są oni mobilni. Żeby uzyskać nową siłę roboczą, deweloper musi podnieść stawki robotnikom mieszkającym w B, zachęcając ich do zmiany bieżącego pracodawcy. Oznaczałoby to spadek rentowności, a może nawet i straty firmy. Aby nie dopuścić do takiej sytuacji, deweloper musi ograniczyć inne koszty, takie jak materiały. Wykorzystuje zatem fakt, że 1 PLN A = 0,5 PLN B, co oznacza, że średnio biorąc materiały może kupić dwa razy taniej w gminie A. Dopóki koszt przewozu zakupionych materiałów z A do B nie będzie zbyt wysoki, deweloper będzie mógł zredukować nadwyżkowe koszty. Większy eksport, tj. zakupy dewelopera, nieco poprawi PKB w gminie A, redukując częściowo tamtejsze bezrobocie.

Pozostaje jeszcze niewyjaśniona kwestia przyczyny deprecjacji waluty w gminie A. Co może być przyczyną tego spadku zaraz po spadku popytu? Ponieważ mamy tu na myśli bardziej popyt wewnętrzny, to nieprzekonujące jest twierdzenie, że kurs walutowy wynika ze spadku importu do A (wtedy sprawa byłaby oczywista - A otrzymuje pieniądze w PLN B, więc musi zamienić na PLN A, a więc przy braku importu nie ma też tej zamiany, tj. popytu na PLN A). Odpowiedzią jest za to obniżona stopa procentowa.

Spadek aktywności gospodarczej powoduje, że banki muszą obniżyć stopę procentową, ponieważ popyt na pieniądz spada. Z Wicksellowskiego punktu widzenia stopa proc. musi po prostu zrównać się z naturalną stopą procentową (zob. Dodatek nr 2  w tym artykule).

Na stopę procentową można by spojrzeć jeszcze inaczej. Charakterystyczną cechą danego obszaru walutowego jest to, że każdy taki obszar posiada własny bank centralny, który kontroluje ilość pieniądza w danej walucie. Jeżeli bank centralny oczekuje recesji, obniża stopę procentową, "wyręczając" niejako tym samym banki komercyjne w ustalaniu stóp procentowych.

Trzeba w tym miejscu zauważyć, że podejście od strony kontrolnej banku centralnego może prowadzić do błędnych wniosków. W prowadzonej analizie nie pozwalam na jakiekolwiek zaburzenie mechanizmu rynkowego. Samo obniżenie stopy procentowej do poziomu naturalnego nie spowoduje żadnej poprawy PKB. Sprawia jedynie, że stopa się dopasowuje do danego poziomu gospodarki (jest to skutek, a nie przyczyna zmian PKB). Natomiast, gdy mówimy o roli banku centralnego, to od razu myślimy właśnie o zaburzeniu tego mechanizmu, tzn. wpływaniu na aktywność gospodarczą. Gdybyśmy szli tym tropem w naszym przykładzie, to pozostałe w gminie A bezrobocie można zlikwidować zmniejszając jeszcze bardziej stopę procentową, tzn. poniżej stopy naturalnej. 

W każdym razie niższa stopa proc. zmniejsza popyt na lokaty i obligacje w danej walucie. Stąd region B zmniejsza skłonność do zamiany PLN B na PLN A. I dlatego kurs walutowy PLN A spada.

W analizie Mundella kluczowym czynnikiem jest mobilność pracowników. Im mniejsza mobilność, tym bardziej prawdopodobne, że różne waluty będą korzystne dla kraju. Gdybyśmy mieli oddalone gminy albo duże obszary kraju, mobilność musiałaby być mniejsza, ale wtedy dla większego kraju mogłoby się okazać, że lepsze są różne waluty. 

Krugman [2] zwraca jednak jeszcze uwagę na argument P. Kenena [1], kiedy to wspólny obszar walutowy może okazać się korzystniejszy. Pisałem już o tym, że w nowoczesnej gospodarce bank centralny kontroluje częściowo gospodarkę całego obszaru walutowego. Innym organem centralnym jest rząd lub samorząd. W tym wypadku nie możemy utożsamiać rządu z państwem, ponieważ mówimy o dowolnym obszarze.  Rząd pobiera podatki i płaci z nich zasiłki dla bezrobotnych oraz wydaje na różne programy pomocowe, szkolenia i inne inwestycje. Tak więc jeśli gmina A wpada w recesję, to płaci mniejsze podatki niż gmina B. Jednocześnie jednak gmina A korzysta z zasiłków i innych wydatków rządowych, które płyną dzięki wyższym podatkom z gminy B. 

Zastąpmy gminy krajami. Kraje mają niezależne od siebie rządy, ale waluty mogą mieć wspólne. Bez wspólnego rządu każdy kraj będzie na łasce i niełasce wspólnego banku centralnego. Jeżeli kraj A wpadł w recesję, a B w ożywienie, to co teraz miałby zrobić bank centralny? Musi obniżyć stopy, aby pomóc krajowi A, ale zaszkodzi B, w którym wywoła inflację. Rząd B musiałby obciąć wydatki, aby zneutralizować inflację. Kraj A z kolei być może wolałby mocniej obniżyć stopy, czego jednak bank centralny już nie robi, ponieważ musi ważyć skutki dla wszystkich krajów. A jeśli kraj A nie może na to liczyć, powinien zwiększać wydatki, co będzie utrudnione z powodu recesji. W sumie nie ma żadnej wartości dodanej przebywanie w takim układzie. Oczywiście jeżeli kraje przyjęłyby własne waluty, to musiałyby wziąć pod uwagę koszty transakcyjne ich wymiany. W czasach Internetu koszty te jednak powinny spadać z uwagi na łatwą komunikację obu stron wymiany.

Gdyby jednak istniał wspólny rząd, to mógłby lepiej regulować przepływy pieniężne między krajami. Zapomnijmy o banku centralnym, żeby nie komplikować analizy. Kraj B przechodzi dobrą koniunkturę, kraj A złą. Kraj B wpłaca więc większe podatki, a A mniejsze, ale dostaje zastrzyk zasiłków dzięki podatkom pochodzącym z B. Wtedy oczywiście B niczego nie zyskuje, ale pamiętajmy, że gdy koniunktura się odwróci i w A nastąpi ożywienie, a w B depresja, to B zyska. Widzimy więc, że argument Kenena opiera się na finansowej koncepcji hedgingu albo dywersyfikacji ryzyka.

Strefa euro stanowi próbę wprowadzenia w życie teorii Mundella i Kenena. Próba ta nie może zakończyć się sukcesem, o czym mówi Krugman. Mobilność siły roboczej z jednego kraju UE do drugiego nie jest wystarczająca, aby podpierać się tym czynnikiem do wdrażania euro. Podobnie integracja fiskalna jest daleka od doskonałości, właśnie dlatego, że nie ma prawdziwego rządu UE. Żeby idea euro miała sens, kraje musiałyby w dużym zakresie zrezygnować z własnej polityki podatkowej i socjalnej, a polegać na decyzjach Komisji Europejskiej będącej odpowiednikiem rządu. Krugman wprost nazywa strefę euro błędem.

Jak pisałem, nie jestem specjalistą, więc nie wiem czy euro jest błędem czy nie. Na pewno przy obecnej strukturze UE strefa euro nie jest optymalna. Dla porównania USA wydają się optymalnym obszarem walutowym, co potwierdza Krugman. Chodzi tu zarówno o większą mobilność ludności z jednego stanu do drugiego, jak i większą integrację socjalno-fiskalną.

Ostatnie lata raczej wskazują na dezintegrację UE. W ciągu 15 lat, od 2008 r., eurodolar spadł o 30% i jeśli nie nastąpi poprawa pod tym względem, to spodziewałbym się jego stopniowego upadku.    


Literatura:

[1] Kenen, P.. 1969. The Theory of Optimum Currency Areas: An Eclectic View. In Monetary Problems of the International Economy, edited by R. Mundell and A. Swoboda, 41–60. Chicago: University of Chicago Press.

[2] Krugman, P., 2013. Revenge of the optimum currency area. NBER Macroeconomics Annual, 27 (1) (2013), pp. 439-448

[3] Mundell, R.. 1961. A Theory of Optimum Currency Areas. American Economic Review 51 (4): 657–65.

niedziela, 25 czerwca 2023

Jak sprawdzić stabilność każdego współczynnika osobno (język R)?

Kontynuuję temat analizy stabilności w języku R. Jest to kontynuacja serii 1, 2, 3, 4, 5. Dotychczas pokazywałem sposób testowania stabilności całego modelu. Był to dobry sposób, aby sprawdzić czy i w którym miejscu zmienił się wyraz wolny lub współczynnik kierunkowy, czyli wpływ zmiennej objaśniającej. Jeżeli jednak model zawierał zarówno wyraz wolny jak i współczynnik kierunkowy, a chciałem ocenić co dokładnie i kiedy się zmieniło, zaczynały się problemy. Dotychczas radziłem sobie w ten sposób, że testowałem stałość średniej poprzez test stacjonarności zmiennych i jeśli została ona zachowana, dopiero testowałem stabilność modelu. Jeśli model okazał się niestabilny, wnioskowałem, że współczynnik kierunkowy się zmienił, ponieważ wyraz wolny pełnił rolę średniej zmiennej objaśnianej (która z testu stacjonarności okazała się stała). To rozwiązanie jednak odpada, jeśli model zawiera wiele zmiennych objaśniających - nie będę wiedział, która zmienna "zmieniła" współczynnik.

Będą potrzebne 2 pakiety: strucchange oraz forecast. Jednak postanowiłem włączyć do swoich analiz trzeci pakiet, xts. 

Czasem może się też przydać zmiana opcji programu, ponieważ wolę unikać notacji naukowej. Często dostajemy małe liczby, np. 0.0001. Normalnie program wyświetli 1e-04. Żeby dostać zapis dziesiętny, moża użyć funkcji options(scipen = 1). Jednak jeszcze mniejsza liczba, jak 0.00001 znów zostanie zapisana naukowo. Żeby przywrócić zapis dziesiętny, musimy zwiększyć scipen do 2. Itd. Domyślne ustawienia to scipen = 0. Czyli jeśli część dziesiętna ma powyżej 3 zer, to zwiększamy scipen o odpowiednio dużą liczbę. Wystarczy, że ustawimy scipen = 100. Trzeba zauważyć, że nie chodzi tu tylko o część dziesiętną, ale ogólnie o liczbę cyfr. Dla scipen = 100 dopiero po wystąpieniu powyżej 103 cyfr w całej liczbie dostaniemy notację naukową. 

Dodatkowym problemem jest nie rzadko wyświetlanie wielu cyfr w części dziesiętnej. Chodzi tu o format, a nie zaokrąglenie. Aby kontrolować ten format stosujemy parametr digits w funkcji options. Domyślnie program używa digits = 7, co oznacza, że maksymalnie wyświetli się 7 cyfr. Np. liczbę 555.00001 pokaże jako 555. Aby zwiększyć precyzję, musimy ustawić digits = 8 lub więcej. W naszym przypadku lepiej będzie na odwrót - zmniejszymy precyzję, ponieważ nasze liczby to zazwyczaj ułamki dziesiętne i nie potrzebujemy aż 7 cyfr po kropce czy przecinku. Dla obecnych potrzeb ustawię digits = 4. 

Dla przypomnienia kody dzielę zazwyczaj na 3 części:

1. Przedwstęp - najbardziej wstępne ustawienia, nie trzeba nic zmieniać

2. Wstęp - ustawianie ścieżek i modelu, tutaj wprowadzamy własne dane

3. Część główna - dotyczy głównego tematu, zazwyczaj nie trzeba nic zmieniać.

Z powodu takiego podziału przykłady zawierają jedynie Wstęp i Część główną.

# Przedwstęp

#precyzja i sposób zapisu liczb: 

options(scipen = 100, digits = 4)

#pakiety

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

  install.packages("strucchange")

  library("strucchange")

}

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

  install.packages("forecast")

  library("forecast")

}

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

  install.packages("xts")

  library("xts") 

}

Pakiet strucchange jest bardzo bogaty w testy stabilności. Dotychczas stosowałem testy bazujące na statystyce F. Pakiet jednak oferuje także drugą grupę testów - fluktuacji procesów empirycznych (The Empirical Process, EFP). Pierwsza grupa może wskazać co najwyżej momenty zmiany struktury modelu. Druga grupa w pewnym sensie jest bardziej ogólna, bo pokazuje jak dochodzi do tej zmiany i co się dzieje po tej zmianie - punkt zostaje zastąpiony przedziałem w czasie. Ma to duże znaczenie dla zrozumienia metodologii i interpretacji tego typu testów. W teście F dzielono próbę na dwie (różne) części i sprawdzano czy gdzieś w modelu empirycznym następuje zmiana. Tutaj natomiast sprawdza się zachowanie modelu w różnych przedziałach próby na tle procesu czysto losowego. W stabilnym modelu empiryczna ocena parametru zawsze losowo fluktuuje wokół teoretycznego parametru. Jeśli jej wariancja staje się za duża, znaczy, że proces przestaje być losowy. Aby określić czy mamy do czynienia ciągle z fluktuacją losową, wykorzystuje się centralne twierdzenie graniczne funkcjonałów. Zgodnie z tym twierdzeniem suma wielu zmiennych niezależnych ze stałą średnią i wariancją będzie w (dodatkowym) wymiarze czasowym błądzeniem przypadkowym. Dlatego w teście sumowane są kolejne odchylenia od parametru teoretycznego w danym przedziale próbnym i sprawdza się czy to zsumowane odchylenie tworzy błądzenie przypadkowe. Jeśli fluktuacja jest za duża i wychodzi poza błądzenie przypadkowe, to znaczy, że wartość prawdziwego parametru się zmieniła. Pełen zakres czasu wynosi (0, 1, ... t, ... T). Kolejne przedziały testowania są budowane na 2 sposoby, (a) od 0 do t albo (b) od t-h do t+h. Przedział (a) dotyczy skumulowanych sum reszt standardowych (ang. cumulative sums of standardized residuals, CUSUM) oraz estymacji rekurencyjnych (recursive estimates, RE), natomiast (b) ruchomych sum reszt (moving sums of residuals, MOSUM) oraz estymacji ruchomych (moving estimates, ME). W przypadku RE i ME odchylenia dotyczą nie tyle sum, co średnich w danym przedziale. Wszystkie te testy zostały powiązane w jedną klasę tzw. uogólnionego testu fluktuacji (The Generalized Fluctuation Test [1]), a ich procesy w uogólniony empiryczny proces fluktuacji (Generalized Empirical Fluctuation Process, GEFP).
Z moich obserwacji wynika, że nie ma najlepszego testu fluktuacji. Dla mniejszych prób ME najlepiej radzi sobie z odkrywaniem nagłych zmian w symulowanych modelach. Żeby mieć pełen obraz należy wykonać wszystkie testy i dopiero wtedy wyciągać wnioski.
Są różne typy testu CUSUM. Standardowy test tego rodzaju opisałem już kiedyś w tym art., gdy swoje badania wykonywałem w gretlu. Nas interesują tylko te, które umożliwiają analizę poszczególnych parametrów. Są to Score-CUSUM, Score-MOSUM, ME, RE i GEFP.
Pierwsze próby wykonam na symulacjach, aby sprawdzić jak test sobie radzi. Symulowana próba składa się z 1500 danych podzielona na 3 identyczne części. Najpierw zrobię przykład stabilnego ARMA(1,1).

Przykład 1. Brak zmian: ARMA(1, 1)

#Wstęp
set.seed(1980)

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

sim123 = ts(c(sim1, sim2, sim3))

#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)])
sim.model <- sim ~ sim_ar1 + sim_ma1


# Część główna

plot(sim)



Są dwie funkcje do użytku: efp() i gefp(). 

# EFP: CUSUM, MOSUM, ME, RE

# jest wiele podtypów testu CUSUM, związane metodą estymacji: OLS, Recursive CUSUM, Score-CUSUM. Tylko ten ostatni pozwala przeanalizować wszystkie współczynniki na wykresie. Podobnie jest z MOSUM.

#Score based CUSUM
# Sumarycznie:
cus = efp(sim.model, type="Score-CUSUM")
plot(cus, functional="max")


# Osobno dla każdego parametru:

cus = efp(sim.model, type="Score-CUSUM")
plot(cus, functional=NULL)



Aby zobaczyć zmiany w każdym parametrze, trzeba użyć functional=NULL.
Zwróćmy uwagę, że testujemy tutaj także wariancję modelu. 

#Score based MOSUM
mos = efp(sim.model, type="Score-MOSUM")
plot(mos, functional="max")



MOSUM wskazał błędnie zmianę struktury modelu - została przebita linia czerwona. Zobaczmy, który czynnik wywołał taki błąd:

plot(mos, functional=NULL)



A więc test błędnie wskazał na zmianę struktury w wariancji modelu. 

#Moving estimate
me <- efp(sim.model, type="ME")
plot(me, functional="max")



plot(me, functional = NULL)



#recursive estimate
par(mfrow=c(3,1), mar=c(2.5,5,1.5,3))
re <- efp(sim.model, type="RE")
plot(re, functional="max")


plot(re, functional = NULL)


#GEFP
par(mfrow=c(3,1), mar=c(2.5,5,1.5,3))
gcus = gefp(sim.model)
plot(gcus, aggregate=TRUE)



plot(gcus, aggregate=FALSE)


Wszystkie testy, za wyjątkiem Score-MOSUM, dały prawidłową odpowiedź.


Przykład 2. Zmiana wyrazu wolnego ARMA(1, 1)

#Wstęp

set.seed(1980)

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

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

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

sim123 = ts(c(sim1, sim2, sim3))

#tworzenie modelu
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)])

plot(sim)



#Część główna
# CUSUM, MOSUM, ME, RE, GEFP

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

#Score based CUSUM
cus = efp(sim.model, type="Score-CUSUM")
plot(cus, functional="max")




Przebicie linii czerwonej dowodzi, że CUSUM prawidłowo rozpoznał dwukrotną zmianę struktury. Sprawdźmy, który czynnik wywołał zmianę.

plot(cus, functional=NULL)

Faktycznie wyraz wolny dwukrotnie wystaje poza losowy obszar, więc CUSUM prawidłowo rozpoznaje niestabilność.


#Score based MOSUM
mos = efp(sim.model, type="Score-MOSUM")
plot(mos, functional="max")



MOSUM - prawidłowo, zróbmy analizę poszczególnych współczynników:

plot(mos, functional=NULL)


MOSUM zupełnie inaczej przedstawia sytuację od CUSUM. Wg MOSUM pierwsza część danych (sim1) jest "normalna", druga (sim2) wskazuje na zmienioną strukturę w całym zakresie, a trzecia (sim3) powrót do "normalności". Normalność wynika z tego, że więcej czasu (2/3) zajmuje w por. do "nienormalności" (1/3).

#Moving estimate
me <- efp(sim.model, type="ME")
plot(me, functional="max")



W agregatach prawidłowo.

plot(me, functional = NULL)



W dezagregatach ME wskazuje poprawnie, że zmienił się wyraz wolny, ale błędnie, że dodatkowo zmienił się wpływ danych z poprzedniego okresu. Jest to zastanawiające, bo 1500 obserwacji to duża próba i nie powinno być takich błędów.

#recursive estimate
re <- efp(sim.model, type="RE")
plot(re, functional="max")


re <- efp(sim.model, type="RE")
plot(re, functional = NULL)


RE także informuje, że nastąpiła dwukrotnie zmiana na plus i na minus tylko dla wyrazu wolnego - OK.


#GEFP
gcus = gefp(sim.model)
plot(gcus, aggregate=TRUE)



plot(gcus, aggregate=FALSE)


GEFP pokazał to samo co CUSUM i RE - prawidłowo.


W omówionym przykładzie najgorzej poradził sobie ME.  


Przykład 3. Zmiana z ARMA(1, 1) na ARMA(0, 1) i powrót do ARMA(1, 1)

#Wstęp

set.seed(1980)

sim1 = arima.sim(list(order=c(1,0,1), ar=0.6, ma=0.3), n=500, mean=1)
sim2 = arima.sim(list(order=c(1,0,1), ar=0.000001, ma=0.3), n=500, mean=1) 
sim3 = arima.sim(list(order=c(1,0,1), ar=0.6, ma=0.3), n=500, mean=1)
sim123 = ts(c(sim1, sim2, sim3))

#tworzenie modelu
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)])

sim.model <- sim ~ sim_ar1 + sim_ma1

plot(sim)



#Część główna
# CUSUM, MOSUM, ME, RE, GEFP
#Score based CUSUM
par(mfrow=c(3,1), mar=c(2.5,5,1.5,3))
cus = efp(sim.model, type="Score-CUSUM")
plot(cus, functional="max")



plot(cus, functional=NULL)


CUSUM wprawdzie ogólnie wskazał dwa momenty zmiany, ale osobno raz pomylił wyraz wolny ze zmianą parametru AR(1), zaś dla drugiej zmiany pokazał prawidłowo. 

#Score based MOSUM
mos = efp(sim.model, type="Score-MOSUM")
plot(mos, functional="max")


plot(mos, functional=NULL)

MOSUM jeszcze gorzej od CUSUM - wskazał zmiany wyraz wolnego i MA(1). Część AR(1) tylko sygnalizowała zmianę.

#Moving estimate
me <- efp(sim.model, type="ME")
plot(me, functional="max")



 plot(me, functional = NULL)




ME tym razem pokazał prawidłowo zmianę par. AR(1), a wyraz wolny nie przebił linii czerwonej. 

#recursive estimate
re <- efp(sim.model, type="RE")
plot(re, functional="max")




plot(re, functional = NULL)



RE przynosi podobny efekt do CUSUM - błędnie informuje o zmianie wyrazu wolnego. 

#GEFP
gcus = gefp(sim.model)
plot(gcus, aggregate=TRUE)



plot(gcus, aggregate=FALSE)

W przypadku GEFP to samo co dla poprzednika.

Podsumowując ten przykład, ME dla odmiany okazał się tu najlepszy. 


Przykład 4. Zmiana części MA

Ponieważ proces MA(1) można wyrazić w postaci procesu AR(T), gdzie T to cały zakres próby (zob. ten wpis), a więc pozostaje częściowo zależny od AR(1), zmiana zachowania MA jest problematyczna w testach. Aby uzyskać efekt, obniżę 0,3 do -0,4. Prawidłowe połączenie prób wymaga wtedy trudniejszego zapisu, który analizowałem tutaj. Wtedy zapis wstępny będzie następujący:

#Wstęp

set.seed(1980)

dlugosc_proby = 500

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

sim2 = arima.sim(list(order=c(1,0,1), ar=0.6, ma=-0.4), n=dlugosc_proby, mean=1, start.innov=c(sim1[dlugosc_proby-1], sim1[dlugosc_proby]), n.start=2)

sim12 = c(sim1, sim2)

sim3 = arima.sim(list(order=c(1,0,1), ar=0.6, ma=0.3), n=dlugosc_proby, mean=1, start.innov=c(sim12[dlugosc_proby-1], sim12[dlugosc_proby]), n.start=2)

sim123 = ts(c(sim12, sim3))

#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)])

sim.model <- sim ~ sim_ar1 + sim_ma1

plot(sim)




Od teraz kod będzie zawierał tylko wykres dla każdego współczynnika osobno.

#Część główna
# CUSUM, MOSUM, ME, RE, GEFP

#Score based CUSUM
par(mfrow=c(3,1), mar=c(2.5,5,1.5,3))
cus = efp(sim.model, type="Score-CUSUM")
plot(cus, functional=NULL)



CUSUM błędnie wskazuje na zmianę każdego współczynnika.

#Score based MOSUM
par(mfrow=c(3,1), mar=c(2.5,5,1.5,3))
mos = efp(sim.model, type="Score-MOSUM")
plot(mos, functional=NULL)


MOSUM wskazał małą (dwukrotną) zmianę wyrazu wolnego i sporą zmianę w MA. Kiedy spojrzymy na wykres, to zauważymy, że faktycznie wydaje się, jakby średnia dla sim2 spadła. W rzeczywistości nie spadła tylko utrzymuje się na poziomie 1,4, a to właśnie sim1 i sim3 mają za wysokie średnie: dla sim1 3,12, sim3 3,41. Nie ma tu błędu, a po prostu silne autokorelacje podwyższają średnią empiryczną. MOSUM obiektywnie więc prawidłowo informuje nas o zmianach na wykresie, bo główną przyczynę utraty stabilności dostrzega w części MA. Część AR nie przebiła linii czerwonej, a wyraz wolny ledwo się poza nią wychylił.


#Moving estimate
par(mfrow=c(3,1), mar=c(2.5,5,1.5,3))
me <- efp(sim.model, type="ME")
plot(me, functional = NULL)


ME z kolei błędnie widzi zmianę w AR, a niewiele zmian w MA.


#recursive estimate
par(mfrow=c(3,1), mar=c(2.5,5,1.5,3))
re <- efp(sim.model, type="RE")
plot(re, functional = NULL)



RE kompletnie zawodzi w tym przypadku.


#GEFP
par(mfrow=c(3,1), mar=c(2.5,5,1.5,3))
gcus = gefp(sim.model)
plot(gcus, aggregate=FALSE)


Podobnie jak CUSUM - błędnie.

W przykładzie 4 najlepszym okazał się MOSUM. I teraz idziemy z żywymi przykładami.


Przykład 5. WIG tygodniowy w korelacji z S&P 500. 
Okres: 02.01.2000 - 18.06.2023 (1225 obserwacji). Źródło: stooq.pl

Dla przypomnienia, ponieważ działam na Windowsie, używanie ścieżek w R nastręcza problemów, stąd kod przekształcający "\" w "/".
 

# Wstęp
#zamieniam "\" na "/"

sciezka = r"(C:\R\testy)"
sciezka = gsub("\\", "/", sciezka, fixed=T)

# albo
# sciezka = gsub("\\\\", "/", sciezka)
# ustawiamy folder roboczy

setwd(sciezka)

# pliki - źródło stooq.pl
nazwa1 = "^spx_w.csv"
nazwa2 = "wig_w.csv"

nazwa <- c(nazwa1, nazwa2)
dane <- list()
stopa <- list()


for (i in 1:2) {

if (grepl(";", readLines(nazwa[i], n=1))) {
plik = read.csv(nazwa[i], sep=";", dec=",")
} else if (grepl(",", readLines(nazwa[i], n=1))) {
plik = read.csv(nazwa[i], sep=",", dec=".")
} else {
stop("Unknown separator in the file: ", nazwa[i])
}

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"))
daty = as.Date(daty)
#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)]
}

# inne dla pliku:

st = TRUE
k = 5

zmienna = as.numeric(plik[,k])
dane[[i]] <- na.omit(data.frame(daty, zmienna))
dane[[i]] <- as.xts(dane[[i]])

# warunek czy stopy zwrotu
if (st==TRUE) {
stopa[[i]] = round(diff(as.numeric(dane[[i]][,1]))/as.numeric(dane[[i]][,1][-length(dane[[i]][,1])]), 4)
stopa[[i]] = as.xts(na.omit(data.frame(daty[-1], stopa[[i]])))
}

dane_merged <- merge.xts(dane[[1]], dane[[2]], join="inner")
stopa_merged <- merge.xts(stopa[[1]], stopa[[2]], join="inner")

colnames(dane_merged) <- c('SP500', 'WIG')
colnames(stopa_merged) <- c('SP500', 'WIG')

sp500 = dane_merged$SP500
wig = dane_merged$WIG
stopa_sp500 = stopa_merged$SP500
stopa_sp500_1 = stopa_sp500[-length(stopa_sp500)]
stopa_wig = stopa_merged$WIG

cor(wig, sp500)
cor(stopa_wig, stopa_sp500)

Korelacja między WIG a S&P500 wynosi 0,72, natomiast ich stopy zwrotu korelują na poziomie 0,51. Regresję liniową przy pomocy funkcji lm [linear model] zbudujemy na stopach zwrotu.

stopa.model <- lm(stopa_wig ~ stopa_sp500)
summary(stopa.model)


Model jest oczywiście istotny (wpływ SP500 wynosi 58%), natomiast ciekawe, że wyraz wolny jest nieistotny - średnia niewarunkowa stopa wychodzi zerowa.


# Część główna
# cusum, mosum, ME, RE, GEFP

stopa.model <- stopa_wig ~ stopa_sp500

#Score based CUSUM
cus = efp(stopa.model, type="Score-CUSUM")
plot(cus, functional=NULL)


Wg CUSUM zmieniła się wariancja modelu.


#Score based MOSUM
mos = efp(stopa.model, type="Score-MOSUM")
plot(mos)
plot(mos, functional=NULL)

MOSUM nie wskazuje na zmiany.

#Moving estimate
me <- efp(stopa.model, type="ME")
plot(me, functional = NULL)



ME sugeruje, że korelacja z S&P500 zmieniła się na chwilę. Wykres agregujący pozwoli pokazać ten moment na WIG. W tym miejscu muszę przyznać, że zrobiłem "błąd", bo operowałem na pakiecie xts, sądząc, że będzie łatwiej przy konstrukcji danych i wykresach. Jednakże pakiet strucchange działa (obecnie) formalnie na obiektach ts i jeśli jest to xts, to nie widać dat na wykresach. Drugim problemem jest odpowiednie dopasowanie osi x do wykresu. Aby pokazać każdy rok, trzeba użyć dodatkowego kodu, który objaśniłem na stronie o parametrach graficznych w R.

par(mfrow=c(2,1), mar=c(1,5,2,3), oma=c(2,0,0,0))
plot(me, functional="max")
plot(x=daty, y=as.ts(as.zoo(wig)), type="l", ylab="WIG", xaxt="n")
datySkr = seq(from=daty[1], to=daty[length(daty)], length.out=round(length(daty)/52))
axis(side=1, cex.axis=1, cex.lab=1, at=datySkr, labels=format(datySkr,"%Y"), las=2)




Z wykresów wynika, że tygodniowa kowariancja z amerykańskim indeksem wzrosła w latach 2005-2006. Obecnie wpływ USA na WIG wynosi niecałe 60% zgodnie z modelem regresji  i nic nie wskazuje, by miało się to zmienić w najbliższym czasie.

Możemy też wnioskować, że słaba kondycja naszej giełdy na tle USA wynikała nie ze spadku korelacji, ale średniej (niewarunkowej) stopy zwrotu. To znaczy średnia tygodniowa pozostała stabilna - jest praktycznie zerowa - i znajduje się w obszarze normalności, ale "losowo" spadła poniżej zera (Intercept < 0). Dzieje się to od 2008 r. Ostatnie tygodnie to jej poprawa i jeżeli iść "logiką" regresji do średniej to powinniśmy oczekiwać wychylenia na plus.

P. S. Pozostałych dwóch testów: RE i GEFP nie pokazuję już, bo wyniki są takie same jak przy CUSUM.