fitSSM: 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

model

Model object of class SSModel.

inits

Initial values for optim.

updatefn

User defined function which updates the model given the parameters. Must be of form updatefn(pars, model, ...), where ... correspond to optional additional arguments. Function should return the original model with updated parameters. See details for description of the default updatefn.

checkfn

Optional function of form checkfn(model) for checking the validity of the model. Should return TRUE if the model is valid, and FALSE otherwise. See details.

update_args

Optional list containing additional arguments to updatefn.

...

Further arguments for functions optim and logLik.SSModel, such as nsim = 1000, marginal = TRUE, and method = "BFGS".

Value

A list with elements

optim.out

Output from function optim.

model

Model with estimated parameters.

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 $\sigma$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.

Examples

Run this code

# NOT RUN {
# 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)
# }