Nothing
##########################################################################
# These functions are
# Copyright (C) 2014-2020 V. Miranda & T. Yee
# Auckland University of Technology & University of Auckland
# All rights reserved.
#
# Links renamed on Jan-2019 conforming with VGAM_1.1-0
# Supports the Wald, score, and lrt tests (20180209)
ARMA.EIM.G2 <- function(y, arCoe, maCoe, estRes, sdError,
var.arg = TRUE, order = c(1, 1),
nodrift = FALSE, addRidge = 0.001) {
if (length(order) !=2 )
stop("Currently, this function only handles ARMA (p, q) models.")
Order <- order; rm(order)
arOrd <- Order[1]; maOrd <- Order[2]
y <- cbind(y); nn <- nrow(y)
arCoe <- matrix(arCoe, nrow = nn, ncol = arOrd)
# R matrix
RMat <- diag(c(sdError^2))
dRMat <- cbind( if(nodrift) NULL else rep(0, nn),
if (var.arg) 1 else 2 * sdError,
matrix(0, nn, arOrd + maOrd))
# dm(t) matrix
arLags <- WN.lags(y = cbind(y), lags = arOrd,
to.complete = 1/(1 - arCoe[1, ] + 0.01)^2)
maLags <- WN.lags(y = cbind(estRes), lags = maOrd,
to.complete = 1/(1 - maCoe[1, ] + 0.01)^2)
dmt <- cbind( if (nodrift) NULL else rep(1, nn),
rep(0, nn),
WN.lags(y = cbind(y), lags = arOrd,
to.complete = 1/(1 - arCoe[1, ] + 0.01)^2),
WN.lags(y = cbind(estRes), lags = maOrd,
to.complete = 1/(1 - maCoe[1, ] + 0.01)^2))
##
M <- 2 + arOrd + maOrd # ARMA model
try.comb <- iam(NA, NA, M = M - nodrift, both = TRUE)
try.comb <- cbind(try.comb$row.index, try.comb$col.index)
### This is equation A.6, Porat and Friedlander (1986)
finMat <- apply(try.comb, 1, function(x){
kl <- x
invRMat <- diag(solve(RMat))
dRthel <- c(dRMat[, kl[1]])
dRthek <- c(dRMat[, kl[2]])
term.1 <- (0.5) * invRMat * dRthel * invRMat * dRthek
term.2 <- c(dmt[, kl[1]]) * invRMat * c(dmt[, kl[2]])
term.1 + term.2
})
finMat[, 1:M] <- (1 + addRidge) * finMat[, 1:M]
finMat
}
# ARMApq density
dARMA <- function(x, mean = 0, sd = 1, log = FALSE)
dnorm(x = x, mean = mean, sd = sd, log = log)
ARMAXff <-
function(order = c(1, 1),
zero = c(if (nodrift) NULL else "drift",
"ARcoeff", "MAcoeff"),
xLag = 0,
type.EIM = c("exact", "approximate")[1],
var.arg = TRUE,
nodrift = FALSE,
noChecks = FALSE,
ldrift = "identitylink",
lsd = "loglink",
lvar = "loglink",
lARcoeff = "identitylink",
lMAcoeff = "identitylink",
idrift = NULL, # Must be a vector of length NOS
isd = NULL, # Must be a vector of length NOS
ivar = NULL, # Must be a vector of length NOS
iARcoeff = NULL, # Must be a vectors
iMAcoeff = NULL) {
if ( !Is.Numeric(order, length.arg = 2, isInteger = TRUE ) )
stop("Invalid 'order'. Should be a vector (p, q) of integers.")
if (xLag < 0)
stop("Wrong input for argument 'xLag'")
ARord <- order[1]
MAord <- order[2]; rm(order)
nOrder <- ARord + MAord
realM1 <- 2 + ARord + MAord - nodrift
if ( ARord == 0 || MAord == 0 )
stop("Only ARMAX(p, q) models, p, q > 0, handled." ,
" Please, refer to ARXff() or MAXff().")
if (!is.logical(noChecks) && !is.logical (var.arg) &&
!is.logical(nodrift))
stop("'noChecks', 'nodrift', and 'var.arg' must be logical")
type.likelihood <- "exact"
type.EIM <- match.arg(type.EIM, c("exact", "approximate"))[1]
ldrift <- as.list(substitute(ldrift))
edrift <- link2list(ldrift)
ldrift <- attr(edrift, "function.name")
lvar <- as.list(substitute(lvar))
evar <- link2list(lvar)
lvar <- attr(evar, "function.name")
lsd <- as.list(substitute(lsd))
esd <- link2list(lsd)
lsd <- attr(esd, "function.name")
lARcoeff <- as.list(substitute(lARcoeff))
eARcoeff <- link2list(lARcoeff)
lARcoeff <- attr(eARcoeff, "function.name")
lMAcoeff <- as.list(substitute(lMAcoeff))
eMAcoeff <- link2list(lMAcoeff)
lMAcoeff <- attr(eMAcoeff, "function.name")
pre.blurb1 <- paste("ARcoeff", 1:ARord, sep = "")
pre.blurb2 <- paste("MAcoeff", 1:MAord, sep = "")
blurb.vec <-
c(if (nodrift) NULL else
namesof("drift",
link = ldrift, earg = edrift, tag = FALSE),
if (var.arg)
namesof("noiseVar", link = lvar, earg = evar, tag = FALSE)
else
namesof("noiseSD", link = lsd, earg = esd, tag = FALSE) )
for (jj in 1:ARord)
blurb.vec <- c( blurb.vec, namesof(pre.blurb1[jj],
link = lARcoeff, earg = eARcoeff))
for (kk in 1:MAord) # tag = FALSE by default.. not needed
blurb.vec <- c( blurb.vec, namesof(pre.blurb2[kk],
link = lMAcoeff, earg = eMAcoeff))
blurb.vec2 <- character(0)
for (kk in 1:realM1) {
if (kk != realM1)
aux <- paste(blurb.vec[kk], ", ", sep = "")
else
aux <- paste(blurb.vec[kk], ".", sep = "")
blurb.vec2 <- c(blurb.vec2, aux)
}
preAR <- if (ARord == 1) "(ARcoeff1) * Y[t - 1] + "
preMA <- if (MAord == 1) "(MAcoeff1) * e[t - 1] "
alo <-
new("vgltsmff",
blurb = c("VGLTMs: ", nOrder + 1 + !nodrift,
"-parameter ARMAX model of order",
"(", ARord, ",", MAord, ").",
"\n", "\n",
"LINKS: ",
blurb.vec2, "\n",
"Model: Y[t] = ",
if (nodrift) NULL else "drift + ",
" B^T * X_t + ",
if (ARord == 1) preAR else
c("(ARcoeff1) * Y[t - 1] +... + (ARcoeff", ARord,")",
" * Y[t -", ARord,"] + "),
if (MAord == 1) preMA else
c(" (MAcoeff1) * e[t - 1] + ... + (MAcoeff", MAord,")",
" * e[t - ", MAord,"]") , "+ e[t].", "\n"),
#"E[Yt] = drift / (1 - ARcoeff1 - ... - ARcoeff",
#ARord, ")"),
constraints = eval(substitute(expression({
M1 <- .realM1
constraints <- cm.zero.VGAM(constraints, x = x,
zero = .zero , M = M,
# Old : predictors.names...
predictors.names = parameters.names,
M1 = M1)
}),
list( .zero = zero , .realM1 = realM1 ) )),
first = eval(substitute(expression({
if ((NCOL(x) == 1) && (colnames(x) == "(Intercept)") && ( .xLag > 0))
stop("Set xLag = 0. No covariates entered.")
if ( .xLag ) {
x.matrix <- x
for (jj in 1:(NCOL(x) - 1)) {
temp1 <- WN.lags(cbind(x[, 1 + jj]), .xLag )
temp2 <- paste(colnames(x)[1 + jj], "Lag", sep = "")
colnames(temp1) <- paste(temp2, 1:( .xLag), sep = "")
x.matrix<- cbind(x.matrix, temp1)
}
rm(temp1, temp2)
list.names <- vector("list", NCOL(x.matrix))
names(list.names) <- colnames(x.matrix)
for (ii in 1:(NCOL(x.matrix)))
list.names[[ii]] <- ii
attr(x.matrix, "assign") <- list.names
x <- x.matrix
}
}), list( .xLag = xLag ))),
infos = eval(substitute(function(...) {
list(M1 = .realM1 ,
Order = .nOrder ,
Q1 = 1,
ldrift = .ldrift ,
lsd = .lsd ,
lvar = .lvar ,
lARcoeff = .lARcoeff ,
lMAcoeff = .lMAcoeff ,
edrift = .edrift ,
esd = .esd ,
evar = .evar ,
eARcoeff = .eARcoeff ,
eMAcoeff = .eMAcoeff ,
ARord = .ARord ,
MAord = .MAord ,
zero = .zero ,
nodrift = .nodrift ,
type.EIM = .type.EIM ,
expected = TRUE,
multipleResponse = TRUE )
}, list( .realM1 = realM1 , .nOrder = nOrder,
.ldrift = ldrift , .lsd = lsd , .lvar = lvar ,
.lARcoeff = lARcoeff , .lMAcoeff = lMAcoeff ,
.edrift = edrift , .esd = esd , .evar = evar ,
.eARcoeff = eARcoeff , .eMAcoeff = eMAcoeff ,
.ARord = ARord , .MAord = MAord ,
.nodrift = nodrift , .type.EIM = type.EIM ,
.zero = zero ))),
initialize = eval(substitute(expression({
M1 <- .realM1
mycheck <- w.y.check(w = w, y = y,
Is.positive = FALSE,
ncol.w.max = Inf,
ncol.y.max = Inf,
out.wy = TRUE,
colsyperw = 1,
maximize = TRUE)
w <- mycheck$w
y <- extra$y <- mycheck$y
NOS <- ncol(mycheck$y)
M <- NOS * M1; n <- nrow(y)
nOrder <- .nOrder
if (length( .iARcoeff ) || length( .iMAcoeff ))
quick.check.coeffs(arc = .iARcoeff, mac = .iMAcoeff,
NOS = NOS, arOrd = .ARord, maOrd = .MAord)
myIdrift <- rep( .idrift , NOS)[1:NOS]
myIvar <- rep( .ivar , NOS)[1:NOS]
myIsd <- rep( .isd , NOS)[1:NOS]
iniar <- if (!length( .iARcoeff )) NULL else
matrix( .iARcoeff , NOS, .ARord , byrow = TRUE)
inima <- if (!length( .iMAcoeff )) NULL else
matrix( .iMAcoeff , NOS, .MAord , byrow = TRUE)
drift.n <- if (NOS == 1) "drift.mean" else
paste("drift.mean", 1:NOS, sep = "")
varsd.n <- if ( .var.arg ) if (NOS == 1) "noiseVar" else
paste("noiseVar", 1:NOS, sep="") else
if (NOS == 1) "noiseSD" else
paste("noiseSD", 1:NOS, sep = "")
ar.n <- ma.n <- character(0)
for (jj in 1:( .ARord ))
ar.n <- if (NOS == 1) c(ar.n, paste("ARcoeff", jj, sep = "")) else
c(ar.n, paste(paste("ARcoeff", jj, sep = ""), 1:NOS, sep = ""))
for (jj in 1:( .MAord ))
ma.n <- if (NOS == 1) c(ma.n, paste("MAcoeff", jj, sep = "")) else
c(ma.n, paste(paste("MAcoeff", jj, sep = ""), 1:NOS, sep = ""))
parameters.names <- c(if ( .nodrift ) NULL else drift.n,
varsd.n, ar.n, ma.n)
predictors.names <- c( if( .nodrift ) NULL else
namesof(drift.n, link = .ldrift , earg = .edrift, tag = FALSE),
namesof(varsd.n, link = if ( .var.arg ) .lvar else .lsd ,
earg = if ( .var.arg ) .evar else .esd ),
namesof(ar.n, link = .lARcoeff , earg = .eARcoeff , tag = FALSE),
namesof(ma.n, link = .lMAcoeff , earg = .eMAcoeff , tag = FALSE))
predictors.names <- predictors.names[interleave.VGAM(M1*NOS, M1 = M1)]
if (!length(etastart)) {
init.dr <- matrix(0.0, nrow = n, ncol = NOS)
init.sig <- matrix(0.0, nrow = n, ncol = NOS)
init.AR <- array( 0.0, dim = c( n, NOS , .ARord ))
init.MA <- array( 0.0, dim = c( n, NOS , .MAord ))
y.sc <- y# scale(y, center = TRUE, scale = FALSE)
#extra$y.sc <- scale(y, center = TRUE, scale = FALSE)
extra$res <- matrix(NA_real_, nrow = n, ncol = NOS)
res2.mat <- matrix(NA_real_, nrow = n, ncol = NOS)
for (rsp in 1:NOS) {
to.fit <- cbind(y.sc[, rsp, drop = FALSE],
WN.lags(y = cbind(y.sc[, rsp, drop = FALSE]),
lags = .ARord + .MAord + 5))
to.fit <- lsfit(x = to.fit[, -1, drop = FALSE],
y = to.fit[, 1, drop = FALSE],
intercept = FALSE)
extra$res <- res2.mat[, rsp] <- residuals(to.fit)
# Updating residuals not allowed.
#res2.mat[, rsp] <- residuals(to.fit)
# Initial values.
if (TRUE) {
to.fit <- cbind(y.sc[, rsp, drop = FALSE],
WN.lags(y = cbind(y.sc[, rsp, drop = FALSE]),
lags = .ARord ),
WN.lags(y = res2.mat[, rsp, drop = FALSE],
lags = .MAord ))
to.fit <- lsfit(x = to.fit[, -1, drop = FALSE],
y = to.fit[, 1, drop = FALSE],
intercept = FALSE)
}
initCoe <- coef(to.fit)
extra$res <- res2.mat[, rsp] <- residuals(to.fit)
just.hlp <-matrix(if (length( .iARcoeff )) iniar[ rsp, ] else
initCoe[1:( .ARord )], n, .ARord , byrow = TRUE)
init.AR[, rsp, ] <- just.hlp
just.hlp <-matrix(if (length( .iMAcoeff )) inima[ rsp, ] else
initCoe[-c(1:( .ARord ))], n, .MAord , byrow = TRUE)
init.MA[, rsp, ] <- just.hlp
init.dr[, rsp] <- if (length( .idrift )) myIdrift[rsp] else
mean(y[, rsp]) * (1 - sum(initCoe[1:( .ARord )]))
init.sig[, rsp] <- if( .var.arg ) var(res2.mat[, rsp]) else
sd(res2.mat[, rsp])
}
#res2.mat[c(1: max( .ARord, .MAord )), 1] <- median(res2.mat[, 1])
extra$res <- res2.mat
etastart <- cbind(if ( .nodrift ) NULL else init.dr ,
init.sig,
matrix(init.AR, n, NOS * ( .ARord )),
matrix(init.MA, n, NOS * ( .MAord )))
etastart <- cbind(if ( .nodrift ) NULL else
theta2eta(init.dr, .ldrift , .edrift ),
theta2eta(init.sig, if ( .var.arg ) .lvar else .lsd ,
if ( .var.arg ) .evar else .esd ),
theta2eta(matrix(init.AR, n, NOS * ( .ARord )),
.lARcoeff , .eARcoeff ),
theta2eta(matrix(init.MA, n, NOS * ( .MAord ) ),
.lMAcoeff , .eMAcoeff))
etastart <- etastart[, interleave.VGAM(M1 * NOS, M1 = M1)]
etastart
} # End of !etastart
}), list( .ldrift = ldrift , .lvar = lvar , .lsd = lsd ,
.lARcoeff = lARcoeff , .lMAcoeff = lMAcoeff ,
.edrift = edrift , .evar = evar , .esd = esd ,
.eARcoeff = eARcoeff , .eMAcoeff = eMAcoeff ,
.idrift = idrift , .ivar = ivar , .isd = isd ,
.iARcoeff = iARcoeff , .iMAcoeff = iMAcoeff ,
.var.arg = var.arg , .realM1 = realM1 ,
.ARord = ARord , .MAord = MAord , .nOrder = nOrder ,
.nodrift = nodrift ))),
linkinv = eval(substitute(function(eta, extra = NULL){
M1 <- .realM1; n <- nrow(eta)
NOS <- ncol(eta)/M1; NNro <- 1
drifts <- if ( .nodrift ) matrix(0.0, nrow = n, ncol = NOS) else
eta2theta(eta[, M1*(1:NOS) - .nOrder - 1, drop = FALSE],
link = .ldrift , earg = .edrift )
sigWN <- eta2theta(eta[, M1*(1:NOS) - .nOrder, drop = FALSE],
link = if ( .var.arg ) .lvar else .lsd ,
earg = if ( .var.arg ) .evar else .esd )
ars <- array(NA_real_, dim = c(n, NOS, .ARord ))
mas <- array(NA_real_, dim = c(n, NOS, .MAord ))
for (ai in 1:( .ARord ))
ars[, , ai] <- eta2theta(eta[, M1 * (1:NOS) -
( .ARord + .MAord ) + ai, drop = FALSE],
link = .lARcoeff , earg = .eARcoeff )
for ( mi in 1:( .MAord ))
mas[, , mi] <- eta2theta(eta[, M1 * (1:NOS) -
( .MAord ) + mi , drop = FALSE],
link = .lMAcoeff , earg = .eMAcoeff )
y.est <- matrix(NA_real_, n , NOS)
for (ii in 1:NOS) {
y.est[, ii] <- drifts[, ii] + rowSums(cbind(ars[, ii, ] *
WN.lags(y = cbind(extra$y[, ii]), lags = .ARord ),
mas[, ii, ] * WN.lags(y = cbind(extra$res[, ii]),
lags = .MAord ))) # + extra$res[, ii]
y.est[c(1:NNro), ii] <- (NNro + 0.10) *
y.est[c(2:(NNro + 1)), ii]
}
y.est
}, list( .ldrift = ldrift , .lvar = lvar , .lsd = lsd ,
.lARcoeff = lARcoeff , .lMAcoeff = lMAcoeff ,
.edrift = edrift , .evar = evar , .esd = esd ,
.eARcoeff = eARcoeff , .eMAcoeff = eMAcoeff ,
.idrift = idrift , .ivar = ivar , .isd = isd ,
.iARcoeff = iARcoeff , .iMAcoeff = iMAcoeff ,
.var.arg = var.arg , .realM1 = realM1 ,
.nodrift = nodrift , .nOrder = nOrder ,
.ARord = ARord , .MAord = MAord ))),
last = eval(substitute(expression({
M1 <- .realM1
res.mat <- extra$res
drifts <- if ( !(.nodrift) )
eta2theta(eta[, M1*(1:NOS) - .nOrder - 1 , drop = FALSE],
link = .ldrift , earg = .edrift ) else
if (.nodrift ) matrix(0.0, nrow = n, ncol = NOS)
sigWN <- if ( .var.arg )
eta2theta(eta[, M1*(1:NOS) - .nOrder , drop = FALSE],
link = .lvar , earg = .evar ) else
eta2theta(eta[ , M1*(1:NOS) - .nOrder , drop = FALSE],
link = .lsd , earg = .esd )
ars <- array(NA, dim = c(n, NOS, .ARord ))
mas <- array(NA, dim = c(n, NOS, .MAord ))
for (ai in 1:.ARord)
ars[, , ai ] <-
eta2theta(eta[, M1*(1:NOS) - (.MAord + .ARord) + ai,
drop = FALSE],
link = .lARcoeff , earg = .eARcoeff )
for ( mi in 1:.MAord )
mas[, , mi ] <-
eta2theta(eta[, M1*(1:NOS) - ( .MAord ) + mi , drop = FALSE],
link = .lMAcoeff , earg = .eMAcoeff )
prov.names <- c(if ( .nodrift ) NULL else drift.n , varsd.n)
for (kk in 1:.ARord )
prov.names <- c(prov.names, ar.n[[kk]])
for (kk in 1:.MAord )
prov.names <- c(prov.names, ma.n[[kk]])
misc$link <- rep( .ldrift , times = M1 * NOS )
names(misc$link) <-
prov.names[interleave.VGAM(M1 * NOS , M1 = M1)]
misc$earg <- vector("list", length = M1 * NOS)
names(misc$earg) <-
prov.names[interleave.VGAM(M1 * NOS , M1 = M1)]
for (jj in 1:NOS) {
if (!(.nodrift ))
misc$link[ M1 * jj - .nOrder - 1 ] <- .ldrift
misc$link[M1 * jj - .nOrder] <- if ( .var.arg ) .lvar else .lsd
if (!(.nodrift ))
misc$earg[[ M1 * jj - .nOrder - 1 ]] <- .edrift
misc$earg[[M1 * jj - .nOrder]] <- if ( .var.arg ) .evar else .esd
for (kk in 1:.ARord ) {
misc$link[ M1 * jj - .nOrder + kk ] <- .lARcoeff
misc$earg[[ M1 * jj - .nOrder + kk ]] <- .eARcoeff
}
for( ll in 1:.MAord ) {
misc$link[ M1 * jj - .MAord + ll ] <- .lMAcoeff
misc$earg[[ M1 * jj - .MAord + ll ]] <- .eMAcoeff
}
}
extra$res[c(1: max( .ARord, .MAord )), 1] <- median(res.mat)
misc$noChecks <- .noChecks
misc$expected <- TRUE
misc$process <- "ARMA"
misc$var.arg <- .var.arg
misc$nomean <- .nodrift
misc$Order <- c( .ARord , .MAord )
misc$M1 <- .realM1
misc$NOS <- NOS
misc$ARord <- .ARord
misc$MAord <- .MAord
misc$flagAR <- FALSE
misc$residuals <- extra$res
misc$theta.names <- parameters.names
misc$multipleResponses <- TRUE
misc$control$wzepsilon <- control$wzepsilon
fit$prior.weights <- w # 20180315
# Checks removed on March/2/2016. Now, refer to 'summary'.
if (!( .noChecks )) {
M <- M1 * NOS
FlagArma <- grep("(Intercept)", names(fit$coefficients))
if (length(FlagArma) == M) {
coefArma <- fit$coefficients
coefMa <- coefAr <- numeric(0)
for ( kk in 1:NOS ) {
myAux2 <- coefArma[ M1 * kk - .nOrder + 1:( .ARord )]
myAux3 <- coefArma[ M1 * kk - .MAord + 1:( .MAord )]
coefAr <- c(coefAr, myAux2)
coefMa <- c(coefMa, myAux3)
}
coefAr <- eta2theta(coefAr, link = .lARcoeff , earg = .eARcoeff)
coefMa <- eta2theta(coefMa, link = .lMAcoeff , earg = .eMAcoeff)
# 2016/March/01. Roots printed out by 'summary(vglm)'. '.
checkAR <- checkTS.ffs(thetaEst = coefAr,
tsclass = "AR",
NofS = NOS,
chOrder = .ARord ,
pRoots = FALSE,
retmod = TRUE)
checkMA <- checkTS.ffs(thetaEst = coefMa,
tsclass = "MA",
NofS = NOS,
chOrder = .MAord ,
pRoots = FALSE,
retmod = TRUE)
flag1 <- (all(checkAR > 0) && any(checkAR < 1 + 5e-3))
flag1 <- ((all(checkMA > 0) && any(checkMA < 1 + 5e-3)) || flag1)
if (!flag1)
cat("\nChecks on stationarity / invertibility successfully",
"performed. \nNo roots liying inside the unit circle.",
"\nFurther details within the 'summary' output.\n") else
cat("\nChecks on stationarity / invertibility successfully",
"performed. \nNOTE: Some roots lie inside the unit ",
"circle. \nFurther details within the 'summary' ",
"output.\n")
} else {
cat("\nApparently, some constraints have been set.",
"\nDetails on stationarity / invertibility",
"checks \nwhitin the 'summary' output.")
}
} # End of !( .noChecks )
}), list( .ldrift = ldrift , .lvar = lvar , .lsd = lsd ,
.lARcoeff = lARcoeff , .lMAcoeff = lMAcoeff ,
.edrift = edrift , .evar = evar , .esd = esd ,
.eARcoeff = eARcoeff , .eMAcoeff = eMAcoeff ,
.var.arg = var.arg , .realM1 = realM1 ,
.ARord = ARord , .MAord = MAord ,
.nodrift = nodrift , .noChecks = noChecks ,
.nOrder = nOrder,
.type.likelihood = type.likelihood ))),
loglikelihood =
eval(substitute(function(mu, y, w, residuals = FALSE,
eta, extra = NULL, summation = TRUE) {
#NOS <- extra$NOS
y <- extra$y
M1 <- .realM1; n <- nrow(eta)
NOS <- ncol(eta)/M1
eim.ap <- ( .type.EIM == "approximate")
drifts <- if ( .nodrift ) matrix(0.0, nrow = n, ncol = NOS) else
eta2theta(eta[, M1*(1:NOS) - .nOrder - 1, drop = FALSE],
link = .ldrift , earg = .edrift )
if ( .var.arg ) {
sigWN <- sqrt(eta2theta(eta[, M1*(1:NOS) - .nOrder,
drop = FALSE], .lvar , .evar ))
} else {
sigWN <- eta2theta(eta[, M1*(1:NOS) - .nOrder,
drop = FALSE], .lsd , .esd )
}
ars <- array(NA_real_, dim = c(n, NOS, .ARord ))
mas <- array(NA_real_, dim = c(n, NOS, .MAord ))
for (ai in 1:( .ARord ))
ars[, , ai] <- eta2theta(eta[, M1 * (1:NOS) -
( .ARord + .MAord ) + ai, drop = FALSE],
link = .lARcoeff , earg = .eARcoeff )
for ( mi in 1:( .MAord ))
mas[, , mi] <- eta2theta(eta[, M1 * (1:NOS) -
( .MAord ) + mi , drop = FALSE],
link = .lMAcoeff , earg = .eMAcoeff )
y.est <- matrix(NA_real_, n , NOS)
for (ii in 1:NOS)
y.est[, ii] <- drifts[, ii] + rowSums(cbind(ars[, ii, ] *
WN.lags(y = cbind(extra$y[, ii]), lags = .ARord ),
mas[, ii, ] * WN.lags(y = cbind(extra$res[, ii]),
lags = .MAord ))) + extra$res
if (residuals) {
stop("Loglikelihood not implemented yet to",
" handle residuals.")
} else {
loglik.terms <- c(w) * dARMA(x = y, mean = y.est,
sd = sigWN, log = TRUE)
}
if ( summation )
loglik.terms <- sum(loglik.terms)
#if (!eim.ap) sum(loglik.terms) else
#round(sum(loglik.terms), 4)
loglik.terms
}, list( .ldrift = ldrift , .lsd = lsd , .lvar = lvar ,
.lARcoeff = lARcoeff , .lMAcoeff = lMAcoeff ,
.edrift = edrift , .esd = esd , .evar = evar ,
.eARcoeff = eARcoeff , .eMAcoeff = eMAcoeff ,
.nOrder = nOrder , .ARord = ARord , .MAord = MAord ,
.type.likelihood = type.likelihood ,
.nodrift = nodrift , .type.EIM = type.EIM ,
.var.arg = var.arg , .realM1 = realM1 ))),
validparams = eval(substitute(function(eta, y, extra = NULL) {
M1 <- .realM1; n <- nrow(eta)
NOS <- ncol(eta)/M1
drifts <- if ( .nodrift ) matrix(0.0, nrow = n, ncol = NOS) else
eta2theta(eta[, M1*(1:NOS) - .nOrder - 1, drop = FALSE],
link = .ldrift , earg = .edrift )
sigWN <- eta2theta(eta[, M1*(1:NOS) - .nOrder, drop = FALSE],
link = if ( .var.arg ) .lvar else .lsd ,
earg = if ( .var.arg ) .evar else .esd )
ars <- array(NA_real_, dim = c(n, NOS, .ARord ))
mas <- array(NA_real_, dim = c(n, NOS, .MAord ))
for (ai in 1:( .ARord ))
ars[, , ai] <- eta2theta(eta[, M1 * (1:NOS) -
( .ARord + .MAord ) + ai, drop = FALSE],
link = .lARcoeff , earg = .eARcoeff )
for ( mi in 1:( .MAord ))
mas[, , mi] <- eta2theta(eta[, M1 * (1:NOS) -
( .MAord ) + mi , drop = FALSE],
link = .lMAcoeff , earg = .eMAcoeff )
okay1 <- (all(is.finite(sigWN)) && all(0 < sigWN) &&
all(is.finite(ars)) && all(is.finite(mas)) &&
all(is.finite(drifts)))
okay1
}, list( .ldrift = ldrift , .lsd = lsd , .lvar = lvar ,
.lARcoeff = lARcoeff , .lMAcoeff = lMAcoeff ,
.edrift = edrift , .esd = esd , .evar = evar ,
.eARcoeff = eARcoeff , .eMAcoeff = eMAcoeff ,
.nOrder = nOrder , .ARord = ARord , .MAord = MAord ,
.type.likelihood = type.likelihood ,
.nodrift = nodrift ,
.var.arg = var.arg , .realM1 = realM1 ))),
#vfamily = c("vgtsff", "vgltsff-class", "ARMAff"),
vfamily = c("ARMAXff", "ARMAvgltsmff"),
deriv = eval(substitute(expression({
M1 <- .realM1; n <- nrow(eta)
NOS <- ncol(eta)/M1
M <- M1 * NOS
drifts <- if ( .nodrift ) matrix(0.0, nrow = n, ncol = NOS) else
eta2theta(eta[, M1*(1:NOS) - .nOrder - 1, drop = FALSE],
link = .ldrift , earg = .edrift )
if ( .var.arg ) {
sigWN <- sqrt(eta2theta(eta[, M1*(1:NOS) - .nOrder,
drop = FALSE], .lvar , .evar ))
} else {
sigWN <- eta2theta(eta[, M1*(1:NOS) - .nOrder,
drop = FALSE], .lsd , .esd )
}
ars <- array(NA_real_, dim = c(n, NOS, .ARord ))
mas <- array(NA_real_, dim = c(n, NOS, .MAord ))
for (ai in 1:( .ARord ))
ars[, , ai] <- eta2theta(eta[, M1 * (1:NOS) -
( .ARord + .MAord ) + ai, drop = FALSE],
link = .lARcoeff , earg = .eARcoeff )
for ( mi in 1:( .MAord ))
mas[, , mi] <- eta2theta(eta[, M1 * (1:NOS) -
( .MAord ) + mi , drop = FALSE],
link = .lMAcoeff , earg = .eMAcoeff )
y.est <- matrix(NA_real_, n , NOS)
for (ii in 1:NOS)
y.est[, ii] <- drifts[, ii] + rowSums(cbind(ars[, ii, ] *
WN.lags(y = cbind(extra$y[, ii]), lags = .ARord ),
mas[, ii, ] *
WN.lags(y = cbind(extra$res[, ii]), lags = .MAord )))
dl.drfmean <- (y - y.est) / sigWN^2
if ( .var.arg ) {
dl.dvar <- (y - y.est)^2 / (2 * sigWN^4) - 1 / (2 * sigWN^2)
} else {
dl.dsd <- (y - y.est)^2 / sigWN^3 - 1 / sigWN
}
dl.dThe <- array(NA_real_, dim = c(n, NOS, .ARord ))
dl.dPhi <- array(NA_real_, dim = c(n, NOS, .MAord ))
for (ii in 1:NOS)
dl.dThe[, ii, ] <- WN.lags(y = cbind(y[, ii]), lags = .ARord ) *
matrix((y[, ii]- y.est[, ii]), n , .ARord ) / sigWN[, ii]^2
for (ii in 1:NOS)
dl.dPhi[, ii, ] <- WN.lags(y = cbind(extra$res[, ii]),
lags = .MAord ) *
matrix((y[, ii] - y.est[, ii]), n , .MAord ) / sigWN[, ii]^2
ddrif.deta <- dtheta.deta(drifts , .ldrift , .edrift)
dsewn.deta <- if (.var.arg) dtheta.deta(sigWN^2, .lvar , .evar) else
dtheta.deta(sigWN, .lsd , .esd)
dThe.deta <- dtheta.deta(ars , .lARcoeff , .eARcoeff)
dPhi.deta <- dtheta.deta(mas , .lMAcoeff , .eMAcoeff)
myderiv <- c(w) * cbind(if ( .nodrift ) NULL else
dl.drfmean * ddrif.deta,
if ( .var.arg ) dl.dvar * dsewn.deta else
dl.dsd * dsewn.deta,
matrix( dl.dThe * dThe.deta, n, NOS * ( .ARord )),
matrix( dl.dPhi * dPhi.deta, n, NOS * ( .MAord )))
myderiv <- myderiv[, interleave.VGAM(( .realM1 )*NOS , M1 = .realM1)]
myderiv
}), list( .ldrift = ldrift , .lsd = lsd , .lvar = lvar ,
.lARcoeff = lARcoeff , .lMAcoeff = lMAcoeff ,
.edrift = edrift , .esd = esd , .evar = evar ,
.eARcoeff = eARcoeff , .eMAcoeff = eMAcoeff ,
.ARord = ARord , .MAord = MAord , .nOrder = nOrder ,
.var.arg = var.arg , .realM1 = realM1,
.nodrift = nodrift ))),
weight = eval(substitute(expression({
if (( .type.EIM ) == "exact" ) {
final.eim <- numeric(0)
for (ii in 1:NOS) {
exact.eim <- ARMA.EIM.G2(y = cbind(y[, ii]),
arCoe = cbind(ars[, ii, ]),
maCoe = cbind(mas[, ii, ]),
estRes = cbind(extra$res[, ii]),
sdError = cbind(sigWN[, ii]),
var.arg = .var.arg ,
order = c(.ARord , .MAord ),
nodrift = .nodrift )
comb.wz <- combVGAMextra(1:M1)
a.bind <- cbind(if ( .nodrift ) NULL else ddrif.deta[, ii],
dsewn.deta[, ii],
dThe.deta[, ii, ],
dPhi.deta[, ii, ])
dthdeta <- apply(comb.wz, 1, function(x) {
a.bind[, x[1]] * a.bind[, x[2]]
})
final.eim <- cbind(final.eim, c(w[, ii]) * exact.eim * dthdeta)
}
final.eim <-
final.eim[, interleave.VGAM( NOS * M1 * (M1 + 1) / 2, M1 = NOS)]
final.eim <- array(c(final.eim), dim = c(n, NOS, M1 * (M1 + 1)/ 2))
wz <- arwz2wz(final.eim, M = M1 * NOS, M1 = M1)
} else {
ned2l.dsmn <- 1 / sigWN^2
ned2l.dvarSD <- if ( .var.arg ) 1 / (2 * sigWN^4) else 2 / sigWN^2
gammas.y <- apply(y, 2, function(x) {
cross.gammas(x = x, lags = 0)
})
ned2l.dthe <- matrix(gammas.y, n, NOS * ( .ARord ), byrow = TRUE) /
matrix(sigWN^2, n, NOS * ( .ARord ))
dThe.deta <- matrix(dThe.deta, n, NOS * ( .ARord ))
ned2l.dphi <- matrix(1.0, n, NOS * ( .MAord ))
dPhi.deta <- matrix(dPhi.deta, n , NOS * ( .MAord ))
wz <- c(w) * cbind( if ( .nodrift ) NULL else
ned2l.dsmn * ddrif.deta^2, ned2l.dvarSD * dsewn.deta^2,
ned2l.dthe * dThe.deta^2 , ned2l.dphi * dPhi.deta^2)
wz <- wz[, interleave.VGAM( NOS * M1 , M1 = M1)]
}
wz[which(wz[, 1:(NOS * M1)] < .Machine$double.eps)] <-
.Machine$double.base^(0.75)
wz
}), list( .var.arg = var.arg , .nOrder = nOrder,
.ARord = ARord , .MAord = MAord ,
.realM1 = realM1 , .nodrift = nodrift ,
.idrift = idrift ,
.iARcoeff = iARcoeff , .iMAcoeff = iMAcoeff ,
.ivar = ivar, .isd = isd ,
.type.EIM = type.EIM )))
) # End of alo
alo
} # End of ARMAXff
## ARMAXff control function.
ARMAXff.control <- function(save.weights = TRUE,
summary.HDEtest = FALSE,...) {
#criterion <- "loglikelihood"
list(save.weights = save.weights,
summary.HDEtest = summary.HDEtest,...)
}
quick.check.coeffs <- function(arc, mac, NOS, arOrd, maOrd) {
if (length(arc) && (length(arc) != arOrd * NOS))
stop("Conflicting number of initial values for the AR component.")
if ( length(mac) && ( length(mac) != maOrd * NOS) )
stop("Conflicting number of initial values for the MA component.")
if (length(arc)) {
arc <- matrix(arc, NOS, arOrd, byrow = TRUE)
initChecks <- apply(arc, 1, function(x) {
checkTS.VGAMextra(thetaEst = x,
tsclass = "AR",
NofS = 1,
chOrder = arOrd,
pRoots = FALSE,
retmod = TRUE)
})
if(any( initChecks < 1 + 1e-8))
warning("Initial values for the AR coefficients for response ",
" do not belong to a stationary model.")
}
if (length(mac)) {
mac <- matrix(mac, NOS, maOrd, byrow = TRUE)
initChecks <- apply(mac, 1, function(x) {
checkTS.VGAMextra(thetaEst = x, tsclass = "MA",
NofS = 1, chOrder = maOrd,
pRoots = FALSE, retmod = TRUE)
})
if(any( initChecks < 1 + 1e-8))
warning("Initial values for the MA coefficients for response ",
" do not belong to an invertible model.")
}
invisible(NULL)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.