##############################################################
## Copula estimation and goodness-of-fit
##############################################################

library(rugarch)
library(xts)
library(copula)
library(qrmdata)
library(qrmtools)

## Select the data 
data("SP500_const") # load the constituents data of the S&P 500
stocks <- c("INTC", "QCOM", "GOOGL", "AAPL", "MSFT") # Intel, Qualcomm, Google, Apple, Microsoft
time <- c("2007-01-03", "2009-12-31") # time period
S <- SP500_const[paste0(time, collapse = "/"), stocks] # data

## Build negative log-returns
X <- -returns(S) # -log-returns
plot.zoo(X, main = "Negative log-returns")

### Fitting marginal ARMA(1,1)-GARCH(1,1) models with standardized t residuals
## componentwise fitting of univariate ARMA-GARCH processes via qrmtools

uspec <- rep(list(ugarchspec(distribution.model = "std")), ncol(X))
fit.ARMA.GARCH <- fit_ARMA_GARCH(X, ugarchspec.list = uspec)
fits <- fit.ARMA.GARCH$fit # fitted models
fits

# extract standardized residuals
Z <- as.matrix(do.call(merge, lapply(fits, residuals, standardize = TRUE))) 
colnames(Z) <- colnames(S)
Z[1:5, ]

# vector of estimated df
nu.mar<-c(fits[[1]]@fit$coef[["shape"]], fits[[2]]@fit$coef[["shape"]], 
  fits[[3]]@fit$coef[["shape"]], fits[[4]]@fit$coef[["shape"]], 
  fits[[5]]@fit$coef[["shape"]])
nu.mar

n <- nrow(X) # sample size
d <- ncol(X) # dimension

### Fitting copula models ##########################################################
## Compute pseudo-observations from the standardized t residuals

U <- pobs(Z)  # componentwise scaled ranks 
pairs2(U, cex = 0.4)
cor.test(U[,1],U[,2], method="kendall")


## Fitting a Gumbel copula
# Estimation of Copula Parameters via the MPLE using the 
#   available pseudo-observations 
fit.gc <- fitCopula(gumbelCopula(dim = d),
        data = U, method = "mpl")  
summary(fit.gc)
fit.gc@estimate # estimated copula parameter
gc <- fit.gc@copula # fitted copula

## Fitting a t copula

# For the t family a convenient estimation technique
# combines method-of-moments estimation based on Kendall’s tau 
# for the estimation of the correlation matrix and 
# maximum pseudo-likelihood estimation for the estimation of the 
# number of dof by Mashal and Zeevi (2002) 

fit.tc <- fitCopula(tCopula(dim = d, dispstr = "un"),
                    data = U, method = "itau.mpl")
summary(fit.tc)
nu <- tail(fit.tc@estimate, n = 1) # estimated degrees of freedom nu
fit.tc@estimate[-11] # vector of d(d-1)/2=10 corr. coeff. 
P <- p2P(fit.tc@estimate[-11]) # estimated correlation matrix

tc <- fit.tc@copula # fitted copula


### Goodness-of-fit ##########################################################

## parametric bootstrap based on the maximum pseudo-likelihood estimator.

set.seed(270) # for reproducibility
N <- 1000

### A test of exchangeability 
exchTest(U, N = N)

## Check the Gumbel copula
gof.gc <- gofCopula(gc, x = U, N = N) 
# warning also because the copula doesn't fit well here
gof.gc 
# p-value=0.0004995

## Check the t copula:
t.copula <- tCopula(dim = d, dispstr = "un", df.fixed = TRUE)
gofCopula(t.copula, x=U, N=N, simulation="mult")  

# Check the model graphically via the Rosenblatt transform
# The Rosenblatt transform of a d-dimensional random vector U ∼ C 
# is the d-dimensional random vector R_C(U)∼ Pi, where P is the independence copula. 
# This can be used for a graphical goodness-of-fit test

U.Rsnbl <- cCopula(U, copula = tc)
pairs2(U.Rsnbl, cex = 0.4)

#Simulate from the chosen t copula with parameter from 'mpl'
Cop<-rCopula(n, copula=tc)
dim(Cop)
pairs2(Cop, cex = 0.4)


###Predict the aggregated loss and VaR_0.99 ##########################
## Predict VaR of the aggregated loss nonparametrically

B <- 200
m <- ceiling(n/20) # length of the simulates paths
X.lst <- lapply(1:B, function(b){
  ## 1) Simulate from the fitted copula
  U <- rCopula(m, copula = tc)
  ## 2) Quantile-transform to standardized t distributions (for ugarchsim())
  Zt <- sapply(1:d, function(j) sqrt((nu.mar[j]-2)/nu.mar[j]) * qt(U[,j], df = nu.mar[j]))
  ## 3) Use these t innovations to sample from the time series
  sim <- lapply(1:d, function(j)
    ugarchsim(fits[[j]], n.sim = m, m.sim = 1,
              custom.dist = list(name = "sample",
              distfit = Zt[,j, drop = FALSE])))
  ## 4) Extract simulated series
  sapply(sim, function(x) fitted(x)) 
})

XSim <- sapply(X.lst, rowSums) # simulated aggregated losses:  
dim(XSim) #(m, B)-matrix

#XS <- rowSums(X) # aggregated loss; n-vector
XSim.M<- rowMeans(XSim) # predicted aggregated loss; m-vector
XSim.CI <- apply(XSim, 1, function(x) 
  quantile(x, probs = c(0.025, 0.975))) # CIs; (2, m)-matrix

alpha <- 0.99 # confidence level
# VaR_alpha; m-vector
VaR <- apply(XSim, 1, function(x) quantile(x, probs = alpha)) 

XSim.M
XSim.CI
VaR