1
I have BEKKK11 function of MTS package
function (rt, include.mean = T, cond.dist = "normal", ini.estimates = NULL)
{
if (!is.matrix(rt))
rt = as.matrix(rt)
nT = dim(rt)[1]
k = dim(rt)[2]
if (k > 3) {
cat("Program Note: Dimension is limited to 3",
"\n")
k = 3
}
RTN <- rt[, 1:k]
mu = apply(RTN, 2, mean)
Cov1 <- cov(RTN)
m1 = chol(Cov1)
S = 1e-06
S1 = -0.5
if (k == 2) {
if (is.null(ini.estimates)) {
A1 = matrix(c(0.1, 0.02, 0.02, 0.1), k, k)
B1 = matrix(c(0.9, 0.051, 0.01, 0.9), k, k)
}
else {
ist = 0
if (include.mean) {
mu = ini.estimates[1:2]
ist = 2
}
m1[1, 1] = ini.estimates[ist + 1]
m1[1, 2] = ini.estimates[ist + 2]
m1[2, 2] = ini.estimates[ist + 3]
ist = ist + 3
A1 = matrix(ini.estimates[(ist + 1):(ist + 4)], k,
k)
ist = ist + 4
B1 = matrix(ini.estimates[(ist + 1):(ist + 4)], k,
k)
}
if (include.mean) {
par = c(mu1 = mu[1], mu2 = mu[2], A011 = m1[1, 1],
A021 = m1[1, 2], A022 = m1[2, 2], A11 = A1[1,
1], A21 = A1[2, 1], A12 = A1[1, 2], A22 = A1[2,
2], B11 = B1[1, 1], B21 = B1[2, 1], B12 = B1[1,
2], B22 = B1[2, 2])
c1 = c(mu1 = -10 * abs(mu[1]), mu2 = -10 * abs(mu[2]),
A011 = m1[1, 1] * 0.2, A021 = m1[1, 2] * 0.2,
A022 = m1[2, 2] * 0.2, A11 = S, A21 = S1, A12 = S1,
A22 = S, B11 = S, B21 = S1, B12 = S1, B22 = S)
c2 = c(mu1 = 10 * abs(mu[1]), mu2 = 10 * abs(mu[2]),
A011 = m1[1, 1] * 1.1, A021 = m1[1, 2] * 1.1,
A022 = m1[2, 2] * 1.1, A11 = 1 - S, A21 = -S1,
A12 = -S1, A22 = 1 - S, B11 = 1 - S, B21 = -S1,
B12 = -S1, B22 = 1 - S)
}
else {
par = c(A011 = m1[1, 1], A021 = m1[1, 2], A022 = m1[2,
2], A11 = A1[1, 1], A21 = A1[2, 1], A12 = A1[1,
2], A22 = A1[2, 2], B11 = B1[1, 1], B21 = B1[2,
1], B12 = B1[1, 2], B22 = B1[2, 2])
c1 = c(A011 = m1[1, 1] * 0.2, A021 = m1[1, 2] * 0.2,
A022 = m1[2, 2] * 0.2, A11 = S, A21 = S1, A12 = S1,
A22 = S, B11 = S, B21 = S1, B12 = S1, B22 = S)
c2 = c(A011 = m1[1, 1] * 1.1, A021 = m1[1, 2] * 1.1,
A022 = m1[2, 2] * 1.1, A11 = 1 - S, A21 = -S1,
A12 = -S1, A22 = 1 - S, B11 = 1 - S, B21 = -S1,
B12 = -S1, B22 = 1 - S)
}
}
if (k == 3) {
if (is.null(ini.estimates)) {
A1 = matrix(c(0.1, 0.02, 0.02, 0.02, 0.1, 0.02, 0.02,
0.02, 0.1), k, k)
B1 = matrix(c(0.8, 0.02, 0.02, 0.02, 0.8, 0.02, 0.02,
0.02, 0.8), k, k)
}
else {
ist = 0
if (include.mean) {
mu = ini.estimates[1:k]
ist = 3
}
m1[1, 1] = ini.estimates[ist + 1]
m1[1, 2] = ini.estimates[ist + 2]
m1[1, 3] = ini.estimates[ist + 3]
m1[2, 2] = ini.estimates[ist + 4]
m1[2, 3] = ini.estimates[ist + 5]
m1[3, 3] = ini.estimates[ist + 6]
ist = ist + 6
A1 = matrix(ini.estimates[(ist + 1):(ist + 9)], k,
k)
ist = ist + 9
B1 = matrix(ini.estimates[(ist + 1):(ist + 9)], k,
k)
}
if (include.mean) {
par = c(mu1 = mu[1], mu2 = mu[2], mu3 = mu[3], A011 = m1[1,
1], A021 = m1[1, 2], A031 = m1[1, 3], A022 = m1[2,
2], A032 = m1[2, 3], A033 = m1[3, 3], A11 = A1[1,
1], A21 = A1[2, 1], A31 = A1[3, 1], A12 = A1[1,
2], A22 = A1[2, 2], A32 = A1[3, 2], A13 = A1[1,
3], A23 = A1[2, 3], A33 = A1[3, 3], B11 = B1[1,
1], B21 = B1[2, 1], B31 = B1[3, 1], B12 = B1[1,
2], B22 = B1[2, 2], B32 = B1[3, 2], B13 = B1[1,
3], B23 = B1[2, 3], B33 = B1[3, 3])
c1 = c(mu1 = -10 * abs(mu[1]), mu2 = -10 * abs(mu[2]),
mu3 = -10 * abs(mu[3]), A011 = m1[1, 1] * 0.2,
A021 = m1[1, 2] * 0.2, A031 = m1[1, 3] * 0.2,
A022 = m1[2, 2] * 0.2, A032 = m1[2, 3] * 0.2,
A033 = m1[3, 3] * 0.2, A11 = S, A21 = S1, A31 = S1,
A12 = S1, A22 = S, A32 = S1, A13 = S1, A23 = S1,
A33 = S, B11 = S, B21 = S1, B31 = S1, B12 = S1,
B22 = S, B32 = S1, B13 = S1, B23 = S1, B33 = S)
c2 = c(mu1 = 10 * abs(mu[1]), mu2 = 10 * abs(mu[2]),
mu3 = 10 * abs(mu[3]), A011 = m1[1, 1] * 1.1,
A021 = m1[1, 2] * 1.1, A031 = m1[1, 3] * 1.1,
A022 = m1[2, 2] * 1.1, A032 = m1[2, 3] * 1.1,
A033 = m1[3, 3] * 1.1, A11 = 1 - S, A21 = -S1,
A31 = -S1, A12 = -S1, A22 = 1 - S, A32 = -S1,
A13 = -S1, A23 = -S1, A33 = 1 - S, B11 = 1 -
S, B21 = -S1, B31 = -S1, B12 = -S1, B22 = 1 -
S, B32 = -S1, B13 = -S1, B23 = -S1, B33 = 1 -
S)
}
else {
par = c(A011 = m1[1, 1], A021 = m1[1, 2], A031 = m1[1,
3], A022 = m1[2, 2], A032 = m1[2, 3], A033 = m1[3,
3], A11 = A1[1, 1], A21 = A1[2, 1], A31 = A1[3,
1], A12 = A1[1, 2], A22 = A1[2, 2], A32 = A1[3,
2], A13 = A1[1, 3], A23 = A1[2, 3], A33 = A1[3,
3], B11 = B1[1, 1], B21 = B1[2, 1], B31 = B1[3,
1], B12 = B1[1, 2], B22 = B1[2, 2], B32 = B1[3,
2], B13 = B1[1, 3], B23 = B1[2, 3], B33 = B1[3,
3])
c1 = c(A011 = m1[1, 1] * 0.2, A021 = m1[1, 2] * 0.2,
A031 = m1[1, 3] * 0.2, A022 = m1[2, 2] * 0.2,
A032 = m1[2, 3] * 0.2, A033 = m1[3, 3] * 0.2,
A11 = S, A21 = S1, A31 = S1, A12 = S1, A22 = S,
A32 = S1, A13 = S1, A23 = S1, A33 = S, B11 = S,
B21 = S1, B31 = S1, B12 = S1, B22 = S, B32 = S1,
B13 = S1, B23 = S1, B33 = S)
c2 = c(A011 = m1[1, 1] * 1.1, A021 = m1[1, 2] * 1.1,
A031 = m1[1, 3] * 1.1, A022 = m1[2, 2] * 1.1,
A032 = m1[2, 3] * 1.1, A033 = m1[3, 3] * 1.1,
A11 = 1 - S, A21 = -S1, A31 = -S1, A12 = -S1,
A22 = 1 - S, A32 = -S1, A13 = -S1, A23 = -S1,
A33 = 1 - S, B11 = 1 - S, B21 = -S1, B31 = -S1,
B12 = -S1, B22 = 1 - S, B32 = -S1, B13 = -S1,
B23 = -S1, B33 = 1 - S)
}
}
cat("Initial estimates: ", par, "\n")
cat("Lower limits: ", c1, "\n")
cat("Upper limits: ", c2, "\n")
A1at <- function(x) {
x = as.matrix(x)
Prod1 = NULL
for (i in 1:nrow(x)) {
Prod1 = rbind(Prod1, kronecker(x[i, ], x[i, ]))
}
Prod1
}
mlikeG <- function(par, RTN = RTN, include.mean = include.mean) {
nT = dim(RTN)[1]
k = dim(RTN)[2]
Cov1 = cov(RTN)
if (k == 2) {
if (include.mean) {
mu1 = par[1]
mu2 = par[2]
A011 = par[3]
A021 = par[4]
A022 = par[5]
A11 = par[6]
A21 = par[7]
A12 = par[8]
A22 = par[9]
B11 = par[10]
B21 = par[11]
B12 = par[12]
B22 = par[13]
}
else {
mu1 = 0
mu2 = 0
A011 = par[1]
A021 = par[2]
A022 = par[3]
A11 = par[4]
A21 = par[5]
A12 = par[6]
A22 = par[7]
B11 = par[8]
B21 = par[9]
B12 = par[10]
B22 = par[11]
}
A0 = matrix(c(A011, A021, 0, A022), 2, 2)
A0A0t = A0 %*% t(A0)
A1 = matrix(c(A11, A21, A12, A22), 2, 2)
B1 = matrix(c(B11, B21, B12, B22), 2, 2)
resi = cbind(RTN[, 1] - mu1, RTN[, 2] - mu2)
Sig = matrix(c(Cov1), 1, 4)
res = resi %*% t(A1)
ArchP = A1at(res)
}
if (k == 3) {
if (include.mean) {
mu1 = par[1]
mu2 = par[2]
mu3 = par[3]
A011 = par[4]
A021 = par[5]
A031 = par[6]
A022 = par[7]
A032 = par[8]
A033 = par[9]
A11 = par[10]
A21 = par[11]
A31 = par[12]
A12 = par[13]
A22 = par[14]
A32 = par[15]
A13 = par[16]
A23 = par[17]
A33 = par[18]
B11 = par[19]
B21 = par[20]
B31 = par[21]
B12 = par[22]
B22 = par[23]
B32 = par[24]
B13 = par[25]
B23 = par[26]
B33 = par[27]
}
else {
mu1 = 0
mu2 = 0
mu3 = 0
A011 = par[1]
A021 = par[2]
A031 = par[3]
A022 = par[4]
A032 = par[5]
A033 = par[6]
A11 = par[7]
A21 = par[8]
A31 = par[9]
A12 = par[10]
A22 = par[11]
A32 = par[12]
A13 = par[13]
A23 = par[14]
A33 = par[15]
B11 = par[16]
B21 = par[17]
B31 = par[18]
B12 = par[19]
B22 = par[20]
B32 = par[21]
B13 = par[22]
B23 = par[23]
B33 = par[24]
}
A0 = matrix(c(A011, A021, A031, 0, A022, A032, 0,
0, A033), k, k)
A0A0t = A0 %*% t(A0)
A1 = matrix(c(A11, A21, A31, A12, A22, A32, A13,
A23, A33), k, k)
B1 = matrix(c(B11, B21, B31, B12, B22, B32, B13,
B23, B33), k, k)
resi = cbind(RTN[, 1] - mu1, RTN[, 2] - mu2, RTN[,
3] - mu3)
Sig = matrix(c(Cov1), 1, k^2)
res = resi %*% t(A1)
ArchP = A1at(res)
}
llike = 0
for (t in 2:nT) {
Sigt = A0A0t + matrix(ArchP[t - 1, ], k, k) + B1 %*%
matrix(Sig[t - 1, ], k, k) %*% t(B1)
Sigt = (Sigt + t(Sigt))/2
Sig = rbind(Sig, c(Sigt))
d1 = dmvnorm(resi[t, ], mean = rep(0, k), Sigt, log = TRUE)
llike = llike - d1
}
llike
}
fit = nlminb(start = par, objective = mlikeG, RTN = RTN,
include.mean = include.mean, lower = c1, upper = c2)
va = mlikeG(fit$par, RTN = RTN, include.mean = include.mean)
cat("Log likelihood function: ", -va, "\n")
epsilon = 0.00025 * fit$par
if (k == 3)
epsilon = 5e-04 * fit$par
npar = length(par)
Hessian = matrix(0, ncol = npar, nrow = npar)
for (i in 1:npar) {
for (j in 1:npar) {
x1 = x2 = x3 = x4 = fit$par
x1[i] = x1[i] + epsilon[i]
x1[j] = x1[j] + epsilon[j]
x2[i] = x2[i] + epsilon[i]
x2[j] = x2[j] - epsilon[j]
x3[i] = x3[i] - epsilon[i]
x3[j] = x3[j] + epsilon[j]
x4[i] = x4[i] - epsilon[i]
x4[j] = x4[j] - epsilon[j]
Hessian[i, j] = (mlikeG(x1, RTN = RTN, include.mean = include.mean) -
mlikeG(x2, RTN = RTN, include.mean = include.mean) -
mlikeG(x3, RTN = RTN, include.mean = include.mean) +
mlikeG(x4, RTN = RTN, include.mean = include.mean))/(4 *
epsilon[i] * epsilon[j])
}
}
est = fit$par
se.coef = sqrt(diag(solve(Hessian)))
tval = fit$par/se.coef
matcoef = cbind(fit$par, se.coef, tval, 2 * (1 - pnorm(abs(tval))))
dimnames(matcoef) = list(names(tval), c(" Estimate",
" Std. Error", " t value", "Pr(>|t|)"))
cat("\nCoefficient(s):\n")
printCoefmat(matcoef, digits = 6, signif.stars = TRUE)
BEKK11vol <- function(rtn, est, include.mean = T) {
if (!is.matrix(rtn))
rtn = as.matrix(rtn)
nT = dim(rtn)[1]
k = dim(rtn)[2]
if (k > 3) {
k = 3
rtn = rtn[, 1:k]
}
mu1 = 0
mu2 = 0
mu3 = 0
icnt = 0
if (k == 2) {
if (include.mean) {
mu1 = est[1]
mu2 = est[2]
icnt = k
}
A0 = matrix(c(est[icnt + 1], est[icnt + 2], 0, est[icnt +
3]), 2, 2)
A0A0t = A0 %*% A0
A1 = matrix(c(est[(icnt + 4):(icnt + 7)]), 2, 2)
B1 = matrix(c(est[(icnt + 8):(icnt + 11)]), 2, 2)
resi = cbind(rtn[, 1] - mu1, rtn[, 2] - mu2)
Cov1 = cov(rtn)
Sig = matrix(c(Cov1), 1, 4)
res = resi %*% t(A1)
ArchP = A1at(res)
}
if (k == 3) {
if (include.mean) {
mu1 = est[1]
mu2 = est[2]
mu3 = est[3]
icnt = k
}
A0 = matrix(c(est[icnt + 1], est[icnt + 2], est[icnt +
3], 0, est[icnt + 4], est[icnt + 5], 0, 0, est[icnt +
6]), k, k)
A0A0t = A0 %*% t(A0)
icnt = icnt + 6
A1 = matrix(c(est[(icnt + 1):(icnt + 9)]), k, k)
icnt = icnt + 9
B1 = matrix(c(est[(icnt + 1):(icnt + 9)]), k, k)
resi = cbind(rtn[, 1] - mu1, rtn[, 2] - mu2, rtn[,
3] - mu3)
Cov1 = cov(rtn)
Sig = matrix(c(Cov1), 1, k^2)
res = resi %*% t(A1)
ArchP = A1at(res)
}
for (t in 2:nT) {
Sigt = A0A0t + matrix(ArchP[t - 1, ], k, k) + B1 %*%
matrix(Sig[t - 1, ], k, k) %*% t(B1)
Sigt = (Sigt + t(Sigt))/2
Sig = rbind(Sig, c(Sigt))
}
BEKK11vol <- list(volmtx = Sig)
}
m2 = BEKK11vol(RTN, est, include.mean = include.mean)
Sigma.t = m2$volmtx
BEKK11 <- list(data = rt, estimates = est, HessianMtx = Hessian,
Sigma.t = Sigma.t, include.mean = include.mean)
}
<bytecode: 0x00000000894c05c0>
<environment: namespace:MTS>
However, to calculate the Hedge Optimal Ratio of Garch BEKK, I would need to make the following calculation:
Where ΔS _ t, ΔF _ t is the return price Spot and Future, and Ω_t-1 is conditional variance and covariance matrix. But I can’t calculate this matrix. How could I do?