5 Multivariate Modellierung

Zunächst beziehen wir Gold- und Silberpreise der letzten 10 Jahre:

##          Date Goldpreis Silberpreis
## 1: 2010-01-04   777.470     11.9319
## 2: 2010-01-05   779.385     12.1845
## 3: 2010-01-06   786.579     12.4669
## 4: 2010-01-07   789.060     12.6195
## 5: 2010-01-08   786.837     12.6713
## 6: 2010-01-11   794.351     12.9618

5.1 VAR Modell

Um die bedingten Erwartungswerte zu modellieren, schätzen wir ein VAR Modell für die beiden Zeitreihen.

## 
## Call:
## ar(x = cbind(diff_gold, diff_silber), aic = FALSE, order.max = 1)
## 
## $ar
## , , 1
## 
##             diff_gold diff_silber
## diff_gold    0.045592     -0.3608
## diff_silber  0.009376     -0.2505
## 
## 
## $var.pred
##             diff_gold diff_silber
## diff_gold     104.746      1.7984
## diff_silber     1.798      0.1334

Insbesondere die beiden Elemente, die nicht auf der Diagonalen der AR-Parameter liegen sind interessant. \(\hat{\phi}_{1, 2} = -0.3608\) spricht dafür, dass der gestrige Silberpreis negativ mit dem heutigen Goldpreis korreliert, während \(\hat{\phi}_{2, 1} = -0.3608\) es sich anders herum unkorreliert verhält. Die geschätzen bedingten Erwartungswerte sehen folgendermaßen aus.

5.2 EWMA und DCC Modell

Nun wollen wir die bedingten Volatitlitäten bzw. Korrelationen durch ein EWMA und ein DCC Modell schätzen.

#zusätzliche Funktionen zum Schätzen des Modells

#######################################
#### multivariate Ljung--Box tests ####
#######################################
"mq" = function(x, lag = 1, df.adj = 0){
  # Compute multivariate Ljung-Box test statistics
  #
  x = as.matrix(x)
  nr = dim(x)[1]
  nc = dim(x)[2]
  g0 = var(x)
  ginv = solve(g0)
  qm = 0.0
  df = 0
  out = matrix(0,lag,4)
  for (i in 1:lag){
    x1 = x[(i+1):nr,]
    x2 = x[1:(nr-i),]
    g = cov(x1,x2)
    g = g*(nr-i-1)/(nr-1)
    h = t(g)%*%ginv%*%g%*%ginv
    qm = qm+nr*nr*sum(diag(h))/(nr-i)
    df = df+nc*nc
    pv = 1-pchisq(qm,df-df.adj)
    #  print(c(i,qm,pv))
    out[i,] = c(i,round(qm,2),round(df-df.adj),round(pv,3))
  }
  output = as.data.frame(matrix(out[lag,],1,4))
  names(output) = c("K", "Q(K)", "d.f.", "p-value")                         
  print(output)
}

mLjungBox = mq

"mq0" = function(A, m){
  N = dim(A)[1]
  k = dim(A)[2]
  temp = numeric(m)
  pval = numeric(m)
  out = as.data.frame(matrix(0,m+1,4))
  names(out) = c("K", "Q(K)", "d.f.", "p-value")
  out[,1] = 0:m
  
  Q.temp = N*sum(cor(A)[lower.tri(cor(A), diag = FALSE)]^2)
  out[1,2] = Q.temp 
  df.temp = k*(k-1)/2
  out[1,3] = df.temp
  out[1,4] = 1-pchisq(Q.temp, df.temp)
  
  for(j in 1:m){
    ccf = cor(A[-(1:j),],A[-((N-j+1):N),])
    Q.temp = Q.temp + N*(N+2)*sum(ccf^2)/(N-j)
    out[(j+1),2] = Q.temp 
    df.temp = df.temp+ k^2
    out[(j+1),3] = df.temp
    out[(j+1),4] = 1-pchisq(Q.temp, df.temp)
  }  
  round(out,3)
}


#EWMA Modell
"nllik.ewma" = function(lambda, innov){
  clambda = 1-lambda
  Sigma.hat = var(innov)
  invSigma.hat = chol2inv(chol(Sigma.hat)) # invSigma.hat = solve(Sigma.hat)
  detSigma.hat = det(Sigma.hat)
  llik = -0.5*log(detSigma.hat) - 0.5*crossprod(innov[1,],invSigma.hat)%*%innov[1,]
  llik = llik -0.5*log(detSigma.hat) - 0.5*crossprod(innov[2,],invSigma.hat)%*%innov[2,]
  n = dim(innov)[1]
  for(i in 3:n){
    atm1 = innov[(i-1),]
    at = innov[i,]
    denom = 1 - lambda^(i-1)
    #    Sigma.hat = clambda * tcrossprod(atm1) + lambda * Sigma.hat #approx
    Sigma.hat = (clambda/denom) * tcrossprod(atm1) + (lambda*(1-lambda^(i-2))/denom) * Sigma.hat #exact 
    invSigma.hat = chol2inv(chol(Sigma.hat)) # invSigma.hat = solve(Sigma.hat)
    detSigma.hat = det(Sigma.hat)
    llik = llik - 0.5*(log(detSigma.hat) + crossprod(at,invSigma.hat)%*%at)
  }
  nllik = -llik; nllik
}

"est.ewma" = function(lambda.0, innov){
  out = optim(lambda.0, nllik.ewma, lower = 0.001, upper = 0.999, innov = innov, method = "L-BFGS-B", hessian = TRUE)
  lambda.hat = out$par
  lambda.hat.se = 1/sqrt(out$hessian)
  list(lambda.hat = lambda.hat, lambda.hat.se = lambda.hat.se)
}

"sigma.ewma" = function(lambda, innov){
  clambda = 1-lambda
  Sigma.hat = var(innov)
  n = dim(innov)[1]
  d = dim(innov)[2]
  Sigma.t = array(0, c(d,d,n))
  Sigma.t[,,1:2] = Sigma.hat
  for(i in 3:n){
    atm1 = innov[(i-1),]
    at = innov[i,]
    denom = 1 - lambda^(i-1)
    #    Sigma.hat = clambda * tcrossprod(atm1) + lambda * Sigma.hat #approx
    Sigma.t[,,i] = (clambda/denom) * tcrossprod(atm1) + (lambda*(1-lambda^(i-2))/denom) * Sigma.t[,,(i-1)] #exact 
  }
  Sigma.t
}


###################################
#### DCCe estimation functions ####
###################################
"llik.DCCe" = function(theta, innov){
  # llik for the correlation component
  E = innov
  n = dim(E)[1]
  d = dim(E)[2]
  
  lambda = theta[1]
  ####  alpha = theta[1]; beta = theta[2];  omega = (1-alpha-beta)
  
  S = cor(E)
  Q = S
  Q.temp = diag(d)
  
  llik = 0
  for(t in 2:n){
    Q = (1-lambda) * E[t-1,] %*% t(E[t-1,]) + lambda * Q
    ####        Q = S * omega + alpha * E[t-1,] %*% t(E[t-1,]) + beta * Q
    diag(Q.temp) = 1/sqrt(diag(Q))
    R = (Q.temp) %*% Q %*% (Q.temp)
    llik = llik + log(det(R)) + t(E[t,])%*%solve(R)%*%E[t,] - t(E[t,])%*%E[t,]
  }
  llik/2    
}


"est.DCCe" = function(theta, innov){
  out = optim(theta, llik.DCCe, lower = 0.001, upper = 0.999, innov = innov, method = "L-BFGS-B", hessian = FALSE)
  out$par
}

"sigma.DCCe" = function(theta, innov){
  Y = innov
  n = dim(Y)[1]
  d = dim(Y)[2]
  
  D = matrix(0,n,d)
  E = matrix(0,n,d)
  
  for(i in 1:d){
    fit.temp = garchFit(data = Y[,i], include.mean = FALSE, trace = FALSE)
    D[,i] = sqrt(fit.temp@h.t)
    E[,i] = Y[,i]/sqrt(fit.temp@h.t)
  }
  
  lambda = theta[1]
  ####  alpha = theta[1]; beta = theta[2];  omega = (1-alpha-beta)
  
  S = cor(E)
  Q = S
  Q.temp = diag(d)
  
  Sigma.t = array(0, c(d,d,n))
  R.t = array(0, c(d,d,n))
  Sigma.t[,,1] = var(Y)
  R.t[,,1] = S
  
  for(t in 2:n){
    Q = (1-lambda) * E[t-1,] %*% t(E[t-1,]) + lambda * Q
    ####        Q = S * omega + alpha * E[t-1,] %*% t(E[t-1,]) + beta * Q
    diag(Q.temp) = 1/sqrt(diag(Q))
    R.t[,,t] = (Q.temp) %*% Q %*% (Q.temp)
    Sigma.t[,,t] = diag(D[t,]) %*% R.t[,,t] %*% diag(D[t,])
  }
  
  list(Sigma.t = Sigma.t, R.t = R.t)
  
}

"fit.DCCe" = function(theta.0 = 0.9, innov){
  Y = innov
  n = dim(Y)[1]
  d = dim(Y)[2]
  
  D = matrix(0,n,d)
  E = matrix(0,n,d)
  
  garch.omega.alpha.beta = matrix(0,d,3)
  
  for(i in 1:d){
    fit.temp = garchFit(data = Y[,i], include.mean = FALSE, trace = FALSE)
    garch.omega.alpha.beta[i,] = fit.temp@fit$matcoef[,1]
    D[,i] = sqrt(fit.temp@h.t)
    E[,i] = Y[,i]/sqrt(fit.temp@h.t)
  }
  
  DCCe.params = est.DCCe(theta = theta.0, innov = E)
  DCCe.params
  DCCe.Sigma = sigma.DCCe(theta = DCCe.params, innov=Y)
  
  list(Sigma.t = DCCe.Sigma$Sigma.t, R.t = DCCe.Sigma$R.t, coef = garch.omega.alpha.beta, lambda = DCCe.params)
  
}


#weitere Funktionen
"matrix.sqrt" = function(A)
{
  sva <- svd(A)
  if (min(sva$d)>=0)
    Asqrt <- t(sva$v %*% (t(sva$u) * sqrt(sva$d)))
  else
    stop("Matrix square root is not defined")
  return(Asqrt)
}

"matrix.sqrt.inv" = function(A)
{
  sva <- svd(A)
  if (min(sva$d)>=0)
    Asqrt <- t(sva$v %*% (t(sva$u) / sqrt(sva$d)))
  else
    stop("Matrix square root is not defined")
  return(Asqrt)
}

"orthogonalize" <- function(M)
{
  solve(matrix.sqrt(M%*%t(M)))%*%M
}

"pd.vcov.check" = function(Sigma){
  dex = NULL
  if(is.na(dim(Sigma)[3])) {
    sva <- svd(Sigma)
    temp = min(sva$d)   
    if (min(sva$d)<=0) {dex = c(dex, t)}
  }
  else{ 
    N = dim(Sigma)[3] 
    temp = numeric(N)
    for(t in 1:N){
      sva <- svd(Sigma[,,t])
      temp[t] = min(sva$d)
      if (min(sva$d)<=0) {dex = c(dex, t)}
    }   
  }         
  list(index = dex, min.sv = temp)
}
##   K Q(K) d.f. p-value
## 1 1 1.13    4    0.89
##   K Q(K) d.f. p-value
## 1 2 25.1    8   0.001
##   K  Q(K) d.f. p-value
## 1 3 26.79   12   0.008
##   K  Q(K) d.f. p-value
## 1 5 30.88   20   0.057
##    K   Q(K) d.f. p-value
## 1 10 117.45   40       0
##              [,1]       [,2]      [,3]
## [1,] 1.0051878084 0.08628705 0.9096534
## [2,] 0.0001215247 0.06556263 0.9383054
t1 <- list(
  text = "Volatilität Gold",
  xref = "paper",
  yref = "paper",
  yanchor = "bottom",
  xanchor = "center",
  align = "center",
  x = 0.5,
  y = 1,
  showarrow = FALSE
)

t2 <- list(
  text = "Volatilität Silber",
  xref = "paper",
  yref = "paper",
  yanchor = "bottom",
  xanchor = "center",
  align = "center",
  x = 0.5,
  y = 1,
  showarrow = FALSE
)

t3 <- list(
  text = "Kovarianz Gold/Silber",
  xref = "paper",
  yref = "paper",
  yanchor = "bottom",
  xanchor = "center",
  align = "center",
  x = 0.5,
  y = 1,
  showarrow = FALSE
)

t4 <- list(
  text = "Korrelation Gold/Silber",
  xref = "paper",
  yref = "paper",
  yanchor = "bottom",
  xanchor = "center",
  align = "center",
  x = 0.5,
  y = 1,
  showarrow = FALSE
)



p_sigma1 <- plot_ly() %>%
  add_lines(x = metal_prices$Date[-c(1,2)], y = EWMA_fit[1,1,]^(0.5), line = list(color = 'grey'), name = 'EWMA') %>% 
  add_lines(x = metal_prices$Date[-c(1,2)], y = DCC_fit$Sigma.t[1,1,]^(0.5), line = list(color = 'olivedrab'), name = 'DCC') %>% 
  layout(annotations = t1)

p_sigma2 <- plot_ly() %>%
  add_lines(x = metal_prices$Date[-c(1,2)], y = EWMA_fit[2,2,]^(0.5), line = list(color = 'grey'), name = 'EWMA', showlegend = FALSE) %>%
  add_lines(x = metal_prices$Date[-c(1,2)], y = DCC_fit$Sigma.t[2,2,]^(0.5), line = list(color = 'olivedrab'), name = 'DCC', showlegend = FALSE) %>%
  layout(annotations = t2)

p_Cov <- plot_ly() %>%
  add_lines(x = metal_prices$Date[-c(1,2)], y = EWMA_fit[1,2,], line = list(color = 'grey'), name = 'EWMA', showlegend = FALSE) %>% 
  add_lines(x = metal_prices$Date[-c(1,2)], y = DCC_fit$Sigma.t[1,2,], line = list(color = 'olivedrab'), name = 'DCC', showlegend = FALSE) %>%
  layout(annotations = t3)

p_Cor <- plot_ly() %>%
  add_lines(x = metal_prices$Date[-c(1,2)], y = EWMA_fit[1,2,] / sqrt(EWMA_fit[1,1,] * EWMA_fit[2,2,]), line = list(color = 'grey'), name = 'EWMA', showlegend = FALSE) %>% 
  add_lines(x = metal_prices$Date[-c(1,2)], y = DCC_fit$R.t[2,1,], line = list(color = 'olivedrab'), name = 'DCC', showlegend = FALSE) %>%
  layout(annotations = t4)

subplot(p_sigma1, p_sigma2, p_Cov, p_Cor, nrows = 2,
        titleX = TRUE, titleY = TRUE, margin = 0.1) %>% layout(showlegend = TRUE)

Beim Vergleich der beiden Modelle fällt auf, dass die Korrelation im DCC wesentlich weniger stark über die Zeit hinweg variiert. Um zu überprüfen, ob das VAR und das multivariate GARCH Modell die bedingten Erwartungswerte und Kovarianzen in den Daten gut abbilden können, betrachten wir im letzen Schritt die geschätzten Innovationen \(\hat{\boldsymbol{Z}}_t = \hat{\Sigma}_t^{-1/2} \hat{\boldsymbol{\epsilon}}_t\). Diese sollten keinerlei von null verschiedenen (Cross-)Autokorrelationen und bedingten Korrelation aufweisen. Ohne hier auf Details einzugehen, kann hierfür der multivariate Ljung-Box-Test auf die quadrierten Werte \(\hat{\boldsymbol{Z}}_t^2\) und der univariate Ljung-Box-Test auf die Produktzeitreihe \(\hat{Z}_{1,t} \hat{Z}_{2,t}\) angewendet werden.

##   K  Q(K) d.f. p-value
## 1 1 58.81    4       0
##   K  Q(K) d.f. p-value
## 1 2 61.66    8       0
##   K  Q(K) d.f. p-value
## 1 3 62.19   12       0
##   K  Q(K) d.f. p-value
## 1 5 69.91   20       0
##    K   Q(K) d.f. p-value
## 1 10 100.16   40       0
##   K Q(K) d.f. p-value
## 1 1 28.5    4       0
##   K  Q(K) d.f. p-value
## 1 2 28.66    8       0
##   K  Q(K) d.f. p-value
## 1 3 29.41   12   0.003
##   K  Q(K) d.f. p-value
## 1 5 32.96   20   0.034
##    K  Q(K) d.f. p-value
## 1 10 58.94   40   0.027
##   K Q(K) d.f. p-value
## 1 1 0.35    1   0.553
##   K Q(K) d.f. p-value
## 1 2 0.45    2   0.798
##   K Q(K) d.f. p-value
## 1 3 0.99    3   0.804
##   K Q(K) d.f. p-value
## 1 5 1.81    5   0.874
##    K Q(K) d.f. p-value
## 1 10 5.68   10   0.841
##   K Q(K) d.f. p-value
## 1 1 0.43    1   0.512
##   K Q(K) d.f. p-value
## 1 2 0.44    2   0.801
##   K Q(K) d.f. p-value
## 1 3 1.03    3   0.795
##   K Q(K) d.f. p-value
## 1 5 1.05    5   0.958
##    K Q(K) d.f. p-value
## 1 10 3.93   10   0.951

Die Resultate sprechen dafür, dass durch beide Modelle zwar die zeitvariierenden Korrelation, jedoch nicht die bedingten Varianzen abgebildet werden. Ein Blick auf den Verlauf der geschätzten Innovationen bestätigt diesen Eindruck, zeigt jedoch auch, dass dieses Resultat auch durch wenige extreme Realisierungen zustande kommt. Für den Moment müssen wir uns mit diesem Ergebnis zufrieden geben.

5.3 Multivariate Modellierung mit Copulas

Ignorieren wir für einen Moment Abhängigkeiten über den Zeitverlauf und betrachten, wie man Abhängigkeiten mittels Copulas analyisieren kann. In R existiert das copula-Paket, welches in meinen Augen den Standard darstellt. Ich persönlich finde jedoch das VineCopula-Paket nützlicher. In diesem ist nicht nur eine Vielzahl an verschiedenen Copulamodellen implementiert, sondern auch nützliche Funktionen zur Umrechnung von Abhängigkeitsparametern und zur Quantifizierung von Tail-Dependence.

Typischerweise wird die Stärke der Abhängigkeit im Zusammenhang mit Copulas durch Kendall’s Tau \(\rho_{\tau}\) gemessen. Kenall’s Tau ist ein Konkordanzmaß und durch:

\[ \rho_{\tau} = P((X - \tilde{X})(Y - \tilde{Y}) > 0) - P((X - \tilde{X})(Y - \tilde{Y}) < 0) \]

definiert. Etwas intuitiver wird es, wenn wir uns dem Schätzer für Kendall’s Tau:

\[ \hat{\rho}_{\tau} = \binom n 2^{-1} \sum_{1 \leq i < j \leq n} \text{sign} \lbrace (X_i - X_j)(Y_i - Y_j) \rbrace \] ansehen. \(\binom n 2^{-1}\) gibt die Anzahl an Paarkombination im Summanden an, womit wir es mit der durchschnittlichen Anzahl gleichläufiger bzw. gegenläufiger Paarkombinationen zu tun haben. \(\text{sign}\) steht für die Vorzeichenfunktion mit

\[ \text{sign}(x) = \begin{cases} ~~~1 & \text{ falls } x > 0 \\ -1 & \text{ falls x < 0 } \\ ~~~0 & \text{ falls } x = 0 \end{cases} \]

Kendall’s Tau wird verwendet, da es invariant gegenüber Transformationen ist. Dies bedeutet, die gemessene Abhängigkeit wird nicht durch beliebige Veränderungen der ursprünglichen Variablen verändert. Im Copulakontext ist diese Eigenschaft besonders wertvoll, da somit:

\[ \hat{\rho}_{\tau} (\boldsymbol{X}, \boldsymbol{Y}) = \hat{\rho}_{\tau} (F_X(\boldsymbol{X}), F_Y(\boldsymbol{Y})) = \hat{\rho}_{\tau} (\boldsymbol{U}, \boldsymbol{V}) \]

Für die meisten Copulas besteht ein direkter Zusammenhang zwischen Kendall’s Tau und dem Copula-spezifischen Parameter. Gehen wir, z.B. von einem Wert $_{} = 0.60 $ aus und betrachten, wie unterschiedlich bei dieser Abhängigkeitsstärke die Abhängigkeiten verschiedener Copulas aussehen.

## [1] 0.809017
## [1] 3
## [1] 2.5

Gegben die Parameter der jeweiligen Copula lässt sich auch die Tail-Dependence bestimmen, falls diese vorhanden ist.

## $lower
## [1] 0
## 
## $upper
## [1] 0
## $lower
## [1] 0.7937005
## 
## $upper
## [1] 0
## $lower
## [1] 0
## 
## $upper
## [1] 0.6804921

Falls Sie sich fragen, was dies mit den von uns in der realen Welt beobachteten Daten zu tun hat. Stellen Sie sich vor \((\boldsymbol{X}, \boldsymbol{Y})\) sei eine Serie von standardnormal verteilten Renditen. Gegeben die Abhängigkeitsstruktur entspricht einer Gauss Copula, dann simulieren wir die Renditen:

##           [,1]      [,2]
## [1,] 0.9148060 0.9777229
## [2,] 0.2861395 0.5418258
## [3,] 0.6417455 0.6262445
## [4,] 0.7365883 0.4454084
## [5,] 0.6569923 0.7401767
## [6,] 0.4577418 0.6007080

Betrachten wir dies zwischen Gold und Silber.

Schätzen wir das gemäß AIC Kriterium beste Copula Modell:

## Bivariate copula: t (par = 0.52, par2 = 5.35, tau = 0.34)

Der Output verät uns, dass die beobachteten Daten am besten zu einer Student t Copula passen. Diese hat eine geschätzte Korrelation \(\hat{\rho} = 0.52\) was einem Wert für Kendall’s Tau von \(\hat{\rho}_{\tau} = 0.34\). Zudem sprechen die geschätzten Freiheitsgrade \(\hat{\nu} = 5.35\) für relativ hohe untere und obere Tail-Dependence, da allgmein \(\nu \in (0, \infty)\). Je kleiner der Wert, umso höher die Tail-Dependence, hier:

## $lower
## [1] 0.2038857
## 
## $upper
## [1] 0.2038857

In der eben eingenommenen Betrachtung ist die Abhängigkeitsstärke in Form von Kendall’s Tau und der unteren sowie oberen Tail-Dependence konstant über die Zeit. Bei Zeitreihen können diese Parameter auch über die Zeit varriieren. Daher bietet es sich beipspielsweise an, die Korrelation durch ein DCC Modell zu modellieren. Aus dieser Überlegung entsteht ein DCC-Copula Modell, welches z.B. durch das rmgarch-Paket geschätzt werden kann.

## 
## *-------------------------------------------------*
## *                  Copula GARCH Fit               *
## *-------------------------------------------------*
## 
## Distribution     :  mvt
## DCC Order            :  1 1
## Asymmetric           :  FALSE
## No. of Parameters    :  18
## [VAR GARCH DCC UncQ]: [0+14+3+1]
## No. of Series        :  2
## No. of Observations  :  2521
## Log-Likelihood       :  -8601.427
## Av.Log-Likelihood    :  -3.412 
## 
## Optimal Parameters
## ---------------------------------------------------
##                       Estimate  Std. Error   t value Pr(>|t|)
## [diff_gold].mu        0.093524    0.151228  0.618429 0.536292
## [diff_gold].ar1      -0.703868    0.115536 -6.092178 0.000000
## [diff_gold].ma1       0.736753    0.109218  6.745725 0.000000
## [diff_gold].omega     0.957829    0.412338  2.322923 0.020183
## [diff_gold].alpha1    0.060411    0.014362  4.206296 0.000026
## [diff_gold].beta1     0.931759    0.016018 58.168643 0.000000
## [diff_gold].shape     5.207043    0.539327  9.654703 0.000000
## [diff_silber].mu     -0.000687    0.003198 -0.214729 0.829978
## [diff_silber].ar1    -0.040952    0.199616 -0.205153 0.837453
## [diff_silber].ma1    -0.015585    0.199252 -0.078215 0.937657
## [diff_silber].omega   0.000261    0.000163  1.603958 0.108723
## [diff_silber].alpha1  0.067837    0.022970  2.953342 0.003144
## [diff_silber].beta1   0.931163    0.023040 40.415666 0.000000
## [diff_silber].shape   5.043629    0.487430 10.347390 0.000000
## [Joint]dcca1          0.023743    0.011441  2.075223 0.037966
## [Joint]dccb1          0.790065    0.134903  5.856533 0.000000
## [Joint]mshape        10.075991    2.272116  4.434628 0.000009
## 
## Information Criteria
## ---------------------
##                    
## Akaike       6.8373
## Bayes        6.8766
## Shibata      6.8372
## Hannan-Quinn 6.8516
## 
## 
## Elapsed time : 4.537034
#Betrachten wir erneut die zeitvariiernde

t1 <- list(
  text = "Volatilität Gold",
  xref = "paper",
  yref = "paper",
  yanchor = "bottom",
  xanchor = "center",
  align = "center",
  x = 0.5,
  y = 1,
  showarrow = FALSE
)

t2 <- list(
  text = "Volatilität Silber",
  xref = "paper",
  yref = "paper",
  yanchor = "bottom",
  xanchor = "center",
  align = "center",
  x = 0.5,
  y = 1,
  showarrow = FALSE
)

t3 <- list(
  text = "Kovarianz Gold/Silber",
  xref = "paper",
  yref = "paper",
  yanchor = "bottom",
  xanchor = "center",
  align = "center",
  x = 0.5,
  y = 1,
  showarrow = FALSE
)

t4 <- list(
  text = "Korrelation Gold/Silber",
  xref = "paper",
  yref = "paper",
  yanchor = "bottom",
  xanchor = "center",
  align = "center",
  x = 0.5,
  y = 1,
  showarrow = FALSE
)



p_sigma1 <- plot_ly() %>%
  add_lines(x = metal_prices$Date[-c(1)], y = fit@mfit$H[1,1,]^(0.5), line = list(color = 'grey'), name = 'DCC-Copula') %>% 
  layout(annotations = t1)

p_sigma2 <- plot_ly() %>%
  add_lines(x = metal_prices$Date[-c(1)], y = fit@mfit$H[2,2,]^(0.5), line = list(color = 'grey'), name = 'DCC-Copula', showlegend = FALSE) %>%
  layout(annotations = t2)

p_Cov <- plot_ly() %>%
  add_lines(x = metal_prices$Date[-c(1)], y = fit@mfit$H[1,2,], line = list(color = 'grey'), name = 'DCC-Copula', showlegend = FALSE) %>% 
  layout(annotations = t3)

p_Cor <- plot_ly() %>%
  add_lines(x = metal_prices$Date[-c(1)], y = fit@mfit$H[1,2,] / sqrt(fit@mfit$H[1,1,] * fit@mfit$H[2,2,]), line = list(color = 'grey'), name = 'DCC-Copula', showlegend = FALSE) %>% 
  layout(annotations = t4)

subplot(p_sigma1, p_sigma2, p_Cov, p_Cor, nrows = 2,
        titleX = TRUE, titleY = TRUE, margin = 0.1) %>% layout(showlegend = TRUE)