Wybrane pakiety do analizy finansowych szeregów czasowych w R

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")

Wykres liniowy dla spółki Adobe od 03-01-2007 do 28-03-2013.’ [TAF1]

# wykres świecowy:
> chartSeries(ADBE,theme= chartTheme("white"),TA=NULL)

Wykres świecowy dla spółki Adobe od 03-01-2007 do 28-03-2013.’ [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")

Wykres świecowy dla spółki Adobe - marzec 2013.’ [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")

Wykres świecowy dla spółki Adobe - marzec 2013.’ [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()))

Wykresy dla spółki Adobe - rok 2012.’ [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 &lt;- 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")

Dzienne ceny zamknięcia dla spółki TPsa od 18-11-1998 do 31-01-2013.’ [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)

Dzienne stopy zwrotu spółki TPsa od 18-11-1998 do 28-03-2013.’ [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&#39;Agostino Normality Test
 
Test Results:
  STATISTIC:
    Chi2 | Omnibus: 1110.8943
    Z3  | Skewness: -17.923
    Z4  | Kurtosis: 28.1009
  P VALUE:
    Omnibus  Test: &lt; 2.2e-16 
    Skewness Test: &lt; 2.2e-16 
    Kurtosis Test: &lt; 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))

Wykres kwantylowy oraz gęstości - dzienne stopy zwrotu.’ [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&#39;Agostino Normality Test
 
Test Results:
  STATISTIC:
    Chi2 | Omnibus: 269.7419
    Z3  | Skewness: -10.7334
    Z4  | Kurtosis: 12.4313
  P VALUE:
    Omnibus  Test: &lt; 2.2e-16 
    Skewness Test: &lt; 2.2e-16 
    Kurtosis Test: &lt; 2.2e-16
# miesięczne stopy zwrotu:
> rtp3= periodReturn(txts1, type="log", period= "monthly")
> z3= as.numeric(rtp3)
> dagoTest(z3)
 
Title:
 D&#39;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&#39;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))

Wykres kwantylowy oraz gęstości - dzienne stopy zwrotu.’ [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)

Autokorelacja kwadratów stóp zwrotu spółki TPsa.’ [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))

Ceny i stópy zwrotu spółki TPsa.’ [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  =  − 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)

Autokorelacja stóp zwrotu spółki TPsa.’ [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  &lt; 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 &lt; 2.2e-16
 
    Box-Ljung test
 
data:  Squared.Residuals 
X-squared = 0.0024, df = 1, p-value = 0.9607
# logarytm wiarygodności:
> logLik(n)
&#39;log Lik.&#39; 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  &lt; 2e-16 ***
delta  1.123e+00   1.372e-01    8.182 2.22e-16 ***
skew   1.059e+00   2.381e-02   44.450  &lt; 2e-16 ***
shape  5.868e+00   5.436e-01   10.794  &lt; 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  &lt; 2e-16 ***
skew   1.058e+00   2.361e-02   44.796  &lt; 2e-16 ***
shape  5.877e+00   5.499e-01   10.689  &lt; 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))

Wykres reszt oraz kwantylowy.’ [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

Wykresy diagnostyczne.’ [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.

2 thoughts on “Wybrane pakiety do analizy finansowych szeregów czasowych w R”

Dodaj komentarz

Twój adres email nie zostanie opublikowany. Pola, których wypełnienie jest wymagane, są oznaczone symbolem *