Maximum Likelihood Estimation of a State Space Model

Description

Function fitSSM finds the maximum likelihood estimates for unknown parameters of an arbitary state space model, given the user-defined model updating function.

Usage

fitSSM(model, inits, updatefn, checkfn, update_args = NULL, ...)

Arguments

Details

Note that fitSSM actually minimizes -logLik(model), so for example the Hessian matrix returned by hessian = TRUE has an opposite sign than expected.

This function is simple wrapper around optim. For optimal performance in complicated problems, it is more efficient to use problem specific codes with calls to logLik method directly.

In fitSSM, the objective function for optim first updates the model based on the current values of the parameters under optimization, using function updatefn. Then function checkfn is used for checking that the resulting model is valid (the default checkfn checks for non-finite values and overly large (>1e7) values in covariance matrices). If checkfn returns TRUE, the log-likelihood is computed using a call -logLik(model,check.model = FALSE). Otherwise objective function returns value corresponding to .Machine$double.xmax^0.75.

The default updatefn can be used to estimate the values marked as NA in unconstrained time-invariant covariance matrices Q and H. Note that the default updatefn function cannot be used with trigonometric seasonal components as its covariance structure is of form σI, i.e. not all NA's correspond to unique value.

The code for the default updatefn can be found in the examples. As can be seen from the function definition, it is assumed that unconstrained optimization method such as BFGS is used.

Note that for non-Gaussian models derivative-free optimization methods such as Nelder-Mead might be more reliable than methods which use finite difference approximations. This is due to noise caused by the relative stopping criterion used for finding approximating Gaussian model. In most cases this does not seem to cause any problems though.

Value

A list with elements

See Also

logLik, KFAS.

Examples

# Example function for updating covariance matrices H and Q
# (also used as a default function in fitSSM)

updatefn <- function(pars, model){
  if(any(is.na(model$Q))){
    Q <- as.matrix(model$Q[,,1])
    naQd  <- which(is.na(diag(Q)))
    naQnd <- which(upper.tri(Q[naQd,naQd]) & is.na(Q[naQd,naQd]))
    Q[naQd,naQd][lower.tri(Q[naQd,naQd])] <- 0
    diag(Q)[naQd] <- exp(0.5 * pars[1:length(naQd)])
    Q[naQd,naQd][naQnd] <- pars[length(naQd)+1:length(naQnd)]
    model$Q[naQd,naQd,1] <- crossprod(Q[naQd,naQd])
  }
 if(!identical(model$H,'Omitted') && any(is.na(model$H))){#'
   H<-as.matrix(model$H[,,1])
   naHd  <- which(is.na(diag(H)))
   naHnd <- which(upper.tri(H[naHd,naHd]) & is.na(H[naHd,naHd]))
   H[naHd,naHd][lower.tri(H[naHd,naHd])] <- 0
   diag(H)[naHd] <-
     exp(0.5 * pars[length(naQd)+length(naQnd)+1:length(naHd)])
   H[naHd,naHd][naHnd] <-
     pars[length(naQd)+length(naQnd)+length(naHd)+1:length(naHnd)]
   model$H[naHd,naHd,1] <- crossprod(H[naHd,naHd])
   }
 model
}

# Example function for checking the validity of covariance matrices.

checkfn <- function(model){
  #test positive semidefiniteness of H and Q
  inherits(try(ldl(model$H[,,1]),TRUE),'try-error') ||
  inherits(try(ldl(model$Q[,,1]),TRUE),'try-error')
}


model <- SSModel(Nile ~ SSMtrend(1, Q = list(matrix(NA))), H = matrix(NA))

#function for updating the model
update_model <- function(pars, model) {
  model["H"] <- pars[1]
  model["Q"] <- pars[2]
  model
}

#check that variances are non-negative
check_model <- function(model) {
  (model["H"] > 0 && model["Q"] > 0)
}

fit <- fitSSM(inits = rep(var(Nile)/5, 2), model = model,
                 updatefn = update_model, checkfn = check_model)

Results

R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(KFAS)
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/KFAS/fitSSM.Rd_%03d_medium.png", width=480, height=480)
> ### Name: fitSSM
> ### Title: Maximum Likelihood Estimation of a State Space Model
> ### Aliases: fitSSM
> 
> ### ** Examples
> 
> 
> # Example function for updating covariance matrices H and Q
> # (also used as a default function in fitSSM)
> 
> updatefn <- function(pars, model){
+   if(any(is.na(model$Q))){
+     Q <- as.matrix(model$Q[,,1])
+     naQd  <- which(is.na(diag(Q)))
+     naQnd <- which(upper.tri(Q[naQd,naQd]) & is.na(Q[naQd,naQd]))
+     Q[naQd,naQd][lower.tri(Q[naQd,naQd])] <- 0
+     diag(Q)[naQd] <- exp(0.5 * pars[1:length(naQd)])
+     Q[naQd,naQd][naQnd] <- pars[length(naQd)+1:length(naQnd)]
+     model$Q[naQd,naQd,1] <- crossprod(Q[naQd,naQd])
+   }
+  if(!identical(model$H,'Omitted') && any(is.na(model$H))){#'
+    H<-as.matrix(model$H[,,1])
+    naHd  <- which(is.na(diag(H)))
+    naHnd <- which(upper.tri(H[naHd,naHd]) & is.na(H[naHd,naHd]))
+    H[naHd,naHd][lower.tri(H[naHd,naHd])] <- 0
+    diag(H)[naHd] <-
+      exp(0.5 * pars[length(naQd)+length(naQnd)+1:length(naHd)])
+    H[naHd,naHd][naHnd] <-
+      pars[length(naQd)+length(naQnd)+length(naHd)+1:length(naHnd)]
+    model$H[naHd,naHd,1] <- crossprod(H[naHd,naHd])
+    }
+  model
+ }
> 
> # Example function for checking the validity of covariance matrices.
> 
> checkfn <- function(model){
+   #test positive semidefiniteness of H and Q
+   inherits(try(ldl(model$H[,,1]),TRUE),'try-error') ||
+   inherits(try(ldl(model$Q[,,1]),TRUE),'try-error')
+ }
> 
> 
> model <- SSModel(Nile ~ SSMtrend(1, Q = list(matrix(NA))), H = matrix(NA))
> 
> #function for updating the model
> update_model <- function(pars, model) {
+   model["H"] <- pars[1]
+   model["Q"] <- pars[2]
+   model
+ }
> 
> #check that variances are non-negative
> check_model <- function(model) {
+   (model["H"] > 0 && model["Q"] > 0)
+ }
> 
> fit <- fitSSM(inits = rep(var(Nile)/5, 2), model = model,
+                  updatefn = update_model, checkfn = check_model)
> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>