ccgarch源程序代码

function (a, A, B, R, u, model = "diagonal")

{

invR <- solve(R)

nobs <- dim(u)[1]

ndim <- dim(u)[2]

npar.h <- ndim * (2 * ndim + 1)

npar.r <- ndim * (ndim - 1)/2

In <- diag(ndim)

vecIn <- as.vector(In)

h <- vector.garch(u, a, A, B)

sq.invh <- 1/sqrt(h)

norm.z <- u * sq.invh

dhw <- vec.garch.derivative(u, B, h)

for (i in 1:ndim) {

dhw[, ((i - 1) * npar.h + 1):(i * npar.h)] <- sq.invh[,

i] * dhw[, ((i - 1) * npar.h + 1):(i * npar.h)] }

dhw <- t(0.5 * dhw)

if (model == "diagonal") {

ind <- as.vector(rbind(1, In, In))

dhw <- dhw[ind == 1, ]

npar.h <- 3 * ndim

}

vecdR <- vdR(ndim)

D <- matrix(0, npar.h, nobs)

Q <- matrix(0, npar.r, nobs)

for (i in 1:nobs) {

invD <- diag(sq.invh[i, ])

Z <- norm.z[i, ] %o% norm.z[i, ]

K <- Z %*% invR %*% invD

d <- matrix(dhw[, i], npar.h, ndim)

W <- -0.5 * as.vector(2 * invD - (K + t(K)))

W <- W[vecIn == 1]

D[, i] <- d %*% W

Q[, i] <- -0.5 * vecdR %*% as.vector(invR - invR %*% Z %*% invR)

}

rbind(D, Q)

}

function (a, A, B, R, u, model)

{

nobs <- dim(u)[1]

ndim <- dim(u)[2]

npar.h <- ndim * (2 * ndim + 1)

npar.r <- ndim * (ndim - 1)/2

In <- diag(ndim)

vecIn <- as.vector(In)

invR <- solve(R)

h <- vector.garch(u, a, A, B)

sq.h <- sqrt(h)

sq.invh <- 1/sq.h

dhw <- vec.garch.derivative(u, B, h)

for (i in 1:ndim) {

dhw[, ((i - 1) * npar.h + 1):(i * npar.h)] <- sq.invh[,

i] * dhw[, ((i - 1) * npar.h + 1):(i * npar.h)]

}

dhw <- t(0.5 * dhw)

if (model == "diagonal") {

ind <- c(rep(1, ndim), as.vector(In), as.vector(In))

dhw <- dhw[ind == 1, ]

npar.h <- 3 * ndim

}

vecdR <- vdR(ndim)

H11 <- matrix(0, npar.h, npar.h)

H21 <- matrix(0, npar.r, npar.h)

for (i in 1:nobs) {

D <- diag(sq.h[i, ])

invD <- diag(sq.invh[i, ])

invH <- invD %*% invR %*% invD

dh <- matrix(dhw[, i], npar.h, ndim)

H11.temp <- 2 * (invD %x% invD) + invH %x% R + R %x%

invH

H11.temp <- H11.temp[vecIn == 1, vecIn == 1]

H11 <- H11 + dh %*% H11.temp %*% t(dh)

RD <- invR %*% invD

H21.temp <- vecdR %*% (RD %x% In + In %x% RD)

H21.temp <- H21.temp[, vecIn == 1]

H21 <- H21 + H21.temp %*% t(dh)

}

H22 <- nobs * vecdR %*% (invR %x% invR) %*% t(vecdR)

0.5 * cbind(rbind(H11, H21), rbind(t(H21), H22))

}

function (u, B, h, model)

{

nobs <- dim(u)[1]

ndim <- dim(u)[2]

In <- diag(ndim)

ind <- as.vector(rbind(1, In, In))

dhw <- vec.garch.derivative(u, B, h)

inv.h.sq <- 1/h^2

if (model == "diagonal") {

npar.h <- 3 * ndim

}

else {

npar.h <- ndim * (2 * ndim + 1)

}

dlv2 <- matrix(0, nobs, npar.h^2)

for (i in 1:nobs) {

dhwt <- matrix(dhw[i, ], ncol = ndim)

if (model == "diagonal") {

dhwt <- dhwt[ind == 1, ]

}

tmp.dlv2 <- matrix(0, npar.h, npar.h)

for (j in 1:ndim) {

tmp.dlv2 <- tmp.dlv2 + inv.h.sq[i, j] * outer(dhwt[,

j], dhwt[, j])

}

dlv2[i, ] <- as.vector(tmp.dlv2)

}

0.5 * dlv2

}

function (dvar, param)

{

uncR <- cov(dvar)

out <- .Call("dcc_est", dvar, uncR, param[1], param[2]) list(DCC = out[[1]], Q = out[[2]])

}

> dcc.estimation

function (inia, iniA, iniB, ini.dcc, dvar, model, method = "BFGS",

gradient = 1, message = 1)

{

dvar <- as.matrix(dvar)

nobs <- dim(dvar)[1]

ndim <- dim(dvar)[2]

In <- diag(ndim)

if (!is.matrix(iniA) || !is.matrix(iniB)) {

stop("iniA or iniB or both must be matrices")

}

if (length(inia) != ndim || dim(iniA)[2] != ndim || dim(iniB)[2] !=

ndim) {

stop("the dimension is incorrect in inia, iniA or iniB.")

}

first.stage <- dcc.estimation1(dvar = dvar, a = inia, A = iniA,

B = iniB, model = model, method = method)

if (first.stage$convergence != 0) {

cat("***********************************************\n")

cat("* The first stage optimization has failed. *\n")

cat("* See the list variable 'second' for details. *\n")

cat("***********************************************\n")

}

tmp.para <- c(first.stage$par, In[lower.tri(In)])

estimates <- p.mat(tmp.para, model = model, ndim = ndim)

esta <- estimates$a

estA <- estimates$A

estB <- estimates$B

h <- vector.garch(dvar, esta, estA, estB)

std.resid <- dvar/sqrt(h)

second.stage <- dcc.estimation2(std.resid, ini.dcc, gradient = gradient)

if (second.stage$convergence != 0) {

cat("***********************************************\n")

cat("* The second stage optimization has failed. *\n")

cat("* See the list variable 'second' for details. *\n")

cat("***********************************************\n")

}

else if (message != 0) {

cat("****************************************************************\n")

cat("* Estimation has been completed. *\n")

cat("* The outputs are saved in a list with components: *\n")

cat("* out : the estimates and their standard errors *\n")

cat("* loglik : the value of the log-likelihood at the estimates *\n")

cat("* h : a matrix of estimated conditional variances *\n")

cat("* DCC : a matrix of DCC estimates *\n")

cat("* std.resid : a matrix of the standardised residuals *\n")

cat("* first : the results of the first stage estimation *\n")

cat("* second : the results of the second stage estimation *\n")

cat("****************************************************************\n") }

dcc <- dcc.est(std.resid, second.stage$par)$DCC

lf1 <- -0.5 * ndim * log(2 * pi) - rowSums(log(h))

lf2 <- numeric(nobs)

for (i in 1:nobs) {

R <- matrix(dcc[i, ], ndim, ndim)

invR <- solve(R)

lf2[i] <- -0.5 * (log(det(R)) + sum(std.resid[i, ] *

crossprod(invR, std.resid[i, ])))

}

loglik <- sum(lf1 + lf2)

output <- dcc.results(dvar, first.stage$par, second.stage$par,

h, model = model)

list(out = output, loglik = loglik, h = h, DCC = dcc, std.resid = std.resid,

first = first.stage, second = second.stage)

}

function (dvar, a, A, B, model, method = "BFGS")

{

nobs <- dim(dvar)[1]

ndim <- dim(dvar)[2]

if (model == "diagonal") {

init <- sqrt(c(a, diag(A), diag(B)))

}

else {

init <- sqrt(c(a, as.vector(A), as.vector(B)))

}

step1 <- optim(par = init, fn = loglik.dcc1, method = method, control = list(maxit = 10^5, reltol = 1e-05), dvar = dvar,

model = model)

step1

}

> dcc.estimation2

function (dvar, para, gradient = 0)

{

resta <- rbind(c(-1, -1), diag(2))

restb <- c(-1, 0, 0)

if (gradient != 0) {

step2 <- constrOptim(theta = para, f = loglik.dcc2, gr = grad.dcc2, ui = resta, ci = restb, mu = 1e-05, dvar = dvar)

}

else {

step2 <- constrOptim(theta = para, f = loglik.dcc2, gr = NULL,

ui = resta, ci = restb, mu = 1e-05, dvar = dvar)

}

step2

}

> dcc.results

function (u, garch.para, dcc.para, h, model)

{

nobs <- dim(u)[1]

ndim <- dim(u)[2]

In <- diag(ndim)

param <- c(garch.para, In[lower.tri(In)])

estimates <- p.mat(param, model, ndim)

a <- estimates$a

A <- estimates$A

B <- estimates$B

if (model == "diagonal") {

parameters <- c(a, diag(A), diag(B), dcc.para)

}

else {

parameters <- c(a, as.vector(A), as.vector(B), dcc.para) }

d2lvdw <- d2lv(u, B, h, model = model)

dlvdw <- dlv(u, a, A, B, model = model)

O <- dlc(dcc.para, B, u, h, model = model)

dlcdf <- O$dlc

d2lcdf <- O$d2lc

d2lcdfdw <- O$dfdwd2lc

G11 <- matrix(colMeans(d2lvdw), ncol = sqrt(dim(d2lvdw)[2]))

G22 <- matrix(colMeans(d2lcdf), 2, 2)

G21 <- matrix(colMeans(d2lcdfdw), nrow = 2)

O11 <- dlvdw %*% t(dlvdw)/nobs

O22 <- t(dlcdf) %*% dlcdf/nobs

O12 <- dlvdw %*% dlcdf/nobs

G. <- cbind(rbind(G11, G21), rbind(matrix(0, nrow = dim(G11)[1],

ncol = 2), G22))

O. <- cbind(rbind(O11, t(O12)), rbind(O12, O22))

invG. <- solve(G., tol = 1e-50)

V <- invG. %*% O. %*% t(invG.)/nobs

se <- sqrt(diag(V))

results <- rbind(parameters, se)

name.a <- numeric(0)

name.A <- numeric(0)

name.B <- numeric(0)

for (i in 1:ndim) {

name.a <- c(name.a, paste("a", paste(i), sep = ""))

for (j in 1:ndim) {

name.A <- c(name.A, paste("A", paste(paste(j), i,

sep = ""), sep = ""))

name.B <- c(name.B, paste("B", paste(paste(j), i,

sep = ""), sep = ""))

}

}

if (model == "diagonal") {

ind <- as.vector(diag(ndim))

name.A <- name.A[ind == 1]

name.B <- name.B[ind == 1]

}

colnames(results) <- c(name.a, name.A, name.B, "dcc alpha", "dcc beta")

rownames(results) <- c("estimates", "std.err")

results

}

> dcc.sim

function (nobs, a, A, B, R, dcc.para, d.f = Inf, cut = 1000,

model)

{

if (model == "diagonal") {

A <- diag(diag(A))

B <- diag(diag(B))

}

if (stationarity(A, B) > 1) {

stop("A and B matrices do not satisfy the stationarity condition \n") }

nobs <- nobs + cut

Id <- diag(length(a))

inih <- solve(Id - A - B) %*% a

dccpar1 <- dcc.para[1]

dccpar2 <- dcc.para[2]

sim <- .Call("dcc_sim", nobs, a, A, B, inih, R, dccpar1,

dccpar2, d.f)

list(z = sim[[1]][(cut + 1):(nobs), ], std.z = sim[[2]][(cut +

1):(nobs), ], dcc = sim[[3]][(cut + 1):(nobs), ], h = sim[[4]][(cut +

1):(nobs), ], eps = sim[[5]][(cut + 1):(nobs), ])

}

> dlc

function (dcc.para, B, u, h, model)

{

nobs <- dim(h)[1]

ndim <- dim(h)[2]

dccpar1 <- dcc.para[1]

dccpar2 <- dcc.para[2]

hlag <- rbind(colMeans(h), h)

dhw <- vec.garch.derivative(u, B, h)

dhwlag <- rbind(colMeans(dhw), dhw)

inv.sq.h <- 1/sqrt(h)

https://www.360docs.net/doc/7c6691715.html,g <- rbind(colMeans(inv.sq.h), inv.sq.h)

z <- u * inv.sq.h

zlag <- rbind(colMeans(z), z)

Q <- cov(z)

dcc <- dcc.est(z, dcc.para)

P <- dcc$DCC

Plag <- rbind(colMeans(P), P)

Qt <- dcc$Q

Qtlag <- rbind(colMeans(Qt), Qt)

In <- diag(ndim)

ind <- as.vector(rbind(1, In, In))

if (model == "diagonal") {

npar.h <- 3 * ndim

}

else {

npar.h <- ndim * (2 * ndim + 1)

}

Qtilde <- 1/sqrt(Qt[, as.vector(In) == 1])

dvecQ <- matrix(0, nobs, 2 * ndim^2)

dvecP <- matrix(0, nobs, 2 * ndim^2)

dlc <- matrix(0, nobs, 2)

d2lc <- matrix(0, nobs, 4)

const <- -rbind(as.vector(Q), as.vector(Q))

dvecQt <- matrix(0, 2, ndim^2)

dwdvecQt <- matrix(0, npar.h, ndim^2)

dwdvecQ <- matrix(0, nobs, npar.h * ndim^2)

dfdwd2lc <- matrix(0, nobs, 2 * npar.h)

for (i in 1:nobs) {

dvecQt <- const + rbind(as.vector(outer(zlag[i, ], zlag[i, ])), Qtlag[i, ]) + dccpar2 * dvecQt

dvecQ[i, ] <- as.vector(dvecQt)

invQtilde <- diag(1/Qtilde[i, ])

Pt <- matrix(P[i, ], ndim, ndim)

invPt <- solve(Pt)

Z <- Pt %x% invQtilde

tZ <- invQtilde %x% Pt

QPQ <- invQtilde %x% invQtilde - 0.5 * diag(as.vector(invQtilde)) %*% (Z + tZ)

dvecPt <- dvecQt %*% QPQ

dvecPt[, as.vector(In) == 1] <- 0

dvecP[i, ] <- as.vector(dvecPt)

zz <- outer(z[i, ], z[i, ])

dlc[i, ] <- -0.5 * as.vector(dvecPt %*% as.vector(invPt -

invPt %*% zz %*% invPt))

d2lc[i, ] <- 0.5 * as.vector(dvecPt %*% (invPt %x% invPt) %*%

t(dvecPt))

dhwtlag <- matrix(dhwlag[i, ], ncol = ndim)

dhwt <- matrix(dhw[i, ], ncol = ndim)

if (model == "diagonal") {

dhwtlag <- dhwtlag[ind == 1, ]

dhwt <- dhwt[ind == 1, ]

}

dwdvecVtlag <- matrix(0, npar.h, ndim^2)

dwdvecVt <- matrix(0, npar.h, ndim^2)

dwdvecVtlag[, as.vector(In) == 1] <- dhwtlag

dwdvecVt[, as.vector(In) == 1] <- dhwt

Dlag <- diag(hlag[i, ])

Dt <- diag(sqrt(h[i, ]))

tmpDIn <- diag(Dt %x% In + In %x% Dt)

tmpDInlag <- diag(Dlag %x% In + In %x% Dlag)

dwdvecDlag <- dwdvecVtlag %*% diag(1/tmpDInlag)

dwdvecD <- dwdvecVtlag %*% diag(1/tmpDIn)

invD <- diag(inv.sq.h[i, ])

Ptlag <- matrix(Plag[i, ], ndim, ndim)

invDlag <- diag(https://www.360docs.net/doc/7c6691715.html,g[i, ])

ZD <- Ptlag %x% invDlag

DZ <- invDlag %x% Ptlag

dwdvecQt <- -0.5 * dccpar1 * dwdvecDlag %*% (ZD + DZ) +

dccpar2 * dwdvecQt

dwdvecQ[i, ] <- as.vector(dwdvecQt)

dwdvecPt <- t(dwdvecQt %*% QPQ)

dfdwd2lc[i, ] <- as.vector(0.5 * dvecPt %*% (invPt %x%

invPt) %*% dwdvecPt + 0.5 * dvecPt %*% ((invPt %*%

invD) %x% In + In %x% (invPt %*% invD)) %*% t(dwdvecD))

}

list(dlc = dlc, dvecP = dvecP, dvecQ = dvecQ, d2lc = d2lc, dfdwd2lc = dfdwd2lc) }

function (u, a, A, B, model)

{

ndim <- dim(u)[2]

dl <- analytical.grad(a, A, B, diag(ndim), u, model = model) dl.row <- dim(dl)[1] - ndim * (ndim - 1)/2

dl[1:dl.row, ]

}

function (par, dvar, model)

{

ndim <- dim(dvar)[2]

In <- diag(ndim)

par <- c(par, In[lower.tri(In)])

para <- p.mat(par, model, ndim)

dl <- analytical.grad(para$a, para$A, para$B, In, dvar, model = model) dl.row <- dim(dl)[1] - ndim * (ndim - 1)/2

rowSums(dl[1:dl.row, ])

}

> eccc.estimation

function (a, A, B, R, dvar, model, method = "BFGS")

{

dvar <- as.matrix(dvar)

nobs <- dim(dvar)[1]

ndim <- dim(dvar)[2]

if (!is.matrix(A) || !is.matrix(B)) {

stop("A or B or both must be matrices")

}

if (model == "diagonal") {

init <- c(sqrt(a), diag(sqrt(A)), diag(sqrt(B)), R[lower.tri(R)])

}

else {

init <- c(sqrt(a), as.vector(sqrt(A)), as.vector(sqrt(B)),

R[lower.tri(R)])

}

results <- optim(par = init, fn = loglik.eccc, method = method,

dvar = dvar, model = model, control = list(maxit = 10^5,

reltol = 1e-15))

if (results$convergence != 0) {

cat("***********************************************************\n")

cat("* Optimization is FAILED. *\n")

stop("* Fine tuning is required. *\n")

cat("***********************************************************\n") }

estimates <- p.mat(results$par, model = model, ndim = ndim)

esta <- estimates$a

estA <- estimates$A

estB <- estimates$B

estR <- estimates$R

h <- vector.garch(dvar, esta, estA, estB)

std.resid <- dvar/sqrt(h)

grad <- analytical.grad(a = esta, A = estA, B = estB, R = estR,

u = dvar, model = model)

grad <- grad %*% t(grad)/nobs

H <- analytical.Hessian(a = esta, A = estA, B = estB, R = estR,

u = dvar, model = model)/nobs

invH <- solve(H)

Avar <- (invH %*% grad %*% invH)/nobs

rob.se <- sqrt(diag(Avar))

out.se <- sqrt(diag(solve(grad)/nobs))

H.se <- sqrt(diag(invH/nobs))

std.error <- rbind(H.se, out.se, rob.se)

rownames(std.error) <- c("Inv. Hessian", "Outer Prod.", "Robust")

name.a <- numeric(0)

name.A <- numeric(0)

name.B <- numeric(0)

name.R <- numeric(0)

for (i in 1:ndim) {

name.a <- c(name.a, paste("a", paste(i), sep = ""))

for (j in 1:ndim) {

name.A <- c(name.A, paste("A", paste(paste(j), i,

sep = ""), sep = ""))

name.B <- c(name.B, paste("B", paste(paste(j), i,

sep = ""), sep = ""))

name.R <- c(name.R, paste("R", paste(paste(j), i,

sep = ""), sep = ""))

}

}

names(esta) <- name.a

if (model == "diagonal") {

ind <- as.vector(diag(ndim))

vecA <- diag(estA)

name.A <- name.A[ind == 1]

vecB <- diag(estB)

name.B <- name.B[ind == 1]

}

else {

vecA <- as.vector(estA)

vecB <- as.vector(estB)

}

names(vecA) <- name.A

names(vecB) <- name.B

vecR <- estR[lower.tri(estR)]

name.R <- matrix(name.R, ndim, ndim)

name.R <- name.R[lower.tri(name.R)]

names(vecR) <- name.R

colnames(std.error) <- c(name.a, name.A, name.B, name.R)

para.estimates <- c(esta, vecA, vecB, vecR)

list(out = rbind(para.estimates, std.error), h = h, std.resid = std.resid, opt = results, para.mat = estimates)

}

function (nobs, a, A, B, R, d.f = Inf, cut = 1000, model)

{

if (model == "diagonal") {

A <- diag(diag(A))

B <- diag(diag(B))

}

if (stationarity(A, B) > 1) {

stop("A and B matrices do not satisfy the stationarity condition. \n") }

nobs <- nobs + cut

Id <- diag(length(a))

inih <- solve(Id - A - B) %*% a

sim <- .Call("eccc_sim", nobs, a, A, B, R, inih, d.f)

list(h = sim[[1]][(cut + 1):(nobs), ], eps = sim[[2]][(cut +

1):(nobs), ])

}

function (A, B, R)

{

AB <- A + B

AA <- A %x% A

G <- AB %x% AB + 2 * AA * (diag(as.vector(R))^2)

max(Mod(eigen(G)$values))

}

function (a, A, B, dcc.para, dvar, d = 1e-05, model)

{

if (model == "diagonal") {

param <- c(a, diag(A), diag(B), dcc.para) }

else {

param <- c(a, as.vector(A), as.vector(B), dcc.para) }

nobs <- dim(dvar)[1]

ndim <- dim(dvar)[2]

npara <- length(param)

Id <- d * diag(npara)

param.d <- matrix(param, npara, npara) + Id

lf <- loglik.dcc(param, dvar, model)

lf.d <- matrix(0, nobs, npara)

for (i in 1:npara) {

lf.d[, i] <- loglik.dcc(param.d[, i], dvar, model) }

(lf.d - lf)/d

}