## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
has_bgev <- requireNamespace("bgev", quietly = TRUE)
knitr::opts_chunk$set(eval = has_bgev)
library(GEVStableGarch)

## ----check-bgev, eval = !has_bgev, echo = FALSE, results = "asis"-------------
cat("> The `bgev` package is not installed, so the code below is shown but not run.\n")

## ----build--------------------------------------------------------------------
# bgev_mean <- function(xi, delta)
#   integrate(function(p) bgev::qbgev(p, 0, 1, xi, delta), 0, 1,
#             subdivisions = 400L, rel.tol = 1e-6)$value
# 
# bgev_innov <- gs_dist(
#   name = "bgev",
#   log_density = function(z, par) {
#     d <- bgev::dbgev(z, 0, 1, par[["bgev_xi"]], par[["bgev_delta"]])
#     ifelse(d > 0, log(d), -Inf)
#   },
#   random   = function(n, par) bgev::rbgev(n, 0, 1, par[["bgev_xi"]], par[["bgev_delta"]]),
#   cdf      = function(q, par) bgev::pbgev(q, 0, 1, par[["bgev_xi"]], par[["bgev_delta"]]),
#   quantile = function(p, par) bgev::qbgev(p, 0, 1, par[["bgev_xi"]], par[["bgev_delta"]]),
#   par_names = c("bgev_xi", "bgev_delta"),
#   lower = c(bgev_xi = 0.02, bgev_delta = 0.5),
#   upper = c(bgev_xi = 0.45, bgev_delta = 5),
#   start = c(bgev_xi = 0.20, bgev_delta = 2),
#   default_delta = 2,                                   # GARCH mode: variance recursion
#   max_power = function(par) (par[["bgev_delta"]] + 1) / par[["bgev_xi"]],
#   mean = function(par) bgev_mean(par[["bgev_xi"]], par[["bgev_delta"]]))

## ----check--------------------------------------------------------------------
# gs_check_dist(bgev_innov)

## ----fit----------------------------------------------------------------------
# spec  <- gs_spec(garch = c(1, 1), dist = bgev_innov)
# model <- gs_model(spec, mu = 0, omega = 0.1, alpha = 0.1, beta = 0.6,
#                   dist_par = c(bgev_xi = 0.1, bgev_delta = 4))
# gs_stationarity(model)
# 
# x   <- gs_sim(model, n = 3000, seed = 7)$y
# fit <- gs_fit(x, spec)
# fit

## ----forecast-----------------------------------------------------------------
# predict(fit, n_ahead = 5, level = c(0.01, 0.05), n_sim = 5000, seed = 1)

## ----pit----------------------------------------------------------------------
# u <- residuals(fit, type = "pit")
# round(quantile(u, c(0.1, 0.5, 0.9)), 3)