#' Model Estimation
#'
#' @description Estimates a model given a specification object using
#' maximum likelihood.
#' @param object an object of class \dQuote{tsissm.spec}.
#' @param solver one of either \dQuote{solnp}, \dQuote{nlminb} or \dQuote{optim}.
#' The latter uses the L-BFGS-B algorithm from the lbfgsb3c package. For option
#' \dQuote{autodiff}, valid solvers are \dQuote{nlminb} and \dQuote{nloptr}.
#' @param control solver control parameters.
#' @param autodiff whether to use automatic differentiation for estimation.
#' This makes use of the tsissmad package.
#' @param ... only additional argument which can be passed when choosing autodiff
#' is that of \dQuote{use_hessian} which tells the solver to make use of second
#' derivatives.
#' @details The maximum likelihood estimation for this model is described in the
#' online book of tsmodels.
#' @return An object of class \dQuote{tsissm.estimate}
#' @aliases estimate
#' @method estimate tsissm.spec
#' @rdname estimate
#' @export
#'
#'
estimate.tsissm.spec <- function(object, solver = "nlminb", control = list(trace = 0), autodiff = TRUE, ...)
{
# initialize parameters
if (object$seasonal$include_seasonal & object$seasonal$seasonal_type == "regular") {
if (autodiff) {
autodiff <- FALSE
warning("\nautodiff only currently supported for trigonometric seasonality (switching to non autodiff)")
}
}
estimate <- NULL
tic <- Sys.time()
solver_env = new.env(hash = TRUE)
assign("issm_llh", 1, envir = solver_env)
object$solver_env <- solver_env
pars <- object$parmatrix[estimate == 1]$initial
assign("issm_llh", 1, envir = solver_env)
if (autodiff) {
opt <- estimate_ad(object, solver = solver, control = control, ...)
object$xseed <- opt$xseed
f <- iss_filter(opt$pars, obj = object)
parameters <- NULL
object$target$y <- object$transform$transform(object$target$y_orig, f$parmatrix[parameters == "lambda"]$optimal)
f$spec <- object
f$opt <- opt$solver_out
f$opt$elapsed <- difftime(Sys.time(), tic, units = "mins")
f$hessian <- opt$hessian
f$autodiff <- TRUE
class(f) <- "tsissm.estimate"
} else {
if (solver == "nlminb") {
if (is.null(control$iter.max)) control$iter.max <- 5000
if (is.null(control$eval.max)) control$eval.max <- 5000
control2 <- list()
control2$maxit <- 100
control2$trace <- 0
opt <- optim(par = pars, control = control2, method = "Nelder-Mead", fn = loglik_fun, obj = object)
opt <- nlminb(start = opt$par, control = control, objective = loglik_fun, lower = object$parmatrix[estimate == 1]$lower,
upper = object$parmatrix[estimate == 1]$upper, obj = object)
pars <- opt$par
} else if (solver == "optim") {
if (is.null(control$maxit)) control$maxit <- 5000
control$parscale <- object$parmatrix[estimate == 1]$scale
pars <- transform_pars(pars, lower = object$parmatrix[estimate == 1]$lower, upper = object$parmatrix[estimate == 1]$upper, inverse = TRUE)
opt <- optim(par = pars, control = control, method = "BFGS", fn = loglik_fun_unc, obj = object)
pars <- opt$par
pars <- transform_pars(pars, object$parmatrix[estimate == 1]$lower, object$parmatrix[estimate == 1]$upper)
} else if (solver == "solnp") {
control2 <- list()
control2$maxit <- 100
control2$trace <- 0
opt <- optim(par = pars, control = control2, method = "Nelder-Mead", fn = loglik_fun, obj = object)
opt <- solnp(pars = opt$par, fun = loglik_fun, LB = object$parmatrix[estimate == 1]$lower,
UB = object$parmatrix[estimate == 1]$upper, obj = object, control = control)
pars <- opt$pars
}
object$xseed <- get("xseed", solver_env)
if (is.null(object$xseed)) {
if (solver == "optim" | solver == "solnp") {
opt <- nlminb(start = pars, control = control, objective = loglik_fun, lower = object$parmatrix[estimate == 1]$lower,
upper = object$parmatrix[estimate == 1]$upper, obj = object)
pars <- opt$par
} else {
control2 <- list()
control2$maxit <- 5000
control2$parscale <- object$parmatrix[estimate == 1]$scale
opt <- optim(par = pars, control = control2, method = "Nelder-Mead", fn = loglik_fun, obj = object)
pars <- opt$par
}
if (is.null(object$xseed)) {
stop("\nproblem estimating model...exiting")
}
}
f <- iss_filter(pars, obj = object)
parameters <- NULL
object$target$y <- object$transform$transform(object$target$y_orig, f$parmatrix[parameters == "lambda"]$optimal)
f$spec <- object
f$opt <- opt
f$opt$elapsed <- difftime(Sys.time(), tic, units = "mins")
f$autodiff <- FALSE
class(f) <- "tsissm.estimate"
}
return(f)
}
iss_filter <- function(pars, obj)
{
estimate <- NULL
obj$parmatrix[estimate == 1, "initial"] <- pars
pars <- obj$parmatrix$initial
names(pars) <- obj$parmatrix$parameters
Mnames <- na.omit(obj$S$pars)
S <- obj$S
S[which(!is.na(pars)),"values"] <- pars[Mnames]
##################################################################
parnames <- names(pars)
mdim = obj$dims
f <- issfilter(f0_ = S[list("F0")]$values,
f1_ = S[list("F1")]$values,
f2_ = S[list("F2")]$values,
w_ = as.numeric(S[list("w")]$values),
g_ = as.numeric(S[list("g")]$values),
y_ = as.numeric(obj$target$y),
lambda_ = as.numeric(pars["lambda"]),
xreg_ = S[list("xreg")]$values,
kappa_ = S[list("kappa")]$values,
mdim = mdim, xseed_ = obj$xseed, good_ = as.numeric(obj$good))
L <- list()
L$fitted <- obj$transform$inverse(f$fitted[-1], lambda = pars["lambda"])
L$residuals <- obj$target$y_orig - L$fitted
L$states <- f$states[-1,,drop = FALSE]
L$xseed <- f$xseed
L$loglik <- f$loglik
L$condition <- f$condition
L$w <- f$w
L$g <- f$g
L$F <- f$F
L$D <- f$D
p <- obj$parmatrix
colnames(p)[2] <- "optimal"
return(list(model = L, parmatrix = p))
}
loglik_fun <- function(pars, obj)
{
estimate <- NULL
pp <- NULL
pp <<- pars
if (any(is.na(pars))) return(Inf)
#print(pars)
solver_env <- obj$solver_env
obj$parmatrix[estimate == 1, "initial"] <- pars
pars <- obj$parmatrix$initial
names(pars) <- obj$parmatrix$parameters
Mnames <- na.omit(obj$S$pars)
S <- obj$S
S[which(!is.na(pars)),"values"] <- pars[Mnames]
##################################################################
parnames <- names(pars)
if (sum(obj$arma$order) > 0) {
if (arma_conditions(pars, parnames)) {
llh <- get("issm_llh", solver_env) + 0.25*(abs(get("issm_llh", solver_env)))
return(llh)
}
}
mdim = obj$dims
f <- issestimation(f0_ = S[list("F0")]$values,
f1_ = S[list("F1")]$values,
f2_ = S[list("F2")]$values,
w_ = as.numeric(S[list("w")]$values),
g_ = as.numeric(S[list("g")]$values),
y_ = as.numeric(obj$target$y),
lambda_ = as.numeric(pars["lambda"]),
xreg_ = S[list("xreg")]$values,
kappa_ = S[list("kappa")]$values,
mdim = mdim, good_ = as.numeric(obj$good))
if (!is.finite(f$loglik) | f$condition == 1 ) {
llh <- get("issm_llh", solver_env) + 0.25*(abs(get("issm_llh", solver_env)))
} else{
llh <- f$loglik
assign("issm_llh", llh, envir = solver_env)
}
assign("xseed", f$xseed, envir = solver_env)
return(llh)
}
loglik_fun_unc <- function(pars, obj)
{
estimate <- NULL
pars <- transform_pars(pars, obj$parmatrix[estimate == 1]$lower, obj$parmatrix[estimate == 1]$upper)
pp <- NULL
pp <<- pars
if (any(is.na(pars))) return(Inf)
#print(pars)
solver_env <- obj$solver_env
obj$parmatrix[estimate == 1, "initial"] <- pars
pars <- obj$parmatrix$initial
names(pars) <- obj$parmatrix$parameters
Mnames <- na.omit(obj$S$pars)
S <- obj$S
S[which(!is.na(pars)),"values"] <- pars[Mnames]
##################################################################
parnames <- names(pars)
if (sum(obj$arma$order) > 0) {
if (arma_conditions(pars, parnames)) {
llh <- get("issm_llh", solver_env) + 0.25*(abs(get("issm_llh", solver_env)))
return(llh)
}
}
mdim = obj$dims
f <- issestimation(f0_ = S[list("F0")]$values,
f1_ = S[list("F1")]$values,
f2_ = S[list("F2")]$values,
w_ = as.numeric(S[list("w")]$values),
g_ = as.numeric(S[list("g")]$values),
y_ = as.numeric(obj$target$y),
lambda_ = as.numeric(pars["lambda"]),
xreg_ = S[list("xreg")]$values,
kappa_ = S[list("kappa")]$values,
mdim = mdim, good_ = as.numeric(obj$good))
if (!is.finite(f$loglik) | f$condition == 1 ) {
llh <- get("issm_llh", solver_env) + 0.25*(abs(get("issm_llh", solver_env)))
} else{
llh <- f$loglik
assign("issm_llh", llh, envir = solver_env)
}
assign("xseed", f$xseed, envir = solver_env)
return(llh)
}
arma_conditions <- function(pars, parnames) {
if (any(grepl("theta",parnames))) {
ar <- pars[grepl("theta",parnames)]
if (all(ar == 0)) {
} else {
if (min(Mod(polyroot(c(1, -ar)))) < 1.01) {
return(TRUE)
}
}
}
if (any(grepl("psi",parnames))) {
ma <- pars[grepl("psi",parnames)]
if (all(ma == 0)) {
} else {
if (min(Mod(polyroot(c(1, ma)))) < 1.01) {
return(TRUE)
}
}
}
return(FALSE)
}
transform_pars <- function(pars, lower, upper, inverse = FALSE)
{
if (!inverse) {
ans <- lower + (upper - lower)/(1 + exp(-1 * pars))
} else{
ans <- -1 * log(-(upper - pars)/(-pars + lower))
}
return(ans)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.