From 5bc974f948607221baf0cffc0f55862aa5a2bfda Mon Sep 17 00:00:00 2001
From: "jouni.helske@jyu.fi" <jouni.helske@jyu.fi>
Date: Wed, 16 Apr 2014 07:15:21 +0300
Subject: [PATCH] Multiple minor corrections

Dimensionality check corrections to model components, degenerate model
checks for loglik.SSModel, corrected subset assignment function and
coef.KFS.
---
 ChangeLog                |  6 +++++
 NAMESPACE                |  1 +
 R/SSMarima.R             |  4 +--
 R/SSMcycle.R             | 13 +++++++---
 R/SSMregression.R        |  4 +--
 R/SSMseasonal.R          |  4 +--
 R/SSMtrend.R             |  9 ++++---
 R/coef.KFS.R             | 10 --------
 R/extract.SSModel.R      | 17 ++++++++++++-
 R/logLik.SSModel.R       | 13 ++++++----
 R/subset.SSModel.R       |  3 ++-
 man/Extract.SSModel.Rd   | 20 +++++++++++++--
 man/KFS.Rd               | 24 +++++++++++-------
 man/SSModel.Rd           |  6 +++--
 man/approxSSM.Rd         | 19 +++++++-------
 man/importanceSSM.Rd     |  4 +--
 man/logLik.SSModel.Rd    | 11 +++++---
 man/predict.SSModel.Rd   |  3 ++-
 src/glogliku.f95         | 22 ++++++++--------
 tests/arimatests.R       | 33 ++++++++++++++++++++++++
 tests/structTest.R       | 55 +++++++++++++++++++++++++++++++++-------
 tests/testthat/glmfits.R |  2 +-
 22 files changed, 204 insertions(+), 79 deletions(-)
 create mode 100644 tests/arimatests.R

diff --git a/ChangeLog b/ChangeLog
index 565b3ef..96e87ac 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -2,10 +2,15 @@ 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.
@@ -13,6 +18,7 @@ Changes from Version 1.0.3 to 1.0.4:
    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").
diff --git a/NAMESPACE b/NAMESPACE
index fca5f2e..bb3978a 100644
--- a/NAMESPACE
+++ b/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)
diff --git a/R/SSMarima.R b/R/SSMarima.R
index 6738e0c..57c98f9 100644
--- a/R/SSMarima.R
+++ b/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))
diff --git a/R/SSMcycle.R b/R/SSMcycle.R
index 97fd258..31bd559 100644
--- a/R/SSMcycle.R
+++ b/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)
+                        
         }
         
         
diff --git a/R/SSMregression.R b/R/SSMregression.R
index 59dd1e9..cea1562 100644
--- a/R/SSMregression.R
+++ b/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)
   }
diff --git a/R/SSMseasonal.R b/R/SSMseasonal.R
index fb97a08..7954eb4 100644
--- a/R/SSMseasonal.R
+++ b/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))
diff --git a/R/SSMtrend.R b/R/SSMtrend.R
index ebcf844..3e7fcbd 100644
--- a/R/SSMtrend.R
+++ b/R/SSMtrend.R
@@ -78,8 +78,11 @@ 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]]
             }
@@ -87,8 +90,8 @@ SSMtrend <- function(degree = 1, type, Q, index, a1, P1, P1inf, n, ynames) {
             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)
diff --git a/R/coef.KFS.R b/R/coef.KFS.R
index 02829cf..29faed5 100644
--- a/R/coef.KFS.R
+++ b/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)
 }
\ No newline at end of file
diff --git a/R/extract.SSModel.R b/R/extract.SSModel.R
index a4fda56..c944d65 100644
--- a/R/extract.SSModel.R
+++ b/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"))
diff --git a/R/logLik.SSModel.R b/R/logLik.SSModel.R
index 0d7ba18..9f1e1b2 100644
--- a/R/logLik.SSModel.R
+++ b/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, 
diff --git a/R/subset.SSModel.R b/R/subset.SSModel.R
index 8c6097f..705d8e4 100644
--- a/R/subset.SSModel.R
+++ b/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
diff --git a/man/Extract.SSModel.Rd b/man/Extract.SSModel.Rd
index dfa3f15..3fe7a57 100644
--- a/man/Extract.SSModel.Rd
+++ b/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!
 }
 
diff --git a/man/KFS.Rd b/man/KFS.Rd
index bf38e09..480848d 100644
--- a/man/KFS.Rd
+++ b/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
diff --git a/man/SSModel.Rd b/man/SSModel.Rd
index 50d330f..182f9b5 100644
--- a/man/SSModel.Rd
+++ b/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.}
diff --git a/man/approxSSM.Rd b/man/approxSSM.Rd
index 13b50cd..f7daa3b 100644
--- a/man/approxSSM.Rd
+++ b/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}}.
 }
 
diff --git a/man/importanceSSM.Rd b/man/importanceSSM.Rd
index ef2e73a..5128aab 100644
--- a/man/importanceSSM.Rd
+++ b/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}
diff --git a/man/logLik.SSModel.Rd b/man/logLik.SSModel.Rd
index 79464f6..e39ce83 100644
--- a/man/logLik.SSModel.Rd
+++ b/man/logLik.SSModel.Rd
@@ -4,8 +4,8 @@
 \title{Log-likelihood of the State Space Model.}
 \usage{
 \method{logLik}{SSModel}(object, nsim = 0, antithetics = TRUE, theta,
-  check.model = FALSE, transform = c("ldl", "augment"), maxiter = 25,
-  seed, ...)
+  check.model = FALSE, transform = c("ldl", "augment"), maxiter = 50,
+  seed, convtol = 1e-15, ...)
 }
 \arguments{
   \item{object}{State space model of class \code{SSModel}.}
@@ -34,12 +34,17 @@
   details.}
 
   \item{maxiter}{The maximum number of iterations used in
-  linearisation. Default is 25. Only used for non-Gaussian
+  linearisation. Default is 50. Only used for non-Gaussian
   model.}
 
   \item{seed}{The value is used as a seed via set.seed
   function. 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))}.}
+
   \item{...}{Ignored.}
 }
 \value{
diff --git a/man/predict.SSModel.Rd b/man/predict.SSModel.Rd
index ed6c62c..3464903 100644
--- a/man/predict.SSModel.Rd
+++ b/man/predict.SSModel.Rd
@@ -5,7 +5,8 @@
 \usage{
 \method{predict}{SSModel}(object, newdata, interval = c("none", "confidence",
   "prediction"), level = 0.95, type = c("response", "link"),
-  states = NULL, se.fit = FALSE, nsim = 0, prob = TRUE, ...)
+  states = NULL, se.fit = FALSE, nsim = 0, prob = TRUE, maxiter = 50,
+  ...)
 }
 \arguments{
   \item{object}{Object of class \code{SSModel}.}
diff --git a/src/glogliku.f95 b/src/glogliku.f95
index 7a68ab7..92951c3 100644
--- a/src/glogliku.f95
+++ b/src/glogliku.f95
@@ -111,20 +111,20 @@ subroutine glogliku(yt, ymiss, timevar, zt, ht,tt, rt, qt, a1, p1, p1inf,&
         end do diffuse
 
         !non-diffuse filtering begins
-         if(rankp .EQ. 0) then
+        if(rankp .EQ. 0) then
 
-        call dgemv('n',m,m,1.0d0,tt(:,:,(d-1)*timevar(3)+1),m,arec,1,0.0d0,at(:),1)  !at(:,t+1) = matmul(tt,a_rec)
+            call dgemv('n',m,m,1.0d0,tt(:,:,(d-1)*timevar(3)+1),m,arec,1,0.0d0,at(:),1)  !at(:,t+1) = matmul(tt,a_rec)
 
-        call dsymm('r','u',m,m,1.0d0,prec,m,tt(:,:,(d-1)*timevar(3)+1),m,0.0d0,mm,m)
-        call dgemm('n','t',m,m,m,1.0d0,mm,m,tt(:,:,(d-1)*timevar(3)+1),m,0.0d0,pt,m)
+            call dsymm('r','u',m,m,1.0d0,prec,m,tt(:,:,(d-1)*timevar(3)+1),m,0.0d0,mm,m)
+            call dgemm('n','t',m,m,m,1.0d0,mm,m,tt(:,:,(d-1)*timevar(3)+1),m,0.0d0,pt,m)
 
-        call dsymm('r','u',m,r,1.0d0,qt(:,:,(d-1)*timevar(5)+1),r,rt(:,:,(d-1)*timevar(4)+1),m,0.0d0,mr,m)
-        call dgemm('n','t',m,m,r,1.0d0,mr,m,rt(:,:,(d-1)*timevar(4)+1),m,1.0d0,pt,m)
+            call dsymm('r','u',m,r,1.0d0,qt(:,:,(d-1)*timevar(5)+1),r,rt(:,:,(d-1)*timevar(4)+1),m,0.0d0,mr,m)
+            call dgemm('n','t',m,m,r,1.0d0,mr,m,rt(:,:,(d-1)*timevar(4)+1),m,1.0d0,pt,m)
 
 
-        call dcopy(m,at(:),1,arec,1)
-        prec = pt
-end if
+            call dcopy(m,at(:),1,arec,1)
+            prec = pt
+        end if
     end if
 
    
@@ -142,11 +142,11 @@ subroutine glogliku(yt, ymiss, timevar, zt, ht,tt, rt, qt, a1, p1, p1inf,&
             vt = yt(t,1) - ddot(m,zt(1,:,(t-1)*timevar(1)+1),1,arec,1)
             call dsymv('u',m,1.0d0,prec,m,zt(1,:,(t-1)*timevar(1)+1),1,0.0d0,kt(:,1),1)
             ft = ddot(m,zt(1,:,(t-1)*timevar(1)+1),1,kt(:,1),1)+ht(1,1,(t-1)*timevar(2)+1)
-           if (ft .GT. meps) then
+            if (ft .GT. meps) then
                 call daxpy(m,vt/ft,kt(:,1),1,arec,1) !a_rec = a_rec + kt(:,i,t)*vt(:,t)
                 call dsyr('u',m,-1.0d0/ft,kt(:,1),1,prec,m) !p_rec = p_rec - kt*kt'*ft(i,i,t)
                 lik = lik - 0.5d0*(log(ft) + vt**2/ft)-c
-           end if
+            end if
         end if
         call dgemv('n',m,m,1.0d0,tt(:,:,(t-1)*timevar(3)+1),m,arec,1,0.0d0,at(:),1)
         call dsymm('r','u',m,m,1.0d0,prec,m,tt(:,:,(t-1)*timevar(3)+1),m,0.0d0,mm,m)
diff --git a/tests/arimatests.R b/tests/arimatests.R
new file mode 100644
index 0000000..2a8799a
--- /dev/null
+++ b/tests/arimatests.R
@@ -0,0 +1,33 @@
+s <- 12
+phis <- 0.99
+phi1 <- 0.0001
+phi <- c(phi1,rep(0,s-2),phis,-phi1*phis)
+theta <- 0.7
+out <- makeARIMA(phi,theta,NULL,SSinit="Ross")
+min(eigen(out$Pn)$value)
+
+set.seed(1)
+ x <- arima.sim(100,model=list(ar=phi,ma=theta))
+KalmanLike(x,mod=out)
+out<-arima(x,c(13,0,1)
+model<-SSModel(x~-1+SSMarima(ar=phi,ma=theta),H=0)
+logLik(model)
+arima(x,order=c(1,0,1),seasonal=list(period=12,order=c(1,0,0)),include.mean=FALSE,init=c(phi1,theta,phis),method='ML',
+      SSinit="Ros")
+
+Q0ter <- function(phi,theta){
+  p <- length(phi)
+  q <- length(theta)
+  r <- max(p,q+1)
+  T <- matrix(0,r,r)
+  if (p) T[1:p,1] <- phi
+  if (r) T[1:(r-1),2:r] <- diag(r-1)
+  V <- matrix(0,r,r)
+  ttheta <- c(1,theta)
+  V[1:(q+1),1:(q+1)] <- ttheta%x%t(ttheta)
+  V <- matrix(V,ncol=1)
+  S <- diag(r*r)-T%x%T
+  Q0 <- solve(S,V)
+  Q0 <- matrix(Q0,ncol=r)
+}
+microbenchmark(makeARIMA(phi,theta,NULL,SSinit="Ross")$Pn,SSMarima(ar=phi,ma=theta)$Q,Q0ter(phi,theta))
diff --git a/tests/structTest.R b/tests/structTest.R
index ede0e2f..0326905 100644
--- a/tests/structTest.R
+++ b/tests/structTest.R
@@ -1,11 +1,48 @@
-t1<-ts(rnorm(100))
-t2<-ts(rnorm(100))
+context("Structural time series tests")
 
-t<-cbind(t1,t2)
+test_that("SSModel works properly",{
+  set.seed(123)
+  t12<-ts(cbind(t1=rnorm(100),t2=rnorm(100)))
+  t12[sample(size=50,1:200)]<-NA
+  expect_that( 
+    model<-SSModel(t12~SSMcycle(period=10, type='common',Q=2)
+                   +SSMcycle(period=10, type='distinct',P1=diag(c(1,1,2,2)),Q=diag(1:2))
+                   +SSMtrend(2,type="common",Q=diag(c(1,0.5)))
+                   +SSMtrend(2,type="distinct",Q=list(diag(0.1,2),diag(0.01,2)),P1=diag(c(0.1,0.01,0.1,0.01)))
+                   +SSMseasonal(period=4,type="common")
+                   +SSMseasonal(period=4,type="distinct",Q=diag(c(2,3)),P1=diag(c(2,2,2,3,3,3)))
+                   +SSMseasonal(period=5,type="common",sea.type="trig")
+                   +SSMseasonal(period=5,type="distinct",sea.type="trig",Q=diag(c(0.1,0.2)),
+                                P1=diag(rep(c(0.1,0.2),each=4)))
+                   +SSMarima(ar=0.9,ma=0.2)+SSMregression(~-1+t2,index=1,Q=1)
+    ),not(gives_warning()))
+  expect_that(print(model),not(gives_warning()))
+  expect_that(logLik(model),not(gives_warning()))
+  expect_equivalent(logLik(model),-442.270604625564)
+  expect_that(out<-KFS(model,filtering=c("state","mean"),smoothing=c("state","mean","disturbance")),
+              not(gives_warning()))
+  expect_equal(out$d,11)
+  expect_equal(out$j,1)
+})
 
-model1<-SSModel(t~SSMcycle(period=10, type='common')+SSMcycle(period=10, type='distinct')+
-                  SSMtrend(2,type="common")+SSMtrend(2,type="distinct")+
-                  SSMseasonal(period=4,type="common")+SSMseasonal(period=4,type="distinct")+
-                  SSMseasonal(period=4,type="common",sea.type="trig")+
-                  SSMseasonal(period=4,type="distinct",sea.type="trig")
-                )
+test_that("StructTS and KFS give equivalent results",{
+require(graphics)
+trees <- window(treering, start = 0)
+fit <- StructTS(trees, type = "level")
+# Construct model per StructTS:
+model<-SSModel(trees~SSMtrend(degree=1,Q=fit$coef[1],P1=fit$model0$P,a1=fit$model0$a),H=fit$coef[2])
+out<-KFS(model)
+expect_equivalent(c(fitted(fit)),coef(out,filtered=TRUE)[-1])
+expect_equivalent(fit$loglik,logLik(model))
+
+fit2 <- StructTS(trees, type = "level",fixed=c(NA,1e-15))
+model<-SSModel(trees~SSMtrend(degree=1,Q=fit2$coef[1],P1=fit2$model0$P,a1=fit2$model0$a),H=fit2$coef[2])
+expect_less_than(abs(fit2$loglik-logLik(model)),1e-4)
+fit2 <- StructTS(trees, type = "level",fixed=c(1e-15,NA))
+model<-SSModel(trees~SSMtrend(degree=1,Q=fit2$coef[1],P1=fit2$model0$P,a1=fit2$model0$a),H=fit2$coef[2])
+expect_less_than(abs(fit2$loglik-logLik(model)),1e-11)
+
+model<-SSModel(trees~SSMtrend(degree=1,Q=NA),H=NA)
+fit2<-fitSSM(model,inits=c(0,0),method="BFGS")
+expect_equivalent(c(fit2$model$H),fit$coef[2])
+expect_less_than(abs(c(fit2$model$Q)-fit$coef[1]),1e-9)
diff --git a/tests/testthat/glmfits.R b/tests/testthat/glmfits.R
index 4dcb0e9..c1f3dc1 100644
--- a/tests/testthat/glmfits.R
+++ b/tests/testthat/glmfits.R
@@ -161,7 +161,7 @@ test_that("negative binomial GLM fitting works properly",{
     model$P1inf[]<-0
     model$a1[]<-pars[2:15]
     if(estimate)
-      return(-logLik(model,check=TRUE,nsim=0,tol=1e-5))
+      return(-logLik(model,check=TRUE,nsim=0))
     model
   }