Nothing
MARSSresiduals.tT <- function(object, Harvey = FALSE, normalize = FALSE, silent = FALSE, fun.kf = c("MARSSkfas", "MARSSkfss")) {
# These are the residuals and their variance conditioned on all the data
# Harvey=TRUE uses Harvey et al (1998) algorithm to compute these
# Harvey=FALSE uses the straight smoother output
# model.residuals y(t|yT)-Zx(t|yT)-a
# state.residuals x(t|yT)-Bx(t-1|yT)-u
# var.residuals variance of above conditioned on y(1)
# for missing values, Harvey=TRUE returns 0 for var for y_i missing and Harvey=FALSE returns R + Z VtT t(Z)
# Note, I think there is a problem with the Harvey algorithm when the variance of the state residuals (Q)
# is non-diagonal and there are missing values; it can become non-invertible
######################################
# Set up variables
MLEobj <- object
if (missing(fun.kf)) {
fun.kf <- MLEobj$fun.kf
} else {
MLEobj$fun.kf <- fun.kf
# to ensure that MARSSkf, MARSShatyt and fitted() use the spec'd fun
}
if (fun.kf == "MARSSkfas" && Harvey == TRUE) stop("MARSSresiduals.tT: Harvey=TRUE requires the Kalman gain thus MARSSkfss must be used. Pass in fun.kf='MARSSkfss'.\n", call. = FALSE)
model.dims <- attr(MLEobj$marss, "model.dims")
TT <- model.dims[["x"]][2]
m <- model.dims[["x"]][1]
n <- model.dims[["y"]][1]
y <- MLEobj$marss$data
# set up holders
et <- st.et <- mar.st.et <- bchol.st.et <- matrix(0, n + m, TT)
var.et <- array(0, dim = c(n + m, n + m, TT))
msg <- NULL
#### list of time-varying parameters
time.varying <- is.timevarying(MLEobj)
kf <- MARSSkf(MLEobj)
Ey <- MARSShatyt(MLEobj)
Rt <- parmat(MLEobj, "R", t = 1)$R # returns matrix
Ht <- parmat(MLEobj, "H", t = 1)$H
Rt <- Ht %*% tcrossprod(Rt, Ht)
Zt <- parmat(MLEobj, "Z", t = 1)$Z
Qtp <- parmat(MLEobj, "Q", t = 2)$Q
Gtp <- parmat(MLEobj, "G", t = 2)$G
Qtp <- Gtp %*% tcrossprod(Qtp, Gtp)
Btp <- parmat(MLEobj, "B", t = 2)$B
utp <- parmat(MLEobj, "U", t = 2)$U
######################################
######################################
# Compute residuals via Holmes algorithm
if (!Harvey) {
# model.et will be 0 where no data E(y)-modeled(y)
model.et <- Ey$ytT - fitted(MLEobj, type = "ytT", output = "matrix") # model residuals
et[1:n, ] <- model.et
for (t in 1:TT) {
# model residuals
if (time.varying$R) Rt <- parmat(MLEobj, "R", t = t)$R # returns matrix
if (time.varying$H) Ht <- parmat(MLEobj, "H", t = t)$H
if (time.varying$R || time.varying$H) Rt <- Ht %*% tcrossprod(Rt, Ht)
if (time.varying$Z) Zt <- parmat(MLEobj, "Z", t = t)$Z
# model.et defined outside for loop
# compute the variance of the residuals and state.et
St <- Ey$yxtT[, , t] - tcrossprod(Ey$ytT[, t, drop = FALSE], kf$xtT[, t, drop = FALSE])
tmpvar.et <- Rt - Zt %*% tcrossprod(kf$VtT[, , t], Zt) + tcrossprod(St, Zt) + tcrossprod(Zt, St)
if (t < TT) { # fill in var.et for t (model resid)
if (time.varying$Q) Qtp <- parmat(MLEobj, "Q", t = t + 1)$Q
if (time.varying$G) Gtp <- parmat(MLEobj, "G", t = t + 1)$G
if (time.varying$Q || time.varying$G) Qtp <- Gtp %*% tcrossprod(Qtp, Gtp)
if (time.varying$B) Btp <- parmat(MLEobj, "B", t = t + 1)$B
if (time.varying$U) utp <- parmat(MLEobj, "U", t = t + 1)$U
Sttp <- Ey$yxttpT[, , t] - tcrossprod(Ey$ytT[, t, drop = FALSE], kf$xtT[, t + 1, drop = FALSE])
cov.et <- tcrossprod(Zt, kf$Vtt1T[, , t + 1]) - Zt %*% tcrossprod(kf$VtT[, , t], Btp) - Sttp + tcrossprod(St, Btp)
tmpvar.state.et <- Qtp - kf$VtT[, , t + 1] - Btp %*% tcrossprod(kf$VtT[, , t], Btp) + tcrossprod(kf$Vtt1T[, , t + 1], Btp) + tcrossprod(Btp, kf$Vtt1T[, , t + 1])
et[(n + 1):(n + m), t] <- kf$xtT[, t + 1] - Btp %*% kf$xtT[, t] - utp
} else {
cov.et <- matrix(0, n, m)
tmpvar.state.et <- matrix(0, m, m)
}
var.et[1:n, , t] <- cbind(tmpvar.et, cov.et)
var.et[(n + 1):(n + m), , t] <- cbind(t(cov.et), tmpvar.state.et)
if (normalize) {
Qpinv <- matrix(0, m, n + m)
Rinv <- matrix(0, n, n + m)
if (t < TT) {
if (time.varying$Q) Qtp <- parmat(MLEobj, "Q", t = t + 1)$Q
Qpinv[, (n + 1):(n + m)] <- psolve(t(pchol(Qtp)))
}
Rinv[, 1:n] <- psolve(t(pchol(Rt)))
RQinv <- rbind(Rinv, Qpinv) # block diag matrix
et[, t] <- RQinv %*% et[, t]
var.et[, , t] <- RQinv %*% tcrossprod(var.et[, , t], RQinv)
}
}
} else {
######################################
# use Harvey algorithm
# Reference page 112-133 in Messy Time Series
# Reference de Jong and Penzer 1998; with model transformed so sigma^2 = 1
# refs in man file
# NOTATION here uses that in Messy Time Series (Harvey et al 1998) with some differences
# MARSS uses Koopman's terminology where w_t = G%*%w, Harvey uses H%*%w
# By definition residual = 0 at time t, where x_0 is defined at t=0 or t=1, **when x_0 is estimated**
# because x_0 is estimated and the max L will be when residual is 0
# x_TT = Bx_{TT-1}+u_{TT-1}+w_{TT-1}, so w_TT cannot be computed (you'd need X_{TT+1}
# so state residual at t=TT is NA
# The state residuals are diff(states) but scaled by t(chol(Q))
rt <- matrix(0, m, TT)
ut <- matrix(0, n, TT)
Nt <- array(0, dim = c(m, m, TT))
Mt <- array(0, dim = c(n, n, TT))
Jt <- matrix(0, m, TT)
vt <- kf$Innov
# If they are time-varying, Q, G and B at t=T will not appear (cancelled out by r_T and N_T = 0).
# Set for t=1. Will update in for loop if time-varying.
pari <- parmat(MLEobj, t = 1)
Qtp <- pari[["Q"]]
Gtp <- pari[["G"]]
Ttp <- pari[["B"]]
Z <- pari[["Z"]] # base, will be modified if missing values
R <- pari[["R"]] # base, will be modified if missing values
Ht <- pari[["H"]]
for (t in seq(TT, 1, -1)) {
# define all the, potential time-varying, parameters
if (t < TT) {
if (time.varying[["Q"]]) Qtp <- parmat(MLEobj, "Q", t = t + 1)$Q
if (time.varying[["G"]]) Gtp <- parmat(MLEobj, "G", t = t + 1)$G
if (time.varying[["B"]]) Ttp <- parmat(MLEobj, "B", t = t + 1)$B
}
# Zt and Rt modified in missing values case below so reset even
# if time-varying
if (time.varying[["Z"]]) {
Zt <- parmat(MLEobj, "Z", t = t)$Z
} else {
Zt <- Z
}
if (time.varying[["R"]]) {
Rt <- parmat(MLEobj, "R", t = t)$R
} else {
Rt <- R
}
if (time.varying[["H"]]) Ht <- parmat(MLEobj, "H", t = t)$H
# implement missing values modifications per Shumway and Stoffer
diag.Rt <- diag(Rt)
Rt[is.na(y[, t]), ] <- 0
Rt[, is.na(y[, t])] <- 0
# diag(Rt)=diag.Rt
Zt[is.na(y[, t]), ] <- 0
# create the m x n+m and n x n+m matrices
Rstar <- matrix(0, n, n + m)
if (normalize) {
Rstar[, 1:n] <- tcrossprod(Ht, pchol(Rt))
} else {
Rstar[, 1:n] <- Ht %*% tcrossprod(Rt, Ht)
}
Qpstar <- matrix(0, m, n + m)
# MARSS uses Koopman's terminology where w_t = G%*%w, Harvey uses H%*%w
if (normalize) {
Qpstar[, (n + 1):(n + m)] <- tcrossprod(Gtp, pchol(Qtp))
} else {
Qpstar[, (n + 1):(n + m)] <- Gtp %*% tcrossprod(Qtp, Gtp)
}
# Don't use kf$Sigma since that has (i,i)=1
Ftinv <- psolve(Zt %*% tcrossprod(kf$Vtt1[, , t], Zt) + Rt)
# Harvey algorithm modified to return non-normalized errors
Kt <- Ttp %*% matrix(kf$Kt[, , t], m, n) # R is dropping the dims so we force it to be mxn
Lt <- Ttp - Kt %*% Zt
Jt <- Qpstar - Kt %*% Rstar
ut[, t] <- Ftinv %*% vt[, t, drop = FALSE] - t(Kt) %*% rt[, t, drop = FALSE]
if (t > 1) rt[, t - 1] <- t(Zt) %*% ut[, t, drop = FALSE] + t(Ttp) %*% rt[, t, drop = FALSE]
# Mt[,,t] = Ftinv + t(Kt)%*%Nt[,,t]%*%Kt #not used
if (t > 1) Nt[, , t - 1] <- t(Zt) %*% Ftinv %*% Zt + t(Lt) %*% Nt[, , t] %*% Lt
et[, t] <- t(Rstar) %*% ut[, t, drop = FALSE] + t(Qpstar) %*% rt[, t, drop = FALSE]
# see deJong and Penzer 1998, page 800, right column, halfway down
var.et[, , t] <- t(Rstar) %*% Ftinv %*% Rstar + t(Jt) %*% Nt[, , t] %*% Jt
}
}
######################################
######################################
# prepare standardized residuals
for (t in 1:TT) {
tmpvar <- sub3D(var.et, t = t)
resids <- et[, t, drop = FALSE]
# don't includ values for resids if there is no residual (no data)
# replace NAs with 0s
is.miss <- c(is.na(y[, t]), rep(FALSE, m))
resids[is.miss] <- 0
tmpvar[abs(tmpvar) < sqrt(.Machine$double.eps)] <- 0
# Marginal
# inverse of diagonal of variance matrix for marginal standardization
# psolve deals with 0s on diagonal
tmpvarinv <- try(psolve(makediag(takediag(tmpvar))), silent = TRUE)
if (inherits(tmpvarinv, "try-error")) {
mar.st.et[, t] <- NA
msg <- c(msg, paste('MARSSresiduals.tT warning: the diagonal matrix of the variance of the residuals at t =", t, "is not invertible. NAs returned for mar.residuals at t =", t, "\n'))
} else { # inverse of the diagonal is ok
mar.st.et[, t] <- sqrt(tmpvarinv) %*% resids
mar.st.et[is.miss, t] <- NA
}
# Block Cholesky
# psolve and pchol deal with 0s on diagonal
tmpchol <- try(pchol(tmpvar[(n + 1):(n + m), (n + 1):(n + m), drop = FALSE]), silent = TRUE)
if (inherits(tmpchol, "try-error")) {
bchol.st.et[(n + 1):(n + m), t] <- NA
msg <- c(msg, paste("MARSSresiduals.tT warning: the chol of the variance of the state residuals at t =", t, "returned errors. NAs returned for bchol.std.residuals at t =", t, ". See MARSSinfo(\"residvarinv\")\n"))
} else {
# chol() returns the upper triangle. We need to lower triangle to t()
tmpcholinv <- try(psolve(t(tmpchol)), silent = TRUE)
if (inherits(tmpcholinv, "try-error")) {
bchol.st.et[(n + 1):(n + m), t] <- NA
msg <- c(msg, paste("MARSSresiduals.tT warning: the variance of the state residuals at t =", t, "is not invertible. NAs returned for bchol.std.residuals at t =", t, ". See MARSSinfo('residvarinv')\n"))
} else {
bchol.st.et[(n + 1):(n + m), t] <- tmpcholinv %*% resids[(n + 1):(n + m), 1, drop = FALSE]
}
}
# Cholesky
# psolve and pchol deal with 0s on diagonal
tmpchol <- try(pchol(tmpvar), silent = TRUE)
if (inherits(tmpchol, "try-error")) {
st.et[, t] <- NA
msg <- c(msg, paste("MARSSresiduals.tT warning: the chol of the variance of the residuals at t =", t, "returned errors. NAs returned for std.residuals at t =", t, ". See MARSSinfo(\"residvarinv\")\n"))
next
}
# chol() returns the upper triangle. We need the lower triangle so t()
tmpcholinv <- try(psolve(t(tmpchol)), silent = TRUE)
if (inherits(tmpcholinv, "try-error")) {
st.et[, t] <- NA
msg <- c(msg, paste("MARSSresiduals.tT warning: the variance of the residuals at t =", t, "is not invertible. NAs returned for std.residuals at t =", t, ". See MARSSinfo('residvarinv')\n"))
next
}
st.et[, t] <- tmpcholinv %*% resids
st.et[is.miss, t] <- NA
}
bchol.st.et[1:n, ] <- st.et[1:n, ] # because the upper right block of the lower tri is 0
######################################
######################################
# NA at last time step
# the state.residual at the last time step is NA because it is x(T+1) - f(x(T)) and T+1 does not exist. For the same reason, the var.residuals at TT will have NAs
et[(n + 1):(n + m), TT] <- NA
var.et[, (n + 1):(n + m), TT] <- NA
var.et[(n + 1):(n + m), , TT] <- NA
st.et[, TT] <- NA
mar.st.et[(n + 1):(n + m), TT] <- NA
bchol.st.et[(n + 1):(n + m), TT] <- NA
if (Harvey == TRUE) {
# Harvey algorithm doesn't calculate var for missing data
for (t in 1:TT) {
is.miss <- c(is.na(y[, t]), rep(FALSE, m))
var.et[is.miss, 1:n, t] <- NA
var.et[1:n, is.miss, t] <- NA
}
}
######################################
######################################
# Variance of missing values conditioned on the data
# et is the expected value of the residuals conditioned on y(1)-the observed data
E.obs.v <- et[1:n, , drop = FALSE]
var.obs.v <- array(0, dim = c(n, n, TT))
for (t in 1:TT) var.obs.v[, , t] <- Ey$OtT[, , t] - tcrossprod(Ey$ytT[, t])
######################################
######################################
# the observed model residuals are data - E(data), so NA for missing data.
model.et <- et[1:n, , drop = FALSE]
model.et[is.na(y)] <- NA
et[1:n, ] <- model.et
######################################
######################################
# add rownames
Y.names <- attr(MLEobj$model, "Y.names")
X.names <- attr(MLEobj$model, "X.names")
rownames(et) <- rownames(st.et) <- rownames(mar.st.et) <- rownames(bchol.st.et) <- rownames(var.et) <- colnames(var.et) <- c(Y.names, X.names)
rownames(E.obs.v) <- Y.names
rownames(var.obs.v) <- colnames(var.obs.v) <- Y.names
######################################
######################################
# output any warnings
if (!is.null(msg) && object[["control"]][["trace"]] >= 0 && !silent) cat("MARSSresiduals.tT reported warnings. See msg element or attribute of returned residuals object.\n")
######################################
return(list(
model.residuals = et[1:n, , drop = FALSE],
state.residuals = et[(n + 1):(n + m), , drop = FALSE],
residuals = et,
var.residuals = var.et,
std.residuals = st.et,
mar.residuals = mar.st.et,
bchol.residuals = bchol.st.et,
E.obs.residuals = E.obs.v,
var.obs.residuals = var.obs.v,
msg = msg
))
}
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.