Wrong number of parameters in VARXorder #14
Description
Hello,
In the code, please see below, k is the dimension of the output process, and m is the dimension of the exogenous process. Thus the number of parameters should be:
kkp + km(mm+1), with p the autoregressive order, and mm the number of lags of the exogenous process, k the dimension of the output process, and m the dimension of the exogenous process.
In the code, you defined:
ksq = k * k
and later in the code:
npar = p * ksq + k * (mm + 1)
But npar is wrong and the correct number of parameters should be
npar = p * ksq + k * m * (mm + 1)
The code of VARXorder is given below. The number of parameters need to be corrected at some places. If needed, we could send you a corrected code.
---
VARXorder
function (x, exog, maxp = 13, maxm = 3, output = T)
{
x1 = as.matrix(x)
exog = as.matrix(exog)
nT = dim(x1)[1]
k = dim(x1)[2]
ksq = k * k
if (maxp < 1)
maxp = 1
nT1 = dim(exog)[1]
m = dim(exog)[2]
if (nT1 > nT) {
cat("Adjustment made for different nobs:", c(nT,
nT1), "\n")
}
if (nT > nT1) {
cat("Adjustment made for different nobs:", c(nT,
nT1), "\n")
nT = nT1
}
aic = matrix(0, maxp + 1, maxm + 1)
bic = aic
hq = aic
for (mm in 0:maxm) {
isto = mm + 1
y = x1[isto:nT, ]
xm = rep(1, rep(nT - mm))
for (i in 0:mm) {
xm = cbind(xm, exog[(isto - i):(nT - i), ])
}
xm = as.matrix(xm)
xpx = crossprod(xm, xm)
xpxi = solve(xpx)
xpy = t(xm) %% y
beta = xpxi %% xpy
y = y - xm %% beta
s = t(y) %% y/(nT - mm)
enob = nT - mm
c1 = log(det(s))
aic[1, mm + 1] = c1 + 2 * k * mm/enob
bic[1, mm + 1] = c1 + log(enob) * k * m/enob
hq[1, mm + 1] = c1 + 2 * log(log(enob)) * k * m/enob
for (p in 1:maxp) {
ist = max(mm + 1, p + 1)
enob = nT - ist + 1
y = as.matrix(x1[ist:nT, ])
xmtx = rep(1, enob)
for (i in 0:mm) {
xmtx = cbind(xmtx, exog[(ist - i):(nT - i), ])
}
for (j in 1:p) {
xmtx = cbind(xmtx, x1[(ist - j):(nT - j), ])
}
xm1 = as.matrix(xmtx)
xpx = t(xm1) %% xm1
xpxinv = solve(xpx)
xpy = t(xm1) %% y
beta = xpxinv %% xpy
resi = y - xm1 %% beta
sse = (t(resi) %*% resi)/enob
d1 = log(det(sse))
npar = p * ksq + k * (mm + 1)
aic[p + 1, mm + 1] = d1 + (2 * npar)/enob
bic[p + 1, mm + 1] = d1 + (log(enob) * npar)/enob
hq[p + 1, mm + 1] = d1 + (2 * log(log(enob)) * npar)/enob
}
}
ind.min <- function(x) {
if (!is.matrix(x))
x = as.matrix(x)
r = dim(x)[1]
c = dim(x)[2]
xm = min(x)
COntin = TRUE
while (COntin) {
for (j in 1:c) {
idx = c(1:r)[x[, j] == xm]
jj = j
if (length(idx) > 0) {
ii = c(1:r)[x[, jj] == xm][1]
kdx = c(ii, jj)
COntin = FALSE
}
}
}
kdx
}
aicor = ind.min(aic) - 1
bicor = ind.min(bic) - 1
hqor = ind.min(hq) - 1
if (output) {
cat("selected order(p,s): aic = ", aicor, "\n")
cat("selected order(p,s): bic = ", bicor, "\n")
cat("selected order(p,s): hq = ", hqor, "\n")
}
VARXorder <- list(aic = aic, aicor = aicor, bic = bic, bicor = bicor,
hq = hq, hqor = hqor)
}