# R/CCC.R In ccgarch2: Conditional Correlation GARCH Models

#### Documented in estimateCCCprint.cccprint.summary.cccsummary.ccc

```####################################
loglik.ccc.t <- function(param, data, model, mode){
if(is.zoo(data)) data <- as.matrix(data)
nobs <- dim(data)[1]
ndim <- dim(data)[2]
mu <- matrix(param[1:ndim], nobs, ndim, byrow = TRUE)    # constant in the mean
data <- data - mu
param <- param[-(1:ndim)]

para.mat <- p.mat(param, model, ndim)
# check if R is positive definite. If not, making R positive definite
R <- make.pd(para.mat\$R)

h <- vgarch(para.mat\$a, para.mat\$A, para.mat\$B, data)    # a call to vgarch function
z <- data/sqrt(h)
invR <- solve(R)
lf <- -0.5*ndim*log(2*pi) - 0.5*rowSums(log(h)) - 0.5*log(det(R)) - 0.5*rowSums((z%*%invR)*z)

lf
} else {
sum(lf)
}
}

loglik.ccc <- function(param, data, model){
if(is.zoo(data)) data <- as.matrix(data)
nobs <- dim(data)[1]
ndim <- dim(data)[2]
mu <- matrix(param[1:ndim], nobs, ndim, byrow = TRUE)    # constant in the mean
data <- data - mu
param <- param[-(1:ndim)]

para.mat <- p.mat(param, model, ndim)
# check if R is positive definite. If not, making R positive definite
R <- make.pd(para.mat\$R)

h <- vgarch(para.mat\$a, para.mat\$A, para.mat\$B, data)    # a call to vgarch function
z <- data/sqrt(h)
invR <- solve(R)
lf <- -0.5*nobs*ndim*log(2*pi) - 0.5*sum(log(h)) - 0.5*nobs*log(det(R)) - 0.5*sum((z%*%invR)*z)
-lf
}

# making a symmetric matrix positive definite. From r-help 2003.12.27
# except normalization.
make.pd <- function(x, tol=1e-6) {
eig <- eigen(x, symmetric=TRUE)
rtol <- tol * eig\$values[1]
if(min(eig\$values) < rtol) {
vals <- eig\$values
vals[vals < rtol] <- rtol
srev <- eig\$vectors %*% (vals * t(eig\$vectors))
dimnames(srev) <- dimnames(x)
srev <- diag(1/sqrt(diag(srev)))%*%srev%*%diag(1/sqrt(diag(srev)))
return(srev)
} else {return(x)}
}

# Inequality constraints for the stationarity of the vector GARCH equation in the (E)CCC GARCH
inEQccc <- function(param, data, model){
ndim <- ncol(data)
param <- param[-(1:ndim)]   # removing constants in the mean
pmat <- p.mat(param, model, ndim)
ret <- max(Mod(eigen(pmat\$A + pmat\$B)\$values))     # common to diagonal/extended
ret
}

#*****************************************************************************************************************
estimateCCC <- function(inia = NULL, iniA = NULL, iniB = NULL, data, model="diagonal", ...){
nobs <- dim(data)[1]
ndim <- dim(data)[2]

if(is.zoo(data)){
d.ind <- index(data)
} else {
d.ind <- 1:nobs
}

data <- as.matrix(data)
mu <- colMeans(data)
R <- cor(data)

# setting upper and lower bounds for the constraints
if(model == "diagonal"){
LB <- c(rep(-Inf, ndim), rep(0, 3*ndim), rep(-1, ndim*(ndim-1)/2))
} else {
LB <- c(rep(-Inf, ndim), rep(0, ndim + 2*ndim^2), rep(-1, ndim*(ndim-1)/2))
}
ineqLB <- 0
ineqUB <- 1

if(!is.null(inia) & !is.null(iniA) & !is.null(iniB)){ # when the initial values are supplied
if(model=="diagonal"){
init <- c(mu, inia, diag(iniA), diag(iniB), R[lower.tri(R)])
} else {
init <- c(mu, inia, as.vector(iniA), as.vector(iniB), R[lower.tri(R)])
}
tryCatch(
suppressWarnings(
results <- solnp(pars = init, fun = loglik.ccc,
ineqfun = inEQccc, ineqLB = ineqLB, ineqUB = ineqUB,
LB = LB,
control = list(trace=0),
data = data, model = model
)
),
error = function(e) conditionMessage(e),
finally=cat("Bad initial values. Try again without initial values.")
)
} else { # when initial values are not supplied
cat("Initial values are not supplied. Random values are used.")
# first stage optimisation
conv <- 1
ntry <- 0
inia <- diag(cov(data))     # initial values for the constants in the GARCH part
while(conv != 0){ # repeating the first stage optimization until it gives a successful convergence
if(model != "diagonal"){
ret <- 2
while(ret > 1){
ub <- runif(ndim, min=0.0001, max=0.04)
iniA <- matrix(runif(ndim^2, min=0, max=ub[sample(1:ndim, 1)]), ndim, ndim)
iniB <- matrix(runif(ndim^2, min=-0.04, max=ub[sample(1:ndim, 1)]), ndim, ndim)
diag(iniA) <- round(runif(ndim, min=0.04, max=0.05), 4)
diag(iniB) <- round(runif(ndim, min=0.8, max=0.9), 4)
init <- c(mu, inia, as.vector(iniA), as.vector(iniB), R[lower.tri(R)])
ret <- max(Mod(eigen(iniA + iniB)\$values))     # common to diagonal/extended
}
} else {
iniA <- diag(round(runif(ndim, min=0.04, max=0.05), 4))
iniB <- diag(round(runif(ndim, min=0.8, max=0.9), 4))
init <- c(mu, inia, diag(iniA), diag(iniB), R[lower.tri(R)])
}
suppressWarnings(
results <- solnp(pars = init, fun = loglik.ccc,
ineqfun = inEQccc, ineqLB = ineqLB, ineqUB = ineqUB,
LB = LB,
control = list(trace=0),
data = data, model = model
)
)
conv <- results\$convergence
ntry <- ntry + 1
}
}

mu <- matrix(results\$pars[1:ndim], nobs, ndim, byrow = TRUE)
eps <- data - mu

param <- results\$pars[-(1:ndim)]
estimates <- p.mat(param, model=model, ndim=ndim)
estimates\$R <- make.pd(estimates\$R)

# computing conditional variance and std. residuals
h <- vgarch(estimates\$a, estimates\$A, estimates\$B, eps)    # estimating conditional variances
std.resid <- eps/sqrt(h)                    # std. residuals

# A character vector/matrix for naming parameters
name.id <- as.character(1:ndim)

if(is.null(colnames(data))){
colnames(std.resid) <- paste("Series", name.id, sep="")
colnames(h) <- paste("Series", name.id, sep="")
colnames(data) <- paste("Series", name.id, sep="")
} else {
colnames(std.resid) <- colnames(data)
colnames(h) <- colnames(data)
}

output <- list(
results = results,
model = model,
method = "SQP by Rsolnp package",
initial = list(a=inia, A=iniA, B=iniB),
data = zoo(data, d.ind),
estimates = estimates,
h = zoo(h, d.ind),              # conditional variances (not volatility)
z = zoo(std.resid, d.ind)       # standardized residuals
)

class(output) <- "ccc"
return(output)
}

# functions for summarizing output
summary.ccc <- function(object, ...){
cat("Summarizing outcomes. This takes a while.")
ndim <- ncol(object\$data)
nobs <- nrow(object\$data)
object\$nobs <- nobs

param <- object\$results\$pars
object\$mu <- param[1:ndim]             # mean estimates (constant)
names(object\$mu) <- paste("mu", 1:ndim, sep="")
param <- param[-(1:ndim)]

# re-arranging parameter vector into a list with paramerer matrices
para.mat <- p.mat(param, object\$model, ndim)
para.mat\$R <- make.pd(para.mat\$R)

# A character vector/matrix for naming parameters
name.id <- as.character(1:ndim)
namev <- diag(0, ndim, ndim)
for(i in 1:ndim){
for(j in 1:ndim){
namev[i, j] <- paste(name.id[i], name.id[j], sep="")
}
}
# naming parameters
if(object\$model=="diagonal"){
vecA <- diag(para.mat\$A)
vecB <- diag(para.mat\$B)
names(para.mat\$a) <- paste("a", 1:ndim, sep="")
names(vecA) <- paste("A", diag(namev), sep="")
names(vecB) <- paste("B", diag(namev), sep="")
} else {
vecA <- as.vector(para.mat\$A)
vecB <- as.vector(para.mat\$B)
names(para.mat\$a) <-paste("a", 1:ndim, sep="")
names(vecA) <-paste("A", namev, sep="")
names(vecB) <-paste("B", namev, sep="")
}
object\$R <- para.mat\$R[lower.tri(para.mat\$R)]
names(object\$R) <- paste("R", namev[lower.tri(namev)], sep="")
object\$garch.par <- c(para.mat\$a, vecA, vecB)  # estimates for conditional variance

# computing Jacobian and Hessian for the mean and GARCH part
all.par <- c(object\$mu, object\$garch.par, object\$R)
npar <- length(all.par)       # the number of total parameters
ja <- jacobian(func=loglik.ccc.t, x=all.par,
data=object\$data, model=object\$model, mode="gradient")  # using jacobian() in numDeriv
J <- crossprod(ja)  # information matrix
H <- hessian(func = loglik.ccc.t, x=all.par,
data=object\$data, model=object\$model, mode="hessian", method.args=list(eps=1e-5, d=1e-5))
invH <- solve(H)
S <- invH%*%J%*%invH
se <- sqrt(diag(S))

coef <- c(object\$mu, object\$garch.par, object\$R)
zstat <- all.par/se
pval <- round(2*pnorm(-abs(zstat)), 6)
coef <- cbind(all.par, se, zstat, pval)
colnames(coef) <- c("Estimate", "Std. Error", "z value", "Pr(>|z|)")

object\$coef <- coef

object\$convergence <- object\$results\$convergence    # convergence status
names(object\$convergence) <- "1st"
object\$counts <- c(object\$results\$outer.iter, object\$results\$nfuneval)   # the number of iterations

object\$logLik <- -tail(object\$results\$values, 1)  # the value of the log-like at the estimate

# Information Criteria
object\$AIC <- -2*object\$logLik + 2*npar
object\$BIC <- -2*object\$logLik + log(nobs)*npar
object\$CAIC <- -2*object\$logLik + (1+log(nobs))*npar

class(object) <- "summary.ccc"
return(object)
}

logLik.ccc <- function(object, ...){
LL <- -tail(object\$results\$values, 1)
cat("Log-likelihood at the estimates: ", formatC(LL, digits = 10), sep = "\n")
}

print.summary.ccc <- function(x, digits = max(3, getOption("digits") - 1), ...){
cat("Constant Conditional Correlation GARCH Model", "\n")
cat("Conditional variance equation:", x\$model, "\n")
cat("\nCoefficients:", "\n")
printCoefmat(x\$coef, digits = 4, dig.tst = 4)    # printCoefmat() is defined in stats package

cat("\nNumber of Obs.:", formatC(x\$nobs, digits = 0), "\n")
cat("Log-likelihood:", formatC(x\$logLik, format="f", digits = 5), "\n\n")

cat("Information Criteria:", "\n")
cat("  AIC:", formatC(x\$AIC, format="f", digits = 3), "\n")
cat("  BIC:", formatC(x\$BIC, format="f", digits = 3), "\n")
cat(" CAIC:", formatC(x\$CAIC, format="f", digits = 3), "\n")

cat("\nOptimization method:", x\$method, "\n")
cat("Convergence :", x\$convergence, "\n")
cat("Iterations  :", "\n")
cat("  Major iterations :", x\$counts[1], "\n")
cat(" Func. evaluations :", x\$counts[2], "\n")

invisible(x)
}

print.ccc <- function(x, ...){
cat("Estimated model:", x\$model, "\n")
cat("Use summary() to see the estimates and related statistics", "\n")
invisible(x)
}
```

## Try the ccgarch2 package in your browser

Any scripts or data that you put into this service are public.

ccgarch2 documentation built on May 31, 2017, 4:23 a.m.