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