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.06679233
#> mu1 -0.02391857 -0.02392656 -0.02391860
#> mu2 0.07727041 0.07728411 0.07727031
#> scale 1.75585289 1.75541862 1.75585346
#> shape -0.18130928 -0.18112018 -0.18130938
#>
#> 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
#> Warning: convergence not reached in optimisation
if (FALSE) {
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
#> Date Period Level Quant L U
#> 1 1920-01-01 5 95% 35.56868 34.40283 36.70862
#> 4 1920-01-01 10 95% 36.50728 35.10215 37.62305
#> 7 1920-01-01 20 95% 37.29521 35.64052 38.38084
#> 10 1920-01-01 50 95% 38.17370 36.19522 39.26107
#> 13 1920-01-01 100 95% 38.74129 36.52633 40.04037
#> 16 1920-01-01 150 95% 39.04057 36.69178 40.50197
#> 19 1920-01-01 200 95% 39.23967 36.79821 40.86567
#> 22 1920-01-01 500 95% 39.80800 37.08048 42.02431
#> 25 1920-01-01 1000 95% 40.17897 37.24080 42.80505
#> 2 1970-01-01 5 95% 34.37272 33.61483 35.08379
#> 5 1970-01-01 10 95% 35.31132 34.43237 36.12135
#> 8 1970-01-01 20 95% 36.09925 34.96321 36.88258
#> 11 1970-01-01 50 95% 36.97774 35.57076 38.05900
#> 14 1970-01-01 100 95% 37.54533 35.87251 38.84861
#> 17 1970-01-01 150 95% 37.84461 36.01993 39.24529
#> 20 1970-01-01 200 95% 38.04371 36.11343 39.52153
#> 23 1970-01-01 500 95% 38.61204 36.38916 40.46521
#> 26 1970-01-01 1000 95% 38.98300 36.55485 41.27273
#> 3 2020-01-01 5 95% 37.04023 35.99739 37.99837
#> 6 2020-01-01 10 95% 37.97883 36.80519 38.98023
#> 9 2020-01-01 20 95% 38.76676 37.45656 39.85428
#> 12 2020-01-01 50 95% 39.64525 38.18986 40.97480
#> 15 2020-01-01 100 95% 40.21284 38.42936 41.79168
#> 18 2020-01-01 150 95% 40.51212 38.53751 42.23072
#> 21 2020-01-01 200 95% 40.71122 38.61795 42.53176
#> 24 2020-01-01 500 95% 41.27955 38.90991 43.43648
#> 27 2020-01-01 1000 95% 41.65052 39.08464 44.06995