dynrq.R
dynrq <- function (formula, tau = 0.5, data, subset, weights, na.action, method = "br",
contrasts = NULL, start = NULL, end = NULL, ...)
{
require(zoo)
Zenv <- new.env(parent = environment(formula))
assign("dynformula", function(x) structure(x, class = unique(c("dynformula",
oldClass(x)))), envir = Zenv)
assign("L", function(x, k = 1) {
if (length(k) > 1) {
rval <- lapply(k, function(i) lag(x, k = -i))
rval <- if (inherits(x, "ts"))
do.call("ts.intersect", rval)
else do.call("merge", c(rval, list(all = FALSE)))
colnames(rval) <- k
}
else {
rval <- lag(x, k = -k)
}
return(rval)
}, envir = Zenv)
assign("d", function(x, lag = 1) diff(x, lag = lag), envir = Zenv)
assign("season", function(x, ref = NULL) {
freq <- frequency(x)
stopifnot(freq > 1 && identical(all.equal(freq, round(freq)),
TRUE))
freq <- ofreq <- round(freq)
freq <- if (freq == 12)
month.abb
else if (freq == 4)
paste("Q", 1:4, sep = "")
else 1:freq
rval <- factor(coredata(cycle(x)), labels = freq)
if (!is.null(ref))
rval <- relevel(rval, ref = ref)
rval <- zoo(rval, index(x), ofreq)
return(rval)
}, envir = Zenv)
assign("trend", function(x, scale = TRUE) {
freq <- ofreq <- if (inherits(x, "ts"))
frequency(x)
else attr(x, "frequency")
if (is.null(freq) | !scale)
freq <- 1
stopifnot(freq >= 1 && identical(all.equal(freq, round(freq)),
TRUE))
freq <- round(freq)
rval <- zoo(seq_along(index(x))/freq, index(x), frequency = ofreq)
return(rval)
}, envir = Zenv)
assign("harmon", function(x, order = 1) {
freq <- frequency(x)
stopifnot(freq > 1 && identical(all.equal(freq, round(freq)),
TRUE))
freq <- round(freq)
order <- round(order)
stopifnot(order <= freq/2)
rval <- outer(2 * pi * index(x), 1:order)
rval <- cbind(apply(rval, 2, cos), apply(rval, 2, sin))
colnames(rval) <- if (order == 1) {
c("cos", "sin")
}
else {
c(paste("cos", 1:order, sep = ""), paste("sin", 1:order,
sep = ""))
}
if ((2 * order) == freq)
rval <- rval[, -(2 * order)]
return(rval)
}, envir = Zenv)
assign("model.frame.dynformula", function(formula, data = NULL,
subset = NULL, na.action = na.omit, drop.unused.levels = FALSE,
xlev = NULL, ...) {
if (is.null(data))
data <- parent.frame()
if (!is.list(data))
data <- as.list(data)
args <- as.list(attr(terms(formula), "variables"))[-1]
args$retclass <- "list"
args$all <- FALSE
formula <- terms(formula)
attr(formula, "predvars") <- as.call(append(merge.zoo,
args))
attr(formula, "predvars")[[1]] <- as.name("merge.zoo")
NextMethod("model.frame", formula = formula)
}, envir = Zenv)
if (missing(data))
data <- Zenv
orig.class <- if (is.data.frame(data) || is.environment(data))
class(eval(attr(terms(formula), "variables")[[2]], data,
Zenv))
else class(data)
cl <- match.call()
mf <- match.call(expand.dots = FALSE)
m <- match(c("formula", "data", "subset", "weights", "na.action"), names(mf), 0)
mf <- mf[c(1, m)]
mf$drop.unused.levels <- TRUE
mf[[1]] <- as.name("model.frame")
mf[[2]] <- as.call(list(as.name("dynformula"), mf[[2]]))
mf <- eval(mf, envir = Zenv)
mfna <- attr(mf, "na.action")
if (length(index(mf[, 1])) > nrow(mf)) {
for (i in 1:NCOL(mf)) attr(mf[, i], "index") <- attr(mf[,
i], "index")[-as.vector(mfna)]
}
is.zoofactor <- function(x) !is.null(attr(x, "oclass")) &&
attr(x, "oclass") == "factor"
for (i in 1:NCOL(mf)) if (is.zoofactor(mf[, i]))
mf[, i] <- coredata(mf[, i])
mf1 <- mf[, 1]
start <- if (is.null(start))
1
else {
if (length(start) > 1)
start <- start[1] + (start[2] - 1)/frequency(mf1)
start <- min(which(index(mf1) >= start))
}
end <- if (is.null(end))
length(mf1)
else {
if (length(end) > 1)
end <- end[1] + (end[2] - 1)/frequency(mf1)
end <- max(which(index(mf1) <= end))
}
if (end < start) {
warning("empty model frame specified")
mf <- head(mf, 0)
mf1 <- head(mf1, 0)
}
else {
mf <- mf[start:end, , drop = FALSE]
mf1 <- mf1[start:end]
if (!is.null(mfna))
attr(mf, "na.action") <- mfna[as.vector(mfna) >=
start & as.vector(mfna) <= end]
}
if ("ts" %in% orig.class && is.regular(mf1, strict = TRUE)) {
for (i in 1:ncol(mf)) if (!is.factor(mf[, i]))
mf[, i] <- zoo:::as.ts.zoo(mf[, i])
}
if (all(orig.class == "numeric")) {
for (i in 1:ncol(mf)) if (!is.factor(mf[, i]))
mf[, i] <- as.vector(mf[, i])
}
rownames(mf) <- index2char(index(mf1), frequency(mf1))
mt <- attr(mf, "terms")
attr(mt, "predvars") <- NULL
attr(mt, "dataClasses") <- NULL
Y <- model.response(mf, "numeric")
weights <- as.vector(model.weights(mf))
if (is.empty.model(mt)) {
X <- NULL
rval <- list(coefficients = numeric(0), residuals = Y,
fitted.values = 0 * Y, weights = weights, rank = 0, df.residual = length(Y))
}
else {
Rho <- function(u, tau) u * (tau - (u < 0))
eps <- .Machine$double.eps^(2/3)
X <- model.matrix(mt, mf, contrasts)
if (length(tau) > 1) {
if (any(tau < -eps) || any(tau > 1 + eps))
stop("invalid tau: taus should be >= 0 and <= 1")
coef <- matrix(0, ncol(X), length(tau))
rho <- rep(0, length(tau))
fitted <- resid <- matrix(0, nrow(X), length(tau))
for (i in 1:length(tau)) {
z <- {
if (length(weights))
rq.wfit(X, Y, tau = tau[i], weights, method, ...)
else rq.fit(X, Y, tau = tau[i], method, ...)
}
coef[, i] <- z$coefficients
resid[, i] <- z$residuals
rho[i] <- sum(Rho(z$residuals, tau[i]))
fitted[, i] <- Y - z$residuals
}
taulabs <- paste("tau=", format(round(tau, 3)))
dimnames(coef) <- list(dimnames(X)[[2]], taulabs)
dimnames(resid) <- list(dimnames(X)[[1]], taulabs)
rval <- z
rval$coefficients <- coef
rval$residuals <- resid
rval$fitted.values <- fitted
class(rval) <- c("dynrqs", "rqs")
}
else {
rval <- {
if(length(weights))
rq.wfit(X, Y, tau = tau, weights, method, ...)
else rq.fit(X, Y, tau = tau, method, ...)
}
dimnames(rval$residuals) <- list(dimnames(X)[[1]], NULL)
rho <- sum(Rho(rval$residuals, tau))
class(rval) <- "rq"
class(rval) <- c("dynrq", "rq")
}
}
rval$na.action <- attr(mf, "na.action")
rval$contrasts <- attr(X, "contrasts")
rval$xlevels <- .getXlevels(mt, mf)
rval$call <- cl
rval$tau <- tau
rval$terms <- mt
rval$model <- mf
rval$index <- index(mf1)
rval$frequency <- frequency(mf1)
rval$residuals <- drop(rval$residuals)
rval$X <- X
rval$y <- Y
rval$rho <- rho
rval$method <- method
rval$fitted.values <- drop(rval$fitted.values)
return(rval)
}
index.dynrq <- time.dynrq <- function(x, ...) {
x$index
}
start.dynrq <- function(x, ...) {
start(x$residuals)
}
end.dynrq <- function(x, ...) {
end(x$residuals)
}
print.dynrq <- function(x, ...) {
rx <- residuals(x)
cat(paste("\nDynamic quantile regression \"", class(rx)[1], "\" data:\n", sep = ""))
cat(paste("Start = ", index2char(index(rx)[1], x$frequency),
", End = ", index2char(index(rx)[length(rx)], x$frequency), "\n", sep = ""))
NextMethod()
}
print.dynrqs <- function(x, ...) {
rx <- residuals(x)
ix <- dimnames(rx)[[1]]
cat(paste("\nDynamic quantile regression \"", class(rx)[1], "\" data:\n", sep = ""))
cat(paste("Start = ", index2char(ix[1], x$frequency),
", End = ", index2char(ix[length(ix)], x$frequency), "\n", sep = ""))
NextMethod()
}
summary.dynrq <- function(object, vcov. = NULL, df = NULL, ...) {
rval <- NextMethod()
rval$frequency <- object$frequency
class(rval) <- c("summary.dynrq", class(rval))
return(rval)
}
print.summary.dynrq <- function(x, ...) {
rx <- residuals(x)
x$residuals <- coredata(x$residuals)
cat(paste("\nDynamic quantile regression \"", class(rx)[1], "\" data:\n", sep = ""))
cat(paste("Start = ", index2char(index(rx)[1], x$frequency),
", End = ", index2char(index(rx)[length(rx)], x$frequency), "\n", sep = ""))
NextMethod()
}
print.summary.dynrqs <- function(x, ...) {
lapply(x, print.summary.dynrq)
}