Bootstrap for a TVGEV model.
A TVGEV object.
Target number of bootstrap samples.
Character. The values "param" and "NP"
can be used to chose between parametric and nonparametric bootstrap,
see Details.
Argument passed to MLE.
Logical. If TRUE the loop over the bootstrap
sample is performed using foreach with %do par%.
Further arguments passed to MLE.
A matrix containing the bootstraped parameter vectors, with one row for each bootstrap sample.
The parametric bootstrap (type = "param") is as follows: for
each bootstrap sample, a response vector
\(\mathbf{y}^\star\) is drawn from the estimated GEV
distribution as described in object; a MLE is performed to
find the bootstraped parameter vector
\(\boldsymbol{\psi}^\star\). For the nonparametric
bootstrap (type = "NP"), the generalised residuals of the
object as returned by residuals.TVGEV are resampled
and are back-transformed to simulated observations; a MLE is then
performed to find the bootstraped parameter vector
\(\boldsymbol{\psi}^\star\) as in the parametric case.
The simulated response does not cover the observations where
the original response was NA.
For some bootstrap samples, the MLE may fail, in
which case the coefficients will be ignored. So the true number
of sampled coefficient vecors will generally be smaller than the
target number as given in R.
example(TVGEV)
#>
#> TVGEV> ## transform a numeric year into a date
#> TVGEV> df <- within(TXMax_Dijon, Date <- as.Date(sprintf("%4d-01-01", Year)))
#>
#> TVGEV> df0 <- subset(df, !is.na(TXMax))
#>
#> TVGEV> ## fit a TVGEV model. Only the location parameter is TV.
#> TVGEV> t1 <- system.time(
#> TVGEV+ res1 <- TVGEV(data = df, response = "TXMax", date = "Date",
#> TVGEV+ design = breaksX(date = Date, breaks = "1970-01-01", degree = 1),
#> TVGEV+ loc = ~ t1 + t1_1970))
#>
#> TVGEV> ## The same using "nloptr" optimisation.
#> TVGEV> t2 <- system.time(
#> TVGEV+ res2 <- TVGEV(data = df, response = "TXMax", date = "Date",
#> TVGEV+ design = breaksX(date = Date, breaks = "1970-01-01", degree = 1),
#> TVGEV+ loc = ~ t1 + t1_1970,
#> TVGEV+ estim = "nloptr",
#> TVGEV+ parTrack = TRUE))
#>
#> TVGEV> ## use extRemes::fevd the required variables need to be added to the data frame
#> TVGEV> ## passed as 'data' argument
#> TVGEV> t0 <- system.time({
#> TVGEV+ df0.evd <- cbind(df0, breaksX(date = df0$Date, breaks = "1970-01-01",
#> TVGEV+ degree = 1));
#> TVGEV+ res0 <- fevd(x = df0.evd$TXMax, data = df0.evd, loc = ~ t1 + t1_1970)
#> TVGEV+ })
#>
#> TVGEV> ## compare estimate and negative log-liks
#> TVGEV> cbind("fevd" = res0$results$par,
#> TVGEV+ "TVGEV_optim" = res1$estimate,
#> TVGEV+ "TVGEV_nloptr" = res2$estimate)
#> fevd TVGEV_optim TVGEV_nloptr
#> mu0 32.06678895 32.06638460 32.06679282
#> mu1 -0.02391857 -0.02392656 -0.02391856
#> mu2 0.07727041 0.07728411 0.07727061
#> scale 1.75585289 1.75541862 1.75585353
#> shape -0.18130928 -0.18112018 -0.18131041
#>
#> TVGEV> cbind("fevd" = res0$results$value,
#> TVGEV+ "VGEV_optim" = res1$negLogLik,
#> TVGEV+ "TVGEV_nloptr" = res2$negLogLik)
#> fevd VGEV_optim TVGEV_nloptr
#> [1,] 177.2014 177.2014 177.2014
#>
#> TVGEV> ## ====================================================================
#> TVGEV> ## use a loop on plausible break years. The fitted models
#> TVGEV> ## are stored within a list
#> TVGEV> ## ====================================================================
#> TVGEV>
#> TVGEV> ## Not run:
#> TVGEV> ##D
#> TVGEV> ##D yearBreaks <- c(1940, 1950, 1955, 1960:2000, 2005, 2010)
#> TVGEV> ##D res <- list()
#> TVGEV> ##D
#> TVGEV> ##D for (ib in seq_along(yearBreaks)) {
#> TVGEV> ##D d <- sprintf("%4d-01-01", yearBreaks[[ib]])
#> TVGEV> ##D floc <- as.formula(sprintf("~ t1 + t1_%4d", yearBreaks[[ib]]))
#> TVGEV> ##D res[[d]] <- TVGEV(data = df, response = "TXMax", date = "Date",
#> TVGEV> ##D design = breaksX(date = Date, breaks = d, degree = 1),
#> TVGEV> ##D loc = floc)
#> TVGEV> ##D }
#> TVGEV> ##D
#> TVGEV> ##D ## [continuing...] ]find the model with maximum likelihood, and plot
#> TVGEV> ##D ## something like a profile likelihood for the break date considered
#> TVGEV> ##D ## as a new parameter. However, the model is not differentiable w.r.t.
#> TVGEV> ##D ## the break!
#> TVGEV> ##D
#> TVGEV> ##D ll <- sapply(res, logLik)
#> TVGEV> ##D plot(yearBreaks, ll, type = "o", pch = 21, col = "orangered",
#> TVGEV> ##D lwd = 2, bg = "gold", xlab = "break", ylab = "log-lik")
#> TVGEV> ##D grid()
#> TVGEV> ##D iMax <- which.max(ll)
#> TVGEV> ##D abline(v = yearBreaks[iMax])
#> TVGEV> ##D abline(h = ll[iMax] - c(0, qchisq(0.95, df = 1) /2),
#> TVGEV> ##D col = "SpringGreen3", lwd = 2)
#> TVGEV> ##D
#> TVGEV> ## End(Not run)
#> TVGEV>
#> TVGEV>
#> TVGEV>
bsb <- bs(res2, R = 50, estim = "nloptr")
#> Warning: convergence not reached in optimisation
if (FALSE) { # \dontrun{
library(parallel)
library(doParallel)
nc <- detectCores()
cl <- makeCluster(nc)
registerDoParallel(cl)
## findings: with 'nloptr', less than 1% of the optimisations
## diverges, while more than 10% diverged with 'optim'.
te <- system.time(bspb <- bs(res2, R = 5000, estim = "nloptr",
parallel = TRUE))
stopCluster(cl)
} # }
## add the 'bs' results to the model object as 'boot' element,
## in order to facilitate the prediction. The number of
## bootstrap replicates in the prediction will correspond
## to the number of replications stored in the model object,
## an no longer to the default
res2$boot <- bsb
predict(res2, confintMethod = "boot")
#> Since 'object' is really time-varying, the Return Levels
#> depend on the date. A default choice of dates is made here.
#> Use the 'newdate' formal to change this.
#> 'object' embeds a bootstrap distribution. Formal arguments
#> intented to be passed to `bs` (if any) will be ignored
#> Period Level Id Date L Quant U
#> 1 5 95% 1 1920-01-01 34.37176 35.56868 36.78309
#> 2 5 95% 2 1970-01-01 33.71856 34.37272 34.94116
#> 3 5 95% 3 2020-01-01 35.83685 37.04025 38.05875
#> 4 10 95% 1 1920-01-01 35.32215 36.50728 37.74905
#> 5 10 95% 2 1970-01-01 34.52955 35.31131 35.89317
#> 6 10 95% 3 2020-01-01 36.75191 37.97884 39.01963
#> 7 20 95% 1 1920-01-01 35.93496 37.29520 38.63042
#> 8 20 95% 2 1970-01-01 35.16045 36.09924 36.84051
#> 9 20 95% 3 2020-01-01 37.53289 38.76677 39.80647
#> 10 50 95% 1 1920-01-01 36.40651 38.17369 39.50825
#> 11 50 95% 2 1970-01-01 35.86145 36.97773 37.83633
#> 12 50 95% 3 2020-01-01 38.16569 39.64526 40.70356
#> 13 100 95% 1 1920-01-01 36.67347 38.74128 40.17186
#> 14 100 95% 2 1970-01-01 36.09604 37.54532 38.63932
#> 15 100 95% 3 2020-01-01 38.45821 40.21285 41.39201
#> 16 150 95% 1 1920-01-01 36.80239 39.04056 40.53249
#> 17 150 95% 2 1970-01-01 36.20348 37.84460 39.08296
#> 18 150 95% 3 2020-01-01 38.60104 40.51213 41.80834
#> 19 200 95% 1 1920-01-01 36.86953 39.23966 40.77621
#> 20 200 95% 2 1970-01-01 36.26991 38.04370 39.38652
#> 21 200 95% 3 2020-01-01 38.69160 40.71123 42.09230
#> 22 500 95% 1 1920-01-01 37.02903 39.80798 41.52190
#> 23 500 95% 2 1970-01-01 36.44974 38.61202 40.19845
#> 24 500 95% 3 2020-01-01 38.93005 41.27955 42.97116
#> 25 1000 95% 1 1920-01-01 37.11605 40.17894 42.05897
#> 26 1000 95% 2 1970-01-01 36.58897 38.98298 40.93577
#> 27 1000 95% 3 2020-01-01 39.06924 41.65051 43.60238