Nothing
ar.egls <-
function (x, R, order.max, na.action = na.fail, series = NULL,
...)
{
R = sync(R)
aic = F
demean = T
intercept = F
n <- nrow(x)
k <- ncol(x)
if (is.null(series))
series <- deparse(substitute(x))
rescale <- TRUE
ists <- is.ts(x)
x <- na.action(as.ts(x))
xfreq <- frequency(x)
if (any(is.na(x)))
stop("NAs in x")
if (ists)
xtsp <- tsp(x)
x <- as.matrix(x)
if (!is.numeric(x))
stop("`x' must be numeric")
n.used <- nrow(x)
nser <- ncol(x)
if (rescale) {
sc <- sqrt(drop(apply(x, 2, var)))
x <- x/rep(sc, rep(n.used, nser))
}
else sc <- rep(1, nser)
order.max <- if (is.null(order.max))
floor(10 * log10(n.used))
else round(order.max)
if (order.max < 0)
stop("order.max must be >= 0")
if (aic)
order.min <- 0
else order.min <- order.max
A <- vector("list", order.max - order.min + 1)
varE <- vector("list", order.max - order.min + 1)
seA <- vector("list", order.max - order.min + 1)
aic <- rep(Inf, order.max - order.min + 1)
det <- function(x) {
prod(diag(qr(x)$qr)) * (-1)^(ncol(x) - 1)
}
if (demean) {
xm <- colMeans(x)
x <- sweep(x, 2, xm)
}
else xm <- rep(0, nser)
for (m in order.min:order.max) {
y <- embed(x, m + 1)
if (intercept) {
if (m > 0)
X <- cbind(rep(1, nrow(y)), y[, (nser + 1):ncol(y)])
else X <- as.matrix(rep(1, nrow(y)))
}
else {
if (m > 0)
X <- y[, (nser + 1):ncol(y)]
else X <- matrix(0, nrow(y), 0)
}
Y <- t(y[, 1:nser])
N <- ncol(Y)
XX <- t(X) %*% X
rank <- qr(XX)$rank
if (rank != nrow(XX)) {
warning(paste("Model order", m))
warning("Singularities in the computation of the projection matrix")
warning(paste("Results are only valid up to model order",
m - 1))
break
}
P <- if (ncol(XX) > 0)
solve(XX)
else XX
A[[m - order.min + 1]] <- Y %*% X %*% P
YH <- A[[m - order.min + 1]] %*% t(X)
E <- (Y - YH)
varE[[m - order.min + 1]] <- E %*% t(E)/N
aic[m - order.min + 1] <- n.used * log(det(varE[[m -
order.min + 1]])) + 2 * nser * (nser * m + intercept)
}
m <- which(aic == min(aic)) + order.min - 1
y <- embed(x, m + 1)
sigmau <- (1/(n - k * order.max - 1)) * E %*% t(E)
for (i in 1:2) {
T1 <- ginv(t(R) %*% (kronecker((XX), ginv(sigmau)) %*%
R))
T2 <- t(R) %*% (kronecker(id(k * order.max), ginv(sigmau)))
gam <- vec(A[[m - order.min + 1]]) + T1 %*% T2 %*% vec(E %*%
(X))
beta <- R %*% gam
AA <- unvec(beta, k, order.max)
YH <- AA %*% t(X)
E <- (Y - YH)
}
sig <- E %*% t(E)/N
varA <- R %*% ginv(t(R) %*% (kronecker((XX/n), ginv(sig))) %*%
R) %*% t(R)/n
seA[[m - order.min + 1]] <- if (ncol(varA) > 0)
sqrt(diag(varA))
else numeric(0)
if (intercept) {
xint <- AA[, 1]
ar <- AA[, -1]
if (m > 0)
X <- cbind(rep(1, nrow(y)), y[, (nser + 1):ncol(y)])
else X <- as.matrix(rep(1, nrow(y)))
}
else {
if (m > 0)
X <- y[, (nser + 1):ncol(y)]
else X <- matrix(0, nrow(y), 0)
xint <- NULL
ar <- AA
}
Y <- t(y[, 1:nser, drop = FALSE])
YH <- AA %*% t(X)
E <- drop(rbind(matrix(NA, m, nser), t(Y - YH)))
aic <- aic - min(aic)
names(aic) <- order.min:order.max
dim(ar) <- c(nser, nser, m)
ar <- aperm(ar, c(3, 1, 2))
ses <- seA[[m - order.min + 1]]
if (intercept) {
sem <- ses[1:nser]
ses <- ses[-(1:nser)]
}
else sem <- rep(0, nser)
dim(ses) <- c(nser, nser, m)
ses <- aperm(ses, c(3, 1, 2))
var.pred <- varE[[m - order.min + 1]]
if (nser > 1) {
snames <- colnames(x)
dimnames(ses) <- dimnames(ar) <- list(seq(length = m),
snames, snames)
dimnames(var.pred) <- list(snames, snames)
names(sem) <- colnames(E)
snames <- colnames(E)
}
if (ists) {
attr(E, "tsp") <- xtsp
attr(E, "class") <- "ts"
}
if (rescale) {
xm <- xm * sc
if (!is.null(xint))
xint <- xint * sc
aa <- outer(sc, 1/sc)
if (nser > 1 && m > 0)
for (i in 1:m) ar[i, , ] <- ar[i, , ] * aa
var.pred <- var.pred * outer(sc, sc)
E <- E * rep(sc, rep(NROW(E), nser))
sem <- sem * sc
if (m > 0)
for (i in 1:m) ses[i, , ] <- ses[i, , ] * aa
}
res <- list(order = m, ar = ar, var.pred = var.pred, x.mean = xm,
x.intercept = xint, aic = aic, n.used = n.used, order.max = order.max,
partialacf = NULL, resid = E, method = "EGLS", series = series,
frequency = xfreq, call = match.call(), asy.se.coef = list(x.mean = sem,
ar = drop(ses)))
class(res) <- "ar"
res
}
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.