How to calculate the variance and covariance matrix of the BEKK11 function of the MTS package

Asked

Viewed 69 times

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:

inserir a descrição da imagem aqui

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?

No answers

Browser other questions tagged

You are not signed in. Login or sign up in order to post.