Nothing
########################################################################################################################
# MARSSinnovationsboot.r
# This bootstraps multivariate time series data using Stoffer and Wall algorithm
# It creates bootstrap data via sampling from the standardized innovations matrix
# In the MARSS code, this is referred to as the nonparametric bootstrap. Strictly speaking, it is not nonparametric.
########################################################################################################################
MARSSinnovationsboot <- function(MLEobj, nboot = 1000, minIndx = 3) {
if (any(is.na(MLEobj$marss$data))) {
stop("Stopped in MARSSinnovationsboot() because this algorithm resamples from the innovations and doesn't allow missing values.\n", call. = FALSE)
}
# The following need to be present in model: parameter estimates + Sigma, Kt, Innov; the latter 3 are part of $kf
if (any(is.null(MLEobj[["par"]]), is.null(MLEobj[["marss"]][["data"]]))) {
stop("Stopped in MARSSinnovationsboot(). The passed in marssMLE object is missing par or model$data elements.\n", call. = FALSE)
}
if (is.null(MLEobj[["kf"]]) || !is.array(MLEobj[["kf"]][["Innov"]]) || !is.array(MLEobj[["kf"]][["Sigma"]])) {
kf <- MARSSkfss(MLEobj)
} else {
kf <- MLEobj[["kf"]]
}
## Rename things for code readability
MODELobj <- MLEobj[["marss"]]
par.dims <- attr(MODELobj, "model.dims")
TT <- par.dims[["data"]][2]
m <- par.dims[["x"]][1]
n <- par.dims[["y"]][1]
f <- MODELobj[["fixed"]]
d <- MODELobj[["free"]]
pari <- parmat(MLEobj, t = 1)
#### make a list of time-varying parameters
time.varying <- list()
for (elem in attr(MODELobj, "par.names")) {
if ((dim(d[[elem]])[3] == 1) & (dim(f[[elem]])[3] == 1)) {
time.varying[[elem]] <- FALSE
} else {
time.varying[[elem]] <- TRUE
} # not time-varying
}
##### Set holders for output
newData <- matrix(NA, n, TT) # store as written for state-space eqn with time across columns
newStates <- matrix(NA, m, TT + 1) # store as written for state-space eqn with time across columns
boot.data <- array(NA, dim = c(par.dims[["data"]], nboot))
boot.states <- array(NA, dim = c(par.dims[["x"]], nboot))
# Calculate the sqrt of sigma matrices, so they don't have to be computed 5000+ times
sigma.Sqrt <- array(0, dim = c(n, n, TT))
BKS <- array(0, dim = c(m, n, TT))
for (i in 1:TT) {
std.innovations <- stdInnov(kf$Sigma, kf$Innov) # standardized innovations; time across rows
eig <- eigen(kf$Sigma[, , i]) # calculate inv sqrt(sigma[1])
sigma.Sqrt[, , i] <- eig$vectors %*% tcrossprod(makediag(sqrt(eig$values)), eig$vectors)
if (time.varying$B & i > 1) pari$B <- parmat(MLEobj, "B", t = i)$B
BKS[, , i] <- pari$B %*% kf$Kt[, , i] %*% sigma.Sqrt[, , i] # this is m x n
}
for (i in 1:nboot) {
# Then the bootstrap algorithm proceeds
# Stoffer & Wall suggest not sampling from innovations 1-3 (not much data)
minIndx <- ifelse(TT > 5, minIndx, 1)
samp <- sample(seq(minIndx + 1, TT), size = (TT - minIndx), replace = TRUE)
e <- as.matrix(std.innovations)
e[, 1:minIndx] <- as.matrix(std.innovations[, 1:minIndx])
e[, (minIndx + 1):TT] <- as.matrix(std.innovations[, samp])
# newStates is a[] in the writeup by EH
newStates[, 1] <- pari$x0
# newStates[,2] = B %*% x0 + U + B %*% Kt[,,1] %*% sigma.Sqrt[,,1] %*% e[,1]
for (t in 1:TT) {
if (time.varying$B & i > 1) pari$B <- parmat(MLEobj, "B", t = i)$B
if (time.varying$U & i > 1) pari$U <- parmat(MLEobj, "U", t = i)$U
if (time.varying$Z & i > 1) pari$Z <- parmat(MLEobj, "Z", t = i)$Z
if (time.varying$A & i > 1) pari$A <- parmat(MLEobj, "A", t = i)$A
newStates[, t + 1] <- pari$B %*% newStates[, t, drop = FALSE] + pari$U + BKS[, , t] %*% e[, t, drop = FALSE]
newData[, t] <- pari$Z %*% newStates[, t, drop = FALSE] + pari$A + sigma.Sqrt[, , t] %*% e[, t, drop = FALSE]
}
newStates <- newStates[, 2:(TT + 1)]
boot.data[, , i] <- newData
boot.states[, , i] <- newStates
# reset newStates to its original dim
newStates <- matrix(NA, m, TT + 1)
}
return(list(boot.states = boot.states, boot.data = boot.data, marss = MODELobj, nboot = nboot))
}
######################################################################################################################
# stdInnov
######################################################################################################################
stdInnov <- function(SIGMA, INNOV) {
# This function added by EW Nov 3, 2008
# SIGMA is covariance matrix, E are original innovations
TT <- dim(INNOV)[2]
n <- dim(INNOV)[1]
SI <- matrix(0, n, TT)
for (i in 1:TT) {
a <- SIGMA[, , i]
a.eig <- eigen(a)
a.sqrt <- a.eig$vectors %*% makediag((sqrt(a.eig$values))) %*% solve(a.eig$vectors)
SI[, i] <- chol2inv(chol(a.sqrt)) %*% INNOV[, i] # eqn 1, p359 S&S
}
SI[which(INNOV == 0)] <- 0 # make sure things that should be 0 are 0
return(SI)
}
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.