View source: R/logLik.SSModel.R
logLik.SSModel | R Documentation |
Function logLik.SSmodel
computes the log-likelihood value of a state
space model.
## S3 method for class 'SSModel'
logLik(
object,
marginal = FALSE,
nsim = 0,
antithetics = TRUE,
theta,
check.model = TRUE,
transform = c("ldl", "augment"),
maxiter = 50,
seed,
convtol = 1e-08,
expected = FALSE,
H_tol = 1e+15,
transform_tol,
...
)
object |
State space model of class |
marginal |
Logical. Compute marginal instead of diffuse likelihood (see
details). Default is |
nsim |
Number of independent samples used in estimating the log-likelihood of the non-Gaussian state space model. Default is 0, which gives good starting value for optimization. Only used for non-Gaussian model. |
antithetics |
Logical. If TRUE, two antithetic variables are used in simulations, one for location and another for scale. Default is TRUE. Only used for non-Gaussian model. |
theta |
Initial values for conditional mode theta. Only used for non-Gaussian model. |
check.model |
Logical. If TRUE, function |
transform |
How to transform the model in case of non-diagonal
covariance matrix |
maxiter |
The maximum number of iterations used in linearisation. Default is 50. Only used for non-Gaussian model. |
seed |
The value is used as a seed via |
convtol |
Tolerance parameter for convergence checks for Gaussian approximation. Only used for non-Gaussian model. |
expected |
Logical value defining the approximation of H_t in case of Gamma
and negative binomial distribution. Default is |
H_tol |
Tolerance parameter for check |
transform_tol |
Tolerance parameter for LDL decomposition in case of a
non-diagonal H and |
... |
Ignored. |
As KFAS is based on diffuse initialization, the log-likelihood is also diffuse,
which coincides with restricted likelihood (REML) in an appropriate (mixed)
models. However, in typical REML estimation constant term log|X'X|
is
omitted from the log-likelihood formula. Similar term is also missing in
diffuse log-likelihood formulations of state space models, but unlike in simpler
linear models this term is not necessarily constant. Therefore omitting this
term can lead to suboptimal results in model estimation if there is unknown
parameters in diffuse parts of Zt or Tt (Francke et al. 2011). Therefore
so called marginal log-likelihood (diffuse likelihood + extra term) is
recommended. See also Gurka (2006) for model comparison in mixed model
settings using REML with and without the additional (constant) term.
The marginal log-likelihood can be computed by setting marginal = TRUE
.
Note that for non-Gaussian models with importance sampling 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.
Log-likelihood of the model.
Francke, M. K., Koopman, S. J. and De Vos, A. F. (2010),
Likelihood functions for state space models with diffuse initial
conditions. Journal of Time Series Analysis, 31: 407–414.
Gurka, M. J (2006), Selecting the Best Linear Mixed Model Under REML. The
American Statistician, Vol. 60.
Casals, J., Sotoca, S., Jerez, M. (2014), Minimally conditioned likelihood for a nonstationary state space model. Mathematics and Computers in Simulation, Vol. 100.
# Example of estimating AR model with covariates, and how to deal with possible
# non-stationarity in optimization.
set.seed(1)
x <- rnorm(100)
y <- 2 * x + arima.sim(n = 100, model = list(ar = c(0.5, -0.3)))
model<- SSModel(y ~ SSMarima(ar = c(0.5, -0.3), Q = 1) + x, H = 0)
obj <- function(pars, model, estimate = TRUE) {
#guard against stationarity
armamod <- try(SSMarima(ar = artransform(pars[1:2]), Q = exp(pars[3])), silent = TRUE)
if(class(armamod) == "try-error") {
return(Inf)
} else {
# use advanced subsetting method for SSModels, see ?`[.SSModel`
model["T", states = "arima"] <- armamod$T
model["Q", eta = "arima"] <- armamod$Q
model["P1", states = "arima"] <- armamod$P1
if(estimate) {
-logLik(model)
} else {
model
}
}
}
fit <- optim(p = c(0.5,-0.5,1), fn = obj, model = model, method ="BFGS")
model <- obj(fit$par, model, FALSE)
model$T
model$Q
coef(KFS(model), last = TRUE)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.