1#' Transform Multivariate State Space Model for Sequential Processing 2#' 3#' \code{transformSSM} transforms the general multivariate Gaussian state space model 4#' to form suitable for sequential processing. 5#' 6#' @details As all the functions in KFAS use univariate approach i.e. sequential processing, 7#' the covariance matrix \eqn{H_t}{H[t]} of the observation equation needs to be 8#' either diagonal or zero matrix. Function \code{transformSSM} performs either 9#' the LDL decomposition of \eqn{H_t}{H[t]}, or augments the state vector with 10#' the disturbances of the observation equation. 11#' 12#' In case of a LDL decomposition, the new \eqn{H_t}{H[t]} contains the diagonal part of the 13#' decomposition, whereas observations \eqn{y_t}{y[t]} and system matrices \eqn{Z_t}{Z[t]} are 14#' multiplied with the inverse of \eqn{L_t}{L[t]}. Note that although the state estimates and 15#' their error covariances obtained by Kalman filtering and smoothing are identical with those 16#' obtained from ordinary multivariate filtering, the one-step-ahead errors 17#' \eqn{v_t}{v[t]} and their variances \eqn{F_t}{F[t]} do differ. The typical 18#' multivariate versions can be obtained from output of \code{\link{KFS}} 19#' using \code{\link{mvInnovations}} function. 20#' 21#' In case of augmentation of the state vector, some care is needed interpreting the 22#' subsequent filtering/smoothing results: For example the \code{muhat} from the output of \code{KFS} 23#' now contains also the smoothed observational level noise as that is part of the state vector. 24#' 25#' 26#' @export 27#' @param object State space model object from function \code{\link{SSModel}}. 28#' @param type Option \code{"ldl"} performs LDL decomposition for covariance matrix \eqn{H_t}{H[t]}, 29#' and multiplies the observation equation with the \eqn{L_t^{-1}}{L[t]^-1}, so \eqn{\epsilon_t^* 30#' \sim N(0,D_t)}{\epsilon[t]* ~ N(0,D[t])}. Option \code{"augment"} adds 31#' \eqn{\epsilon_t}{\epsilon[t]} to the state vector, so \eqn{Q_t}{Q[t]} becomes block diagonal 32#' with blocks \eqn{Q_t}{Q[t]} and \eqn{H_t}{H[t]}. 33#' @param tol Tolerance parameter for LDL decomposition (see \code{\link{ldl}}). Default is 34#' \code{max(100, max(abs(apply(object$H, 3, diag)))) * .Machine$double.eps}. 35#' @return \item{model}{Transformed model.} 36transformSSM <- function(object, type = c("ldl", "augment"), tol) { 37 38 if (any(object$distribution != "gaussian")) 39 stop("Nothing to transform as matrix H is not defined for non-gaussian model.") 40 is.SSModel(object, return.logical = FALSE) 41 type <- match.arg(type, choices = c("ldl", "augment")) 42 p <- attr(object, "p") 43 n <- attr(object, "n") 44 m <- attr(object, "m") 45 r <- attr(object, "k") 46 tv <- attr(object,"tv") 47 tvh <- tv[2] 48 if (type == "ldl") { 49 if (p > 1) { #do nothing for univariate series 50 yt <- t(object$y) 51 ymiss <- is.na(yt) 52 tv[1] <- max(tv[1], tv[2]) 53 if (sum(ymiss) > 0) { 54 # find unique combinations of missing observations 55 # even though H (and Z) becomes time varying in case of partial missingness, 56 # no need to compute Cholesky for all t 57 positions <- unique(ymiss, MARGIN = 2) 58 nh <- dim(positions)[2] 59 tv[1:2] <- 1 60 Z <- array(object$Z, dim = c(p, m, n)) 61 } else { 62 Z <- array(object$Z, dim = c(p, m, (n - 1) * tv[1] + 1)) 63 positions <- rep(FALSE, p) 64 nh <- 1 65 } 66 positions <- as.matrix(positions) 67 H <- array(object$H, c(p, p, n)) 68 if (tvh) { 69 # if H was already time varying, compute cholesky decompositions for all t 70 nh <- n 71 hchol <- 1:n 72 uniqs <- 1:n 73 } else { 74 hchol <- rep(0, n) 75 uniqs <- numeric(nh) 76 for (i in 1:nh) { 77 nhn <- which(colSums(ymiss == positions[, i]) == p) 78 hchol[nhn] <- i 79 uniqs[i] <- nhn[1] 80 } 81 } 82 ichols <- H[, , uniqs, drop = FALSE] 83 ydims <- as.integer(colSums(!ymiss)) 84 yobs <- array(1:p, c(p, n)) 85 if (sum(ymiss) > 0) #positions of missing observations 86 for (i in 1:n) { 87 if (ydims[i] != p && ydims[i] != 0) 88 yobs[1:ydims[i], i] <- yobs[!ymiss[, i], i] 89 if (ydims[i] < p) 90 yobs[(ydims[i] + 1):p, i] <- NA 91 } 92 unidim <- ydims[uniqs] 93 hobs <- yobs[, uniqs, drop = FALSE] 94 storage.mode(yobs) <- storage.mode(hobs) <- storage.mode(hchol) <- "integer" 95 if(missing(tol)) tol <- max(100, max(abs(apply(object$H, 3, diag)))) * .Machine$double.eps 96 # compute the transformations for y, Z and H 97 out <- .Fortran(fldlssm, NAOK = TRUE, yt = yt, ydims = ydims, yobs = yobs, 98 tv = as.integer(tv), Zt = Z, p = p, m = m, 99 n = n, ichols = ichols, nh = as.integer(nh), hchol = hchol, 100 unidim = as.integer(unidim), info = as.integer(0), hobs = hobs, 101 tol = tol) 102 if(out$info!=0){ 103 stop(switch(as.character(out$info), 104 "1" = "LDL decomposition of H failed.", 105 "2" = "Computing the inverse of L failed." 106 )) 107 } 108 H <- array(0, c(p, p, ((n - 1) * tv[2] + 1))) 109 for (i in 1:((n - 1) * tv[2] + 1)) #construct diagonal H 110 diag(H[, , i]) <- diag(out$ichols[, , out$hchol[i]]) 111 attry <- attributes(object$y) 112 object$y <- t(out$yt) 113 attributes(object$y) <- attry 114 object$Z <- out$Z 115 object$H <- H 116 attr(object, "tv") <- as.integer(tv) 117 } 118 } else { #augmentation, add new states for epsilon disturbances and set H=0 119 T <- array(object$T, dim = c(m, m, (n - 1) * tv[3] + 1)) 120 R <- array(object$R, dim = c(m, r, (n - 1) * tv[4] + 1)) 121 Q <- array(object$Q, dim = c(r, r, (n - 1) * tv[5] + 1)) 122 H <- array(object$H, c(p, p, (n - 1) * tv[2] + 1)) 123 Z <- array(object$Z, dim = c(p, m, (n - 1) * tv[1] + 1)) 124 r2 <- r + p 125 m2 <- m + p 126 tv[5] <- max(tv[c(2, 5)]) 127 Qt2 <- array(0, c(r2, r2, 1 + (n - 1) * tv[5])) 128 Qt2[1:r, 1:r, ] <- Q 129 if (tv[2]) { 130 Qt2[(r + 1):r2, (r + 1):r2, -n] <- H[, , -1] 131 } else Qt2[(r + 1):r2, (r + 1):r2, ] <- H 132 Zt2 <- array(0, c(p, m2, (n - 1) * tv[1] + 1)) 133 Zt2[1:p, 1:m, ] <- Z 134 Zt2[1:p, (m + 1):m2, ] <- diag(p) 135 Tt2 <- array(0, c(m2, m2, (n - 1) * tv[3] + 1)) 136 Tt2[1:m, 1:m, ] <- T 137 Rt2 <- array(0, c(m2, r + p, (n - 1) * tv[4] + 1)) 138 Rt2[1:m, 1:r, ] <- R 139 Rt2[(m + 1):m2, (r + 1):r2, ] <- diag(p) 140 P12 <- P1inf2 <- matrix(0, m2, m2) 141 P12[1:m, 1:m] <- object$P1 142 P1inf2[1:m, 1:m] <- object$P1inf 143 P12[(m + 1):m2, (m + 1):m2] <- H[, , 1] 144 a12 <- matrix(0, m2, 1) 145 a12[1:m, ] <- object$a1 146 object$Z <- Zt2 147 object$H <- array(0, c(p, p, 1)) 148 object$T <- Tt2 149 object$R <- Rt2 150 object$Q <- Qt2 151 if (p == 1) { 152 rownames(a12) <- c(rownames(object$a1), "eps") 153 } else { 154 rownames(a12) <- c(rownames(object$a1),paste0(rep("eps.", p), 1:p)) 155 } 156 object$a1 <- a12 157 object$P1 <- P12 158 object$P1inf <- P1inf2 159 attr(object, "m") <- m2 160 attr(object, "k") <- r2 161 attr(object, "tv") <- as.integer(tv) 162 163 } 164 invisible(object) 165} 166