Jakiś czas temu pisałem o zestawie dokumentów przygotowanych przez Krzyśka Trajkowskiego. Zastanawiałem się później, w jaki sposób na dostępność materiałów wpływa format np. pdf czy html. Można się zastanawiać, można sprawdzić. Gdy Krzysiek przygotował kolejny dokument skonwertowałem go do html’a i … dodałem do listy gościnnych wpisów.
Konwersja została wykonana z użyciem konwertera pandoc i była 'prawie’ bezbolesna, tzn te same źródła, które były użyte do kompilacji do pdf’a, po drobnych zmianach posłużyły do kompilacji do html’a. Jedna z niewielu zmian to zastąpienie środowiska lstlisting tagiem <pre lang=”rsplus”>
Tak więc, dziś będzie gościnny wpis Krzysztofa Trajkowskiego, dotyczący wybranych narzędzi dostępnych w R pozwalających na analizę pewnych szeregów czasowych w finansach.
Zapraszamy!
Wersja pdf znajduje się tutaj, wersja HTML poniżej.
Wybrane pakiety do analizy finansowych szeregów czasowych w R
Krzysztof Trajkowski
Pozyskiwanie danych
Do pozyskania aktualnych danych giełdowych (obiekt OHLC) wykorzystamy funkcję getSymbols{quantmod}. Poniżej przykład wczytania danych giełdowych spółki Adobe Systems Inc. z datą końcową 2013-03-31.
> library(quantmod) > getSymbols("ADBE",src="yahoo",to="2013-03-31") [1] "ADBE" > head(ADBE)[,-6] ADBE.Open ADBE.High ADBE.Low ADBE.Close ADBE.Volume 2007-01-03 40.72 41.32 38.89 39.92 7126000 2007-01-04 39.88 41.00 39.43 40.82 4503700 2007-01-05 40.78 40.90 40.12 40.62 2730200 2007-01-08 40.41 40.97 40.13 40.45 5234800 2007-01-09 40.50 40.51 39.38 39.63 5672900 2007-01-10 39.15 39.59 38.73 39.22 5652600 |
Aby korzystać z szeregu funkcji do przetwarzania danych dostępnych w pakiecie quantmod (np. pozyskiwanie informacji o szeregu, agregacja danych) zalecane jest aby szereg czasowy był obiektem klasy xts{xts} lub zoo{zoo}.
> class(ADBE) [1] "xts" "zoo" |
Poniżej kilka przykładów funkcji do pozyskiwania informacji o szeregu ADBE.
# data początkowa: > start(ADBE) [1] "2007-01-03" # data końcowa: > end(ADBE) [1] "2013-03-28" |
# ilość lat: > nyears(ADBE) [1] 7 # ilość kwartałów: > nquarters(ADBE) [1] 25 # ilość miesięcy: > nmonths(ADBE) [1] 75 # ilość tygodni: > nweeks(ADBE) [1] 326 # ilość dni: > ndays(ADBE) [1] 1570 |
Możemy również wyodrębnić dane z interesującego nas okresu.
# pierwsze 5 dni: > first(ADBE, "5 days")[,-6] ADBE.Open ADBE.High ADBE.Low ADBE.Close ADBE.Volume 2007-01-03 40.72 41.32 38.89 39.92 7126000 2007-01-04 39.88 41.00 39.43 40.82 4503700 2007-01-05 40.78 40.90 40.12 40.62 2730200 2007-01-08 40.41 40.97 40.13 40.45 5234800 2007-01-09 40.50 40.51 39.38 39.63 5672900 # ostatnie 5 dni: > last(ADBE, "5 days")[,-6] ADBE.Open ADBE.High ADBE.Low ADBE.Close ADBE.Volume 2013-03-22 42.17 43.23 42.11 42.97 5595200 2013-03-25 42.90 43.42 42.41 42.49 5296600 2013-03-26 42.66 42.95 42.46 42.75 2125400 2013-03-27 42.60 42.88 42.38 42.66 1707000 2013-03-28 42.54 43.59 42.42 43.52 4075200 |
Zamiast opcji days możemy stosować także inne warianty np. years, quarters, months i weeks.
# wyodrębnienie danych dla dowolnie wybranego okresu: > adbe.o= ADBE["2010::2010-01-08"] > adbe.o[,-6] ADBE.Open ADBE.High ADBE.Low ADBE.Close ADBE.Volume 2010-01-04 36.65 37.30 36.65 37.09 4710200 2010-01-05 37.04 37.80 36.87 37.70 7108800 2010-01-06 37.33 37.74 37.20 37.62 5336400 2010-01-07 37.41 37.59 36.81 36.89 5576700 2010-01-08 36.75 36.97 36.34 36.69 5429200 |
# agregacja danych dziennych do tygodniowych: > adbe.weekly= to.weekly(ADBE) > adbe.weekly[1:5,-6] ADBE.Open ADBE.High ADBE.Low ADBE.Close ADBE.Volume 2007-01-05 40.72 41.32 38.89 40.62 14359900 2007-01-12 40.41 40.97 38.73 39.96 24601800 2007-01-19 40.25 40.75 37.50 38.40 22010400 2007-01-26 38.44 39.97 37.20 39.26 28408600 2007-02-02 39.40 39.75 38.28 38.97 21605100 |
Pozostałe dostępne opcje to: to.yearly, to.quarterly, to.monthly i to.daily.
Na podstawie wcześniej wczytanych danych tzn. ADBE dotyczących spółki Adobe obliczymy logarytmiczne stopy zwrotu z wykorzystaniem cen zamknięcia.
# dzienne stopy zwrotu: > r.daily= periodReturn(ADBE,type="log",period="daily") # tygodniowe stopy zwrotu: > r.weekly= periodReturn(ADBE,type="log",period="weekly") # miesięczne stopy zwrotu: > r.monthly= periodReturn(ADBE,type="log",period="monthly") # kwartalne stopy zwrotu: > r.quartely= periodReturn(ADBE,type="log",period="quarterly") # roczne stopy zwrotu: > r.yearly= periodReturn(ADBE,type="log",period="yearly") # wszystkie razem: > r.s= allReturns(ADBE,type="log") |
Teraz wczytamy kilka obiektów OHLC dla głównych indeksów giełdowych: niemiecki DAX, francuski CAC40, brytyjski FTSE, japoński NIKKEY oraz amerykański NASDAQ.
> indeks= c("^GDAXI","^FCHI","^FTSE","^N225","^IXIC") > getSymbols(indeks, src="yahoo", from="2012-01-01",to="2012-12-31") [1] "GDAXI" "FCHI" "FTSE" "N225" "IXIC" |
Oprócz danych dotyczących cen akcji lub indeksów giełdowych możemy także wczytywać dane dotyczące np. kursów walutowych, cen metali kolorowych i surowców oraz dane ekonomiczne różnych krajów.
Osoby zajmujące się analizą techniczną z pewnością zinteresuje możliwość wykonania np. wykresów świecowych. Warto podkreślić, że dla dużego zakresu danych taki wykres upodabnia się do wykresu liniowego.
# wykres liniowy: > chartSeries(AAPL,theme= chartTheme("white"),TA=NULL,type="line") |
[TAF1]
# wykres świecowy: > chartSeries(ADBE,theme= chartTheme("white"),TA=NULL) |
[TAF2]
Świece są bardziej widoczne gdy wykonamy wykres dla krótszego okresu np. dla stycznia 2013 roku.
# wykres świecowy: > chartSeries(ADBE,theme=chartTheme("white"), TA=NULL, subset="2013-01") |
[TAF3]
Jak można łatwo zauważyć wykres świecowy dostarcza nam o wiele więcej informacji niż wykres liniowy. Jest to szczególnie widoczne gdy obserwujemy dane o małym zakresie czasowym (rys. [TAF3]). Każda świeca zawiera informacje zawarte w szeregu czasowym typu OHLC tj. o cenie otwarcia, maksymalnej, minimalnej i zamknięcia. Z kolei jej kolor jest uzależniony od tego czy cena otwarcia była większa od ceny zamknięcia – kolor pomarańczowy czy też nie – kolor zielony.
Gdy w kodzie pominiemy opcję TA=NULL dodatkowo pojawi się wykres słupkowy obrazujący wielkość obrotów w danym okresie czasu.
# wykres świecowy + wykres wielkości obrotów: > chartSeries(AAPL,theme = chartTheme("white"),subset="2013-03") |
[TAF4]
Zwróćmy także uwagę na pojawiające się wartości liczbowe w górnym lewym rogu każdego okienka wykresu (rys. [TAF4]). Na wykresie świecowym oraz słupkowym jest podana informacja odpowiednio o cenie zamknięcia i wielkości obrotów w ostatnim okresie tj. dla 28-03-2013.
> last(ADBE)[,-6] ADBE.Open ADBE.High ADBE.Low ADBE.Close ADBE.Volume 2013-03-28 42.54 43.59 42.42 43.52 4075200 |
Ciekawą opcją jest również możliwość dodawania wszystkich funkcji dostępnych w pakiecie TTR. Na rysunku [TAF5] przykład wykresu świecowego oraz wstęga Bollingera. Natomiast w dolnej części wykresy MACAD i słupkowy wielkości obrotów.
> chartSeries(ADBE,theme = chartTheme("white"),subset="2012", TA=c(addBBands(),addMACD(),addVo())) |
[TAF5]
Budowa szeregu czasowego
W tym podrozdziale pokażemy jak narysować wykres szeregu czasowego – pakiet ggplot2 oraz jak zbudować szereg czasowy – pakiet xts. Zastosowanie pakietu xts do budowy szeregu czasowego umożliwi modyfikację danych za pomocą wielu funkcji dostępnych w bibliotece quantmod. Dane które będziemy wykorzystywać dotyczą cen akcji spółki Telekomunikacja Polska s.a. i zostały pozyskane ze strony www.stooq.pl w formie pliku .csv.
# wczytujemy plik tps_d.csv: > t= read.csv("/home/.../tps_d.csv",h=T) > head(t) Data Otwarcie Najwyzszy Najnizszy Zamkniecie Wolumen 1 1998-11-18 9.8466 10.0806 9.6193 9.6193 9029239 2 1998-11-19 9.5607 9.5607 9.0927 9.3853 5982341 3 1998-11-20 9.5607 9.7882 9.5607 9.7882 4538351 4 1998-11-23 9.9052 10.1392 9.9052 9.9638 4411099 5 1998-11-24 9.9638 9.9638 9.5607 9.6193 2754558 6 1998-11-25 9.6193 9.7297 9.5022 9.6712 1472786 # wykres cen zamknięcia: > library(ggplot2) > t$Data <- as.Date( t$Data, "%Y-%m-%d") > ggplot( data= t, aes(x=Data, y=Zamkniecie,colour=Zamkniecie))+ geom_line()+ scale_colour_gradient(low = "red")+ labs(colour = "cena")+ ylab("cena zamknięcia")+xlab("data") |
[cana1]
Aby uzyskać logarytmiczne stopy zwrotu wykorzystamy pakiet xts oraz quantmod.
> library(xts) # szereg czasowy w formaie xts: > txts= xts(t[,-1], order.by=t[,1]) # kilka pierwszych wierszy: > head(txts) Otwarcie Najwyzszy Najnizszy Zamkniecie Wolumen 1998-11-18 9.8466 10.0806 9.6193 9.6193 9029239 1998-11-19 9.5607 9.5607 9.0927 9.3853 5982341 1998-11-20 9.5607 9.7882 9.5607 9.7882 4538351 1998-11-23 9.9052 10.1392 9.9052 9.9638 4411099 1998-11-24 9.9638 9.9638 9.5607 9.6193 2754558 1998-11-25 9.6193 9.7297 9.5022 9.6712 1472786 # szereg czasowy w formaie xts - tylko ceny zamkniecia: > txts1= txts[,4] # szereg czasowy - logarytmiczna stopy zwrotu: > library(quantmod) > rtp1= periodReturn(txts1, type="log", period= "daily") > rtp1[1:3,] daily.returns 1998-11-18 0.00000000 1998-11-19 -0.02462686 1998-11-20 0.04203294 # wykresy: > data1= data.frame(czas=index(rtp1), rtp1) > library(scales) > plot1= ggplot(data=data1, aes(x=czas, y=daily.returns))+ geom_line(aes(colour = as.numeric(rtp1)),size=0.5)+ scale_colour_gradient(low = "red3",labels = percent)+ scale_y_continuous(labels = percent)+ ylab("stopa zwrotu")+xlab("czas")+labs(colour = "zwrot") > plot2= ggplot(data=data1, aes(x=czas, y=daily.returns))+ geom_point(aes(colour = as.numeric(rtp1)),size=1.5)+ scale_colour_gradient(low = "red3",labels = percent)+ scale_y_continuous(labels = percent)+ ylab("stopa zwrotu")+xlab("czas")+labs(colour = "zwrot") > library(gridExtra) > grid.arrange(plot2,plot1,nrow=2) |
[return1]
Właściwości finansowych szeregów czasowych
Rozkład prawdopodobieństwa
Każdy rozkład stopy zwrotu (np. z akcji) charakteryzuje się pewnymi cechami które mają wpływ na kształt rozkładu prawdopodobieństwa zmiennej. Do tych cech można zaliczyć:
-
skośność — miara asymetrii gęstości prawdopodobieństwa w stosunku do rozkładu normalnego,
-
leptokurtoza — wyższy szczyt gęstości prawdopodobieństwa w porównaniu z rozkładem normalnym,
-
grube ogony — gęstość prawdopodobieństwa na końcach rozkładu jest większa niż w rozkładzie normalnym.
Sprawdźmy zatem czy rzeczywiście zmienna rtp1 nie ma rozkładu normalnego.
# dzienne stopy zwrotu: > z1=as.numeric(rtp1) > library(fBasics) > dagoTest(z1) Title: D'Agostino Normality Test Test Results: STATISTIC: Chi2 | Omnibus: 1110.8943 Z3 | Skewness: -17.923 Z4 | Kurtosis: 28.1009 P VALUE: Omnibus Test: < 2.2e-16 Skewness Test: < 2.2e-16 Kurtosis Test: < 2.2e-16 |
Na podstawie powyższych obliczeń dochodzimy do wniosku, że rozkład stóp zwrotu nie ma rozkładu normalnego. Wniosek ten potwierdzają także dwa wykresy: kwantylowy oraz gęstości.
# parametry linii: > y= quantile(z1, c(0.25, 0.75)) > x= qnorm(c(0.25, 0.75),mean(z1),sd(z1)) > slope= diff(y)/diff(x) > int= y[1L] - slope * x[1L] # wykres kwantylowy: > QQ1= ggplot(NULL, aes(sample= z1))+ stat_qq(distribution= qnorm, dparams= list(mean=mean(z1),sd=sd(z1)), colour="red3")+ geom_abline(slope= slope, intercept= int,size=1) # punkty funkcji gęstości: > D1= with(density(z1),data.frame(x,y)) # wykres gęstości: > DE1=ggplot(data= D1, aes(x= x, y= y)) + geom_line(lwd=0.2,col="red3")+ geom_area(aes(fill="red"), alpha=0.5)+ stat_function(fun= dnorm, arg= list(mean= mean(z1),sd= sd(z1)), colour= "black")+ theme(legend.position = "none") # rozmieszczenie wykresów: > library(grid) > pushViewport(viewport(layout= grid.layout(1, 2))) > print(QQ1, vp= viewport(layout.pos.row= 1, layout.pos.col= 1)) > print(DE1, vp= viewport(layout.pos.row= 1, layout.pos.col= 2)) |
[QD1]
Połączymy teraz nasze dane w większe jednostki czasowe (agregacja szeregów czasowych) i sprawdzimy jak to przekształcenie szeregu wpłynie na rozkład prawdopodobieństwa. W pierwszej kolejności dokonamy agregacji danych dziennych do tygodniowych, następnie do miesięcznych i kwartalnych.
# tygodniowe stopy zwrotu: > rtp2= periodReturn(txts1, type="log", period= "weekly") > z2= as.numeric(rtp2) > dagoTest(z2) Title: D'Agostino Normality Test Test Results: STATISTIC: Chi2 | Omnibus: 269.7419 Z3 | Skewness: -10.7334 Z4 | Kurtosis: 12.4313 P VALUE: Omnibus Test: < 2.2e-16 Skewness Test: < 2.2e-16 Kurtosis Test: < 2.2e-16 |
# miesięczne stopy zwrotu: > rtp3= periodReturn(txts1, type="log", period= "monthly") > z3= as.numeric(rtp3) > dagoTest(z3) Title: D'Agostino Normality Test Test Results: STATISTIC: Chi2 | Omnibus: 42.6827 Z3 | Skewness: -3.1011 Z4 | Kurtosis: 5.7503 P VALUE: Omnibus Test: 5.39e-10 Skewness Test: 0.001928 Kurtosis Test: 8.909e-09 |
# kwartalne stopy zwrotu: > rtp4= periodReturn(txts1, type="log", period= "quarterly") > z4= as.numeric(rtp4) > dagoTest(z4) Title: D'Agostino Normality Test Test Results: STATISTIC: Chi2 | Omnibus: 9.5917 Z3 | Skewness: -2.3388 Z4 | Kurtosis: 2.0302 P VALUE: Omnibus Test: 0.008264 Skewness Test: 0.01935 Kurtosis Test: 0.04233 |
Na podstawie powyższych obliczeń możemy wysnuć wniosek, że im większa agregacja zmiennej (stopa zwrotu) tym bardziej rozkład prawdopodobieństwa upodabnia się do rozkładu normalnego.
W dalszej części tego opracowania będziemy analizować dzienne stopy zwrotu z akcji spółki Telekomunikacja Polska s.a. Zatem interesująca może się okazać odpowiedź na pytanie: jakiemu rozkładowi prawdopodobieństwa może podlegać zmienna rtp1.
# szacowanie parametrów rozkładu skośnego t-Studenta: > library(fGarch) > fit1= sstdFit(z1)$estimate # test zgodności z rozkładem skośnym t-Studenta: > library(ADGofTest) > ad.test(z1,psstd,fit1[1],fit1[2],fit1[3],fit1[4]) Anderson-Darling GoF Test data: z1 and psstd AD = 1.5565, p-value = 0.1635 alternative hypothesis: NA |
Na podstawie testu Andersona-Darlinga możemy stwierdzić, że badana stopa zwrotu ma rozkład skośny t-Studenta. Również wykresy (rys. [QD2]) potwierdzają ten wniosek.
# wykres kwantylowy: > QQ2= ggplot(NULL, aes(sample= z1))+ stat_qq(dist= qsstd, dparam= fit1, col= "red3",alpha=0.5)+ geom_abline(intercept= 0, slope= 1) # wykres gęstości: > DE2= ggplot(data= D1, aes(x= x, y= y)) + geom_line(lwd=0.2,colour="red")+ geom_area(aes(fill="red"), alpha=0.5)+ stat_function(fun= dsstd, arg= list(fit1[1], fit1[2], fit1[3], fit1[4]), colour= "black")+ theme(legend.position = "none") # rozmieszczenie wykresów: > pushViewport(viewport(layout= grid.layout(1, 2))) > print(QQ2, vp= viewport(layout.pos.row= 1, layout.pos.col= 1)) > print(DE2, vp= viewport(layout.pos.row= 1, layout.pos.col= 2)) |
[QD2]
Niejednorodność wariancji
Efekt skupiania zmienności (volatility clustering) to kolejna własność stóp zwrotu. Do oceny czy w badanym szeregu występuje zjawisko grupowania wariancji wykorzystamy test Engle’a.
> library(FinTS) > ArchTest(rtp1,lag=1) ARCH LM-test; Null hypothesis: no ARCH effects data: rtp1 Chi-squared = 17.4541, df = 1, p-value = 2.943e-05 |
Zatem na poziomie istotności α = 0, 05 odrzucamy hipotezę zerową która zakłada jednorodność wariancji. Potwierdzeniem efektu ARCH może być również autokorelacja kwadratów stóp zwrotu.
# przygotowanie danych: > lag= 1:length(acf(rtp1^2,plot=F)$lag[-1]) > acf= acf(rtp1^2,plot=F)$acf[-1] > C=1.96/sqrt(length(rtp1^2)) # wykres: > ggplot(NULL,aes(x=lag, y=acf)) + geom_bar(stat="identity", position= "identity")+ geom_hline(yintercept= C, color= "red3", size= 0.5)+ geom_hline(yintercept= -C, color= "red3", size= 0.5) |
[AUT2]
Z problemem grupowania wariancji wiąże się kolejna własność szeregów finansowych tzw. efekt dźwigni (leverage effects). Na rysunku [CS1] można zaobserwować silniejszą reakcję zmienności na duże spadki cen, niż na wzrosty.
> CE1= ggplot(data= t, aes(x=Data, y=Zamkniecie,colour=Zamkniecie))+ geom_line()+ scale_colour_gradient(low= "red")+ labs(colour= "cena")+ ylab("cena zamknięcia")+xlab("data") > ST1= ggplot(data=data1, aes(x=czas, y=daily.returns))+ geom_line(aes(colour = as.numeric(rtp1)),size=0.5)+ scale_colour_gradient(low = "red3",labels = percent)+ scale_y_continuous(labels = percent)+ ylab("stopa zwrotu")+xlab("data")+labs(colour = "zwrot") > pushViewport(viewport(layout= grid.layout(2, 1))) > print(CE1, vp= viewport(layout.pos.row= 1, layout.pos.col= 1)) > print(ST1, vp= viewport(layout.pos.row= 2, layout.pos.col= 1)) |
[CS1]
Efekt długiej pamięci
Na podstawie wartości współczynnika d (ułamkowy pierwiastek jednostkowy) lub H (wykładnik Hursta) możemy ocenić czy w badanym szeregu czasowym występuje efekt długiej pamięci (long memory processes). Warto podkreślić, że między omawianymi parametrami występuje pewien związek który można opisać równaniem:
d=H-0,5
Poniżej przedstawimy własności szeregów czasowych w zależności od wartości parametru d oraz H:
-
szereg antypersystentny (krótka pamięć):
d ∈ ( − 0, 5; 0) lub H ∈ (0; 0, 5) -
szereg losowy (biały szum):
d = 0 lub H = 0, 5 -
szereg persystentny (długa pamięć):
d ∈ (0; 0, 5) lub H ∈ (0, 5; 1)
Estymacja współczynnika integracji ułamkowej d.
> library(FGN) # estymator Whittle: > GetFitFD(rtp1,algorithm="wmle",ciQ=T) $d [1] -0.00618375 $loglikelihood [1] -1037.969 $alpha [1] 1.012367 $algorithm [1] "wmle" $ci [1] -0.03194746 0.02047613 |
Estymacja wykładnika Hursta H.
> library(FGN) # estymator Whittle: > GetFitFGN(rtp1,algorithm="wmle",ciQ=T) $H [1] 0.4975533 $loglikelihood [1] 4544.512 $alpha [1] 1.004893 $algorithm [1] "wmle" $ci [1] 0.4826359 0.5123336 |
Ponieważ oszacowana wartość niecałkowitego pierwiastka jednostkowego jest równa d̂ = − 0. 00618375 oraz wartość wykładnia Hursta wynosi Ĥ = 0. 4975533 to możemy stwierdzić, że w badanej stopie zwrotu nie występuje efekt długiej pamięci. Warto podkreślić, że w sytuacji odwrotnej tzn. gdy w badanym szeregu wystąpi efekt długiej pamięci można zastosować model z warunkową wartością oczekiwaną ARFIMA z wykorzystaniem np. funkcji fracdiff{FGN}. Dodajmy jeszcze, że w przypadku krótkiej pamięci funkcja autokorelacji ACF opada szybko tzn. w tempie wykładniczym, natomiast gdy występuje długa pamięć opadanie funkcji ACF jest wolniejsze tzn. hiperboliczne.
Funkcja autokorelacji stóp zwrotu.
> lag= 1:length(acf(rtp1,plot=F)$lag[-1]) > acf= acf(rtp1,plot=F)$acf[-1] > C=1.96/sqrt(length(rtp1)) > ggplot(NULL,aes(x=lag, y=acf)) + geom_bar(stat="identity", position= "identity")+ geom_hline(yintercept= C, color= "red3", size= 0.5)+ geom_hline(yintercept= -C, color= "red3", size= 0.5) |
[AUT1]
Jak się okazuje (rys. [AUT1]) szereg finansowy jest ciągiem nieskorelowanych zmiennych losowych. Zjawisko autokorelacji w badanym szeregu nie występuje.
Należy także wspomnieć o bardzo bogatej kolekcji metod estymacji współczynników d oraz H dostępnych w pakiecie fArma.
Modele ARCH/GARCH
Do estymacji standardowego modelu GARCH możemy użyć funkcji garch{tseries}. Dla modelu ARCH(1) [Engle, 1982] składnia jest następująca:
> library(tseries) # Engle [1982]: > n= garch(rtp1,order=c(0,1), trace= F) |
Natomiast dla modelu GARCH(1,1) [Bollerslev, 1986]:
# Bollerslev [1986]: > n=garch(rtp1, order=c(1,1) ,trace= F) > summary(n) Call: garch(x = rtp1, order = c(1, 1)) Model: GARCH(1,1) Residuals: Min 1Q Median 3Q Max -18.0356 -0.5605 0.0000 0.5477 4.8106 Coefficient(s): Estimate Std. Error t value Pr(>|t|) a0 4.274e-05 5.640e-06 7.579 3.49e-14 *** a1 6.558e-02 9.059e-03 7.240 4.50e-13 *** b1 8.485e-01 1.934e-02 43.862 < 2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Diagnostic Tests: Jarque Bera Test data: Residuals X-squared = 153390, df = 2, p-value < 2.2e-16 Box-Ljung test data: Squared.Residuals X-squared = 0.0024, df = 1, p-value = 0.9607 |
# logarytm wiarygodności: > logLik(n) 'log Lik.' 8716.982 (df=3) |
Zwróćmy uwagę na to, że funkcja garch{tseries} zawsze zakłada rozkład normalny reszt. Jak się okazało to założenie zostało naruszone – test Jarque-Bera. W związku z tym należy oszacować model, który uwzględni rozkład inny niż normalny np. rozkład skośny t-Studenta.
Modele APARCH
Gdy w szeregu czasowym występuje asymetria lub efekt dźwigni dobrym rozwiązaniem jest zastosowanie modelu GARCH lub APARCH z odpowiednim rozkładem innowacji. Modele tej klasy możemy estymować za pomocą funkcji garchFit{fGarch}. Opcja cond.dist określa jaki rozkład prawdopodobieństwa będzie wykorzystany w modelu: std, sstd, ged, sged, snig, QMLE, snorm. Rozkład norm czyli rozkład normalny, jest ustawiony jako domyślny. Jest też możliwość uwzględnienia modelu ARMA np. ARMA(2,1)-APARCH(1,1). W takim przypadku formuła określająca model będzie miała postać ~arma(2,1)+aparch(1,1).
-
APARCH – Asymmetric Power Autoregressive Conditional Heteroscedasticity:
# model APARCH [Ding, Granger, Engle 1993]: > ap=garchFit(data=rtp1, ~aparch(1,1),trace=F,cond.dist="sstd", include.delta=T, leverage=T) |
Dodajmy jeszcze, że include.delta=T oraz leverage=T są ustawieniami domyślnymi. Oznacza to, że parametry δ oraz γ w modelu APARCH będą oszacowane. W przeciwnym wypadku tzn. dla include.delta=F oraz leverage=F parametry δ oraz γ nie będą szacowane. Wtedy domyślnymi wartościami będą delta=2 i gamma=0. Oczywiście tym parametrom można przypisywać dowolne wartości. Wynika z tego, że w modelu APARCH jest zagnieżdżonych kilka innych modeli z rodziny GARCH. Poniżej wymienimy niektóre z nich:
-
ARCH — Autoregressive Conditional Heteroscedasticity:
# model ARCH [Engle 1982]: > fit1=garchFit(data=rtp1, ~aparch(1,0),trace=F,cond.dist="sstd", include.delta=F, delta=2, leverage=F, gamma=0)
-
GARCH — Generalized Autoregressive Conditional Heteroscedasticity:
# model GARCH [Bollerslev 1986]: > fit2=garchFit(data=rtp1, ~aparch(1,1),trace=F,cond.dist="sstd", include.delta=F, delta=2, leverage=F, gamma=0)
-
TS-GARCH
# model TS-GARCH [Taylor 1986, Schwert 1989]: > fit3=garchFit(data=rtp1, ~aparch(1,1),trace=F,cond.dist="sstd", include.delta=F, delta=1, leverage=F, gamma=0)
-
GJR-GARCH
# model GJR-ARCH [Glosten, Jaganathan, Runkle 1993]: > fit4=garchFit(data=rtp1, ~aparch(1,1),trace=F,cond.dist="sstd", include.delta=F, delta=2, leverage=T)
-
T-GARCH — Threshold GARCH:
# model T-ARCH [Zakoian 1993]: > fit5=garchFit(data=rtp1, ~aparch(1,1),trace=F,cond.dist="sstd", include.delta=F, delta=1, leverage=T)
-
N-ARCH — Nonlinear ARCH:
# model N-ARCH [Higgins, Bera 1992]: > fit6=garchFit(data=rtp1, ~aparch(1,1),trace=F,cond.dist="sstd", include.delta=T, leverage=F, gamma=0)
Poniżej sprawdzimy czy model sstd-N-ARCH(1,1) będzie dobrze dopasowany do badanej zmienej rtp1. Zauważmy, że include.delta=T oraz gamma=0 (dla opcji leverage=F) są wartościami domyślnymi i nie musimy ich uwzględniać w kodzie.
# model N-ARCH [Higgins, Bera 1992]: > fit6=garchFit(data=rtp1, ~aparch(1,1),trace=F,cond.dist= "sstd", include.delta=T, leverage=F, gamma=0) > fit6 Conditional Distribution: sstd Coefficient(s): mu omega alpha1 beta1 delta 0.00031365 0.00017166 0.07911262 0.92780326 1.12274413 skew shape 1.05850381 5.86789967 Std. Errors: based on Hessian Error Analysis: Estimate Std. Error t value Pr(>|t|) mu 3.137e-04 3.003e-04 1.045 0.29619 omega 1.717e-04 6.355e-05 2.701 0.00691 ** alpha1 7.911e-02 1.244e-02 6.360 2.02e-10 *** beta1 9.278e-01 1.259e-02 73.711 < 2e-16 *** delta 1.123e+00 1.372e-01 8.182 2.22e-16 *** skew 1.059e+00 2.381e-02 44.450 < 2e-16 *** shape 5.868e+00 5.436e-01 10.794 < 2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Log Likelihood: 9006.067 normalized: 2.499602 |
Dokonajmy jeszcze estymacji drugiego modelu tzn. sstd-GARCH(1,1), ale z wykorzystaniem formuły ~garch zamiast ~aparch:
# model GARCH [Bollerslev 1986]: > fit2=garchFit(data=rtp1, ~garch(1,1),trace=F,cond.dist= "sstd") > fit2 Conditional Distribution: sstd Coefficient(s): mu omega alpha1 beta1 skew 2.8110e-04 8.1436e-06 8.8715e-02 8.9968e-01 1.0577e+00 shape 5.8772e+00 Std. Errors: based on Hessian Error Analysis: Estimate Std. Error t value Pr(>|t|) mu 2.811e-04 3.016e-04 0.932 0.35127 omega 8.144e-06 2.883e-06 2.824 0.00474 ** alpha1 8.871e-02 1.565e-02 5.668 1.45e-08 *** beta1 8.997e-01 1.783e-02 50.457 < 2e-16 *** skew 1.058e+00 2.361e-02 44.796 < 2e-16 *** shape 5.877e+00 5.499e-01 10.689 < 2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Log Likelihood: 9013.413 normalized: 2.501641 |
Ponieważ logarytm wiarygodności jest większy niż w modelu sstd-N-ARCH to do dalszych analiz wykorzystamy model sstd-GARCH. Poniżej kilka podstawowych testów diagnostycznych dla modelu sstd-GARCH(1,1).
> summary(fit2) ..... Standardised Residuals Tests: Statistic p-Value Jarque-Bera Test R Chi^2 930873.6 0 Shapiro-Wilk Test R W 0.8598556 0 Ljung-Box Test R Q(10) 6.682323 0.7550566 Ljung-Box Test R Q(15) 13.99558 0.5258644 Ljung-Box Test R Q(20) 15.81096 0.7282805 Ljung-Box Test R^2 Q(10) 0.4081282 0.9999975 Ljung-Box Test R^2 Q(15) 0.4747168 1 Ljung-Box Test R^2 Q(20) 0.7100467 1 LM Arch Test R TR^2 0.4319661 0.9999999 Information Criterion Statistics: AIC BIC SIC HQIC -4.999951 -4.989644 -4.999957 -4.996278 |
W omawianym modelu fit2 założyliśmy rozkład skośny t-Studenta (na powyższym wydruku jest badana zgodność reszt tylko z rozkładem normalnym) więc sprawdzimy to założenie przy pomocy testu Andersona-Darlinga.
> R= residuals(fit2, standardize= T) > fit= sstdFit(R)$estimate > library(ADGofTest) > ad.test(R, psstd, fit[1], fit[2], fit[3], fit[4]) Anderson-Darling GoF Test data: R and psstd AD = 1.0609, p-value = 0.3269 alternative hypothesis: NA |
Dodatkowo mamy do wyboru kilkanaście wykresów diagnostycznych.
# wykresy diagnostyczne: > plot(fit2) Make a plot selection (or 0 to exit): 1: Time Series 2: Conditional SD 3: Series with 2 Conditional SD Superimposed 4: ACF of Observations 5: ACF of Squared Observations 6: Cross Correlation 7: Residuals 8: Conditional SDs 9: Standardized Residuals 10: ACF of Standardized Residuals 11: ACF of Squared Standardized Residuals 12: Cross Correlation between r^2 and r 13: QQ-Plot of Standardized Residuals |
Poniżej przedstawimy wykresy o numerach 9 i 13 czyli wykres standaryzowanych reszt oraz wykres kwantylowy dla standaryzowanych reszt.
> par(mfrow=c(1,2), mar=c(2.5,2.5,1.5,1)+0.1,mgp=c(1.5,0.3,0),tcl=0.3) > plot(fit2,which=c(9,13)) |
[ga1]
Inne modele
Dużo większe możliwości od omawianych powyżej bibliotek posiada pakiet rugarch. Oprócz wyżej wymienionych metod mamy możliwość estymacji modeli które uwzględniają np. efekt długiej pamięci. Poniżej kilka przykładów estymacji modeli z rodziny GARCH.
-
I-GARCH –
> library(rugarch) # model I-GARCH [Bollerslev, Engle, 1986]: > spec1= ugarchspec( variance.model= list(model= "iGARCH", garchOrder= c(1,1)), mean.model= list(armaOrder= c(0,0), include.mean= TRUE), distribution.model= "sstd")
-
EWMA
# model EWMA: > spec2= ugarchspec( variance.model=list(model="iGARCH", garchOrder=c(1,1)), mean.model=list(armaOrder=c(0,0), include.mean=TRUE), distribution.model="sstd", fixed.pars=list(omega=0))
-
E-GARCH
# model E-GARCH [Nelson 1991]: > spec3= ugarchspec( variance.model= list(model= "eGARCH", garchOrder= c(1,1)), mean.model= list(armaOrder= c(0,0), include.mean= TRUE), distribution.model= "sstd")
-
GJR-GARCH
# model GJR-GARCH [Glosten, Jaganathan, Runkle 1993]: > spec4= ugarchspec( variance.model= list(model= "gjrGARCH", garchOrder= c(1,1)), mean.model= list(armaOrder= c(0,0), include.mean= TRUE), distribution.model= "sstd")
-
T-ARCH
# model T-ARCH [Zakoian 1993]: > spec5= ugarchspec( variance.model= list(model= "apARCH", garchOrder= c(1,1)), mean.model= list(armaOrder= c(0,0), include.mean= TRUE), distribution.model= "sstd",fixed.pars= list(delta=1))
-
F-GARCH
# model F-GARCH [Hentschel 1995]: > spec6= ugarchspec( variance.model= list(model= "fGARCH", garchOrder= c(1,1), submodel="ALLGARCH"), mean.model= list(armaOrder= c(0,0), include.mean= TRUE), distribution.model= "sstd")
Ponieważ w analizowanym szeregu rtp1 nie występuje efekt długiej pamięci, więc oszacujemy model sstd-GARCH(1,1).
# model GARCH [Bollerslev, 1986]: > spec= ugarchspec( variance.model= list( model= "fGARCH", garchOrder= c(1,1), submodel= "GARCH"), mean.model= list(armaOrder= c(0,0), include.mean= TRUE), distribution.model= "sstd") > m= ugarchfit(spec= spec, data= rtp1);m *---------------------------------* * GARCH Model Fit * *---------------------------------* Conditional Variance Dynamics ----------------------------------- GARCH Model : fGARCH(1,1) fGARCH Sub-Model : GARCH Mean Model : ARFIMA(0,0,0) Distribution : sstd Optimal Parameters ------------------------------------ Estimate Std. Error t value Pr(>|t|) mu 0.000283 0.000302 0.93651 0.349010 omega 0.000008 0.000003 2.77806 0.005469 alpha1 0.088687 0.015849 5.59563 0.000000 beta1 0.899764 0.018122 49.65021 0.000000 skew 1.057767 0.023613 44.79586 0.000000 shape 5.872931 0.547738 10.72216 0.000000 Robust Standard Errors: Estimate Std. Error t value Pr(>|t|) mu 0.000283 0.000273 1.0384 0.29909 omega 0.000008 0.000005 1.6280 0.10352 alpha1 0.088687 0.025558 3.4700 0.00052 beta1 0.899764 0.031983 28.1322 0.00000 skew 1.057767 0.020246 52.2465 0.00000 shape 5.872931 0.848862 6.9186 0.00000 LogLikelihood : 9013.408 Information Criteria ------------------------------------ Akaike -4.9999 Bayes -4.9896 Shibata -5.0000 Hannan-Quinn -4.9963 Q-Statistics on Standardized Residuals ------------------------------------ statistic p-value Lag[1] 0.5233 0.4694 Lag[p+q+1][1] 0.5233 0.4694 Lag[p+q+5][5] 1.9711 0.8531 d.o.f=0 H0 : No serial correlation Q-Statistics on Standardized Squared Residuals ------------------------------------ statistic p-value Lag[1] 0.02531 0.8736 Lag[p+q+1][3] 0.08572 0.7697 Lag[p+q+5][7] 0.24237 0.9986 d.o.f=2 ARCH LM Tests ------------------------------------ Statistic DoF P-Value ARCH Lag[2] 0.02692 2 0.9866 ARCH Lag[5] 0.19748 5 0.9991 ARCH Lag[10] 0.41379 10 1.0000 Nyblom stability test ------------------------------------ Joint Statistic: 5.5481 Individual Statistics: mu 0.0558 omega 2.5267 alpha1 1.9511 beta1 2.5300 skew 0.3074 shape 2.5282 Asymptotic Critical Values (10% 5% 1%) Joint Statistic: 1.49 1.68 2.12 Individual Statistic: 0.35 0.47 0.75 Sign Bias Test ------------------------------------ t-value prob sig Sign Bias 0.3975 0.6910 Negative Sign Bias 0.8507 0.3950 Positive Sign Bias 0.6133 0.5397 Joint Effect 1.4143 0.7022 Adjusted Pearson Goodness-of-Fit Test: ------------------------------------ group statistic p-value(g-1) 1 20 156.2 1.389e-23 2 30 316.4 4.646e-50 3 40 498.1 5.752e-81 4 50 600.6 5.566e-96 |
Jak można zauważyć oprócz wyników estymacji danego modelu zostały podane także wyniki kilku testów diagnostycznych. Zaznaczmy, że wynik dowolnego testu oraz kilka innych informacji możemy uzyskać wpisując odpowiednią komendę.
# parametry modelu: > coef(m) mu omega alpha1 beta1 2.831571e-04 8.109069e-06 8.868727e-02 8.997642e-01 skew shape 1.057767e+00 5.872931e+00 # kryteria informacyjne: > infocriteria(m) Akaike -4.999949 Bayes -4.989642 Shibata -4.999955 Hannan-Quinn -4.996276 # logarytm wiarygodności: > likelihood(m) [1] 9013.408 # zestaw danych uzyskanych na podstawie modelu: > head(as.data.frame(m)) data fitted residuals sigma 1998-11-18 0.00000000 0.0002831571 -0.0002831571 0.02223207 1998-11-19 -0.02462686 0.0002831571 -0.0249100187 0.02128000 1998-11-20 0.04203294 0.0002831571 0.0417497860 0.02169304 1998-11-23 0.01778095 0.0002831571 0.0174977893 0.02420978 1998-11-24 -0.03518703 0.0002831571 -0.0354701851 0.02371975 1998-11-25 0.00538090 0.0002831571 0.0050977428 0.02501841 # test stabilności parametrów: > test1= nyblom(m) # testy obciążenia znaków: > test2= signbias(m) # testy zgodności rozkładu: > test3= gof(m,c(20,30,40,50)) |
Podobnie jak w pakiecie fGarch mamy możliwość wygenerować kilkanaście wykresów diagnostycznych dla modelu oszacowanego z wykorzystaniem biblioteki rugarch.
# wykresy diagnostyczne: > plot(m) Make a plot selection (or 0 to exit): 1: Series with 2 Conditional SD Superimposed 2: Series with 2.5% VaR Limits (with unconditional mean) 3: Conditional SD 4: ACF of Observations 5: ACF of Squared Observations 6: ACF of Absolute Observations 7: Cross Correlation 8: Empirical Density of Standardized Residuals 9: QQ-Plot of Standardized Residuals 10: ACF of Standardized Residuals 11: ACF of Squared Standardized Residuals 12: News-Impact Curve |
> par(mar=c(2.5,2.5,1.5,1.5)+0.1,mgp=c(1.5,0.3,0),tcl=0.3) > plot(m,which="all") # wszystkie wykresy |
[ga1]
Na koniec zaznaczmy, że nie sposób przedstawić wszystkie modele z rodziny GARCH, cały czas pojawiają się nowe metody. Ciekawą propozycją są modele beta-t-EGARCH dostępne w bibliotece betategarch. Z kolei zwolenników metod bayesowskich z pewnością zainteresuje biblioteka bayesGARCH. Wiele interesujących informacji można znaleźć na bologu http://www.unstarched.net/blog/. Są tam prezentowane ciekawe przykłady z wykorzystaniem różnych pakietów środowiska R np. biblioteki rugarch. Bardzo dobry opis możliwości tego pakietu jest dostępny na stronie http://cran.r-project.org/web/packages/rugarch/vignettes/Introduction_to_the_rugarch_package.pdf. Osoby zainteresowane modelami wielowymiarowymi GARCH (M-GARCH) powinny zwrócić uwagę np. na pakiet gogarch uogólnione modele ortogonalne GARCH (GO-GARCH) oraz ccgarch modele korelacji warunkowej GARCH (CC-GARCH) np. DCC-GARCH (dynamic conditional correlations GARCH) i CCC-GARCH (constant conditional correlation GARCH).
Krzysztof Trajkowski.
Konwersja i tak zostawila kwiatka: w jednym miejscu jest LaTeXowa para \begin{}/\end{}
Kawał dobrej roboty gratulacje! 😀