Nothing
R/main_estimate.R
In gasmodel: Generalized Autoregressive Score Models
Defines functions
plot.gas
logLik.gas
residuals.gas
fitted.gas
vcov.gas
coef.gas
summary.gas
print.gas
gas
Documented in
gas
# ESTIMATION OF GAS MODEL FUNCTION
# Estimate GAS Model -----------------------------------------------------------
#' @title Estimate GAS Model
#'
#' @description
#' A versatile function for estimation of generalized autoregressive score (GAS) models of Creal et al. (2013) and Harvey (2013).
#' Model specification allows for various conditional distributions, different parametrizations, exogenous variables, higher score and autoregressive orders, custom and unconditional initial values of time-varying parameters, fixed and bounded values of coefficients, and NA values.
#' Model estimation is performed by the maximum likelihood method and the Hessian matrix.
#' The function can be supplied with any optimization and Hessian functions.
#'
#' @param y A time series. For univariate time series, a numeric vector or a matrix with a single column. For multivariate times series, a numeric matrix with observations in rows.
#' @param x Optional exogenous variables. For a single variable common for all time-varying parameters, a numeric vector. For multiple variables common for all time-varying parameters, a numeric matrix with observations in rows. For individual variables for each time-varying parameter, a list of numeric vectors or matrices in the above form. The number of observation must be equal to the number of observations of \code{y}.
#' @param distr A conditional distribution. See \code{\link[gasmodel:distr]{distr()}} for available distributions.
#' @param param A parametrization of the conditional distribution. If \code{NULL}, default parametrization is used. See \code{\link[gasmodel:distr]{distr()}} for available parametrizations.
#' @param scaling A scaling function for the score. The supported scaling functions are the unit scaling (\code{scaling = "unit"}), the inverse of the Fisher information matrix scaling (\code{scaling = "fisher_inv"}), and the inverse square root of the Fisher information matrix scaling (\code{scaling = "fisher_inv_sqrt"}). The latter two scalings use the Fisher information for the time-varying parameters only. For the full Fisher information matrix for both time-varying and static parameters, there are the \code{"full_fisher_inv"} and \code{"full_fisher_inv_sqrt"} scalings. For the individual Fisher information for each parameter, there are the \code{"diag_fisher_inv"} and \code{"diag_fisher_inv_sqrt"} scalings. Note that when the parametrization is orthogonal (see \code{\link[gasmodel:distr]{distr()}}), there are no differences between these scaling variants.
#' @param regress A specification of the regression and dynamic equation with regard to exogenous variables. The supported specifications are exogenous variables and dynamics within the same equation (\code{regress = "joint"}) and separate equations for exogenous variables and dynamics in the fashion of regression models with dynamic errors (\code{regress = "sep"}). In a stationary model without exogenous variables, the two specifications are equivalent, although with differently parametrized intercept.
#' @param p A score order. For order common for all parameters, a numeric vector of length 1. For individual order for each parameter, a numeric vector of length equal to the number of parameters. Defaults to \code{1L}.
#' @param q An autoregressive order. For order common for all parameters, a numeric vector of length 1. For individual order for each parameter, a numeric vector of length equal to the number of parameters. Defaults to \code{1L}.
#' @param par_static An optional logical vector indicating static parameters. Overrides \code{x}, \code{p}, and \code{q}.
#' @param par_link An optional logical vector indicating whether the logarithmic/logistic link should be applied to restricted parameters in order to obtain unrestricted values. Defaults to applying the logarithmic/logistic link for time-varying parameters and keeping the original link for constant parameters.
#' @param par_init An optional numeric vector of initial values of time-varying parameters. For \code{NA} values or when \code{NULL}, set initial values to unconditional values of time-varying parameters. For example, in the case of GAS(1,1) model with \code{regress = "joint"}, to \code{omega / (1 - phi1)}. Not to be confused with starting values for the optimization \code{coef_start}.
#' @param lik_skip A numeric value specifying the number of skipped observations at the beginning of the time series or after \code{NA} values in the likelihood computation. Defaults to \code{0L}, i.e. the full likelihood. If \code{NULL}, it is selected as \code{max(p,q)}, i.e. the conditional likelihood.
#' @param coef_fix_value,coef_fix_other,coef_fix_special,coef_bound_lower,coef_bound_upper Restrictions on the coefficients, see the section \emph{Restricting the Coefficients} below for more details.
#' @param coef_start,optim_function,optim_arguments,hessian_function,hessian_arguments,print_progress Controls for numerical procedures; see the section \emph{Controlling Numerical Procedures} below for more details.
#'
#' @details
#' The generalized autoregressive score (GAS) models of Creal et al. (2013) and Harvey (2013), also known as dynamic conditional score (DCS) models or score-driven (SD) models, have established themselves as a useful modern framework for time series modeling.
#'
#' The GAS models are observation-driven models allowing for any underlying probability distribution \eqn{p(y_t|f_t)} with any time-varying parameters \eqn{f_t} for time series \eqn{y_t}.
#' They capture the dynamics of time-varying parameters using the autoregressive term and the lagged score, i.e. the gradient of the log-likelihood function.
#' Exogenous variables can also be included.
#' Specifically, time-varying parameters \eqn{f_t} follow the recursion
#' \deqn{f_t = \omega + \sum_{i=1}^M \beta_i x_{ti} + \sum_{j=1}^P \alpha_j S(f_{t-j}) \nabla(y_{t-j}, f_{t-j}) + \sum_{k=1}^Q \varphi_k f_{t-k},}{f_t = \omega + \sum_{i=1}^M \beta_i x_{ti} + \sum_{j=1}^P \alpha_j S(f_{t-j}) ∇(y_{t-j}, f_{t-j}) + \sum_{k=1}^Q \phi_k f_{t-k},}
#' where \eqn{\omega} is the intercept, \eqn{\beta_i} are the regression parameters, \eqn{\alpha_j} are the score parameters, \eqn{\varphi_k}{\phi_k} are the autoregressive parameters, \eqn{x_{ti}} are the exogenous variables, \eqn{S(f_t)} is a scaling function for the score, and \eqn{\nabla(y_t, f_t)}{∇(y_t, f_t)} is the score given by
#' \deqn{\nabla(y_t, f_t) = \frac{\partial \ln p(y_t | f_t)}{\partial f_t}.}{∇(y_t, f_t) = d ln p(y_t | f_t) / d f_t.}
#' In the case of a single time-varying parameter, \eqn{\omega}, \eqn{\beta_i}, \eqn{\alpha_j}, \eqn{\varphi_k}, \eqn{x_{ti}}, \eqn{S(f_t)}, and \eqn{\nabla(y_t, f_t)} are all scalar.
#' In the case of multiple time-varying parameters, \eqn{x_{ti}} are scalar, \eqn{\omega}, \eqn{\beta_i}, and \eqn{\nabla(y_{t - j}, f_{t - j})} are vectors, \eqn{\alpha_j} and \eqn{\varphi_k} are diagonal matrices, and \eqn{S(f_t)} is a square matrix.
#' Alternatively, a different model can be obtained by defining the recursion in the fashion of regression models with dynamic errors as
#' \deqn{f_t = \omega + \sum_{i=1}^M \beta_i x_{ti} + e_{t}, \quad e_t = \sum_{j=1}^P \alpha_j S(f_{t-j}) \nabla(y_{t-j}, f_{t-j}) + \sum_{k=1}^Q \varphi_k e_{t-k}.}{f_t = \omega + \sum_{i=1}^M \beta_i x_{ti} + e_t, e_t = \sum_{j=1}^P \alpha_j S(f_{t-j}) ∇(y_{t-j}, f_{t-j}) + \sum_{k=1}^Q \phi_k e_{t-k}.}
#'
#' The GAS models can be straightforwardly estimated by the maximum likelihood method.
#' For the asymptotic theory regarding the GAS models and maximum likelihood estimation, see Blasques et al. (2014), Blasques et al. (2018), and Blasques et al. (2022).
#'
#' The use of the score for updating time-varying parameters is optimal in an information theoretic sense.
#' For an investigation of the optimality properties of GAS models, see Blasques et al. (2015) and Blasques et al. (2021).
#'
#' Generally, the GAS models perform quite well when compared to alternatives, including parameter-driven models.
#' For a comparison of the GAS models to alternative models, see Koopman et al. (2016) and Blazsek and Licht (2020).
#'
#' The GAS class includes many well-known econometric models, such as the generalized autoregressive conditional heteroskedasticity (GARCH) model of Bollerslev (1986), the autoregressive conditional duration (ACD) model of Engle and Russel (1998), and the Poisson count model of Davis et al. (2003).
#' More recently, a variety of novel score-driven models has been proposed, such as the Beta-t-(E)GARCH model of Harvey and Chakravarty (2008), the discrete price changes model of Koopman et al. (2018), the circular model of Harvey et al. (2024), the bivariate Poisson model of Koopman and Lit (2019), and the ranking model of Holý and Zouhar (2022).
#' For an overview of various GAS models, see Harvey (2022).
#'
#' The extensive GAS literature is listed on \href{https://www.gasmodel.com}{www.gasmodel.com}.
#'
#' @section Restricting the Coefficients:
#'
#' The coefficients can be restricted in several ways during estimation. They can be set to specific values using a numeric vector supplied via \code{coef_fix_value}. \code{NA} values indicate coefficients to be estimated.
#'
#' Furthermore, the fixed coefficients can be expressed as linear combinations of the estimated coefficients using an optional square numeric matrix of multiples of the estimated coefficients, supplied via \code{coef_fix_other}. A coefficient given by a row is fixed on the coefficient given by a column. By this logic, all rows corresponding to the estimated coefficients should contain only \code{NA} values. Likewise, all columns corresponding to the fixed coefficients should also contain only \code{NA} values.
#'
#' For convenience, there are several predefined structures for \code{coef_fix_value} and \code{coef_fix_other}, which can be specified as a character vector supplied via \code{coef_fix_special}. This is useful mainly for multidimensional models. The value \code{"panel_structure"} forces all regression, autoregression, and score coefficients to be the same for all time-varying parameters within their group. The value \code{"zero_sum_intercept"} forces all constant parameters to sum to zero within their group. The value \code{"random_walk"} forces all autoregressive coefficients to be equal to one (should be used with caution due to nonstationarity; \code{par_init} must be specified). Multiple predefined structures can be used together and can also be combined with custom \code{coef_fix_value} and \code{coef_fix_other}.
#'
#' Finally, instead of fixing the coefficients, they can be bounded. Lower bounds can be set using a numeric vector supplied via \code{coef_bound_lower}, while upper bounds can be set in a smilar way via \code{coef_bound_upper}. For example, the autoregressive coefficient can be bounded above by 1.
#'
#' @section Controlling Numerical Procedures:
#'
#' Starting values for the optimization can be supplied as a numeric vector via \code{coef_start}. If not supplied, starting values are selected from a small grid of values.
#'
#' A custom optimizer can be supplied as a function via \code{optim_function}. For suitable wrappers of common R optimization functions, see \code{\link[gasmodel:wrappers_optim]{wrappers_optim}}. Can be set to \code{NULL} if the optimal solution should not be computed, which can be useful if the goal is only to evaluate the fit for the coefficients specified in argument \code{coef_start}. An optional list of arguments to be passed to the optimization function can be supplied via \code{optim_arguments}.
#'
#' A custom Hessian computation can be supplied as a function via \code{hessian_function}. For suitable wrappers of common R Hessian functions, see \code{\link[gasmodel:wrappers_hessian]{wrappers_hessian}}. Can be set to \code{NULL} if the Hessian matrix should not be computed, which can speed up computations when asymptotic inference is not desired. An optional list of arguments to be passed to the Hessian function can be supplied via \code{hessian_arguments}.
#'
#' Setting \code{print_progress = TRUE} enables a detailed computation report that is printed progressively.
#'
#' @return A \code{list} of S3 class \code{gas} with components:
#' \item{data$y}{The time series.}
#' \item{data$x}{The exogenous variables.}
#' \item{model$distr}{The conditional distribution.}
#' \item{model$param}{The parametrization of the conditional distribution.}
#' \item{model$scaling}{The scaling function.}
#' \item{model$regress}{The specification of the regression and dynamic equation.}
#' \item{model$t}{The length of the time series.}
#' \item{model$n}{The dimension of the model.}
#' \item{model$m}{The number of exogenous variables.}
#' \item{model$p}{The score order.}
#' \item{model$q}{The autoregressive order.}
#' \item{model$par_static}{The static parameters.}
#' \item{model$par_link}{The parameters with the logarithmic/logistic links.}
#' \item{model$par_init}{The initial values of the time-varying parameters.}
#' \item{model$lik_skip}{The number of skipped observations at the beginning of the time series or after \code{NA} values in the likelihood computation.}
#' \item{model$coef_fix_value}{The values to which coefficients are fixed.}
#' \item{model$coef_fix_other}{The multiples of the estimated coefficients, which are added to the fixed coefficients.}
#' \item{model$coef_fix_special}{The predefined structures of \code{coef_fix_value} and \code{coef_fix_other}.}
#' \item{model$coef_bound_lower}{The lower bounds on coefficients.}
#' \item{model$coef_bound_upper}{The upper bounds on coefficients.}
#' \item{model$num_obs}{The actual number of observations used in the likelihood.}
#' \item{model$num_coef}{The actual number of estimated coefficients.}
#' \item{control$optim_function}{The optimization function.}
#' \item{control$optim_arguments}{The arguments which are passed to the optimization function.}
#' \item{control$hessian_function}{The Hessian function.}
#' \item{control$hessian_arguments}{The arguments which are passed to the Hessian function.}
#' \item{solution$status_start}{The status of the starting values computation.}
#' \item{solution$theta_start}{The computed starting values.}
#' \item{solution$status_optim}{The status of the optimization computation.}
#' \item{solution$theta_optim}{The computed optimal values.}
#' \item{solution$status_hessian}{The status of the Hessian computation.}
#' \item{solution$theta_hessian}{The computed Hessian.}
#' \item{fit$coef_est}{The estimated coefficients.}
#' \item{fit$coef_vcov}{The estimated variance-covariance matrix.}
#' \item{fit$coef_sd}{The estimated standard deviations.}
#' \item{fit$coef_zstat}{The statistics of the Z-test.}
#' \item{fit$coef_pval}{The p-values of the Z-test.}
#' \item{fit$par_unc}{The unconditional values of time-varying parameters.}
#' \item{fit$par_tv}{The individual values of time-varying parameter.}
#' \item{fit$score_tv}{The individual scores of time-varying parameters.}
#' \item{fit$mean_tv}{The expected values given by the model.}
#' \item{fit$var_tv}{The variances given by the model.}
#' \item{fit$resid_tv}{The residuals of the model.}
#' \item{fit$loglik_tv}{The log-likelihoods for the individual observations.}
#' \item{fit$loglik_sum}{The overall log-likelihood.}
#' \item{fit$aic}{The Akaike information criterion.}
#' \item{fit$bic}{The Bayesian information criterion.}
#'
#' @note
#' Supported generic functions for S3 class \code{gas} include \code{\link[base:summary]{summary()}}, \code{\link[base:plot]{plot()}}, \code{\link[stats:coef]{coef()}}, \code{\link[stats:vcov]{vcov()}}, \code{\link[stats:fitted]{fitted()}}, \code{\link[stats:residuals]{residuals()}}, \code{\link[stats:logLik]{logLik()}}, \code{\link[stats:AIC]{AIC()}}, \code{\link[stats:BIC]{BIC()}}, and \code{\link[stats:confint]{confint()}}.
#'
#' @references
#' Blasques, F., Gorgi, P., Koopman, S. J., and Wintenberger, O. (2018). Feasible Invertibility Conditions and Maximum Likelihood Estimation for Observation-Driven Models. \emph{Electronic Journal of Statistics}, \strong{12}(1), 1019–1052. \doi{10.1214/18-ejs1416}.
#'
#' Blasques, F., Koopman, S. J., and Lucas, A. (2014). Stationarity and Ergodicity of Univariate Generalized Autoregressive Score Processes. \emph{Electronic Journal of Statistics}, \strong{8}(1), 1088–1112. \doi{10.1214/14-ejs924}.
#'
#' Blasques, F., Koopman, S. J., and Lucas, A. (2015). Information-Theoretic Optimality of Observation-Driven Time Series Models for Continuous Responses. \emph{Biometrika}, \strong{102}(2), 325–343. \doi{10.1093/biomet/asu076}.
#'
#' Blasques, F., Lucas, A., and van Vlodrop, A. C. (2021). Finite Sample Optimality of Score-Driven Volatility Models: Some Monte Carlo Evidence. \emph{Econometrics and Statistics}, \strong{19}, 47–57. \doi{10.1016/j.ecosta.2020.03.010}.
#'
#' Blasques, F., van Brummelen, J., Koopman, S. J., and Lucas, A. (2022). Maximum Likelihood Estimation for Score-Driven Models. \emph{Journal of Econometrics}, \strong{227}(2), 325–346. \doi{10.1016/j.jeconom.2021.06.003}.
#'
#' Blazsek, S. and Licht, A. (2020). Dynamic Conditional Score Models: A Review of Their Applications. \emph{Applied Economics}, \strong{52}(11), 1181–1199. \doi{10.1080/00036846.2019.1659498}.
#'
#' Bollerslev, T. (1986). Generalized Autoregressive Conditional Heteroskedasticity. \emph{Journal of Econometrics}, \strong{31}(3), 307–327. \doi{10.1016/0304-4076(86)90063-1}.
#'
#' Creal, D., Koopman, S. J., and Lucas, A. (2013). Generalized Autoregressive Score Models with Applications. \emph{Journal of Applied Econometrics}, \strong{28}(5), 777–795. \doi{10.1002/jae.1279}.
#'
#' Davis, R. A., Dunsmuir, W. T. M., and Street, S. B. (2003). Observation-Driven Models for Poisson Counts. \emph{Biometrika}, \strong{90}(4), 777–790. \doi{10.1093/biomet/90.4.777}.
#'
#' Engle, R. F. and Russell, J. R. (1998). Autoregressive Conditional Duration: A New Model for Irregularly Spaced Transaction Data. \emph{Econometrica}, \strong{66}(5), 1127–1162. \doi{10.2307/2999632}.
#'
#' Harvey, A. C. (2013). \emph{Dynamic Models for Volatility and Heavy Tails: With Applications to Financial and Economic Time Series}. Cambridge University Press. \doi{10.1017/cbo9781139540933}.
#'
#' Harvey, A. C. (2022). Score-Driven Time Series Models. \emph{Annual Review of Statistics and Its Application}, \strong{9}(1), 321–342. \doi{10.1146/annurev-statistics-040120-021023}.
#'
#' Harvey, A. C. and Chakravarty, T. (2008). Beta-t-(E)GARCH. \emph{Cambridge Working Papers in Economics}, CWPE 0840. \doi{10.17863/cam.5286}.
#'
#' Harvey, A., Hurn, S., Palumbo, D., and Thiele, S. (2024). Modelling Circular Time Series. \emph{Journal of Econometrics}, \strong{239}(1), 105450. \doi{10.1016/j.jeconom.2023.02.016}.
#'
#' Holý, V. and Zouhar, J. (2022). Modelling Time-Varying Rankings with Autoregressive and Score-Driven Dynamics. Journal of the Royal Statistical Society: Series C (Applied Statistics), \strong{71}(5). \doi{10.1111/rssc.12584}.
#'
#' Koopman, S. J. and Lit, R. (2019). Forecasting Football Match Results in National League Competitions Using Score-Driven Time Series Models. \emph{International Journal of Forecasting}, \strong{35}(2), 797–809. \doi{10.1016/j.ijforecast.2018.10.011}.
#'
#' Koopman, S. J., Lit, R., Lucas, A., and Opschoor, A. (2018). Dynamic Discrete Copula Models for High-Frequency Stock Price Changes. \emph{Journal of Applied Econometrics}, \strong{33}(7), 966–985. \doi{10.1002/jae.2645}.
#'
#' Koopman, S. J., Lucas, A., and Scharth, M. (2016). Predicting Time-Varying Parameters with Parameter-Driven and Observation-Driven Models. \emph{Review of Economics and Statistics}, \strong{98}(1), 97–110. \doi{10.1162/rest_a_00533}.
#'
#' @seealso
#' \code{\link[gasmodel:distr]{distr()}},
#' \code{\link[gasmodel:gas_bootstrap]{gas_bootstrap()}},
#' \code{\link[gasmodel:gas_filter]{gas_filter()}},
#' \code{\link[gasmodel:gas_forecast]{gas_forecast()}},
#' \code{\link[gasmodel:gas_simulate]{gas_simulate()}},
#' \code{\link[gasmodel:wrappers_optim]{wrappers_optim}},
#' \code{\link[gasmodel:wrappers_hessian]{wrappers_hessian}}
#'
#' @examples
#' \donttest{# Load the Daily Toilet Paper Sales dataset
#' data("toilet_paper_sales")
#' y <- toilet_paper_sales$quantity
#' x <- as.matrix(toilet_paper_sales[3:9])
#'
#' # Estimate GAS model based on the negative binomial distribution
#' est_negbin <- gas(y = y, x = x, distr = "negbin", regress = "sep")
#' est_negbin
#'
#' # Obtain the estimated coefficients
#' coef(est_negbin)
#'
#' # Obtain the estimated variance-covariance matrix
#' vcov(est_negbin)
#'
#' # Obtain the log-likelihood, AIC, and BIC
#' logLik(est_negbin)
#' AIC(est_negbin)
#' BIC(est_negbin)
#'
#' # Obtain the confidence intervals of coefficients
#' confint(est_negbin)
#'
#' # Plot the time-varying parameters
#' plot(est_negbin)}
#'
#' @export
gas <- function(y, x = NULL, distr, param = NULL, scaling = "unit", regress = "joint", p = 1L, q = 1L, par_static = NULL, par_link = NULL, par_init = NULL, lik_skip = 0L, coef_fix_value = NULL, coef_fix_other = NULL, coef_fix_special = NULL, coef_bound_lower = NULL, coef_bound_upper = NULL, coef_start = NULL, optim_function = wrapper_optim_nloptr, optim_arguments = list(opts = list(algorithm = 'NLOPT_LN_NELDERMEAD', xtol_rel = 0, maxeval = 1e6)), hessian_function = wrapper_hessian_stats, hessian_arguments = list(), print_progress = FALSE) {
# Load auxiliary variables:
load <- load_estimate(y = y, x = x, distr = distr, param = param, scaling = scaling, regress = regress, p = p, q = q, par_static = par_static, par_link = par_link, par_init = par_init, lik_skip = lik_skip, coef_fix_value = coef_fix_value, coef_fix_other = coef_fix_other, coef_fix_special = coef_fix_special, coef_bound_lower = coef_bound_lower, coef_bound_upper = coef_bound_upper, coef_start = coef_start, optim_function = optim_function, optim_arguments = optim_arguments, hessian_function = hessian_function, hessian_arguments = hessian_arguments, print_progress = print_progress)
data <- load$data
model <- load$model
control <- load$control
fun <- load$fun
info_distr <- load$info_distr
info_par <- load$info_par
info_coef <- load$info_coef
info_theta <- load$info_theta
comp <- load$comp
solution <- list()
# Compute starting solution:
if (comp$compute_start) {
if (comp$print_progress) { message("Computing a starting solution...") }
comp$result_start <- starting_theta(theta_start = solution$theta_start, theta_bound_lower = comp$theta_bound_lower, theta_bound_upper = comp$theta_bound_upper, data = data, model = model, fun = fun, info_distr = info_distr, info_par = info_par, info_coef = info_coef, info_theta = info_theta, print_progress = comp$print_progress)
solution$status_start <- comp$result_start$status_start
solution$theta_start <- name_vector(comp$result_start$theta_start, info_theta$theta_names)
if (solution$status_start != "success") {
warning("Computation of a starting solution ended with status '", solution$status_start, "'.")
}
} else {
solution$status_start <- "starting_values_supplied"
solution$theta_start <- name_vector(comp$theta_start, info_theta$theta_names)
}
# Compute optimal solution:
if (comp$compute_optim) {
if (comp$print_progress) { message("Computing the optimal solution...") }
comp$result_optim <- be_silent_but_theta(do.call(control$optim_function, args = c(list(obj_fun = likelihood_objective, theta_start = solution$theta_start, theta_bound_lower = comp$theta_bound_lower, theta_bound_upper = comp$theta_bound_upper, est_details = comp$est_details), control$optim_arguments)))
solution$status_optim <- comp$result_optim$status_optim
solution$theta_optim <- name_vector(comp$result_optim$theta_optim, info_theta$theta_names)
if (!(solution$status_optim %in% c("success", "objective_tolerance_reached", "desired_objective_reached", "variables_tolerance_reached"))) {
warning("Computation of the optimal solution ended with status '", solution$status_optim, "'.")
}
} else {
solution$status_optim <- "computation_skipped"
solution$theta_optim <- solution$theta_start
}
# Compute Hessian matrix:
if (comp$compute_hessian) {
if (comp$print_progress) { message("Computing the Hessian matrix...") }
comp$result_hessian <- be_silent_but_theta(do.call(control$hessian_function, args = c(list(obj_fun = likelihood_objective, theta_optim = solution$theta_optim, est_details = comp$est_details), control$hessian_arguments)))
solution$status_hessian <- comp$result_hessian$status_hessian
solution$theta_hessian <- name_matrix(comp$result_hessian$theta_hessian, info_theta$theta_names, info_theta$theta_names)
if (solution$status_hessian != "success") {
warning("Computation of the Hessian matrix ended with status '", solution$status_hessian, "'.")
}
} else {
solution$status_hessian <- "computation_skipped"
solution$theta_hessian <- name_matrix(matrix(NA_real_, nrow = info_theta$theta_num, ncol = info_theta$theta_num), info_theta$theta_names, info_theta$theta_names)
}
# Compute inference:
fit <- list()
fit$coef_est <- name_vector(convert_theta_vector_to_coef_vector(solution$theta_optim, coef_fix_value = model$coef_fix_value, coef_fix_other = model$coef_fix_other), info_coef$coef_names)
comp$eval_tv <- be_silent(likelihood_evaluate(coef = fit$coef_est, data = data, model = model, fun = fun, info_par = info_par, info_coef = info_coef))
model$num_obs <- sum(!is.na(comp$eval_tv$lik))
model$num_coef <- info_theta$theta_num
comp$theta_vcov <- matrix_inv(solution$theta_hessian) / model$num_obs
comp$struc <- convert_coef_vector_to_struc_list(coef_vec = fit$coef_est, m = model$m, p = model$p, q = model$q, par_names = info_par$par_names, par_of_coef_names = info_coef$par_of_coef_names)
fit$coef_vcov <- name_matrix(convert_theta_matrix_to_coef_matrix(comp$theta_vcov, coef_fix_value = model$coef_fix_value, coef_fix_other = model$coef_fix_other), info_coef$coef_names, info_coef$coef_names)
fit$coef_sd <- be_silent(sqrt(diag(fit$coef_vcov)))
fit$coef_zstat <- fit$coef_est / fit$coef_sd
fit$coef_pval <- 2 * stats::pnorm(-abs(fit$coef_zstat))
if (model$regress == "joint") {
fit$par_unc <- sapply(1:info_par$par_num, function(i) { (comp$struc[[i]]$omega + colMeans(data$x[[i]], na.rm = TRUE) %*% comp$struc[[i]]$beta) / (1 - sum(comp$struc[[i]]$phi)) })
} else if (model$regress == "sep") {
fit$par_unc <- sapply(1:info_par$par_num, function(i) { comp$struc[[i]]$omega + colMeans(data$x[[i]], na.rm = TRUE) %*% comp$struc[[i]]$beta })
}
# Format data:
info_data <- info_data(y = data$y, x = data$x)
data$y <- name_matrix(data$y, info_data$index_time, info_data$index_series, drop = c(FALSE, TRUE))
data$x <- name_list_of_matrices(data$x, info_par$par_names, info_data$index_time_list, info_data$index_vars_list, drop = c(FALSE, TRUE), zero = c(FALSE, TRUE))
# Format time-varying quantities:
fit$par_tv <- name_matrix(comp$eval_tv$par, info_data$index_time, info_par$par_names, drop = c(FALSE, TRUE))
fit$score_tv <- name_matrix(comp$eval_tv$score, info_data$index_time, info_par$par_names, drop = c(FALSE, TRUE))
fit$mean_tv <- name_matrix(fun$mean(comp$eval_tv$par), info_data$index_time, info_data$index_series, drop = c(FALSE, TRUE))
fit$var_tv <- name_matrix(convert_varcov_array_to_var_matrix(fun$var(comp$eval_tv$par)), info_data$index_time, info_data$index_series, drop = c(FALSE, TRUE))
fit$resid_tv <- (data$y - fit$mean_tv) / sqrt(fit$var_tv)
fit$loglik_tv <- name_vector(comp$eval_tv$lik, info_data$index_time)
# Determine fit:
fit$loglik_sum <- sum(fit$loglik_tv, na.rm = TRUE)
fit$aic <- 2 * model$num_coef - 2 * fit$loglik_sum
fit$bic <- log(model$num_obs) * model$num_coef - 2 * fit$loglik_sum
if (mean(fit$loglik_tv, na.rm = TRUE) <= -1e100) {
warning("The likelihood function has zero value. The results are not reliable.")
}
# Return results:
report <- list(data = data, model = model, control = control, solution = solution, fit = fit)
class(report) <- "gas"
return(report)
}
# ------------------------------------------------------------------------------
# Print Estimate ---------------------------------------------------------------
#' @export
print.gas <- function(x, ...) {
info_title <- info_title(distr = x$model$distr, param = x$model$param, scaling = x$model$scaling)
cat("GAS Model:", info_title$title, "\n")
cat("\n")
cat("Coefficients:", "\n")
coef_table <- cbind("Estimate" = x$fit$coef_est, "Std. Error" = x$fit$coef_sd, "Z-Test" = x$fit$coef_zstat, "Pr(>|Z|)" = x$fit$coef_pval)
stats::printCoefmat(coef_table)
cat("\n")
cat("Log-Likelihood:", format(x$fit$loglik_sum))
cat(", AIC:", format(x$fit$aic))
cat(", BIC:", format(x$fit$bic))
cat("\n")
invisible(x)
}
# ------------------------------------------------------------------------------
# Summarize Estimate -----------------------------------------------------------
#' @export
summary.gas <- function(object, ...) {
print(object)
cat("\n")
cat("Unconditional Parameters:", "\n")
print(object$fit$par_unc)
cat("\n")
cat("Time-Varying Parameters:", "\n")
print(object$fit$par_tv)
invisible(object)
}
# ------------------------------------------------------------------------------
# Obtain Coefficients ----------------------------------------------------------
#' @export
coef.gas <- function(object, ...) {
coef_est <- object$fit$coef_est
return(coef_est)
}
# ------------------------------------------------------------------------------
# Obtain Variance-Covariance Matrix --------------------------------------------
#' @export
vcov.gas <- function(object, ...) {
coef_vcov <- object$fit$coef_vcov
return(coef_vcov)
}
# ------------------------------------------------------------------------------
# Obtain Fitted Values ---------------------------------------------------------
#' @export
fitted.gas <- function(object, ...) {
mean_tv <- object$fit$mean_tv
return(mean_tv)
}
# ------------------------------------------------------------------------------
# Obtain Residuals -------------------------------------------------------------
#' @export
residuals.gas <- function(object, ...) {
resid <- object$fit$resid
return(resid)
}
# ------------------------------------------------------------------------------
# Obtain Log-Likelihood --------------------------------------------------------
#' @export
logLik.gas <- function(object, ...) {
loglik_sum <- object$fit$loglik_sum
attr(loglik_sum, "nobs") <- object$model$num_obs
attr(loglik_sum, "df") <- object$model$num_coef
class(loglik_sum) <- "logLik"
return(loglik_sum)
}
# ------------------------------------------------------------------------------
# Plot Time-Varying Parameters -------------------------------------------------
#' @importFrom ggplot2 .data
#' @export
plot.gas <- function(x, which = NULL, ...) {
par_static <- x$model$par_static
par_tv <- as.matrix(x$fit$par_tv)
par_unc <- x$fit$par_unc
par_names <- names(par_unc)
par_num <- length(par_unc)
be_silent(ts_index <- as.numeric(rownames(par_tv)))
if (any(is.na(ts_index)) || any(diff(ts_index) <= 0)) {
ts_index <- 1:nrow(par_tv)
}
gg_list <- list()
for (i in which(!par_static)) {
gg_data <- data.frame(index = ts_index, value = par_tv[, i])
gg_fig <- ggplot2::ggplot(gg_data, mapping = ggplot2::aes(.data$index, .data$value)) +
ggplot2::geom_hline(yintercept = par_unc[i], linetype = "dashed") +
ggplot2::geom_line(color = "#800000") +
ggplot2::geom_point(color = "#800000") +
ggplot2::labs(title = paste("Time-Varying Parameter", par_names[i]), x = "Time Index", y = "Parameter Value")
gg_list <- append(gg_list, list(gg_fig))
}
gg_which <- 1:length(gg_list)
if (!is.null(which)) {
gg_which <- gg_which[gg_which %in% which]
}
if (length(gg_which) == 1) {
be_silent(print(gg_list[[gg_which[1]]]))
} else if (length(gg_which) > 1) {
be_silent(print(gg_list[[gg_which[1]]]))
old_par <- grDevices::devAskNewPage(ask = TRUE)
for (i in 2:length(gg_which)) {
be_silent(print(gg_list[[gg_which[i]]]))
}
on.exit(grDevices::devAskNewPage(ask = old_par))
}
invisible(gg_list)
}
# ------------------------------------------------------------------------------
Try the gasmodel package in your browser
Any scripts or data that you put into this service are public.
gasmodel documentation built on Aug. 19, 2025, 1:15 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions plot.gas logLik.gas residuals.gas fitted.gas vcov.gas coef.gas summary.gas print.gas gas
Documented in gas
# ESTIMATION OF GAS MODEL FUNCTION
# Estimate GAS Model -----------------------------------------------------------
#' @title Estimate GAS Model
#'
#' @description
#' A versatile function for estimation of generalized autoregressive score (GAS) models of Creal et al. (2013) and Harvey (2013).
#' Model specification allows for various conditional distributions, different parametrizations, exogenous variables, higher score and autoregressive orders, custom and unconditional initial values of time-varying parameters, fixed and bounded values of coefficients, and NA values.
#' Model estimation is performed by the maximum likelihood method and the Hessian matrix.
#' The function can be supplied with any optimization and Hessian functions.
#'
#' @param y A time series. For univariate time series, a numeric vector or a matrix with a single column. For multivariate times series, a numeric matrix with observations in rows.
#' @param x Optional exogenous variables. For a single variable common for all time-varying parameters, a numeric vector. For multiple variables common for all time-varying parameters, a numeric matrix with observations in rows. For individual variables for each time-varying parameter, a list of numeric vectors or matrices in the above form. The number of observation must be equal to the number of observations of \code{y}.
#' @param distr A conditional distribution. See \code{\link[gasmodel:distr]{distr()}} for available distributions.
#' @param param A parametrization of the conditional distribution. If \code{NULL}, default parametrization is used. See \code{\link[gasmodel:distr]{distr()}} for available parametrizations.
#' @param scaling A scaling function for the score. The supported scaling functions are the unit scaling (\code{scaling = "unit"}), the inverse of the Fisher information matrix scaling (\code{scaling = "fisher_inv"}), and the inverse square root of the Fisher information matrix scaling (\code{scaling = "fisher_inv_sqrt"}). The latter two scalings use the Fisher information for the time-varying parameters only. For the full Fisher information matrix for both time-varying and static parameters, there are the \code{"full_fisher_inv"} and \code{"full_fisher_inv_sqrt"} scalings. For the individual Fisher information for each parameter, there are the \code{"diag_fisher_inv"} and \code{"diag_fisher_inv_sqrt"} scalings. Note that when the parametrization is orthogonal (see \code{\link[gasmodel:distr]{distr()}}), there are no differences between these scaling variants.
#' @param regress A specification of the regression and dynamic equation with regard to exogenous variables. The supported specifications are exogenous variables and dynamics within the same equation (\code{regress = "joint"}) and separate equations for exogenous variables and dynamics in the fashion of regression models with dynamic errors (\code{regress = "sep"}). In a stationary model without exogenous variables, the two specifications are equivalent, although with differently parametrized intercept.
#' @param p A score order. For order common for all parameters, a numeric vector of length 1. For individual order for each parameter, a numeric vector of length equal to the number of parameters. Defaults to \code{1L}.
#' @param q An autoregressive order. For order common for all parameters, a numeric vector of length 1. For individual order for each parameter, a numeric vector of length equal to the number of parameters. Defaults to \code{1L}.
#' @param par_static An optional logical vector indicating static parameters. Overrides \code{x}, \code{p}, and \code{q}.
#' @param par_link An optional logical vector indicating whether the logarithmic/logistic link should be applied to restricted parameters in order to obtain unrestricted values. Defaults to applying the logarithmic/logistic link for time-varying parameters and keeping the original link for constant parameters.
#' @param par_init An optional numeric vector of initial values of time-varying parameters. For \code{NA} values or when \code{NULL}, set initial values to unconditional values of time-varying parameters. For example, in the case of GAS(1,1) model with \code{regress = "joint"}, to \code{omega / (1 - phi1)}. Not to be confused with starting values for the optimization \code{coef_start}.
#' @param lik_skip A numeric value specifying the number of skipped observations at the beginning of the time series or after \code{NA} values in the likelihood computation. Defaults to \code{0L}, i.e. the full likelihood. If \code{NULL}, it is selected as \code{max(p,q)}, i.e. the conditional likelihood.
#' @param coef_fix_value,coef_fix_other,coef_fix_special,coef_bound_lower,coef_bound_upper Restrictions on the coefficients, see the section \emph{Restricting the Coefficients} below for more details.
#' @param coef_start,optim_function,optim_arguments,hessian_function,hessian_arguments,print_progress Controls for numerical procedures; see the section \emph{Controlling Numerical Procedures} below for more details.
#'
#' @details
#' The generalized autoregressive score (GAS) models of Creal et al. (2013) and Harvey (2013), also known as dynamic conditional score (DCS) models or score-driven (SD) models, have established themselves as a useful modern framework for time series modeling.
#'
#' The GAS models are observation-driven models allowing for any underlying probability distribution \eqn{p(y_t|f_t)} with any time-varying parameters \eqn{f_t} for time series \eqn{y_t}.
#' They capture the dynamics of time-varying parameters using the autoregressive term and the lagged score, i.e. the gradient of the log-likelihood function.
#' Exogenous variables can also be included.
#' Specifically, time-varying parameters \eqn{f_t} follow the recursion
#' \deqn{f_t = \omega + \sum_{i=1}^M \beta_i x_{ti} + \sum_{j=1}^P \alpha_j S(f_{t-j}) \nabla(y_{t-j}, f_{t-j}) + \sum_{k=1}^Q \varphi_k f_{t-k},}{f_t = \omega + \sum_{i=1}^M \beta_i x_{ti} + \sum_{j=1}^P \alpha_j S(f_{t-j}) ∇(y_{t-j}, f_{t-j}) + \sum_{k=1}^Q \phi_k f_{t-k},}
#' where \eqn{\omega} is the intercept, \eqn{\beta_i} are the regression parameters, \eqn{\alpha_j} are the score parameters, \eqn{\varphi_k}{\phi_k} are the autoregressive parameters, \eqn{x_{ti}} are the exogenous variables, \eqn{S(f_t)} is a scaling function for the score, and \eqn{\nabla(y_t, f_t)}{∇(y_t, f_t)} is the score given by
#' \deqn{\nabla(y_t, f_t) = \frac{\partial \ln p(y_t | f_t)}{\partial f_t}.}{∇(y_t, f_t) = d ln p(y_t | f_t) / d f_t.}
#' In the case of a single time-varying parameter, \eqn{\omega}, \eqn{\beta_i}, \eqn{\alpha_j}, \eqn{\varphi_k}, \eqn{x_{ti}}, \eqn{S(f_t)}, and \eqn{\nabla(y_t, f_t)} are all scalar.
#' In the case of multiple time-varying parameters, \eqn{x_{ti}} are scalar, \eqn{\omega}, \eqn{\beta_i}, and \eqn{\nabla(y_{t - j}, f_{t - j})} are vectors, \eqn{\alpha_j} and \eqn{\varphi_k} are diagonal matrices, and \eqn{S(f_t)} is a square matrix.
#' Alternatively, a different model can be obtained by defining the recursion in the fashion of regression models with dynamic errors as
#' \deqn{f_t = \omega + \sum_{i=1}^M \beta_i x_{ti} + e_{t}, \quad e_t = \sum_{j=1}^P \alpha_j S(f_{t-j}) \nabla(y_{t-j}, f_{t-j}) + \sum_{k=1}^Q \varphi_k e_{t-k}.}{f_t = \omega + \sum_{i=1}^M \beta_i x_{ti} + e_t, e_t = \sum_{j=1}^P \alpha_j S(f_{t-j}) ∇(y_{t-j}, f_{t-j}) + \sum_{k=1}^Q \phi_k e_{t-k}.}
#'
#' The GAS models can be straightforwardly estimated by the maximum likelihood method.
#' For the asymptotic theory regarding the GAS models and maximum likelihood estimation, see Blasques et al. (2014), Blasques et al. (2018), and Blasques et al. (2022).
#'
#' The use of the score for updating time-varying parameters is optimal in an information theoretic sense.
#' For an investigation of the optimality properties of GAS models, see Blasques et al. (2015) and Blasques et al. (2021).
#'
#' Generally, the GAS models perform quite well when compared to alternatives, including parameter-driven models.
#' For a comparison of the GAS models to alternative models, see Koopman et al. (2016) and Blazsek and Licht (2020).
#'
#' The GAS class includes many well-known econometric models, such as the generalized autoregressive conditional heteroskedasticity (GARCH) model of Bollerslev (1986), the autoregressive conditional duration (ACD) model of Engle and Russel (1998), and the Poisson count model of Davis et al. (2003).
#' More recently, a variety of novel score-driven models has been proposed, such as the Beta-t-(E)GARCH model of Harvey and Chakravarty (2008), the discrete price changes model of Koopman et al. (2018), the circular model of Harvey et al. (2024), the bivariate Poisson model of Koopman and Lit (2019), and the ranking model of Holý and Zouhar (2022).
#' For an overview of various GAS models, see Harvey (2022).
#'
#' The extensive GAS literature is listed on \href{https://www.gasmodel.com}{www.gasmodel.com}.
#'
#' @section Restricting the Coefficients:
#'
#' The coefficients can be restricted in several ways during estimation. They can be set to specific values using a numeric vector supplied via \code{coef_fix_value}. \code{NA} values indicate coefficients to be estimated.
#'
#' Furthermore, the fixed coefficients can be expressed as linear combinations of the estimated coefficients using an optional square numeric matrix of multiples of the estimated coefficients, supplied via \code{coef_fix_other}. A coefficient given by a row is fixed on the coefficient given by a column. By this logic, all rows corresponding to the estimated coefficients should contain only \code{NA} values. Likewise, all columns corresponding to the fixed coefficients should also contain only \code{NA} values.
#'
#' For convenience, there are several predefined structures for \code{coef_fix_value} and \code{coef_fix_other}, which can be specified as a character vector supplied via \code{coef_fix_special}. This is useful mainly for multidimensional models. The value \code{"panel_structure"} forces all regression, autoregression, and score coefficients to be the same for all time-varying parameters within their group. The value \code{"zero_sum_intercept"} forces all constant parameters to sum to zero within their group. The value \code{"random_walk"} forces all autoregressive coefficients to be equal to one (should be used with caution due to nonstationarity; \code{par_init} must be specified). Multiple predefined structures can be used together and can also be combined with custom \code{coef_fix_value} and \code{coef_fix_other}.
#'
#' Finally, instead of fixing the coefficients, they can be bounded. Lower bounds can be set using a numeric vector supplied via \code{coef_bound_lower}, while upper bounds can be set in a smilar way via \code{coef_bound_upper}. For example, the autoregressive coefficient can be bounded above by 1.
#'
#' @section Controlling Numerical Procedures:
#'
#' Starting values for the optimization can be supplied as a numeric vector via \code{coef_start}. If not supplied, starting values are selected from a small grid of values.
#'
#' A custom optimizer can be supplied as a function via \code{optim_function}. For suitable wrappers of common R optimization functions, see \code{\link[gasmodel:wrappers_optim]{wrappers_optim}}. Can be set to \code{NULL} if the optimal solution should not be computed, which can be useful if the goal is only to evaluate the fit for the coefficients specified in argument \code{coef_start}. An optional list of arguments to be passed to the optimization function can be supplied via \code{optim_arguments}.
#'
#' A custom Hessian computation can be supplied as a function via \code{hessian_function}. For suitable wrappers of common R Hessian functions, see \code{\link[gasmodel:wrappers_hessian]{wrappers_hessian}}. Can be set to \code{NULL} if the Hessian matrix should not be computed, which can speed up computations when asymptotic inference is not desired. An optional list of arguments to be passed to the Hessian function can be supplied via \code{hessian_arguments}.
#'
#' Setting \code{print_progress = TRUE} enables a detailed computation report that is printed progressively.
#'
#' @return A \code{list} of S3 class \code{gas} with components:
#' \item{data$y}{The time series.}
#' \item{data$x}{The exogenous variables.}
#' \item{model$distr}{The conditional distribution.}
#' \item{model$param}{The parametrization of the conditional distribution.}
#' \item{model$scaling}{The scaling function.}
#' \item{model$regress}{The specification of the regression and dynamic equation.}
#' \item{model$t}{The length of the time series.}
#' \item{model$n}{The dimension of the model.}
#' \item{model$m}{The number of exogenous variables.}
#' \item{model$p}{The score order.}
#' \item{model$q}{The autoregressive order.}
#' \item{model$par_static}{The static parameters.}
#' \item{model$par_link}{The parameters with the logarithmic/logistic links.}
#' \item{model$par_init}{The initial values of the time-varying parameters.}
#' \item{model$lik_skip}{The number of skipped observations at the beginning of the time series or after \code{NA} values in the likelihood computation.}
#' \item{model$coef_fix_value}{The values to which coefficients are fixed.}
#' \item{model$coef_fix_other}{The multiples of the estimated coefficients, which are added to the fixed coefficients.}
#' \item{model$coef_fix_special}{The predefined structures of \code{coef_fix_value} and \code{coef_fix_other}.}
#' \item{model$coef_bound_lower}{The lower bounds on coefficients.}
#' \item{model$coef_bound_upper}{The upper bounds on coefficients.}
#' \item{model$num_obs}{The actual number of observations used in the likelihood.}
#' \item{model$num_coef}{The actual number of estimated coefficients.}
#' \item{control$optim_function}{The optimization function.}
#' \item{control$optim_arguments}{The arguments which are passed to the optimization function.}
#' \item{control$hessian_function}{The Hessian function.}
#' \item{control$hessian_arguments}{The arguments which are passed to the Hessian function.}
#' \item{solution$status_start}{The status of the starting values computation.}
#' \item{solution$theta_start}{The computed starting values.}
#' \item{solution$status_optim}{The status of the optimization computation.}
#' \item{solution$theta_optim}{The computed optimal values.}
#' \item{solution$status_hessian}{The status of the Hessian computation.}
#' \item{solution$theta_hessian}{The computed Hessian.}
#' \item{fit$coef_est}{The estimated coefficients.}
#' \item{fit$coef_vcov}{The estimated variance-covariance matrix.}
#' \item{fit$coef_sd}{The estimated standard deviations.}
#' \item{fit$coef_zstat}{The statistics of the Z-test.}
#' \item{fit$coef_pval}{The p-values of the Z-test.}
#' \item{fit$par_unc}{The unconditional values of time-varying parameters.}
#' \item{fit$par_tv}{The individual values of time-varying parameter.}
#' \item{fit$score_tv}{The individual scores of time-varying parameters.}
#' \item{fit$mean_tv}{The expected values given by the model.}
#' \item{fit$var_tv}{The variances given by the model.}
#' \item{fit$resid_tv}{The residuals of the model.}
#' \item{fit$loglik_tv}{The log-likelihoods for the individual observations.}
#' \item{fit$loglik_sum}{The overall log-likelihood.}
#' \item{fit$aic}{The Akaike information criterion.}
#' \item{fit$bic}{The Bayesian information criterion.}
#'
#' @note
#' Supported generic functions for S3 class \code{gas} include \code{\link[base:summary]{summary()}}, \code{\link[base:plot]{plot()}}, \code{\link[stats:coef]{coef()}}, \code{\link[stats:vcov]{vcov()}}, \code{\link[stats:fitted]{fitted()}}, \code{\link[stats:residuals]{residuals()}}, \code{\link[stats:logLik]{logLik()}}, \code{\link[stats:AIC]{AIC()}}, \code{\link[stats:BIC]{BIC()}}, and \code{\link[stats:confint]{confint()}}.
#'
#' @references
#' Blasques, F., Gorgi, P., Koopman, S. J., and Wintenberger, O. (2018). Feasible Invertibility Conditions and Maximum Likelihood Estimation for Observation-Driven Models. \emph{Electronic Journal of Statistics}, \strong{12}(1), 1019–1052. \doi{10.1214/18-ejs1416}.
#'
#' Blasques, F., Koopman, S. J., and Lucas, A. (2014). Stationarity and Ergodicity of Univariate Generalized Autoregressive Score Processes. \emph{Electronic Journal of Statistics}, \strong{8}(1), 1088–1112. \doi{10.1214/14-ejs924}.
#'
#' Blasques, F., Koopman, S. J., and Lucas, A. (2015). Information-Theoretic Optimality of Observation-Driven Time Series Models for Continuous Responses. \emph{Biometrika}, \strong{102}(2), 325–343. \doi{10.1093/biomet/asu076}.
#'
#' Blasques, F., Lucas, A., and van Vlodrop, A. C. (2021). Finite Sample Optimality of Score-Driven Volatility Models: Some Monte Carlo Evidence. \emph{Econometrics and Statistics}, \strong{19}, 47–57. \doi{10.1016/j.ecosta.2020.03.010}.
#'
#' Blasques, F., van Brummelen, J., Koopman, S. J., and Lucas, A. (2022). Maximum Likelihood Estimation for Score-Driven Models. \emph{Journal of Econometrics}, \strong{227}(2), 325–346. \doi{10.1016/j.jeconom.2021.06.003}.
#'
#' Blazsek, S. and Licht, A. (2020). Dynamic Conditional Score Models: A Review of Their Applications. \emph{Applied Economics}, \strong{52}(11), 1181–1199. \doi{10.1080/00036846.2019.1659498}.
#'
#' Bollerslev, T. (1986). Generalized Autoregressive Conditional Heteroskedasticity. \emph{Journal of Econometrics}, \strong{31}(3), 307–327. \doi{10.1016/0304-4076(86)90063-1}.
#'
#' Creal, D., Koopman, S. J., and Lucas, A. (2013). Generalized Autoregressive Score Models with Applications. \emph{Journal of Applied Econometrics}, \strong{28}(5), 777–795. \doi{10.1002/jae.1279}.
#'
#' Davis, R. A., Dunsmuir, W. T. M., and Street, S. B. (2003). Observation-Driven Models for Poisson Counts. \emph{Biometrika}, \strong{90}(4), 777–790. \doi{10.1093/biomet/90.4.777}.
#'
#' Engle, R. F. and Russell, J. R. (1998). Autoregressive Conditional Duration: A New Model for Irregularly Spaced Transaction Data. \emph{Econometrica}, \strong{66}(5), 1127–1162. \doi{10.2307/2999632}.
#'
#' Harvey, A. C. (2013). \emph{Dynamic Models for Volatility and Heavy Tails: With Applications to Financial and Economic Time Series}. Cambridge University Press. \doi{10.1017/cbo9781139540933}.
#'
#' Harvey, A. C. (2022). Score-Driven Time Series Models. \emph{Annual Review of Statistics and Its Application}, \strong{9}(1), 321–342. \doi{10.1146/annurev-statistics-040120-021023}.
#'
#' Harvey, A. C. and Chakravarty, T. (2008). Beta-t-(E)GARCH. \emph{Cambridge Working Papers in Economics}, CWPE 0840. \doi{10.17863/cam.5286}.
#'
#' Harvey, A., Hurn, S., Palumbo, D., and Thiele, S. (2024). Modelling Circular Time Series. \emph{Journal of Econometrics}, \strong{239}(1), 105450. \doi{10.1016/j.jeconom.2023.02.016}.
#'
#' Holý, V. and Zouhar, J. (2022). Modelling Time-Varying Rankings with Autoregressive and Score-Driven Dynamics. Journal of the Royal Statistical Society: Series C (Applied Statistics), \strong{71}(5). \doi{10.1111/rssc.12584}.
#'
#' Koopman, S. J. and Lit, R. (2019). Forecasting Football Match Results in National League Competitions Using Score-Driven Time Series Models. \emph{International Journal of Forecasting}, \strong{35}(2), 797–809. \doi{10.1016/j.ijforecast.2018.10.011}.
#'
#' Koopman, S. J., Lit, R., Lucas, A., and Opschoor, A. (2018). Dynamic Discrete Copula Models for High-Frequency Stock Price Changes. \emph{Journal of Applied Econometrics}, \strong{33}(7), 966–985. \doi{10.1002/jae.2645}.
#'
#' Koopman, S. J., Lucas, A., and Scharth, M. (2016). Predicting Time-Varying Parameters with Parameter-Driven and Observation-Driven Models. \emph{Review of Economics and Statistics}, \strong{98}(1), 97–110. \doi{10.1162/rest_a_00533}.
#'
#' @seealso
#' \code{\link[gasmodel:distr]{distr()}},
#' \code{\link[gasmodel:gas_bootstrap]{gas_bootstrap()}},
#' \code{\link[gasmodel:gas_filter]{gas_filter()}},
#' \code{\link[gasmodel:gas_forecast]{gas_forecast()}},
#' \code{\link[gasmodel:gas_simulate]{gas_simulate()}},
#' \code{\link[gasmodel:wrappers_optim]{wrappers_optim}},
#' \code{\link[gasmodel:wrappers_hessian]{wrappers_hessian}}
#'
#' @examples
#' \donttest{# Load the Daily Toilet Paper Sales dataset
#' data("toilet_paper_sales")
#' y <- toilet_paper_sales$quantity
#' x <- as.matrix(toilet_paper_sales[3:9])
#'
#' # Estimate GAS model based on the negative binomial distribution
#' est_negbin <- gas(y = y, x = x, distr = "negbin", regress = "sep")
#' est_negbin
#'
#' # Obtain the estimated coefficients
#' coef(est_negbin)
#'
#' # Obtain the estimated variance-covariance matrix
#' vcov(est_negbin)
#'
#' # Obtain the log-likelihood, AIC, and BIC
#' logLik(est_negbin)
#' AIC(est_negbin)
#' BIC(est_negbin)
#'
#' # Obtain the confidence intervals of coefficients
#' confint(est_negbin)
#'
#' # Plot the time-varying parameters
#' plot(est_negbin)}
#'
#' @export
gas <- function(y, x = NULL, distr, param = NULL, scaling = "unit", regress = "joint", p = 1L, q = 1L, par_static = NULL, par_link = NULL, par_init = NULL, lik_skip = 0L, coef_fix_value = NULL, coef_fix_other = NULL, coef_fix_special = NULL, coef_bound_lower = NULL, coef_bound_upper = NULL, coef_start = NULL, optim_function = wrapper_optim_nloptr, optim_arguments = list(opts = list(algorithm = 'NLOPT_LN_NELDERMEAD', xtol_rel = 0, maxeval = 1e6)), hessian_function = wrapper_hessian_stats, hessian_arguments = list(), print_progress = FALSE) {
# Load auxiliary variables:
load <- load_estimate(y = y, x = x, distr = distr, param = param, scaling = scaling, regress = regress, p = p, q = q, par_static = par_static, par_link = par_link, par_init = par_init, lik_skip = lik_skip, coef_fix_value = coef_fix_value, coef_fix_other = coef_fix_other, coef_fix_special = coef_fix_special, coef_bound_lower = coef_bound_lower, coef_bound_upper = coef_bound_upper, coef_start = coef_start, optim_function = optim_function, optim_arguments = optim_arguments, hessian_function = hessian_function, hessian_arguments = hessian_arguments, print_progress = print_progress)
data <- load$data
model <- load$model
control <- load$control
fun <- load$fun
info_distr <- load$info_distr
info_par <- load$info_par
info_coef <- load$info_coef
info_theta <- load$info_theta
comp <- load$comp
solution <- list()
# Compute starting solution:
if (comp$compute_start) {
if (comp$print_progress) { message("Computing a starting solution...") }
comp$result_start <- starting_theta(theta_start = solution$theta_start, theta_bound_lower = comp$theta_bound_lower, theta_bound_upper = comp$theta_bound_upper, data = data, model = model, fun = fun, info_distr = info_distr, info_par = info_par, info_coef = info_coef, info_theta = info_theta, print_progress = comp$print_progress)
solution$status_start <- comp$result_start$status_start
solution$theta_start <- name_vector(comp$result_start$theta_start, info_theta$theta_names)
if (solution$status_start != "success") {
warning("Computation of a starting solution ended with status '", solution$status_start, "'.")
}
} else {
solution$status_start <- "starting_values_supplied"
solution$theta_start <- name_vector(comp$theta_start, info_theta$theta_names)
}
# Compute optimal solution:
if (comp$compute_optim) {
if (comp$print_progress) { message("Computing the optimal solution...") }
comp$result_optim <- be_silent_but_theta(do.call(control$optim_function, args = c(list(obj_fun = likelihood_objective, theta_start = solution$theta_start, theta_bound_lower = comp$theta_bound_lower, theta_bound_upper = comp$theta_bound_upper, est_details = comp$est_details), control$optim_arguments)))
solution$status_optim <- comp$result_optim$status_optim
solution$theta_optim <- name_vector(comp$result_optim$theta_optim, info_theta$theta_names)
if (!(solution$status_optim %in% c("success", "objective_tolerance_reached", "desired_objective_reached", "variables_tolerance_reached"))) {
warning("Computation of the optimal solution ended with status '", solution$status_optim, "'.")
}
} else {
solution$status_optim <- "computation_skipped"
solution$theta_optim <- solution$theta_start
}
# Compute Hessian matrix:
if (comp$compute_hessian) {
if (comp$print_progress) { message("Computing the Hessian matrix...") }
comp$result_hessian <- be_silent_but_theta(do.call(control$hessian_function, args = c(list(obj_fun = likelihood_objective, theta_optim = solution$theta_optim, est_details = comp$est_details), control$hessian_arguments)))
solution$status_hessian <- comp$result_hessian$status_hessian
solution$theta_hessian <- name_matrix(comp$result_hessian$theta_hessian, info_theta$theta_names, info_theta$theta_names)
if (solution$status_hessian != "success") {
warning("Computation of the Hessian matrix ended with status '", solution$status_hessian, "'.")
}
} else {
solution$status_hessian <- "computation_skipped"
solution$theta_hessian <- name_matrix(matrix(NA_real_, nrow = info_theta$theta_num, ncol = info_theta$theta_num), info_theta$theta_names, info_theta$theta_names)
}
# Compute inference:
fit <- list()
fit$coef_est <- name_vector(convert_theta_vector_to_coef_vector(solution$theta_optim, coef_fix_value = model$coef_fix_value, coef_fix_other = model$coef_fix_other), info_coef$coef_names)
comp$eval_tv <- be_silent(likelihood_evaluate(coef = fit$coef_est, data = data, model = model, fun = fun, info_par = info_par, info_coef = info_coef))
model$num_obs <- sum(!is.na(comp$eval_tv$lik))
model$num_coef <- info_theta$theta_num
comp$theta_vcov <- matrix_inv(solution$theta_hessian) / model$num_obs
comp$struc <- convert_coef_vector_to_struc_list(coef_vec = fit$coef_est, m = model$m, p = model$p, q = model$q, par_names = info_par$par_names, par_of_coef_names = info_coef$par_of_coef_names)
fit$coef_vcov <- name_matrix(convert_theta_matrix_to_coef_matrix(comp$theta_vcov, coef_fix_value = model$coef_fix_value, coef_fix_other = model$coef_fix_other), info_coef$coef_names, info_coef$coef_names)
fit$coef_sd <- be_silent(sqrt(diag(fit$coef_vcov)))
fit$coef_zstat <- fit$coef_est / fit$coef_sd
fit$coef_pval <- 2 * stats::pnorm(-abs(fit$coef_zstat))
if (model$regress == "joint") {
fit$par_unc <- sapply(1:info_par$par_num, function(i) { (comp$struc[[i]]$omega + colMeans(data$x[[i]], na.rm = TRUE) %*% comp$struc[[i]]$beta) / (1 - sum(comp$struc[[i]]$phi)) })
} else if (model$regress == "sep") {
fit$par_unc <- sapply(1:info_par$par_num, function(i) { comp$struc[[i]]$omega + colMeans(data$x[[i]], na.rm = TRUE) %*% comp$struc[[i]]$beta })
}
# Format data:
info_data <- info_data(y = data$y, x = data$x)
data$y <- name_matrix(data$y, info_data$index_time, info_data$index_series, drop = c(FALSE, TRUE))
data$x <- name_list_of_matrices(data$x, info_par$par_names, info_data$index_time_list, info_data$index_vars_list, drop = c(FALSE, TRUE), zero = c(FALSE, TRUE))
# Format time-varying quantities:
fit$par_tv <- name_matrix(comp$eval_tv$par, info_data$index_time, info_par$par_names, drop = c(FALSE, TRUE))
fit$score_tv <- name_matrix(comp$eval_tv$score, info_data$index_time, info_par$par_names, drop = c(FALSE, TRUE))
fit$mean_tv <- name_matrix(fun$mean(comp$eval_tv$par), info_data$index_time, info_data$index_series, drop = c(FALSE, TRUE))
fit$var_tv <- name_matrix(convert_varcov_array_to_var_matrix(fun$var(comp$eval_tv$par)), info_data$index_time, info_data$index_series, drop = c(FALSE, TRUE))
fit$resid_tv <- (data$y - fit$mean_tv) / sqrt(fit$var_tv)
fit$loglik_tv <- name_vector(comp$eval_tv$lik, info_data$index_time)
# Determine fit:
fit$loglik_sum <- sum(fit$loglik_tv, na.rm = TRUE)
fit$aic <- 2 * model$num_coef - 2 * fit$loglik_sum
fit$bic <- log(model$num_obs) * model$num_coef - 2 * fit$loglik_sum
if (mean(fit$loglik_tv, na.rm = TRUE) <= -1e100) {
warning("The likelihood function has zero value. The results are not reliable.")
}
# Return results:
report <- list(data = data, model = model, control = control, solution = solution, fit = fit)
class(report) <- "gas"
return(report)
}
# ------------------------------------------------------------------------------
# Print Estimate ---------------------------------------------------------------
#' @export
print.gas <- function(x, ...) {
info_title <- info_title(distr = x$model$distr, param = x$model$param, scaling = x$model$scaling)
cat("GAS Model:", info_title$title, "\n")
cat("\n")
cat("Coefficients:", "\n")
coef_table <- cbind("Estimate" = x$fit$coef_est, "Std. Error" = x$fit$coef_sd, "Z-Test" = x$fit$coef_zstat, "Pr(>|Z|)" = x$fit$coef_pval)
stats::printCoefmat(coef_table)
cat("\n")
cat("Log-Likelihood:", format(x$fit$loglik_sum))
cat(", AIC:", format(x$fit$aic))
cat(", BIC:", format(x$fit$bic))
cat("\n")
invisible(x)
}
# ------------------------------------------------------------------------------
# Summarize Estimate -----------------------------------------------------------
#' @export
summary.gas <- function(object, ...) {
print(object)
cat("\n")
cat("Unconditional Parameters:", "\n")
print(object$fit$par_unc)
cat("\n")
cat("Time-Varying Parameters:", "\n")
print(object$fit$par_tv)
invisible(object)
}
# ------------------------------------------------------------------------------
# Obtain Coefficients ----------------------------------------------------------
#' @export
coef.gas <- function(object, ...) {
coef_est <- object$fit$coef_est
return(coef_est)
}
# ------------------------------------------------------------------------------
# Obtain Variance-Covariance Matrix --------------------------------------------
#' @export
vcov.gas <- function(object, ...) {
coef_vcov <- object$fit$coef_vcov
return(coef_vcov)
}
# ------------------------------------------------------------------------------
# Obtain Fitted Values ---------------------------------------------------------
#' @export
fitted.gas <- function(object, ...) {
mean_tv <- object$fit$mean_tv
return(mean_tv)
}
# ------------------------------------------------------------------------------
# Obtain Residuals -------------------------------------------------------------
#' @export
residuals.gas <- function(object, ...) {
resid <- object$fit$resid
return(resid)
}
# ------------------------------------------------------------------------------
# Obtain Log-Likelihood --------------------------------------------------------
#' @export
logLik.gas <- function(object, ...) {
loglik_sum <- object$fit$loglik_sum
attr(loglik_sum, "nobs") <- object$model$num_obs
attr(loglik_sum, "df") <- object$model$num_coef
class(loglik_sum) <- "logLik"
return(loglik_sum)
}
# ------------------------------------------------------------------------------
# Plot Time-Varying Parameters -------------------------------------------------
#' @importFrom ggplot2 .data
#' @export
plot.gas <- function(x, which = NULL, ...) {
par_static <- x$model$par_static
par_tv <- as.matrix(x$fit$par_tv)
par_unc <- x$fit$par_unc
par_names <- names(par_unc)
par_num <- length(par_unc)
be_silent(ts_index <- as.numeric(rownames(par_tv)))
if (any(is.na(ts_index)) || any(diff(ts_index) <= 0)) {
ts_index <- 1:nrow(par_tv)
}
gg_list <- list()
for (i in which(!par_static)) {
gg_data <- data.frame(index = ts_index, value = par_tv[, i])
gg_fig <- ggplot2::ggplot(gg_data, mapping = ggplot2::aes(.data$index, .data$value)) +
ggplot2::geom_hline(yintercept = par_unc[i], linetype = "dashed") +
ggplot2::geom_line(color = "#800000") +
ggplot2::geom_point(color = "#800000") +
ggplot2::labs(title = paste("Time-Varying Parameter", par_names[i]), x = "Time Index", y = "Parameter Value")
gg_list <- append(gg_list, list(gg_fig))
}
gg_which <- 1:length(gg_list)
if (!is.null(which)) {
gg_which <- gg_which[gg_which %in% which]
}
if (length(gg_which) == 1) {
be_silent(print(gg_list[[gg_which[1]]]))
} else if (length(gg_which) > 1) {
be_silent(print(gg_list[[gg_which[1]]]))
old_par <- grDevices::devAskNewPage(ask = TRUE)
for (i in 2:length(gg_which)) {
be_silent(print(gg_list[[gg_which[i]]]))
}
on.exit(grDevices::devAskNewPage(ask = old_par))
}
invisible(gg_list)
}
# ------------------------------------------------------------------------------
Try the gasmodel package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.