Multiple minor corrections

Dimensionality check corrections to model components, degenerate model
checks for loglik.SSModel, corrected subset assignment function and
coef.KFS.

ChangeLog

@@ -2,17 +2,23 @@ Changes from Version 1.0.3 to 1.0.4:
  * Corrected output of LogLik method for non-Gaussian models: It now returns Inf only when the 
    approximation algoritm failed  completely (resulting NA), and issues only warning about 
    non-convergence in other cases.
+ * Added checks of degenerate model to LogLik method. If all elements in R, Q and H/u are zero, or 
+   they contain any non-finite values, -Inf is returned.
  * Corrected a bug in approximation algorithm which caused the approximation to fail for seemingly 
    random models.
  * Fixed a bug in SSMcycle which caused error with common components.
  * Fixed a bug in SSMseasonal which caused error in SSModel when using common components.
+ * SSMseasonal with trigonometric seasonal now works properly when period is odd.  
+ * Fixed a bug in coef.KFS which caused function to return smoothed states even with argument 
+   filtered=TRUE if they were present in KFS object.
  * Added argument "maxiter" to predict.SSModel and changed its default value in all functions to 50.
  * Corrected a bug in function ldl which caused the decomposition of semidefinite matrices to fail 
    silently.
  * Changed variable mu to m for mean filtering for non-Gaussian models without simulation just like 
    in other cases.
  * Changed convergence criterion in approximation algorithm from linear predictor based to 
    deviance based.
+ * Properly exported assigment using subset method. See ?subset.SSModel for details.
 
 Changes from Version 1.0.2 to 1.0.3:
  * Changed default filtering option for Gaussian models back to "state" (was previously "none").

NAMESPACE

@@ -12,6 +12,7 @@ S3method(print,SSModel)
 S3method(residuals,KFS)
 S3method(rstandard,KFS)
 S3method(subset,SSModel)
+export("subset<-")
 export(KFS)
 export(SSMarima)
 export(SSMcustom)

R/SSMarima.R

@@ -64,8 +64,8 @@ SSMarima <- function(ar = NULL, ma = NULL, d = 0, Q, stationary = TRUE, index,
         nd <- which(diag(P1inf) == 0)
         temp <- try(solve(a = diag((m - d * p)^2) - matrix(kronecker(T[nd, nd], T[nd, nd]), (m - d * p)^2, (m - d * p)^2), b = c(R[nd, 
             , drop = FALSE] %*% Q %*% t(R[nd, , drop = FALSE]))), TRUE)        
-        if (class(temp) == "try-error" || class(try(chol(matrix(temp, m - d * p, m - d * p)), TRUE)) == "try-error") {            
-            stop("ARIMA part is non-stationary.")
+        if (class(temp) == "try-error") {            
+            stop("ARIMA part is numerically too close to non-stationarity.")
         } else P1[nd, nd] <- temp
     } else diag(P1inf) <- 1
     state_names <- paste0(rep(paste0("arima", 1:m1), p), rep(ynames, each = m1))

R/SSMcycle.R

@@ -73,16 +73,21 @@ SSMcycle <- function(period, type, Q, index, a1, P1, P1inf, n, ynames) {
             tvq <- max(dim(Q)[3] == n, 0, na.rm = TRUE)
 
             Qm <- array(0, c(m, m, tvq * (n - 1) + 1))
-            if (tvq == 1) {
+            if (tvq) {
                 for (i in 1:(tvq * (n - 1) + 1)) Qm[cbind(rep(1:(p * 2), p), rep(1:2, p^2) + rep(0:(p - 1) * 2, each = p * 2), 
                   i)] <- rep(Q[, , i], each = 2)
             } else Qm[cbind(rep(1:(p * 2), p), rep(1:2, p^2) + rep(0:(p - 1) * 2, each = p * 2), 1)] <- rep(Q, each = 2)
         } else {
-            if (length(Q) != 1 && (!identical(dim(Q)[1], dim(Q)[2]) || dim(Q)[1] != 1 || !(max(dim(Q)[3], 1, na.rm = TRUE) %in% 
+            if (length(Q) != 1 && (!identical(dim(Q)[1], dim(Q)[2]) || isTRUE(dim(Q)[1] != 1) || !(max(dim(Q)[3], 1, na.rm = TRUE) %in% 
                 c(1, n)))) 
                 stop("Misspecified Q, argument Q must be a scalar, (1 x 1) matrix, or (1 x 1 x 1)/(1 x 1 x n) array.")
-            Qm <- array(Q, c(m, m, tvq * (n - 1) + 1))
-
+            tvq <- max(dim(Q)[3] == n, 0, na.rm = TRUE)
+            Qm <- array(0, c(m, m, tvq * (n - 1) + 1))
+            if (tvq) {
+              for (i in 1:(tvq * (n - 1) + 1)) 
+                Qm[cbind(1:2, 1:2, i)] <- rep(Q[1, 1, i], 2)
+            } else Qm[cbind(1:2, 1:2, 1)] <- rep(Q, 2)
+
         }
 
 

R/SSMregression.R

@@ -178,7 +178,7 @@ SSMregression <- function(rformula, data, type, Q, index, R, a1, P1, P1inf,  n,
   } else {
     if (length(Q) == 1) 
       Q <- matrix(Q)
-    if (!identical(dim(Q)[1], dim(Q)[2]) || dim(Q)[1] > m || !(max(dim(Q)[3], 1, na.rm = TRUE) %in% c(1, n))) 
+    if (!identical(dim(Q)[1], dim(Q)[2]) || isTRUE(dim(Q)[1] > m) || !(max(dim(Q)[3], 1, na.rm = TRUE) %in% c(1, n))) 
       stop("Misspecified Q, argument Q must be (k x k) matrix, (k x k x 1), or (k x k x n) array where m is the number of disturbance terms.")
     k <- dim(Q)[1]
     tvq <- max(dim(Q)[3] == n, 0, na.rm = TRUE)
@@ -189,7 +189,7 @@ SSMregression <- function(rformula, data, type, Q, index, R, a1, P1, P1inf,  n,
       R <- diag(m)[, 1:k, drop = FALSE]
     } else R <- NULL
   } else {
-    if (!(dim(R)[1]==m) || dim(R)[2] != k || !(max(dim(R)[3], 1, na.rm = TRUE) %in% c(1, n))) 
+    if (isTRUE(!(dim(R)[1]==m)) || isTRUE(dim(R)[2] != k) || !(max(dim(R)[3], 1, na.rm = TRUE) %in% c(1, n))) 
       stop("Misspecified R, argument R must be (m x k) matrix, (m x k x 1), or (m x k x n) array where m is the number of states and k is the number of disturbance terms.")
     tvr <- max(dim(R)[3] == n, 0, na.rm = TRUE)
   }

R/SSMseasonal.R

@@ -40,8 +40,8 @@ SSMseasonal <- function(period, sea.type = c("dummy", "trigonometric"), type, Q,
 
     lambda <- 2 * pi * 1:floor((period - 1)/2)/period
     T_univariate[cbind(1:m1, 1:m1)] <- rep(c(cos(lambda), -1), each = 2, length = m1)
-    T_univariate[cbind(2 * 1:((m1 - 1)/2) - 1, 2 * 1:((m1 - 1)/2))] <- sin(lambda)
-    T_univariate[cbind(2 * 1:((m1 - 1)/2), 2 * 1:((m1 - 1)/2) - 1)] <- -sin(lambda)
+    T_univariate[which((col(T_univariate)-row(T_univariate))==1)[seq(from=1,by=2,length=length(lambda))]] <- sin(lambda)
+    T_univariate[which((col(T_univariate)-row(T_univariate))==-1)[seq(from=1,by=2,length=length(lambda))]] <- -sin(lambda)
   }
   m <- ((p - 1) * (type != 2) + 1) * (period - 1)
   # k <- ((p-1)*(type != 2)+1)*((sea.type=='dummy') + (sea.type=='trigonometric')*(period - 1))

R/SSMtrend.R

@@ -78,17 +78,20 @@ SSMtrend <- function(degree = 1, type, Q, index, a1, P1, P1inf, n, ynames) {
             tvq <- max(unlist(sapply(lapply(Q, dim), "[", 3)) > 1, 0, na.rm = TRUE)
 
             Qm <- array(0, c(m, m, tvq * (n - 1) + 1))
+            if(!is.list(Q))
+              stop("Q must be a list of length degree.")
             for (i in 1:degree) {
-                if (length(Q[[i]]) != 1 && (any(dim(Q[[i]])[1:2] != p) || !(dim(Q[[i]])[3] %in% c(1, n, NA)))) 
+                if ((p>1 && length(Q[[i]]) == 1) || 
+                      (p==1 && length(Q[[i]]) != 1) || (p>1 &&( dim(Q[[i]])[1:2] != p || !(max(1,dim(Q[[i]])[3],na.rm=TRUE) %in% c(1, n)))))
                   stop("Q must be a list of length degree, which contains (p x p) matrices, (p x p x 1), or (p x p x n) arrays, where p is the number of series. ")
                 Qm[seq(from = i, by = degree, length = p), seq(from = i, by = degree, length = p), ] <- Q[[i]]
             }
             k <- dim(Qm)[1]
             R <- diag(k)
 
         } else {
-            if (is.list(Q) || length(Q) != degree && (!identical(dim(Q)[1], dim(Q)[2]) || dim(Q)[1] != degree || !(max(dim(Q)[3], 
-                1, na.rm = TRUE) %in% c(1, n)))) 
+            if (is.list(Q) || (length(Q) != degree && is.null(dim(Q))) || 
+                  (any(dim(Q)[1:2] != degree) || !(max(1,dim(Q)[3],na.rm=TRUE) %in% c(1, n, NA))))
                 stop("Misspecified Q, argument Q must be a vector of length d, (d x d) matrix, or (d x d x 1)/(d x d x n) array where d is the degree of the trend.")
             if (length(Q) == degree) 
                 Q <- diag(drop(Q), degree)

R/coef.KFS.R

@@ -27,14 +27,4 @@ coef.KFS<-function(object,start=NULL,end=NULL,filtered=FALSE,...){
   if(start==end && !is.null(start))
     tsp(tmp) <- class(tmp) <- NULL
   drop(tmp)
-  if(!is.null(object$alphahat)){
-    tmp<-window(object$alphahat,start,end)
-  } else {
-    if(!is.null(object[["a",exact=TRUE]])){
-      tmp<-window(object$a,start,end)
-    } else stop("Input does not contain estimates for states, rerun KFS with state filtering or smoothing.")
-  }
-  if(start==end && !is.null(start))
-    tsp(tmp) <- class(tmp) <- NULL
-  drop(tmp)
 }

R/extract.SSModel.R

@@ -1,6 +1,8 @@
 #' Extract or Replace Parts of a State Space Model
 #'
-#' S3 methods for extracting or replacing parts of objects of class \code{SSModel}.
+#' S3 methods for extracting or replacing parts of objects of class \code{SSModel}. These methods 
+#' ensure that dimensions of system matrices are not altered. \code{[} and \code{subset} and 
+#' corresponding replacement methods are identical methods with different method names.
 #'
 #' @method [<- SSModel
 #' @S3method [<- SSModel
@@ -19,6 +21,19 @@
 #' @param times Numeric. Which time points are chosen.
 #' @param value A value to be assigned to x.
 #' @return A selected subset of the chosen element or a value.
+#' @examples
+#' set.seed(1)
+#' model<-SSModel(rnorm(10)~1)
+#' model["H"]
+#' model["H"]<-10
+#' # H is still an array:
+#' model["H"]
+#' logLik(model)
+#' model$H<-1
+#' # model["H"] throws an error as H is now scalar:
+#' model$H
+#' logLik(model,check.model=TRUE) #with check.model=FALSE (default) R crashes!
+
 `[<-.SSModel` <- function(x, element, states, etas, series, times, value) {
 
     element <- match.arg(arg = element, choices = c("y", "Z", "H", "T", "R", "Q", "a1", "P1", "P1inf", "u"))

R/logLik.SSModel.R

@@ -30,10 +30,10 @@ logLik.SSModel <- function(object, nsim = 0, antithetics = TRUE, theta, check.mo
                            seed, convtol=1e-15,...) {
   if (check.model) {
     if (!is.SSModel(object, na.check = TRUE)) {
-      warning("Not a valid SSModel object. Returning -Inf")
-      return(-Inf)
+      return(-.Machine$double.xmax)
     }
   }
+
   p <- attr(object, "p")
   m <- attr(object, "m")
   k <- attr(object, "k")
@@ -47,6 +47,8 @@ logLik.SSModel <- function(object, nsim = 0, antithetics = TRUE, theta, check.mo
   tv[4] <- dim(object$R)[3] > 1
   tv[5] <- dim(object$Q)[3] > 1
   if (all(object$distribution == "gaussian")) {
+      if(all(c(object$Q,object$H)==0) || all(c(object$R,object$H)==0)|| any(!is.finite(c(object$R,object$Q,object$H)==0)))
+        return(-.Machine$double.xmax^0.75)
     kfout <- NULL
     if (p == 1) {
       kfout <- .Fortran(fglogliku, NAOK = TRUE, object$y, ymiss, as.integer(tv), object$Z, object$H, object$T, object$R, 
@@ -57,7 +59,7 @@ logLik.SSModel <- function(object, nsim = 0, antithetics = TRUE, theta, check.mo
         object <- tryCatch(transformSSM(object, type = match.arg(arg = transform, choices = c("ldl","augment"))), error = function(e) e)
         if (!inherits(object, "SSModel")) {
           warning(object$message)
-          return(-Inf)
+          return(-.Machine$double.xmax^0.75)
         }
         tv[1] <- dim(object$Z)[3] > 1
         tv[2] <- dim(object$H)[3] > 1
@@ -69,7 +71,8 @@ logLik.SSModel <- function(object, nsim = 0, antithetics = TRUE, theta, check.mo
     }
     logLik <- kfout$lik
   } else {
-
+    if(all(c(object$Q,object$u)==0) || all(c(object$R,object$u)==0) || any(!is.finite(c(object$R,object$Q,object$u)==0)))
+      return(-.Machine$double.xmax^0.75)
     if (missing(theta)) {
       theta <- sapply(1:p, function(i)
         switch(object$distribution[i], 
@@ -138,7 +141,7 @@ logLik.SSModel <- function(object, nsim = 0, antithetics = TRUE, theta, check.mo
                     as.integer(sim), nsim2, as.integer(nd), as.integer(length(nd)),diff=double(1))
     if (!is.finite(out$diff)){
       warning("Non-finite difference in approximation algoritm. Returning -Inf.")
-      return(-Inf)
+      return(-.Machine$double.xmax^0.75)
     }
     if(out$maxiter==maxiter){
       warning(paste("Maximum number of iterations reached, 

R/subset.SSModel.R

@@ -63,7 +63,8 @@
     x
 }
 #' @rdname Extract.SSModel
-`subset<-` <- function(x, value) UseMethod("subset<-")
+#' @export
+`subset<-` <- function(x, value,...) UseMethod("subset<-")
 
 #' @method subset SSModel
 #' @S3method subset SSModel

man/Extract.SSModel.Rd

@@ -12,7 +12,7 @@
 
 \method{subset}{SSModel}(x, element, states, etas, series, times) <- value
 
-subset(x) <- value
+subset(x, ...) <- value
 
 \method{subset}{SSModel}(x, element, states, etas, series, times, ...)
 }
@@ -54,6 +54,22 @@ A selected subset of the chosen element or a value.
 }
 \description{
 S3 methods for extracting or replacing parts of objects of
-class \code{SSModel}.
+class \code{SSModel}. These methods ensure that dimensions
+of system matrices are not altered. \code{[} and
+\code{subset} and corresponding replacement methods are
+identical methods with different method names.
+}
+\examples{
+set.seed(1)
+model<-SSModel(rnorm(10)~1)
+model["H"]
+model["H"]<-10
+# H is still an array:
+model["H"]
+logLik(model)
+model$H<-1
+# model["H"] throws an error as H is now scalar:
+model$H
+logLik(model,check.model=TRUE) #with check.model=FALSE (default) R crashes!
 }
 

man/KFS.Rd

@@ -3,13 +3,14 @@
 \title{Kalman Filter and Smoother with Exact Diffuse Initialization for Exponential Family State Space Models}
 \usage{
 KFS(model, filtering, smoothing, simplify = TRUE, transform = c("ldl",
-  "augment"), nsim = 0, theta, maxiter = 100)
+  "augment"), nsim = 0, theta, maxiter = 50, convtol = 1e-15)
 }
 \arguments{
   \item{model}{Object of class \code{SSModel}.}
 
   \item{filtering}{Types of filtering. Possible choices are
-  'state', 'signal', 'mean', and 'none'. Default is 'none'.
+  'state', 'signal', 'mean', and 'none'. Default is 'state'
+  for Gaussian and 'none' for non-Gaussian models.
   Multiple values are allowed. Note that for Gaussian
   models, signal is mean. Note that filtering for
   non-Gaussian models with importance sampling can be very
@@ -43,8 +44,13 @@ KFS(model, filtering, smoothing, simplify = TRUE, transform = c("ldl",
   Only used for non-Gaussian model.}
 
   \item{maxiter}{The maximum number of iterations used in
-  linearisation. Default is 25. Only used for non-Gaussian
+  approximation Default is 50. Only used for non-Gaussian
   model.}
+
+  \item{convtol}{Tolerance parameter for convergence checks
+  for Gaussian approximation.  Iterations are continued
+  until
+  \eqn{tol>abs(dev_{old}-dev_{new})/(abs(dev_{new})+0.1))}.}
 }
 \value{
 What \code{KFS} returns depends on the arguments
@@ -148,12 +154,12 @@ at the time, the elements of prediction error
 \code{Finf}, \code{K} and \code{Kinf} contain only the
 diagonal elemens of the corresponding covariance matrices.
 
-In rare cases of a very long diffuse initialization phase
-with highly correlated states, cumulative rounding errors
-in computing \code{Finf} and \code{Pinf} can sometimes
-cause the diffuse phase end too early. Changing the
-tolerance parameter \code{tol} of the model (see
-\code{\link{SSModel}}) to smaller (or larger) should help.
+In rare cases of a diffuse initialization phase with highly
+correlated states, cumulative rounding errors in computing
+\code{Finf} and \code{Pinf} can sometimes cause the diffuse
+phase end too early. Changing the tolerance parameter
+\code{tol} of the model (see \code{\link{SSModel}}) to
+smaller (or larger) should help.
 
 In case of non-Gaussian models with \code{nsim=0}, the
 smoothed estimates relate the conditional mode of

man/SSModel.Rd

@@ -52,8 +52,10 @@ SSMtrend(degree = 1, type, Q, index, a1, P1, P1inf, n, ynames)
   observations. Default is \code{rep('gaussian',p)}.}
 
   \item{tol}{a tolerance parameter for a diffuse phase.
-  Smallest value of Finf not counted for zero. Defaults to
-  \code{.Machine$double.eps^0.5}.}
+  Smallest value of Finf not counted for zero.  Defaults to
+  \code{.Machine$double.eps^0.5}.  If smoothing gives
+  negative variances for smoothed states, try adjusting
+  this.}
 
   \item{index}{a vector indicating for which series the
   corresponding components are constructed.}

man/approxSSM.Rd

@@ -1,8 +1,8 @@
 \name{approxSSM}
 \alias{approxSSM}
-\title{Linear gaussian Approximation for Non-gaussian State Space Model}
+\title{Linear Gaussian Approximation for Exponential Family State Space Model}
 \usage{
-approxSSM(model, theta, maxiter = 100, tol = 1e-08)
+approxSSM(model, theta, maxiter = 50, tol = 1e-15)
 }
 \arguments{
   \item{model}{A non-Gaussian state space model object of
@@ -11,17 +11,17 @@ approxSSM(model, theta, maxiter = 100, tol = 1e-08)
   \item{theta}{Initial values for conditional mode theta.}
 
   \item{maxiter}{The maximum number of iterations used in
-  linearisation. Default is 100.}
+  approximation Default is 50.}
 
   \item{tol}{Tolerance parameter for convergence checks.
   Iterations are continued until
-  \eqn{tol>sum(abs(theta_{new}-theta_{old})/(abs(theta_{old})+0.1))/(n*p)}.}
+  \eqn{tol>abs(dev_{old}-dev_{new})/(abs(dev_{new})+0.1))}.}
 }
 \value{
 An object which contains the approximating Gaussian state
 space model with following additional components:
-\item{thetahat}{mode of \eqn{p(\theta|y)}.}
-\item{iterations}{number of iterations used.}
+\item{thetahat}{Mode of \eqn{p(\theta|y)}. }
+\item{iterations}{Number of iterations used. }
 }
 \description{
 Function \code{approxSMM} computes the linear Gaussian
@@ -43,14 +43,15 @@ conditional mode of \eqn{\theta=Z\alpha} given the
 observations \eqn{y} as the original non-gaussian model.
 Models also have a same curvature at the mode.
 
-The linearization of the exponential family state space
+The approximation of the exponential family state space
 model is based on iterative weighted least squares method,
 see McCullagh and Nelder (1983) p.31 and Durbin Koopman
 (2012) p. 243.
 }
 \seealso{
-Importance sampling of non-gaussian state space models
+Importance sampling of non-Gaussian state space models
 \code{\link{importanceSSM}}, construct a \code{SSModel}
-object \code{\link{SSModel}}.
+object \code{\link{SSModel}}, and examples in
+\code{\link{KFAS}}.
 }
 

man/importanceSSM.Rd

@@ -4,7 +4,7 @@
 \usage{
 importanceSSM(model, type = c("states", "signals"), filtered = FALSE,
   nsim = 1000, save.model = FALSE, theta, antithetics = FALSE,
-  maxiter = 100)
+  maxiter = 50)
 }
 \arguments{
   \item{model}{Exponential family state space model of
@@ -30,7 +30,7 @@ importanceSSM(model, type = c("states", "signals"), filtered = FALSE,
   another for scale. Default is FALSE.}
 
   \item{maxiter}{Maximum number of iterations used in
-  linearisation. Default is 25.}
+  linearisation. Default is 50.}
 }
 \value{
 A list containing elements \code{samples}, \code{weights}