Przyczynowość w sensie Grangera. Test Grangera w języku R

Normalnie gdy chcemy sprawdzić czy jedna zmienna wpływa na drugą, budujemy model regresji, w którym zmienna objaśniona jest opóźniona o jakiś okres, i testujemy go. Co jednak zrobić, gdy występuje autokorelacja, jak to ma miejsce często w danych ekonomicznych? Dodawanie do takiego  modelu autoregresji jest teoretycznie błędne, bo wiadomo, że w normalnej regresji zmienne objaśniające powinny być niezależne od siebie. 

Z pomocą przychodzi nam test Grangera. Test sprawdza przyczynowość w sensie Grangera. C. Granger zdefiniował przyczynowość w następujący sposób [1]:

Mówimy, że Xt powoduje Yt, jeśli jesteśmy w stanie lepiej przewidzieć Yt przy użyciu wszystkich dostępnych informacji, niż gdyby wykorzystano wszystkie informacje oprócz Xt.

Jest to definicja bardzo ogólna i abstrakcyjna. Zauważmy, że w tej definicji nie ma opóźnienia Y względem X. Dlaczego? Łatwiejsze zrozumienie przyczynowości Grangera w kontekście procesu ARMA podają autorzy w [2]:

X wywołuje Y, jeśli wcześniejsze informacje o X przewidują zachowanie Y lepiej niż same wcześniejsze informacje o Y.

Innymi słowy Xt i Yt są zmiennymi, które wpływają na siebie w teraźniejszości - natychmiast, ale Granger oddziela zbiór informacji o tych zmiennych, który może mieć przesunięcie w czasie lub ściślej mówiąc zbiór z okresu t zawiera informacje z okresu t-1. Taki podział pozwala rozróżnić zwykłą korelację od przyczynowości. Powiedzmy, że mamy rozkład jazdy: autobus X przejeżdża część trasy od godz. 8 do 9, a autobus Y kolejną część trasy od 9 do 10. Gdy X dojeżdża do swojego końca, Y rusza. Jednak X nie jest przyczyną Y, a jedynie z nim koreluje, ponieważ Y może ruszyć niezależnie od tego czy X dojeżdża czy nie. A teraz wyobraźmy sobie, że jazda obydwu ma się odbywać w tym samym czasie. Gdyby X bardzo przyspieszył i dojeżdżał już o 8:45, to Y musiałby ruszyć o 8:45, a nie 9:00, żeby nie zostać uderzonym przez X. W ten sposób X(t) staje się przyczyną Y(t). Ale ocenić to możemy dopiero, gdy mamy informację, że X wcześniej z jakiegoś powodu przyspieszył, np. nie było nikogo na drodze.

W świecie gospodarki sytuacja wygląda podobnie, choć jest dużo bardziej złożona. Produkcja w jednym kraju będzie wpływać na produkcję w innym kraju, bo np. wzrost produkcji w kraju X to wzrost dochodów, który można przeznaczyć na nowe inwestycje w kraju Y. 

Zajmijmy się teraz stroną praktyczną, czyli zastosujmy test Grangera w języku R. Na stronie OECD możemy ściągnąć statystyki produkcji miesięcznej dla krajów OECD. Wybieramy Time -> Select Date Range -> Monthly -> np. od 2001 r. Eksportujemy text file (CSV).  Interesują mnie ogólnie 4 obszary:

1. Polska

2. UE

3. USA

4. całe OECD

Niestety wskaźniki na OECD są słabo opisane i trzeba się trochę dokopać, żeby zrozumieć, co jest czym. Przykładowo produkcji nie można oddzielać od sprzedaży, bo zgodnie z zasadami rachunkowości tylko sprzedana produkcja się liczy do rachunków. Podstawową miarą jest indeks produkcji. Żeby odkryć co wchodzi w skład produkcji, możemy poniżej kliknąć w Structural Analysis (STAN) Databases. Pojawi się tablica ze składem industry. Możemy kliknąć w Industry, żeby obejrzeć cały skład. Okazuje się, że mieści się tu cała gospodarka, w tym usługi.

Jeśli chodzi o same dane, to są tu dwie istotne cechy. Po pierwsze jest to indeks w odniesieniu do określonego roku. Nie jest ważne do jakiego, ale należy zmienną przekształcić w stopę zmian (tutaj m/m). Po drugie nie są to dane surowe, ale poddane obróbce - usunięto z nich sezonowość za pomocą X12-ARIMA i TRAMO-SEATS , stąd oznaczenie sa (ang. seasonal adjustment). To drugie może stanowić podstawę różnicy między danymi OECD a GUS.

No dobrze, na początek mamy dwie zmienne:

- stopy zmian produkcji w Polsce, m/m: stopa_pl 

- stopy zmian produkcji w USA, m/m: stopa_us

Kod:

plot(stopa_pl, lwd=2, col="blue", main = "Industry, change% m/m (OECD)")

lines(stopa_us, lwd=2, col="red")

legend("topleft", legend=c("Poland", "USA"), col=c("blue", "red"), lty = 1, lwd = 2)

Wykonajmy korelację krzyżową między nimi:

ccf(stopa_pl, stopa_us)

Występuje silna autokorelacja bieżąca (0,6), ale też dość mocna korelacja z pospieszeniem zmiennej stopa_pl o 1 okres. Czyli to produkcja USA jest opóźniona o 1 okres w stos. do polskiej. Wydawałoby się, że powinno być na odwrót - że Polska podąża za USA. Tutaj dostajemy coś mało intuicyjnego. Może problemem są autokorelacje? Sprawdźmy PACF (autokorelacje cząstkowe):

pacf(stopa_pl)

 

Okazuje się, że autokorelacji 1 rzędu brak. Interesujące, że występuje ujemna korelacja co drugi miesiąc, ale na razie zostawiamy to. 

Zastosujmy test Grangera. Użyjemy do tego funkcji grangertest w pakiecie lmtest. Pakiet ten nie jest niczym nowym, bo używamy go zawsze do testów homoskedastyczności. Trzeba nadmienić, że test Grangera dotyczy zmiennych stacjonarnych, więc warto na początku wykonać test stacjonarności, ale to pomijam. 

Przetestujmy to co jest mało intuicyjne - czy Polska (x) wpływa na USA (y).

Hipoteza zerowa: Polska nie wpływa na USA.

Hipoteza alternatywna: Polska wpływa na USA.

Jak widzieliśmy w ccf występowało tylko jedno opóźnienie w korelacji, więc przypuszczamy, że w teście wystarczy wskazać order=1: 

grangertest(x=stopa_pl, y=stopa_us, order=1)

Granger causality test

Model 1: stopa_us ~ Lags(stopa_us, 1:1) + Lags(stopa_pl, 1:1)
Model 2: stopa_us ~ Lags(stopa_us, 1:1)
  Res.Df Df    F   Pr(>F)    
1    265                     
2    266 -1 17.5 0.000038 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Hipoteza zerowa zostaje odrzucona. Nasze wątpliwości potwierdzają się - zmiany produkcji w Polsce powodują w sensie Grangera zmiany w produkcji USA. Sprawdźmy w drugą stronę - czy produkcja USA (x) wpływa na produkcję Polski (y).

Hipoteza zerowa: USA nie wpływa po miesiącu na Polskę.

Hipoteza alternatywna: USA wpływa po miesiącu na Polskę.

grangertest(x=stopa_us, y=stopa_pl, order=1)

Efekt:

Granger causality test

Model 1: stopa_pl ~ Lags(stopa_pl, 1:1) + Lags(stopa_us, 1:1)
Model 2: stopa_pl ~ Lags(stopa_pl, 1:1)
  Res.Df Df    F Pr(>F)
1    265               
2    266 -1 0.08   0.78

Nie można odrzucić hipotezy zerowej. Okazuje się, że USA nie wpływa w sensie Grangera na Polskę. To wydaje się pozbawione sensu, ale pamiętajmy o silnej bieżącej korelacji między krajami. Możliwe, że opóźnienia występują wewnątrz miesiąca. 

Pozostaje zagadka, dlaczego Polska wyprzedza produkcję USA. Zwróćmy uwagę, że Polska jest częścią Unii Europejskiej, która z kolei może wpływać na gospodarkę USA. Tak więc Polska może mieć pośredni wpływ na USA. Sprawdźmy tę hipotezę.

Porównajmy ze sobą:

- stopy zmian produkcji w UE, m/m: stopa_eu

- stopy zmian produkcji w USA, m/m: stopa_us

Nie pokazuję już efektu krzyżowej korelacji, bo jest bardzo podobny do tej pierwszej - występuje głównie korelacja 1 rzędu.

Hipoteza zerowa: UE nie wpływa po miesiącu na USA.

Hipoteza alternatywna: UE wpływa po miesiącu na USA.

Wykonujemy test Grangera. Ponieważ obecnie nie ma jeszcze danych dla UE za czerwiec, odejmuję ostatnią wartość dla USA.  

grangertest(stopa_eu, stopa_us[-length(stopa_us)], order=1)

Granger causality test

Model 1: stopa_us[-length(stopa_us)] ~ Lags(stopa_us[-length(stopa_us)], 1:1) + Lags(stopa_eu, 1:1)
Model 2: stopa_us[-length(stopa_us)] ~ Lags(stopa_us[-length(stopa_us)], 1:1)
  Res.Df Df    F        Pr(>F)    
1    264                          
2    265 -1 43.5 0.00000000023 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Hipoteza zerowa zostaje odrzucona, zatem UE istotnie wpływa na USA. Sprawdźmy to samo na odwrót:

Hipoteza zerowa: USA nie wpływa po miesiącu na UE.

Hipoteza alternatywna: USA wpływa po miesiącu na UE.

grangertest(stopa_us[-length(stopa_us)], stopa_eu, order=1)

Granger causality test

Model 1: stopa_eu ~ Lags(stopa_eu, 1:1) + Lags(stopa_us[-length(stopa_us)], 1:1)
Model 2: stopa_eu ~ Lags(stopa_eu, 1:1)
  Res.Df Df    F Pr(>F)  
1    264                 
2    265 -1 3.48  0.063 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Hipoteza o tym, że USA wpływa na UE musi zostać odrzucona na poziomie 5% istotności, ale już na poziomie 10% nie. Można przypuszczać, że jakiś niewielki feedback, tzn. sprzężenie zwrotne między USA a UE także występuje. 

Zaproponowana hipoteza, że Polska wpływa pośrednio na USA poprzez UE jest już częściowo potwierdzona. Żeby w pełni to potwierdzić przetestujmy związek między UE a Polską.

Hipoteza zerowa: UE nie wpływa po miesiącu na Polskę.

Hipoteza alternatywna: UE wpływa po miesiącu na Polskę.

Podobnie jak z USA, odejmuję ostatnią wartość dla Polski i robię test:

grangertest(stopa_eu, stopa_pl[-length(stopa_pl)], order=1)

Granger causality test

Model 1: stopa_pl[-length(stopa_pl)] ~ Lags(stopa_pl[-length(stopa_pl)], 1:1) + Lags(stopa_eu, 1:1)
Model 2: stopa_pl[-length(stopa_pl)] ~ Lags(stopa_pl[-length(stopa_pl)], 1:1)
  Res.Df Df    F   Pr(>F)    
1    264                     
2    265 -1 16.5 0.000065 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Hipoteza zerowa - odrzucona. Produkcja UE wpływa na produkcję Polski. Odwrotnie:

Hipoteza zerowa: Polska nie wpływa po miesiącu na UE.

Hipoteza alternatywna: Polska wpływa po miesiącu na UE.

grangertest(stopa_pl[-length(stopa_pl)], stopa_eu, order=1)

Granger causality test

Model 1: stopa_eu ~ Lags(stopa_eu, 1:1) + Lags(stopa_pl[-length(stopa_pl)], 1:1)
Model 2: stopa_eu ~ Lags(stopa_eu, 1:1)
  Res.Df Df    F Pr(>F)
1    264               
2    265 -1 1.09    0.3

Hipoteza zerowa - zostaje. Produkcja Polski nie wpływa na produkcję UE z opóźnieniem miesięcznym. 

W sumie potwierdziliśmy hipotezę, że Polska wpływa na USA pośrednio przez UE. 

Zagadka rozwiązana, chociaż pozostaje jeszcze jedna ciekawa zależność, mianowicie między całym OECD a USA, UE i Polską. Podobnie jak w przypadku UE, dane dla całego OECD pojawiają się z opóźnieniem, stąd podobnie trzeba odjąć ostatnie dane USA i Polski. Ponieważ pojawi się tu sporo kombinacji, nie będę już opisywał każdej hipotezy, a podam sam kod, efekt wraz z wnioskami:

OECD -> USA:

grangertest(stopa_oecd, stopa_us[-length(stopa_us)], order=1)

OECD powoduje zmiany w gospodarce USA po miesiącu (duża istotność, p = 0).

USA -> OECD:

grangertest(stopa_us[-length(stopa_us)], stopa_oecd, order=1)

USA powoduje zmiany w gospodarce OECD po miesiącu na poziomie 5% istotności (p = 0.022).

OECD -> UE:

grangertest(x=stopa_oecd, y=stopa_eu, order=1)

OECD nie powoduje zmiany w gospodarce UE po miesiącu (p = 0.36).

UE -> OECD:

grangertest(x=stopa_eu, y=stopa_oecd, order=1)

UE powoduje zmiany w gospodarce OECD po miesiącu (p = 0.0027).

OECD -> Polska:

grangertest(stopa_oecd, stopa_pl[-length(stopa_pl)], order=1)

OECD nie powoduje zmiany w gospodarce Polski po miesiącu na poziomie istotności 5% (p = 0.1).

Polska -> OECD:

grangertest(stopa_pl[-length(stopa_pl)], stopa_oecd, order=1)

Polska nie powoduje zmiany w gospodarce OECD po miesiącu (p = 0.86).

Dostajemy interesujący wniosek. Zmiany w UE wpływają po miesiącu na OECD, w tym na Polskę. Polska nie wpływa na OECD, mimo że wpływa pośrednio na USA, a USA stanowi dużą część OECD. Jednocześnie istnieje sprzężenie zwrotne między USA a OECD. W końcu, UE indukuje zmiany w OECD, ale nie na odwrót. 

Wszystkie odkryte zależności w tym artykule ilustruje poniższy schemat:

Strzałka od USA do UE jest przerywana, bo jest to słaba zależność (sprzężenie zwrotne) oraz ma inny kształt, bo omija Polskę, a wskazuje na Unię jako całość.

 

Literatura:

[1] Granger C.W.J. Investigating causal relations by econometric models and cross spectral methods. Econometrica. 1969;37:424–438.

[2] Amornbunchornvej C., Zheleva E., Berger-Wolf T.Y. Variable-lag Granger Causality for Time Series Analysis. Machine Learning (cs.LG); Econometrics (econ.EM); Quantitative Methods, 2019.

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

Cały kod:

if (require(lmtest)==FALSE) {

  install.packages("lmtest")

  library("lmtest")

}

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

  install.packages("zoo")

  library("zoo")

}

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

  install.packages("tidyverse")

  library(tidyverse)

}

#zamieniam "\" na "/"

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

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

# albo

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

# ustawiamy folder roboczy

setwd(sciezka)

fr = 12

nazwaPliku = "MEI_REAL.csv"

if (grepl(";", readLines(nazwaPliku, n=1))) {

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

} else if (grepl(",", readLines(nazwaPliku, n=1))) {

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

}

produkcja = "Production of total industry sa, Index"

kraj1 = "Poland"

kraj2 = "United States"

kraj3 = "European Union – 27 countries (from 01/02/2020)"

kraj4 = "OECD - Total"

kraje = c(kraj1, kraj2, kraj3, kraj4)

moje_formaty = c("%d.%m.%Y", "%Y.%m.%d", "%d-%m-%Y", "%Y-%m-%d", "%d/%m/%Y", "%Y/%m/%d", "%b-%Y", "%B-%Y")

for (kraj in kraje) {

  filtered <- subset(plik, plik$Subject == produkcja & plik$Country == kraj)

  daty <- filtered[["Time"]]

  daty = parse_date_time(daty, orders=moje_formaty)

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

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

  wartosci <- filtered[["Value"]]

  stopy1m <- diff(wartosci)/wartosci[-length(wartosci)]

  if (kraj==kraje[1]) {

    stopa1 <- stats::lag(ts(stopy1m, start = c(rok[1], mc[1]), frequency= fr), -1)

  } else if (kraj==kraje[2]) {

    stopa2 <- stats::lag(ts(stopy1m, start = c(rok[1], mc[1]), frequency= fr), -1)

  } else if (kraj==kraje[3]) {

    stopa3 <- stats::lag(ts(stopy1m, start = c(rok[1], mc[1]), frequency= fr), -1)

  } else if (kraj==kraje[4]) {

    stopa4 <- stats::lag(ts(stopy1m, start = c(rok[1], mc[1]), frequency= fr), -1)

  }

}

stopa_pl = as.numeric(stopa1)

stopa_us = as.numeric(stopa2)

stopa_eu = as.numeric(stopa3)

stopa_oecd = as.numeric(stopa4)

# test grangera

grangertest(x=stopa_pl, y=stopa_us, order=1)

grangertest(x=stopa_us, y=stopa_pl, order=1)

grangertest(x=stopa_us[-length(stopa_us)], y=stopa_eu, order=1)

grangertest(x=stopa_eu, y=stopa_us[-length(stopa_us)], order=1)

grangertest(x=stopa_eu, y=stopa_pl[-length(stopa_pl)], order=1)

grangertest(x=stopa_pl[-length(stopa_pl)], y=stopa_eu, order=1)

grangertest(x=stopa_oecd, y=stopa_us[-length(stopa_us)], order=1)

grangertest(x=stopa_us[-length(stopa_us)], y=stopa_oecd, order=1)

grangertest(x=stopa_oecd, y=stopa_eu, order=1)

grangertest(x=stopa_eu, y=stopa_oecd, order=1)

grangertest(x=stopa_oecd, y=stopa_pl[-length(stopa_pl)], order=1)

grangertest(x=stopa_pl[-length(stopa_pl)], y=stopa_oecd, order=1)

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.

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.