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)
## ARIMAXff control function.
ARIMAXff.control <- function(save.weights = TRUE,
summary.HDEtest = FALSE,
...) {
list(save.weights = save.weights,
summary.HDEtest = summary.HDEtest,...)
}
ARIMAXff <-
function(order = c(1, 1, 0),
zero = c("ARcoeff", "MAcoeff"),
diffCovs = TRUE,
xLag = 0,
include.current = FALSE,
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 (xLag < 0)
stop("Wrong input for argument 'xLag'")
if ( !Is.Numeric(order, length.arg = 3) )
stop("Invalid 'order'. Should be a vector (p, d, q), length 3.")
ARord2 <- order[1]
dOrd <- order[2]
MAord2 <- order[3]; rm(order)
flagAR <- (ARord2 == 0)
flag.2 <- ((ARord2 == 0) && (MAord2 == 0))
if (flag.2)
stop("Seems like a random walk is modelled.",
" Refer to uninormal() for this. ")
# if ((ARord2 == 0) && (MAord2 == 0))
# stop("Wrong input for order. Enter the AR order and/or MA order.")
ARord <- if (ARord2 == 0) 1 else ARord2
MAord <- if (MAord2 == 0) 1 else MAord2
if (ARord2 == 0) {
rem.lat <- grep("ARcoeff", zero)
zero <- zero[-rem.lat]
}
if (MAord2 == 0) {
rem.lat <- grep("MAcoeff", zero)
zero <- zero[-rem.lat]
}
if (!is.logical(diffCovs))
stop("Bad input for argument 'diffCovs'")
if (!is.logical(include.current))
stop("Bad input for argument 'include.current'")
if (!dOrd)
stop("Refer to ARMAXff(), ARXff() or MAXff(). ",
"Only the ARIMA--class handled", "\n ",
"by this family function")
nOrder <- ARord2 + MAord2
realM1 <- 2 + ARord2 + MAord2 - nodrift
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 <- if (flag.2) "approximate" else
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("VGLTSMs: ", nOrder + 1 + !nodrift,
"-parameter ARIMAX model ",
"of order",
"(", ARord2, ",", dOrd, ",", MAord2, ").",
"\n", "\n",
"LINKS: ",
blurb.vec2, "\n",
"Model: ", "\n" , " Y[t] = ", if (nodrift) NULL else
if (!ARord2) "mu + " else "drift + ",
" B^T * X_t + ",
if (ARord2) {
if (ARord2 == 1) "(ARcoeff1) * Y[t - 1] + " else
c("(ARcoeff1) * Y[t - 1] +... + (ARcoeff", ARord2,")",
" * Y[t -", ARord2,"] + ")
},
if (MAord2) {
if (MAord2 == 1) "(MAcoeff1) * e[t - 1]" else
c("(MAcoeff1) * e[t - 1] + ... + (MAcoeff", MAord,")",
" * e[t - ", MAord,"]")
},
if ( !MAord2 ) " e[t] " else
" + e[t].", "\n"),
#if ( !ARord2 ) "E[Yt] = mu" else
#"E[Yt] = drift / (1 - SUM(ARcoeff))"),
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 ) )),
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 ,
dOrd = .dOrd ,
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 , .dOrd = dOrd ,
.nodrift = nodrift , .type.EIM = type.EIM ,
.zero = zero ))),
first = eval(substitute(expression({
extra$y.r <- y
y <- cbind(y)
n <- NROW(y)
NOS <- NCOL(y)
mat.save <- matrix(0, n, .dOrd)
if (NOS != 1)
stop("Currently this VGLSTM family function handles",
"univariate responses only.")
if ((NCOL(x) == 1) && (colnames(x) == "(Intercept)") && ( .xLag > 0))
stop("Set xLag = 0. No covariates entered.")
y.diff <- diff(y, differences = .dOrd )
for (ii in 1:NOS) {
yaux1 <- y[, ii , drop = FALSE]
for (jj in 1:( .dOrd )) {
mat.save[c(1:(n - jj)), jj] <- diff(yaux1)
yaux1 <- mat.save[c(1:(n - jj)), jj, drop = FALSE]
}
}
y <- extra$y <- y.diff
extra$mat.save <- mat.save
difxs <- .diffCovs
if ((NCOL(x) == 1) && (colnames(x) == "(Intercept)") && difxs) {
difxs <- FALSE
warning("No covariates entered. 'diffCovs' set to FALSE")
}
x.int <- x[, -1, drop = FALSE]
if (difxs) {
x.diff <- apply(cbind(x.int), 2, function(x) {
diff(x, differences = .dOrd )
})
colnames(x.diff) <- paste("Diff", colnames(x.int), sep = "")
} else{
x.diff <- x.int[ c(1:NROW(y.diff)), , drop = FALSE]
}
x.matrix <- cbind(1, x.diff)
colnames(x.matrix) <- c("(Intercept)", colnames(x.diff))
x.temp <- temp3 <- NULL
if ( .xLag ) {
for (jj in 1:NCOL(x.diff)) {
temp1 <- WN.lags(cbind(x.diff[, jj]), .xLag )
temp2 <- paste(colnames(x.diff)[jj], "Lag", sep = "")
colnames(temp1) <- paste(temp2, 1:( .xLag ), sep = "")
x.temp <- cbind(x.temp, temp1)
}
rm(temp1, temp2)
}
if ( .inc.c ) {
temp3 <- x[1:NROW(x.matrix), -1, drop = FALSE]
}
x.matrix <- cbind(x.matrix, x.temp, temp3)
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
w <- w[1:NROW(y)]
}), list( .dOrd = dOrd , .xLag = xLag , .diffCovs = diffCovs ,
.inc.c = include.current ))),
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
y.r <- extra$y
dOrd <- .dOrd
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,
if (!( .ARord2 )) NULL else ar.n,
if (!( .MAord2 )) NULL else 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 ),
if (!( .ARord2 )) NULL else
namesof(ar.n, link = .lARcoeff , earg = .eARcoeff , tag = FALSE),
if (!( .MAord2 )) NULL else
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 <- scale(y, center = TRUE, scale = FALSE) # y 180519
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 = y.sc[, rsp, drop = FALSE],
lags = .ARord + .MAord + 5))
to.fit2 <- lm( to.fit[, 1, drop = FALSE] ~
to.fit[, -1, drop = FALSE] - 1)
extra$res <- res2.mat[, rsp] <- cbind(residuals(to.fit2))
if (TRUE) {
to.fit2 <- cbind(y.sc[, rsp, drop = FALSE],
if ( .ARord2 )
WN.lags(y = cbind(y.sc[, rsp, drop = FALSE]),
lags = .ARord ) else NULL,
if ( .MAord2 )
WN.lags(y = res2.mat[, rsp, drop = FALSE],
lags = .MAord ) else NULL)
to.fit2 <- lm( to.fit2[, 1, drop = FALSE] ~
to.fit2[, -1, drop = FALSE] - 1)
}
initCoe <- coef(to.fit2)
extra$res <- res2.mat[, rsp] <- cbind(residuals(to.fit2))
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[-(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])
}
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 ),
if (!( .ARord2 )) NULL else
theta2eta(matrix(init.AR, n, NOS * ( .ARord )),
.lARcoeff , .eARcoeff ),
if (!( .MAord2 )) NULL else
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 , .ARord2 = ARord2 ,
.MAord2 = MAord2 , .dOrd = dOrd ))),
linkinv = eval(substitute(function(eta, extra = NULL){
M1 <- .realM1; n <- nrow(eta)
NOS <- ncol(eta)/M1
mat.save <- extra$mat.save
dOrd < .dOrd; NNro <- 1
nn <- NROW(mat.save)
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(0, dim = c(n, NOS, .ARord ))
mas <- array(0, dim = c(n, NOS, .MAord ))
if ( .ARord2 )
for (ai in 1:( .ARord ))
ars[, , ai] <- eta2theta(eta[, M1 * (1:NOS) -
( .ARord + .MAord ) + ai, drop = FALSE],
link = .lARcoeff , earg = .eARcoeff )
if ( .MAord2 )
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] <- y.tt <- 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]
extra$res <- extra$y - y.est
sum.mat <- cbind(extra$y.r, extra$mat.save[, -c(dOrd), drop = FALSE],
c(y.est, rep(0, dOrd )))
t.rev <- rev(1:NCOL(sum.mat))
for (ii in 1:dOrd) {
yaux1 <- sum.mat[c(1:(nn - dOrd + ii - 1)), t.rev[ii], drop = FALSE]
yaux2 <- sum.mat[c(1:(nn - dOrd + ii - 1)),
t.rev[ii] - 1, drop = FALSE]
yaux3 <- yaux1 + yaux2
yaux3 <- rbind(sum.mat[1, t.rev[ii] - 1, drop = FALSE], yaux3)
sum.mat[ , t.rev[ii] - 1] <- c(yaux3, rep(0, dOrd - ii))
}
y.est <- sum.mat[, 1]
y.est[c(1:NNro)] <- (NNro + 0.10)*y.est[c(2:(NNro + 1))]
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 , .dOrd = dOrd ,
.ARord2 = ARord2 , .MAord2 = MAord2 ))),
last = eval(substitute(expression({
M1 <- .realM1
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 <- mas <- NULL
prov.names <- c(if ( .nodrift ) NULL else drift.n , varsd.n)
if ( .ARord2 ) {
ars <- array(NA, dim = c(n, NOS, .ARord ))
for (ai in 1:( .ARord ))
ars[, , ai ] <-
eta2theta(eta[, M1*(1:NOS) - (.MAord + .ARord) + ai,
drop = FALSE],
link = .lARcoeff , earg = .eARcoeff )
for (kk in 1:( .ARord ))
prov.names <- c(prov.names, ar.n[[kk]])
}
if ( .MAord2 ) {
mas <- array(NA, dim = c(n, NOS, .MAord ))
for ( mi in 1:( .MAord ))
mas[, , mi ] <-
eta2theta(eta[, M1*(1:NOS) - ( .MAord ) + mi , drop = FALSE],
link = .lMAcoeff , earg = .eMAcoeff )
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
if ( .ARord2 )
for (kk in 1:( .ARord )) {
misc$link[ M1 * jj - .nOrder + kk ] <- .lARcoeff
misc$earg[[ M1 * jj - .nOrder + kk ]] <- .eARcoeff
}
if ( .MAord2 )
for( ll in 1:( .MAord )) {
misc$link[ M1 * jj - .MAord + ll ] <- .lMAcoeff
misc$earg[[ M1 * jj - .MAord + ll ]] <- .eMAcoeff
}
}
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 <- .flagAR
misc$residuals <- cbind(c(rep(0, .dOrd ), extra$res))
misc$theta.names <- parameters.names
misc$multipleResponses <- TRUE
misc$control$wzepsilon <- control$wzepsilon
fit$prior.weights <- w #20180315
y <- c(extra$y.r)
# Checks removed on March/2/2016. Now, refer to 'summary'.
if (!( .noChecks ) && !( .flag.2 )) {
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, .flagAR = flagAR ,
.type.likelihood = type.likelihood ,
.ARord2 = ARord2 , .MAord2 = MAord2,
.flag.2 = flag.2 , .dOrd = dOrd ))),
loglikelihood =
eval(substitute(function(mu, y, w, residuals = FALSE,
eta, extra = NULL, summation = TRUE) {
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 )
y.est <- drifts + log(1) * ifelse(M1 > 1, 0.15, 0.15) * extra$res
ars <- mas <- 0
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 )
}
if ( .ARord2 ) {
#ars <- matrix(NA_real_, n, .ARord )
ars <- eta2theta(eta[, (M1 - ( .ARord2 + .MAord2 )
+ 1):(M1 - .MAord2 ), drop = FALSE],
link = .lARcoeff , earg = .eARcoeff )
y.est <- y.est + cbind(rowSums(ars * WN.lags(y = cbind(extra$y),
lags = .ARord2 )))
}
if ( .MAord2 ) {
#ars <- matrix(NA_real_, n, .ARord )
mas <- eta2theta(eta[, (M1 - ( .MAord2 ) + 1):M1,
drop = FALSE],
link = .lMAcoeff , earg = .eMAcoeff )
y.est <- y.est + cbind(rowSums(mas * WN.lags(y = cbind(extra$res),
lags = .MAord2 )))
}
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 ,
.ARord2 = ARord2 , .MAord2 = MAord2 ))),
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(0, dim = c(n, NOS, .ARord ))
mas <- array(0, 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", "ARIMAff"),
vfamily = c("ARIMAXff", "ARMAvgltsmff"),
deriv = eval(substitute(expression({
M1 <- .realM1; n <- nrow(eta)
NOS <- ncol(eta)/M1
M <- M1 * NOS
y <- extra$y
drifts <- if ( .nodrift ) matrix(0.0, nrow = n, ncol = NOS) else
eta2theta(eta[, M1*(1:NOS) - .nOrder - 1, drop = FALSE],
link = .ldrift , earg = .edrift )
y.est <- drifts
ars <- mas <- 0
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 )
}
if ( .ARord2 ) {
#ars <- matrix(NA_real_, n, .ARord )
ars <- eta2theta(eta[, (M1 - ( .ARord2 + .MAord2 )
+ 1):(M1 - .MAord2 ), drop = FALSE],
link = .lARcoeff , earg = .eARcoeff )
y.est <- y.est + cbind(rowSums(ars * WN.lags(y = cbind(extra$y),
lags = .ARord2 )))
}
if ( .MAord2 ) {
#ars <- matrix(NA_real_, n, .ARord )
mas <- eta2theta(eta[, (M1 - ( .MAord2 ) + 1):M1,
drop = FALSE],
link = .lMAcoeff , earg = .eMAcoeff )
y.est <- y.est + cbind(rowSums(mas * WN.lags(y = cbind(extra$res),
lags = .MAord2 )))
}
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,
if (!( .ARord2 )) NULL else
matrix( dl.dThe[, 1, ] * dThe.deta, n, NOS * ( .ARord )),
if (!( .MAord2 )) NULL else
matrix( dl.dPhi[, 1, ] * dPhi.deta, n, NOS * ( .MAord )))
myderiv <- myderiv[, interleave.VGAM(( .realM1 )*NOS , M1 = .realM1)]
colnames(myderiv) <- NULL
# Required at @weights
ars <- array(ars, dim = c(n, NOS, .ARord2 ))
mas <- array(mas, dim = c(n, NOS, .MAord2 ))
dThe.deta <- dtheta.deta(ars , .lARcoeff , .eARcoeff)
dPhi.deta <- dtheta.deta(mas , .lMAcoeff , .eMAcoeff)
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 , .ARord2 = ARord2 ,
.MAord2 = MAord2 ))),
weight = eval(substitute(expression({
if (( .type.EIM ) == "exact" ) {
final.eim <- numeric(0)
for (ii in 1:NOS) {
if ( ( .ARord2 ) && !( .MAord2 ) )
exact.eim <- ARpEIM.G2(y = cbind(y[, ii]),
drift = NULL,
sdError = cbind(sigWN[, ii]),
ARpcoeff = cbind(ars[, ii, ]),
var.arg = .var.arg ,
order = .ARord2 ,
nodrift = .nodrift )
if ( !( .ARord2 ) && ( .MAord2 ) )
exact.eim <- MAqEIM.G2(y = cbind(y[, ii]),
mean = NULL,
sdError = cbind(sigWN[, ii]),
MAqcoeff = cbind(mas[, ii, ]),
var.arg = .var.arg ,
order = .MAord2 ,
nomean = .nodrift )
if ( ( .ARord2 ) && ( .MAord2 ) )
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, 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,
if (!( .ARord2 )) NULL else ned2l.dthe * dThe.deta^2 ,
if (!( .MAord2 )) NULL else ned2l.dphi * dPhi.deta^2)
wz <- wz[, interleave.VGAM( NOS * M1 , M1 = M1)]
}
wz[which(wz[, 1:(NOS * M1)] < .Machine$double.eps)] <-
.Machine$double.eps^(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 ,
.ARord2 = ARord2 , .MAord2 = MAord2 )))
) # End of alo
alo
}
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.