Bug fixes, improved numerical accuracy

See changelog.

ChangeLog

@@ -1,10 +1,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 
+ * Tweaked the underlying algorithms for increased numerical stability of filtering and smoothing 
+   in KFS. Note that it is still possible that exact diffuse initialization fails due to to numerical
+   issues whereas traditional 'big value' approach works and vice versa.
+ * Corrected a bug in residuals.KFS which threw an error when computing recursive residuals without 
+   diffuse initialization.
+ * Corrected output of LogLik method for non-Gaussian models: It now returns -Inf only when the 
+   approximation algoritm failedcompletely (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 
+ * Fixed 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 bug in SSMcycle which resulted erroneus system matrix T in all cases.

DESCRIPTION

@@ -7,7 +7,7 @@ Maintainer: Jouni Helske <jouni.helske@jyu.fi>
 Depends:
     R (>= 3.0.0)
 Imports: stats
-Suggests: MASS, minqa, testthat
+Suggests: MASS, nloptr, testthat
 Description: Package KFAS provides functions for Kalman filtering, smoothing,
     forecasting and simulation of multivariate exponential family state space
     models with exact diffuse initialization when distributions of some or all

R/KFAS-package.R

@@ -188,7 +188,7 @@
 #'                  data=Seatbelts,H=matrix(NA,2,2)) 
 #' 
 #' likfn<-function(pars,model,estimate=TRUE){
-#'   model$H[,,1]<-exp(0.5*pars[1:2])
+#'   diag(model$H[,,1])<-exp(0.5*pars[1:2])
 #'   model$H[1,2,1]<-model$H[2,1,1]<-tanh(pars[3])*prod(sqrt(exp(0.5*pars[1:2]))) 
 #'   model$R[28:29]<-exp(pars[4:5])
 #'   if(estimate) return(-logLik(model))

R/KFS.R

@@ -110,7 +110,7 @@ KFS <- function(model, filtering, smoothing, simplify = TRUE, transform = c("ldl
 
   if(missing(filtering)){
     if(all(model$distribution == "gaussian")){
-    filtering <- "state"
+      filtering <- "state"
     } else filtering <- "none"
   } else{
     filtering <- match.arg(arg = filtering, 
@@ -336,6 +336,7 @@ KFS <- function(model, filtering, smoothing, simplify = TRUE, transform = c("ldl
   } else KFS_transform <- "none"
   filtersignal<-("signal"%in%filtering) || ("mean"%in%filtering)
   storage.mode(ymiss) <- "integer"
+
   filterout <- .Fortran(fkfilter, NAOK = TRUE, model$y, ymiss, as.integer(tv), model$Z, model$H, model$T, 
                         model$R, model$Q, model$a1, P1 = model$P1, model$P1inf, as.integer(p), 
                         as.integer(n), as.integer(m),as.integer(k), d = integer(1), j = integer(1), 
@@ -350,7 +351,7 @@ KFS <- function(model, filtering, smoothing, simplify = TRUE, transform = c("ldl
 
   if (filterout$d == n & filterout$j == p) 
     warning("Model is degenerate, diffuse phase did not end.")
-  if (filterout$d > 0 & m > 1 & !isTRUE(all.equal(min(apply(filterout$Pinf, 3, diag)), 0))) 
+  if (filterout$d > 0 & m > 1 & min(apply(filterout$Pinf, 3, diag))<0) 
     warning("Possible error in diffuse filtering: Negative variances in Pinf, 
             try changing the tolerance parameter tol of the model.")
   if (sum(filterout$Finf > 0) != sum(diag(model$P1inf))) 
@@ -434,6 +435,7 @@ KFS <- function(model, filtering, smoothing, simplify = TRUE, transform = c("ldl
   }
 
   if (!("none"%in%smoothing)) {
+
     smoothout <- .Fortran(fgsmoothall, NAOK = TRUE, ymiss, as.integer(tv), model$Z, model$H, 
                           model$T, model$R, model$Q, as.integer(p), as.integer(n), as.integer(m), 
                           as.integer(k), filterout$d, filterout$j, filterout$a, filterout$P, 
@@ -449,13 +451,14 @@ KFS <- function(model, filtering, smoothing, simplify = TRUE, transform = c("ldl
                           epshat = array(0, dim = c(p, n)), V_eps = array(0, dim = c(p, n)), 
                           etahat = array(0, dim = c(k, n)), V_eta = array(0, dim = c(k, k, n)), 
                           thetahat = array(0, dim = c(p, n)), V_theta = array(0, dim = c(p, p, n)), 
-                          as.integer(KFS_transform=="ldl" && ("signal" %in% smoothing || "mean" %in% smoothing)), 
-{if (KFS_transform=="ldl" && ("signal" %in% smoothing || "mean" %in% smoothing)) out$model$Z else double(1)}, 
+                          as.integer(KFS_transform=="ldl" && ("signal" %in% smoothing || "mean" %in% smoothing)), {if (KFS_transform=="ldl" && ("signal" %in% smoothing || "mean" %in% smoothing)) out$model$Z else double(1)}, 
                           as.integer(dim(out$model$Z)[3]>1), as.integer(KFS_transform != "augment"), 
                           as.integer("state" %in% smoothing), as.integer("disturbance" %in% smoothing), 
                           as.integer(("signal" %in% smoothing || "mean" %in% smoothing)))
 
-
+    if (m > 1 & min(apply(smoothout$V, 3, diag))<0) 
+      warning("Possible error in smoothing: Negative variances in V, 
+            try changing the tolerance parameter tol of the model.")
 
     if ("state" %in% smoothing) {
       out$alphahat <- ts(t(smoothout$alphahat),start=start(model$y),frequency=frequency(model$y))

R/SSMseasonal.R

@@ -17,9 +17,9 @@ SSMseasonal <- function(period, sea.type = c("dummy", "trigonometric"), type, Q,
   }
   sea.type <- match.arg(arg = sea.type)
 
-  if (!(length(period) == 1 & period > 1 & abs(period - round(period)) == 0)) 
-    stop("Period of the seasonal component must be integer larger than 1.")
-
+  if (!(length(period) == 1 & period > 1)) 
+    stop("Period of the seasonal component must be larger than 1.")
+  period <- floor(period)
   m1 <- period - 1
 
   Z_univariate <- matrix(0, 1, m1)

R/logLik.SSModel.R

@@ -140,7 +140,7 @@ logLik.SSModel <- function(object, nsim = 0, antithetics = TRUE, theta, check.mo
                     convtol, as.integer(nnd), as.integer(nsim), epsplus, etaplus, aplus1, c2, object$tol, info = integer(1), as.integer(antithetics), 
                     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.")
+      warning("Non-finite difference in approximation algoritm.")
       return(-.Machine$double.xmax^0.75)
     }
     if(out$maxiter==maxiter){

R/residuals.KFS.R

@@ -47,8 +47,10 @@ residuals.KFS <- function(object, type = c("recursive","deviance", "pearson", "r
              if (is.null(object[["a",exact=TRUE]])) 
                stop("KFS object needs to contain filtered estimates of states. ")
              series <- object$v
+             if(object$d>0){
              series[1:(object$d - 1), ] <- NA
              series[object$d, 1:object$j] <- NA
+             }
              series
            },
            response = {

man/predict.SSModel.Rd

@@ -3,8 +3,8 @@
 \alias{predict.SSModel}
 \title{State Space Model Predictions}
 \usage{
-\method{predict}{SSModel}(object, newdata, interval = c("none", "confidence",
-  "prediction"), level = 0.95, type = c("response", "link"),
+\method{predict}{SSModel}(object, newdata, n.ahead, interval = c("none",
+  "confidence", "prediction"), level = 0.95, type = c("response", "link"),
   states = NULL, se.fit = FALSE, nsim = 0, prob = TRUE, maxiter = 50,
   ...)
 }
@@ -13,7 +13,15 @@
 
   \item{newdata}{A compatible \code{SSModel} object to be
   added in the end of the old object for which the
-  predictions are required.}
+  predictions are required. If omitted, predictions are
+  either for the whole data (fitted values), or if argument
+  \code{n.ahead} is given, \code{n.ahead} time steps
+  ahead.}
+
+  \item{n.ahead}{Number of steps ahead at which to predict.
+  Only used if \code{newdata} is omitted. Note that when
+  using \code{n.ahead}, object cannot contain time varying
+  system matrices.}
 
   \item{interval}{Type of interval calculation.}
 
@@ -70,4 +78,20 @@ distribution \eqn{p(y|\theta^i)}.
 If no simulations are used, the standard errors in response
 scale are computed using delta method.
 }
+\examples{
+\dontrun{
+set.seed(1)
+x<-runif(n=100,min=1,max=3)
+y<-rpois(n=100,lambda=exp(-1+x))
+model<-SSModel(y~x,distribution="poisson")
+xnew<-seq(0.5,3.5,by=0.1)
+newdata<-SSModel(rep(NA,length(xnew))~xnew,distribution="poisson")
+pred<-predict(model,newdata=newdata,interval="prediction",level=0.9,nsim=1000)
+plot(x=x,y=y,pch=19,ylim=c(0,25),xlim=c(0.5,3.5))
+matlines(x=xnew,y=pred,col=c(2,2,2),lty=c(1,2,2),type="l")
+
+model<-SSModel(Nile~SSMtrend(1,Q=1469),H=15099)
+pred<-predict(model,n.ahead=10,interval="prediction",level=0.9)
+}
+}
 

src/approx.f95

@@ -34,7 +34,8 @@ subroutine approx(yt, ymiss, timevar, zt, tt, rtv, ht, qt, a1, p1,p1inf, p,n,m,r
     !compute rqr
     tv = max(timevar(4),timevar(5))
     do i=1, (n-1)*tv+1
-        call dsymm('r','u',m,r,1.0d0,qt(:,:,(i-1)*timevar(5)+1),r,rtv(:,:,(i-1)*timevar(4)+1),m,0.0d0,mr,m)
+        call dgemm('n','n',m,r,r,1.0d0,rtv(:,:,(i-1)*timevar(4)+1),m,&
+        qt(:,:,(i-1)*timevar(5)+1),r,0.0d0,mr,m)
         call dgemm('n','t',m,m,r,1.0d0,mr,m,rtv(:,:,(i-1)*timevar(4)+1),m,0.0d0,rqr(:,:,i),m)
     end do
 
@@ -68,18 +69,13 @@ subroutine approx(yt, ymiss, timevar, zt, tt, rtv, ht, qt, a1, p1,p1inf, p,n,m,r
             case(4)
                 do i=1,n
                     if(ymiss(i,j).EQ.0) then
-                        !ht(j,j,i) = exp(theta(i,j))/(u(i,j)*yt(i,j))
-                        !ytilde(i,j) = theta(i,j)+1.0d0-exp(theta(i,j))/yt(i,j)
-                        ht(j,j,i) =1.0d0/u(i,j) !1.0d0
+                        ht(j,j,i) =1.0d0/u(i,j)
                         ytilde(i,j) = theta(i,j)+yt(i,j)/exp(theta(i,j))-1.0d0
                     end if
                 end do
             case(5)
                 do i=1,n
                     if(ymiss(i,j).EQ.0) then
-                        !ht(j,j,i) = (exp(theta(i,j))+u(i,j))**2/(u(i,j)*exp(theta(i,j))*(yt(i,j)+u(i,j)))
-                        !ytilde(i,j) = theta(i,j) + ht(j,j,i)*u(i,j)*(yt(i,j)-exp(theta(i,j)))/(u(i,j)+exp(theta(i,j)))
-                        !ht(j,j,i) = u(i,j)*exp(theta(i,j))/(u(i,j)+exp(theta(i,j)))
                         ht(j,j,i) = (1.0d0/u(i,j)+1.0d0/exp(theta(i,j)))
                         ytilde(i,j) = theta(i,j)+yt(i,j)/exp(theta(i,j))-1.0d0
                     end if
@@ -89,14 +85,12 @@ subroutine approx(yt, ymiss, timevar, zt, tt, rtv, ht, qt, a1, p1,p1inf, p,n,m,r
 
 
     muhat = theta
-        call mu(dist,u,n,p,muhat)
-        call deviance(yt,muhat,u,ymiss,n,p,dist,devold)
+    call mu(dist,u,n,p,muhat)
+    call deviance(yt,muhat,u,ymiss,n,p,dist,devold)
 
     do while(diff > convtol .AND. k < maxiter)
 
         k=k+1
-
-
         rankp2 = rankp
 
         call kfstheta(ytilde, ymiss, timevar, zt, ht,tt, rtv,qt,rqr, a1, p1, p1inf, &
@@ -121,18 +115,12 @@ subroutine approx(yt, ymiss, timevar, zt, tt, rtv, ht, qt, a1, p1,p1inf, p,n,m,r
                 case(4)
                     do i=1,n
                         if(ymiss(i,j).EQ.0) then
-                               ! ht(j,j,i) = exp(theta(i,j))/(u(i,j)*yt(i,j))
-                               ! ytilde(i,j) = theta(i,j)+1.0d0-exp(theta(i,j))/yt(i,j)
-                               !ht(j,j,i) = 1.0d0
                             ytilde(i,j) = theta(i,j)+yt(i,j)/exp(theta(i,j))-1.0d0
                         end if
                     end do
                 case(5)
                     do i=1,n
                         if(ymiss(i,j).EQ.0) then
-                                !ht(j,j,i) = (exp(theta(i,j))+u(i,j))**2/(u(i,j)*exp(theta(i,j))*(yt(i,j)+u(i,j)))
-                                !ytilde(i,j) = theta(i,j) + ht(j,j,i)*u(i,j)*(yt(i,j)-exp(theta(i,j)))/(u(i,j)+exp(theta(i,j)))
-                            !ht(j,j,i) = u(i,j)*exp(theta(i,j))/(u(i,j)+exp(theta(i,j)))
                             ht(j,j,i) = (1.0d0/u(i,j)+1.0d0/exp(theta(i,j)))
                             ytilde(i,j) = theta(i,j)+yt(i,j)/exp(theta(i,j))-1.0d0
                         end if