4 ARMA-GARCH-Modell
Zunächst beziehen wir für ein paar Unternehmen die Daten:
#Wir verwenden hierzu das tidyquant Paket, welches die Daten über Yahoo Finance bezieht
rm(list = ls())
library(tidyquant)
library(data.table)
library(plotly)
#lade yahoo finance ticker für den DAX
yh_tickers <- read.csv('/Users/ralfkellner/Library/Mobile Documents/com~apple~CloudDocs/Lehre Saarbrücken/QuantMeth/Codes/Markdown/Data/dax_yahoo_tickers.csv')
#yahoo finance tickers
tickers <- as.character(yh_tickers$Ticker)
#Wir beziehen manuell einen Datensatz mit den Daten der letzten 10 Jahre
our_data <- tq_get(tickers, from = "2010-02-01", to = "2020-02-03")
#transform to data table
df <- data.table(our_data)4.1 ARMA Modell
Im ersten Schritt schätzen wir ein ARMA-Modell für den bedingten Mittelwert. Hier beispielsweise für die Bayer AG Aktie. Die Log-Returns weisen leichte Autokorrelation auf.
temp_df <- df[symbol == 'BAYN.DE', ]
log_r <- diff(log(temp_df[['close']]))
auto <- acf(log_r, plot = FALSE)
p <- plot_ly() %>%
add_bars(x = as.vector(auto$lag), y = as.vector(auto$acf), width = 0.1, marker = list(color = 'black'), name = 'Autokorrelation') %>%
add_lines(x = c(0, 35), y = c(1.96/sqrt(auto$n.used), 1.96/sqrt(auto$n.used)), line = list(dash = 'dash', color = 'blue'), name = 'Kritischer Wert') %>%
add_lines(x = c(0, 35), y = c(-1.96/sqrt(auto$n.used), -1.96/sqrt(auto$n.used)), line = list(dash = 'dash', color = 'blue'), name = 'Kritischer Wert', showlegend = FALSE) %>%
layout(
title = 'Autokorrelation Bayer AG - Log Returns (2010/01 - 2020/02)',
xaxis = list(title = 'Lag')
)
pSchätzen des Modells mittels der Auto-Arima Funktion des Forecast-Pakets:
library(forecast)
#Wir verwenden die arima Funktion zum schätzen des Modells
fit_arma <- auto.arima(log_r)
#Die geschätzten Parameter des Modells sind:
print(fit_arma$coef)## ar1 ar2 ma1 ma2
## 0.5037188 -0.8359622 -0.5178449 0.8756247auto <- acf(fit_arma$residuals, plot = FALSE)
p <- plot_ly() %>%
add_bars(x = as.vector(auto$lag), y = as.vector(auto$acf), width = 0.1, marker = list(color = 'black'), name = 'Autokorrelation') %>%
add_lines(x = c(0, 35), y = c(1.96/sqrt(auto$n.used), 1.96/sqrt(auto$n.used)), line = list(dash = 'dash', color = 'blue'), name = 'Kritischer Wert') %>%
add_lines(x = c(0, 35), y = c(-1.96/sqrt(auto$n.used), -1.96/sqrt(auto$n.used)), line = list(dash = 'dash', color = 'blue'), name = 'Kritischer Wert', showlegend = FALSE) %>%
layout(
title = 'Autokorrelation Bayer AG - Residuals (2010/01 - 2020/02)',
xaxis = list(title = 'Lag')
)
pDie Autokorrelation wird geringer, ist jedoch nicht vollständig heraus gefiltert. Vergleichen wir den Verlauf der Log-Renditen mit einem simulierten Verlauf des von uns geschätzten Modells, wobei wir bei der Simulation konstante Varianz annehmen und diese durch die empirische Stichprobenvarianz geschätzt wird.
#Simuliere Daten mittels sim.arima()-Funktion
sim <- arima.sim(model = list(ar = fit_arma$coef[1:2], ma = fit_arma$coef[3:4]), n = length(log_r),
rand.gen = function(n){rnorm(n, 0, sd(log_r))})
p <- plot_ly() %>%
add_lines(x = temp_df$date[-c(1)], y = log_r, name = 'reale Daten', line = list(color = 'grey')) %>%
add_lines(x = temp_df$date[-c(1)], y = sim, name = 'simulierte Daten', line = list(color = 'olivedrab')) %>%
layout(
title = 'Log-Renditen Bayer - (2010/01 - 2020/02)',
xaxis = list(title = 'Datum')
)
pMan kann erkennen, dass die realen Daten über die Zeit unterschiedlich stark variieren, dies ist bei dem simulierten Prozess nicht so, da konstante Varianz angenommen wird. Um das reale Verhalten besser abbilden zu können, bedarf es eines Modells, welches über die Zeit varrierende bedingte Varianz abbildet. Hierzu eignet sich ein GARCH Modell.
4.2 GARCH Modell
Um ein GARCH Modell schätzen zu können, existieren zwei gute Pakete: fGarch und rugarch. Wir verwenden im Kurs das rugarch Paket.
library(rugarch)
#Zunächst müssen wir ein Modell definieren
spec <- ugarchspec(variance.model = list(model = 'sGARCH', garchOrder = c(1,1)),
mean.model = (list(armaOrder = c(0,0), include.mean = TRUE)),
distribution.model = 'norm')
#Wie fitten dieses Modell für die Residuen des vorherigen ARMA-Prozesses
fit_garch <- ugarchfit(spec, data = fit_arma$residuals)
# Der Zusammenfassung des Modells zeigt uns, dass für die betrachteten Lags keine Autokorrelation mehr in den Residuen und
#in den quadrierten Residuen statistisch nachgewiesen werden kann.
print(fit_garch)##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : sGARCH(1,1)
## Mean Model : ARFIMA(0,0,0)
## Distribution : norm
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 0.000385 0.000304 1.2683 0.204679
## omega 0.000002 0.000001 2.8126 0.004914
## alpha1 0.025376 0.002283 11.1165 0.000000
## beta1 0.965696 0.002235 431.9950 0.000000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.000385 0.000269 1.4317 0.15224
## omega 0.000002 0.000004 0.5652 0.57194
## alpha1 0.025376 0.002152 11.7910 0.00000
## beta1 0.965696 0.004536 212.9078 0.00000
##
## LogLikelihood : 6899.757
##
## Information Criteria
## ------------------------------------
##
## Akaike -5.4362
## Bayes -5.4269
## Shibata -5.4362
## Hannan-Quinn -5.4328
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.2553 0.6134
## Lag[2*(p+q)+(p+q)-1][2] 1.2389 0.4268
## Lag[4*(p+q)+(p+q)-1][5] 2.6583 0.4730
## d.o.f=0
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 1.432 0.2315
## Lag[2*(p+q)+(p+q)-1][5] 2.933 0.4193
## Lag[4*(p+q)+(p+q)-1][9] 3.768 0.6284
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 1.173 0.500 2.000 0.2787
## ARCH Lag[5] 1.390 1.440 1.667 0.6217
## ARCH Lag[7] 1.721 2.315 1.543 0.7761
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 64.1711
## Individual Statistics:
## mu 0.21348
## omega 4.02070
## alpha1 0.05813
## beta1 0.08770
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.07 1.24 1.6
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 0.25987 0.794982
## Negative Sign Bias 2.83781 0.004579 ***
## Positive Sign Bias 0.01921 0.984672
## Joint Effect 10.89584 0.012303 **
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 78.55 3.297e-09
## 2 30 90.26 3.314e-08
## 3 40 99.86 3.061e-07
## 4 50 131.84 1.599e-09
##
##
## Elapsed time : 0.1493371In Zukunft wollen wir die Modelle jedoch nicht getrennt, sondern auf einmal schätzen, daher definieren wir das ARMA-Modell in der spech-Funktion mit:
spec <- ugarchspec(variance.model = list(model = 'sGARCH', garchOrder = c(1,1)),
mean.model = (list(armaOrder = c(1,1), include.mean = TRUE)),
distribution.model = 'norm')
#Wie fitten dieses Modell für die Residuen des vorherigen ARMA-Prozesses
fit_arma_garch <- ugarchfit(spec, data = log_r)
# Der Zusammenfassung des Modells zeigt uns, dass für die betrachteten Lags keine Autokorrelation mehr in den Residuen und
#in den quadrierten Residuen statistisch nachgewiesen werden kann.
print(fit_arma_garch)##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : sGARCH(1,1)
## Mean Model : ARFIMA(1,0,1)
## Distribution : norm
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 0.000389 0.000302 1.28998 0.197057
## ar1 0.138573 1.269844 0.10913 0.913103
## ma1 -0.146164 1.267634 -0.11530 0.908204
## omega 0.000002 0.000001 2.80758 0.004992
## alpha1 0.025497 0.002283 11.16966 0.000000
## beta1 0.965605 0.002239 431.27722 0.000000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.000389 0.000274 1.42187 0.15507
## ar1 0.138573 0.616288 0.22485 0.82210
## ma1 -0.146164 0.615506 -0.23747 0.81229
## omega 0.000002 0.000004 0.56107 0.57475
## alpha1 0.025497 0.002132 11.95942 0.00000
## beta1 0.965605 0.004628 208.64768 0.00000
##
## LogLikelihood : 6895.02
##
## Information Criteria
## ------------------------------------
##
## Akaike -5.4308
## Bayes -5.4170
## Shibata -5.4309
## Hannan-Quinn -5.4258
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.06022 0.80614
## Lag[2*(p+q)+(p+q)-1][5] 4.45182 0.01914
## Lag[4*(p+q)+(p+q)-1][9] 9.08153 0.02216
## d.o.f=2
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 1.344 0.2464
## Lag[2*(p+q)+(p+q)-1][5] 3.184 0.3744
## Lag[4*(p+q)+(p+q)-1][9] 4.095 0.5729
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 1.404 0.500 2.000 0.2360
## ARCH Lag[5] 1.564 1.440 1.667 0.5757
## ARCH Lag[7] 1.897 2.315 1.543 0.7390
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 59.9724
## Individual Statistics:
## mu 0.22984
## ar1 0.61151
## ma1 0.61292
## omega 3.94708
## alpha1 0.05738
## beta1 0.08730
##
## 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.3902 0.696406
## Negative Sign Bias 2.8799 0.004012 ***
## Positive Sign Bias 0.2074 0.835698
## Joint Effect 11.0049 0.011699 **
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 72.89 3.019e-08
## 2 30 84.23 2.735e-07
## 3 40 100.59 2.425e-07
## 4 50 123.29 2.469e-08
##
##
## Elapsed time : 0.186178Aktuell verwendet unser Modell die Normalverteilung, wir können anhand eines pp-Plots oder eine qq-Plots graphisch überprüfen, wie gut die Verteilung an die standardisierten Residuen angepasst ist:
#Gemäß unserer Notation sind die Residuen des Modells: epsilon
eps_t <- fit_arma_garch@fit$residuals
#und die standardisierten Residuen können wir bestimmen über:
z_t <- fit_arma_garch@fit$residuals / fit_arma_garch@fit$sigma
#alternativ könnte man hierfür auch die generische residuals Funktion verwenden
head(eps_t)## [1] 0.007441601 -0.016563614 -0.017833540 -0.032697328 0.014918029
## [6] 0.005662821## [,1]
## 1970-01-02 01:00:00 0.007441601
## 1970-01-03 01:00:00 -0.016563614
## 1970-01-04 01:00:00 -0.017833540
## 1970-01-05 01:00:00 -0.032697328
## 1970-01-06 01:00:00 0.014918029
## 1970-01-07 01:00:00 0.005662821## [1] 0.4495746 -1.0109968 -1.0880929 -1.9901346 0.8754849 0.3334056## [,1]
## 1970-01-02 01:00:00 0.4495746
## 1970-01-03 01:00:00 -1.0109968
## 1970-01-04 01:00:00 -1.0880929
## 1970-01-05 01:00:00 -1.9901346
## 1970-01-06 01:00:00 0.8754849
## 1970-01-07 01:00:00 0.3334056#Betrachten wir nun einen pp-Plot unter Normalverteilungsannahme
p_theor <- pdist(distribution = 'norm', sort(z_t), mu = 0, sigma = 1)
p_emp <- rank(sort(z_t))/length(z_t)
pp_plot <- plot_ly() %>%
add_markers(x = p_theor, y = p_emp, name = 'pp-plot', marker = list(color = 'grey', symbol = 'x')) %>%
add_lines(x = c(0, 1), y = c(0, 1), line = list(color = 'black'), showlegend = FALSE) %>%
layout(
xaxis = list(title = 'Theoretische Verteilung'),
yaxis = list(title = 'Empirische Verteilung')
)
pp_plotEs scheint, als würde die Verteilungsannahme nicht ganz passen. Wir probieren es mit einer t Verteilung:
spec <- ugarchspec(variance.model = list(model = 'sGARCH', garchOrder = c(1,1)),
mean.model = (list(armaOrder = c(1,1), include.mean = TRUE)),
distribution.model = 'std')
#Wie fitten dieses Modell für die Residuen des vorherigen ARMA-Prozesses
fit_arma_garch <- ugarchfit(spec, data = log_r)
# Der Zusammenfassung des Modells zeigt uns, dass für die betrachteten Lags keine Autokorrelation mehr in den Residuen und
#in den quadrierten Residuen statistisch nachgewiesen werden kann.
print(fit_arma_garch)##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : sGARCH(1,1)
## Mean Model : ARFIMA(1,0,1)
## Distribution : std
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 0.000329 0.000265 1.2425 0.21406
## ar1 0.087649 0.421395 0.2080 0.83523
## ma1 -0.116053 0.419871 -0.2764 0.78224
## omega 0.000003 0.000002 1.4806 0.13871
## alpha1 0.037309 0.006773 5.5083 0.00000
## beta1 0.953325 0.007739 123.1804 0.00000
## shape 5.305393 0.617771 8.5880 0.00000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.000329 0.000233 1.41322 0.157591
## ar1 0.087649 0.273180 0.32085 0.748326
## ma1 -0.116053 0.274983 -0.42204 0.672997
## omega 0.000003 0.000006 0.43299 0.665023
## alpha1 0.037309 0.016469 2.26546 0.023484
## beta1 0.953325 0.021578 44.18117 0.000000
## shape 5.305393 1.096019 4.84060 0.000001
##
## LogLikelihood : 6995.859
##
## Information Criteria
## ------------------------------------
##
## Akaike -5.5095
## Bayes -5.4934
## Shibata -5.5096
## Hannan-Quinn -5.5037
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 1.365 0.2426214
## Lag[2*(p+q)+(p+q)-1][5] 5.617 0.0003524
## Lag[4*(p+q)+(p+q)-1][9] 10.314 0.0067548
## d.o.f=2
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.5572 0.4554
## Lag[2*(p+q)+(p+q)-1][5] 1.0413 0.8502
## Lag[4*(p+q)+(p+q)-1][9] 1.9490 0.9110
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 0.4679 0.500 2.000 0.4939
## ARCH Lag[5] 0.5313 1.440 1.667 0.8743
## ARCH Lag[7] 1.0688 2.315 1.543 0.9021
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 70.2842
## Individual Statistics:
## mu 0.2049
## ar1 0.4532
## ma1 0.4611
## omega 4.8532
## alpha1 0.3259
## beta1 0.2049
## shape 0.7011
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.69 1.9 2.35
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 0.08166 0.93492
## Negative Sign Bias 2.16623 0.03039 **
## Positive Sign Bias 0.41994 0.67457
## Joint Effect 7.94180 0.04723 **
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 18.36 0.4987
## 2 30 22.03 0.8188
## 3 40 39.35 0.4542
## 4 50 32.08 0.9706
##
##
## Elapsed time : 0.348285#Gemäß unserer Notation sind die Residuen des Modells: epsilon
eps_t <- fit_arma_garch@fit$residuals
#und die standardisierten Residuen können wir bestimmen über:
z_t <- fit_arma_garch@fit$residuals / fit_arma_garch@fit$sigma
#alternativ könnte man hierfür auch die generische residuals Funktion verwenden
head(eps_t)## [1] 0.007501502 -0.016347125 -0.018101387 -0.033038935 0.014260217
## [6] 0.005965336## [,1]
## 1970-01-02 01:00:00 0.007501502
## 1970-01-03 01:00:00 -0.016347125
## 1970-01-04 01:00:00 -0.018101387
## 1970-01-05 01:00:00 -0.033038935
## 1970-01-06 01:00:00 0.014260217
## 1970-01-07 01:00:00 0.005965336## [1] 0.4530395 -1.0017666 -1.1085354 -2.0135571 0.8230380 0.3463844## [,1]
## 1970-01-02 01:00:00 0.4530395
## 1970-01-03 01:00:00 -1.0017666
## 1970-01-04 01:00:00 -1.1085354
## 1970-01-05 01:00:00 -2.0135571
## 1970-01-06 01:00:00 0.8230380
## 1970-01-07 01:00:00 0.3463844#Betrachten wir nun einen pp-Plot unter Normalverteilungsannahme
p_theor <- pdist(distribution = 'std', sort(z_t), mu = 0, sigma = 1, shape = fit_arma_garch@fit$coef['shape'])
p_emp <- rank(sort(z_t))/length(z_t)
pp_plot <- plot_ly() %>%
add_markers(x = p_theor, y = p_emp, name = 'pp-plot', marker = list(color = 'grey', symbol = 'x')) %>%
add_lines(x = c(0, 1), y = c(0, 1), line = list(color = 'black'), showlegend = FALSE) %>%
layout(
xaxis = list(title = 'Theoretische Verteilung'),
yaxis = list(title = 'Empirische Verteilung')
)
pp_plotDies sieht besser aus! Alternativ kann man natürlich gemäß dem AIC oder anderen Maßen das beste Modell finden.
4.2.1 Prognose von Kursen
Nun wollen wir Kurse und Quantile für den Kurs prognostizieren. Wenn wir direkt den Kurs betrachten möchten, so werden wir feststellen, dass der Kurs an sich nicht stationär sein wird. Betrachten wir beispielsweise direkt den Schlusskurs der Bayer AG.
library(tseries)
t <- S_t <- df[symbol == 'BAYN.DE', ][['date']]
S_t <- df[symbol == 'BAYN.DE', ]$close
p <- plot_ly() %>%
add_lines(x = t, y = S_t, line = list(color = 'grey')) %>%
layout(
title = 'Bayer AG - Aktienkurs',
xaxis = list(title = 'Zeit'),
yaxis = list(title = 'Tagesschlusskurs')
)
p## [1] "Beim KPSS Test wird in der Nullhypothese Stationarität angenommen"##
## KPSS Test for Level Stationarity
##
## data: S_t
## KPSS Level = 10.968, Truncation lag parameter = 8, p-value = 0.01## [1] "Gemäß dem KPSS Test wird das Vorliegen von Stationarität bei einem Signifikanzniveau von 5% abgelehnt"## [1] "Beim ADF Test wird in der Nullhypothese Nicht-Stationarität angenommen"##
## Augmented Dickey-Fuller Test
##
## data: S_t
## Dickey-Fuller = -0.99586, Lag order = 13, p-value = 0.9394
## alternative hypothesis: stationary## [1] "Gemäß dem ADF Test wird die das Vorliegen von Nicht-Stationarität bei einem Signifikanzniveau von 5% nicht abgelehnt"Um die Zeitreihenmodelle verwenden zu können verwenden wir die erste Differenz des Kurses, dies kann alternativ zur (Log-)Rendite gemacht werden und vereinfacht eine direkte Prognose für den Kurs in diesem Fall. Gegeben die erste Differenz: \(\Delta_{S_t} = S_t - S_{t - 1}\), kann der Kurs für den nächsten Tag beispielsweise durch \(S_{t + 1} = S_t + \hat{\Delta_{S_{t + 1}}}\) mit z.B. \(\hat{\Delta_{S_{t + 1}}} = \hat{\mu_t + 1}\) des ARMA Prozesses von \(\Delta_{S_t}\) prognostiziert werden.
#Betrachten wir zunächst die ersten Differenzen
S_diff <- diff(S_t)
#und testen erneut auf Stationarität
kpss.test(S_diff)##
## KPSS Test for Level Stationarity
##
## data: S_diff
## KPSS Level = 0.21531, Truncation lag parameter = 8, p-value = 0.1## [1] "Gemäß dem KPSS Test wird die das Vorliegen von Stationarität bei einem Signifikanzniveau von 5% nicht abgelehnt"##
## Augmented Dickey-Fuller Test
##
## data: S_diff
## Dickey-Fuller = -14.996, Lag order = 13, p-value = 0.01
## alternative hypothesis: stationary## [1] "Gemäß dem ADF Test wird die das Vorliegen von Nicht-Stationarität bei einem Signifikanzniveau von 5% abgelehnt"p <- plot_ly() %>%
add_lines(x = t[-c(1)], y = S_diff, line = list(color = 'grey')) %>%
layout(
title = 'Bayer AG',
xaxis = list(title = 'Zeit'),
yaxis = list(title = 'Erste Differenz - Schlusskurse')
)
pIst ein Prozess nach dem ersten differenzieren stationär spricht man von einem integrierten Prozess des Grades 1. Wie vorher für die Log-Renditen schätzen wir ein ARMA(1,1)-GARCH(1,1) Modell für \(\Delta_{S_t}\) mit einer t Verteilung
spec <- ugarchspec(variance.model = list(model = 'sGARCH', garchOrder = c(1,1)),
mean.model = (list(armaOrder = c(1,1), include.mean = TRUE)),
distribution.model = 'std')
#Wie fitten dieses Modell für die Residuen des vorherigen ARMA-Prozesses
fit_arma_garch <- ugarchfit(spec, data = S_diff)
# Der Zusammenfassung des Modells zeigt uns, dass für die betrachteten Lags keine Autokorrelation mehr in den Residuen und
#in den quadrierten Residuen statistisch nachgewiesen werden kann.
print(fit_arma_garch)##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : sGARCH(1,1)
## Mean Model : ARFIMA(1,0,1)
## Distribution : std
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 0.023754 0.019603 1.21175 0.225607
## ar1 0.061213 0.380327 0.16095 0.872133
## ma1 -0.090476 0.379144 -0.23863 0.811391
## omega 0.006511 0.002831 2.29966 0.021467
## alpha1 0.035849 0.003455 10.37688 0.000000
## beta1 0.962195 0.002378 404.62820 0.000000
## shape 5.053684 0.522296 9.67590 0.000000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.023754 0.017081 1.39066 0.164329
## ar1 0.061213 0.231801 0.26408 0.791720
## ma1 -0.090476 0.233556 -0.38738 0.698472
## omega 0.006511 0.002804 2.32222 0.020221
## alpha1 0.035849 0.002860 12.53537 0.000000
## beta1 0.962195 0.000832 1156.51842 0.000000
## shape 5.053684 0.565221 8.94107 0.000000
##
## LogLikelihood : -4056.459
##
## Information Criteria
## ------------------------------------
##
## Akaike 3.2034
## Bayes 3.2195
## Shibata 3.2033
## Hannan-Quinn 3.2092
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.9827 0.3215277
## Lag[2*(p+q)+(p+q)-1][5] 5.4162 0.0007512
## Lag[4*(p+q)+(p+q)-1][9] 10.0288 0.0089634
## d.o.f=2
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.4506 0.5020
## Lag[2*(p+q)+(p+q)-1][5] 0.7362 0.9159
## Lag[4*(p+q)+(p+q)-1][9] 1.8178 0.9256
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 0.4736 0.500 2.000 0.4913
## ARCH Lag[5] 0.4749 1.440 1.667 0.8910
## ARCH Lag[7] 1.2971 2.315 1.543 0.8610
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 2.213
## Individual Statistics:
## mu 0.1339
## ar1 0.4887
## ma1 0.4991
## omega 0.1694
## alpha1 0.2951
## beta1 0.2229
## shape 0.4962
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.69 1.9 2.35
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 0.310 0.75656
## Negative Sign Bias 1.752 0.07996 *
## Positive Sign Bias 0.617 0.53729
## Joint Effect 4.872 0.18140
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 16.94 0.5941
## 2 30 32.04 0.3182
## 3 40 36.54 0.5825
## 4 50 50.01 0.4330
##
##
## Elapsed time : 0.2867072#Gemäß unserer Notation sind die Residuen des Modells: epsilon
eps_t <- fit_arma_garch@fit$residuals
#und die standardisierten Residuen können wir bestimmen über:
z_t <- fit_arma_garch@fit$residuals / fit_arma_garch@fit$sigma
#alternativ könnte man hierfür auch die generische residuals Funktion verwenden
head(eps_t)## [1] 0.3600468 -0.8054174 -0.8782777 -1.5565610 0.6523390 0.2729481## [,1]
## 1970-01-02 01:00:00 0.3600468
## 1970-01-03 01:00:00 -0.8054174
## 1970-01-04 01:00:00 -0.8782777
## 1970-01-05 01:00:00 -1.5565610
## 1970-01-06 01:00:00 0.6523390
## 1970-01-07 01:00:00 0.2729481## [1] 0.2618479 -0.5953219 -0.6562839 -1.1741715 0.4883557 0.2070026## [,1]
## 1970-01-02 01:00:00 0.2618479
## 1970-01-03 01:00:00 -0.5953219
## 1970-01-04 01:00:00 -0.6562839
## 1970-01-05 01:00:00 -1.1741715
## 1970-01-06 01:00:00 0.4883557
## 1970-01-07 01:00:00 0.2070026#Betrachten wir nun einen pp-Plot unter Normalverteilungsannahme
p_theor <- pdist(distribution = 'std', sort(z_t), mu = 0, sigma = 1, shape = fit_arma_garch@fit$coef['shape'])
p_emp <- rank(sort(z_t))/length(z_t)
pp_plot <- plot_ly() %>%
add_markers(x = p_theor, y = p_emp, name = 'pp-plot', marker = list(color = 'grey', symbol = 'x')) %>%
add_lines(x = c(0, 1), y = c(0, 1), line = list(color = 'black'), showlegend = FALSE) %>%
layout(
xaxis = list(title = 'Theoretische Verteilung'),
yaxis = list(title = 'Empirische Verteilung')
)
pp_plotZunächst können wir den sogenanten in-sample-fit betrachten, um zu sehen, wie gut das Modell beim schätzen an die Daten angepasst wurde.
#die vom Modell geschätzten Werte der bedingten Erwartungswerte über die Zeit sind:
in_fit <- fit_arma_garch@fit$fitted.values
#hierzu gehörige
p <- plot_ly() %>%
add_lines(x = t[-c(1)], y = S_diff, line = list(color = 'grey'), name = 'Erste Differenzen') %>%
add_lines(x = t[-c(1)], y = in_fit, line = list(color = 'purple'), name = 'In-Sample-Fit-Mean') %>%
add_lines(x = t[-c(1)], y = as.vector(quantile(fit_arma_garch, 0.01)),
line = list(color = 'olivedrab', dash = 'dash'), name = 'In-Sample-(0.01)-Quantile') %>%
add_lines(x = t[-c(1)], y = as.vector(quantile(fit_arma_garch, 0.99)),
line = list(color = 'olivedrab', dash = 'dot'), name = 'In-Sample-(0.99)-Quantile') %>%
layout(
title = 'Bayer AG',
xaxis = list(title = 'Zeit'),
yaxis = list(title = 'Erste Differenz - Schlusskurse')
)
pInteressanter ist es jeoch, ob wir für unser Modell auch neu hinzu kommende Daten gut prognostizieren können. Um dies zu evaluieren, wählen wir einen Rolling-Window-Ansatz. Im Beispiel erzeugen wir Prognosen für die Werte der letzten 250 Handelstage. Um einen Prognose für einen Tag zu erstellen, werden jeweils die vorherigen 1.000 Beobachtungen verwendet, um ein Modell zu schätzen und mit diesem Modell die Prognose zu erzeugen.
modelroll <- ugarchroll(spec = spec, S_diff, n.ahead = 1, forecast.length = 250,
refit.every = 1, refit.window = 'moving', window.size = 1000)
p <- plot_ly() %>%
add_lines(x = tail(t, n = 250), y = tail(S_diff, n = 250), line = list(color = 'grey'), name = 'Erste Differenzen') %>%
add_lines(x = tail(t, n = 250), y = modelroll@forecast$density$Mu,
line = list(color = 'purple'), name = 'Prognose-Erwartungswert') %>%
add_lines(x = tail(t, n = 250), y = modelroll@forecast$VaR$`alpha(1%)`,
line = list(color = 'olivedrab', dash = 'dash'), name = 'Prognose-(0.01)-Quantile') %>%
add_lines(x = tail(t, n = 250), y = modelroll@forecast$VaR$`alpha(5%)`,
line = list(color = 'olivedrab', dash = 'dash'), name = 'Prognose-(0.05)-Quantile') %>%
layout(
title = 'Bayer AG',
xaxis = list(title = 'Zeit'),
yaxis = list(title = 'Erste Differenz - Schlusskurse')
)
pp <- plot_ly() %>%
add_lines(x = tail(t, n = 250), y = tail(S_t, n = 250), line = list(color = 'grey'), name = 'Aktienkurs') %>%
add_lines(x = tail(t, n = 250), y = tail(S_t, 251)[-c(251)] + modelroll@forecast$density$Mu,
line = list(color = 'purple'), name = 'Prognose-Erwartungswert') %>%
add_lines(x = tail(t, n = 250), y = tail(S_t, 251)[-c(251)] + modelroll@forecast$VaR$`alpha(1%)`,
line = list(color = 'olivedrab', dash = 'dash'), name = 'Prognose-(0.01)-Quantile') %>%
add_lines(x = tail(t, n = 250), y = tail(S_t, 251)[-c(251)] + modelroll@forecast$VaR$`alpha(5%)`,
line = list(color = 'olivedrab', dash = 'dash'), name = 'Prognose-(0.05)-Quantile') %>%
layout(
title = 'Bayer AG',
xaxis = list(title = 'Zeit'),
yaxis = list(title = 'Schlusskurse')
)
pWenn unser Modell gut die Quantile bzw. den Value-at-Risk prognostiziert, sollte die Unterschreitungsfrequenz der realisierten Daten unterhalb des prognostizierten Quantils der zugehörigen Quantilsunterschreitungswahrscheinlichkeit entsprechen. Erwarten würden wir beispielsweise bei 250 Handelstagen, dass die realisierten Daten an 2.5 Tagen unter dem prognostizierten Quantil liegen.
## [1] "Unterscheitungsfrequenz für das 1%-Quantil"## [1] 0.004## [1] "Unterscheitungsfrequenz für das 5%-Quantil"## [1] 0.036Unser Modell scheint für den betrachteten Zeitraum etwas zu konservativ zu sein.
4.2.2 Weitere Variablen
Bei der Prognose können neben den eigenen vorherigen Werten weitere Variablen verwendet werden. Eine Möglichkeit bietet sich hierfür durch mutlivariate ARMA bzw. GARCH Modelle, die als VAR und MGARCH bezeichnet werden. Alternativ, wenn es wirklich nur auf die Prognose einer Zeitreihe ankommt, können auch multiple Regressionsmodelle mit Residuen, die einem GARCH Proess bzw. einem ARMA-GARCH Prozess folgen. Oftmals steckt jedoch in zeitverzögerten Variablen nicht allzu viel Information für die Zukunft. Betrachten wir beispielsweise die Korrelationen zwischen den ersten Differenzen der Volkswagen und der BMW Aktie inklusive zeitlichen Lags bis zu drei Tagen.
create_lag_data <- function(S, lag = 3, ticker = 'S'){
mat <- cbind(S)
col_names <- c(paste(ticker, 't'))
for(t in 1:lag){
temp_vec <- mat[,t][-c(1)]
mat <- cbind(mat[c(1:(dim(mat)[1] - 1)), ], temp_vec)
col_names <- c(paste(ticker, 't',t, sep = ""), col_names)
}
colnames(mat) <- col_names
return(data.frame(mat))
}
S_vw <- df[symbol == 'VOW3.DE']$close
S_bmw <- df[symbol == 'BMW.DE']$close
vw_data <- create_lag_data(diff(S_vw, 1), ticker = 'VW')
bmw_data <- create_lag_data(diff(S_bmw, 1), ticker = 'BMW')
cor(cbind(vw_data, bmw_data))## VWt3 VWt2 VWt1 VW.t BMWt3
## VWt3 1.000000000 0.10331613 -0.020384429 0.022359685 0.727678282
## VWt2 0.103316126 1.00000000 0.103964727 -0.020150194 0.067772857
## VWt1 -0.020384429 0.10396473 1.000000000 0.104378284 -0.002705644
## VW.t 0.022359685 -0.02015019 0.104378284 1.000000000 0.028644901
## BMWt3 0.727678282 0.06777286 -0.002705644 0.028644901 1.000000000
## BMWt2 0.078404973 0.72763301 0.067901059 -0.002641594 0.066446854
## BMWt1 -0.007443127 0.07867584 0.727579122 0.068082509 0.004888820
## BMW.t -0.002689478 -0.00703974 0.079479304 0.727615053 0.002320848
## BMWt2 BMWt1 BMW.t
## VWt3 0.078404973 -0.007443127 -0.002689478
## VWt2 0.727633013 0.078675840 -0.007039740
## VWt1 0.067901059 0.727579122 0.079479304
## VW.t -0.002641594 0.068082509 0.727615053
## BMWt3 0.066446854 0.004888820 0.002320848
## BMWt2 1.000000000 0.066511161 0.005002647
## BMWt1 0.066511161 1.000000000 0.066828668
## BMW.t 0.005002647 0.066828668 1.000000000Wir können erkennen, dass die kontemporäre Korrelation relativ hoch ist, jedoch die Korrelation der zeitverzögerten Variablen zu den aktuellen relativ gering ist. Dort wo kein Zusammenhang messbar ist, können Prognosen auch nicht mehr als den Zufall vorhersagen.