The Copula GARCH Model

Marius Hofert

2016-05-03

require(copula)
require(rugarch)
print(sessionInfo(), # for reproducibility
      locale=FALSE)
## R version 3.2.5 Patched (2016-04-18 r70508)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Debian GNU/Linux 8 (jessie)
## 
## attached base packages:
##  [1] parallel  grid      stats4    tools     stats     graphics  grDevices
##  [8] utils     datasets  methods   base     
## 
## other attached packages:
## [1] rugarch_1.3-6   gsl_1.9-10.1    lattice_0.20-33 bbmle_1.0.18   
## [5] copula_0.999-15
## 
## loaded via a namespace (and not attached):
##  [1] pcaPP_1.9-60                Rcpp_0.12.4                
##  [3] nloptr_1.0.4                formatR_1.3                
##  [5] ks_1.10.3                   xts_0.9.874                
##  [7] digest_0.6.9                evaluate_0.9               
##  [9] pspline_1.0-17              Matrix_1.2-6               
## [11] multicool_0.1-9             yaml_2.1.13                
## [13] expm_0.999-0                mvtnorm_1.0-5              
## [15] Runuran_0.23.0              stringr_1.0.0              
## [17] knitr_1.12.3                DistributionUtils_0.5-2    
## [19] ADGofTest_0.3               rgl_0.95.1447              
## [21] rmarkdown_0.9.6             Rsolnp_1.16                
## [23] magrittr_1.5                GeneralizedHyperbolic_0.8-3
## [25] SkewHyperbolic_0.3-4        codetools_0.2-14           
## [27] htmltools_0.3.5             MASS_7.3-45                
## [29] stabledist_0.7-1            spd_2.0-1                  
## [31] misc3d_0.8-4                numDeriv_2014.2-1          
## [33] KernSmooth_2.23-15          stringi_1.0-1              
## [35] truncnorm_1.0-7             zoo_1.7-13
set.seed(271)

1 Simulate data

First, we simulate the innovation distribution.

n <- 200 # sample size
d <- 2 # dimension
nu <- 3 # degrees of freedom for t
tau <- 0.5 # Kendall's tau
th <- iTau(ellipCopula("t", df=nu), tau) # corresponding parameter
cop <- ellipCopula("t", param=th, dim=d, df=nu) # define copula object
U <- rCopula(n, cop) # sample the copula
Z <- qnorm(U) # adjust margins
Now we simulate two ARMA(1,1)-GARCH(1,1) processes with these copula-dependent innovations. To this end, recall that an ARMA(\(p_1\),\(q_1\))-GARCH(\(p_2\),\(q_2\)) model is given by
## Set parameters
fixed.p <- list(mu  = 1,
                ar1 = 0.5,
                ma1 = 0.3,
                omega = 2, # alpha_0 (conditional variance intercept)
                alpha1= 0.4,
                beta1 = 0.2)
meanModel <- list(armaOrder=c(1,1))
varModel <- list(model="sGARCH", garchOrder=c(1,1)) # standard GARCH
uspec <- ugarchspec(varModel, mean.model=meanModel,
                    fixed.pars=fixed.p,
                    distribution.model="norm") # conditional innovation density (or use, e.g., "std")
## Note: ugarchpath(): simulate from a spec; ugarchsim(): simulate from a fitted object
X <- ugarchpath(uspec,
                n.sim= n, # simulated path length
                m.sim= d, # number of paths to simulate
                custom.dist=list(name="sample", distfit=Z)) # passing sample (n x d)-matrix
str(X, max.level=2) # => @path$sigmaSim, $seriesSim, $residSim
## Formal class 'uGARCHpath' [package "rugarch"] with 3 slots
##   ..@ path :List of 3
##   ..@ model:List of 10
##   ..@ seed : int 313548882
matplot(X@path$sigmaSim,  type="l") # plot of sigma's (conditional standard deviations)

matplot(X@path$seriesSim, type="l") # plot of X's

matplot(X@path$residSim,  type="l") # plot of Z's

par(pty="s")
plot(pobs(X@path$residSim)) # plot of Z's pseudo-observations => seem fine

2 Fitting procedure based on the simulated data

We now show how to fit an ARMA(1,1)-GARCH(1,1) process to X (remove ‘fixed.pars’ from specification to be able to fit).

uspec <- ugarchspec(varModel, mean.model=meanModel,
                    distribution.model="norm")
fit <- apply(X@path$seriesSim, 2, function(x) ugarchfit(uspec, x))
str(fit, max.level=3)
## List of 2
##  $ :Formal class 'uGARCHfit' [package "rugarch"] with 2 slots
##   .. ..@ fit  :List of 25
##   .. ..@ model:List of 11
##  $ :Formal class 'uGARCHfit' [package "rugarch"] with 2 slots
##   .. ..@ fit  :List of 25
##   .. ..@ model:List of 11
str(fit[[1]], max.level=2) # for first time series
## Formal class 'uGARCHfit' [package "rugarch"] with 2 slots
##   ..@ fit  :List of 25
##   ..@ model:List of 11
stopifnot(identical(fit[[1]]@fit$residuals, residuals(fit[[1]]@fit))) # => the same

Check the residuals.

Z. <- sapply(fit, function(x) x@fit$residuals)
U. <- pobs(Z.)
par(pty="s")
plot(U., # plot of Z's pseudo-observations => seem fine
     xlab=expression(italic(hat(U)[1])),
     ylab=expression(italic(hat(U)[2])))

Fit a t copula to the residuals Z.

fitcop <- fitCopula(ellipCopula("t", dim=2), data=U., method="mpl")
rbind(est = fitcop@estimate, true = c(th, nu)) # hat{rho}, hat{nu}; close to th, nu
##           [,1]     [,2]
## est  0.6766895 2.696942
## true 0.7071068 3.000000

3 Simulate from the fitted time series model

Simulate from the fitted copula model.

U.. <- rCopula(n, fitcop@copula)
Z.. <- qnorm(U..)
X..sim <- lapply(1:d, function(j)
                 ugarchsim(fit[[j]], n.sim=n, m.sim=1,
                           custom.dist=list(name="sample",
                           distfit=Z..[,j, drop=FALSE]))@simulation)
str(X..sim, max.level=3)
## List of 2
##  $ :List of 3
##   ..$ sigmaSim : num [1:200, 1] 1.51 1.49 1.88 1.65 2.37 ...
##   ..$ seriesSim: num [1:200, 1] 1.74 3.35 1.88 3.99 -2.05 ...
##   ..$ residSim : num [1:200, 1] 0.244 1.657 -0.996 2.63 -5.533 ...
##  $ :List of 3
##   ..$ sigmaSim : num [1:200, 1] 2.11 2.11 2.11 2.11 2.11 ...
##   ..$ seriesSim: num [1:200, 1] 1.204 1.792 0.835 2.944 -2.498 ...
##   ..$ residSim : num [1:200, 1] 0.25 0.643 -0.708 2.277 -5.099 ...
X..Z <- sapply(X..sim, `[[`, "residSim")
X.. <- sapply(X..sim, `[[`, "seriesSim")
plot(X..Z, main="residSim"); abline(h=0,v=0, lty=2, col=adjustcolor(1, .5))

plot(X.., main="seriesSim"); abline(h=0,v=0, lty=2, col=adjustcolor(1, .5))

matplot(pobs(X..), type="l")

par(pty="s")
plot(pobs(X..), main="pobs(series..)"); rect(0,0,1,1, border=adjustcolor(1, 1/2))