Expectation or Moments of a TVGEV Model Object

Source: R/statsTVGEV.R

Expectation or moments of a TVGEV model.

# S3 method for TVGEV
mean(x, date = NULL, psi = NULL, ...)

# S3 method for TVGEV
moment(x, which = "variance", date = NULL, psi = NULL, ...)

Arguments

x

An TVGEV object.

date

An object that can be coerced to the "Date" class. By default the date in x is used.

psi

Vector of parameters. By default the vector of the parameters for the model x is used.

...

Not used yet.

which

Description of the moment. For now only the value "variance" is accepted, but a description for the 3-rd and 4-th order moments should be accepted in the future.

Value

An object with class "bts" representing a block time series with its value set to the expectation or the moment of the GEV distribution.

See also

quantile.TVGEV to compute quantiles of the time-varying distribution.

Examples

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> 
autoplot(mean(res1))