0
The function BEKK11 of package MTS generates directly on screen because of the functions:
cat("Initial estimates: ", par, "\n")
cat("Lower limits: ", c1, "\n")
cat("Upper limits: ", c2, "\n")
cat("Log likelihood function: ", -va, "\n") e
printCoefmat(matcoef, digits = 6, signif.stars = TRUE)
How do these functions not directly generate on R, but without losing the functions? I exclude them, however, in the end only generate the coefficient and I would like it to generate everything.
complete code
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: 0x00000000661fee38>
<environment: namespace:MTS>
Thanks Rui, just a doubt, I keep deleting those functions "cat" and "printCoefmat"?
– Everton Toledo
@Evertontoled If you want to see the results while the function runs, keep them that does no harm.
– Rui Barradas
Gee Rui, I don’t have the mastery over R. Do you have the whole function as it would be? If it’s not uncomfortable. It’s not working for me, I think I’m including wrong in function.
– Everton Toledo
@Evertontoled The two answer instructions (edited, left to close the last parenthesis) must be the last function instructions. That’s it. It’s only by these instructions right at the end of the function.
– Rui Barradas
@Evertontoledo Look now.
– Rui Barradas
It worked out! Thanks Rui
– Everton Toledo
Just one more question, Rui, if you can help me. I could include signif.codes?
– Everton Toledo
@Evertontoled In function this is done with the instruction
printCoefmat(result$Everton$coef, digits = 6, signif.stars = TRUE).– Rui Barradas
Actually I deleted printCoefmat to generate PDF in Latex
– Everton Toledo