Simulate a `pgpTList` Object. — simulate.pgpTList • dailymet

Simulate from a `pgpTList` Object. For each of the thresholds of the object, a collection of n random drawings of the exeedance events and the related values are provided.

Usage

# S3 method for pgpTList
simulate(
  object,
  nsim = 1,
  seed = NULL,
  newdata,
  lastFullYear = FALSE,
  tailOnly = TRUE,
  tau = NULL,
  trace = 0,
  how = c("list", "vector", "data.frame"),
  ...
)

Arguments

object: A pgpTList object representing a list of Poisson-GP fitted models
nsim: Number of simulation wanted.
seed: Not used yet.
newdata: An optional object that giving the "new" dates for which the simulation will be done. It can simply be an object with class "Date" or an object with class "dailyMet". Depending on the class of newdata, the variables required for the prediction (such as sine waves) will be recomputed or not. When newdata is not given or is NULL, object$data is used.
lastFullYear: Logical, used only when newdata is not provided or is NULL. When TRUE, only the last full year in the data used when fitting the object will be used.
tailOnly: Logical. If TRUE the simulation will be only for tail events. More precisely only for each of the threshold $u_i(t)$ corresponding to a non-exceedance probability $\tau_i$ the maximum $v_i:=\max_t u_i(t)$ of the the thresholds over the predicted period is computed, and the simulated events correspond to exceedances over $v_i$.
tau: If provided, the simulation is only made for the provided value of tau.
trace: Integer level of verbosity.
how: A technical argument telling how the big data frame on output is to be constructed by concatenation. The "list" method (default and recommended) makes a list containing a large number of data frames, one by simulation and then makes use of the very efficient rbindlist function. With "vector", once concatenates vectors, and with "data.frame" one rbinds data fames at each simulation. The tow later options are by far slower when n is large, say n > 100.
...: Not used yet.

Value

An object with class "summary.pgpTList", inheriting from the "data.frame" class. This is a data frame in long format containing the simulated events with the corresponding simulated values of the meteorological variable.

Details

The thresholds and hence the exceedance events are related to a probability tau of non-exceedance which is quite high, say 0.95 or more. Note that even if the GP model for the excesses over the threshold is perfectly adequate, the simulated exceedance events come at a rate which is smaller than the real rate because the simulation is based on declustered exceedance events. So for a given date, if a large number n of simulations is used, the frequency of the exceedance will be smaller than the real one.

Note

In the returned data frame, the (factor) column Sim gives the simulation number. It can happen that a simulation generates no event, especially if tau is close to 1 and if the simulation period is short.

Caution

the returned data frame can be very large since n simulations are made for each value of tau, and each simulation can will typically give several dozens of events.

Examples

RqU <- rqTList(dailyMet = Rennes, tau = c(0.94, 0.95, 0.96, 0.97, 0.98))
Pgp1 <- pgpTList(dailyMet = Rennes, thresholds = RqU, declust = TRUE,
                 fitLambda = TRUE, logLambda.fun = ~YearNum - 1)
#> o Using meteorological variable : "TX"
#> o Adding new variables 
#> o Using K = 3 and the following phases
#> sinjPhi1 sinjPhi2 sinjPhi3 
#>   105.94     8.53    84.21 
#> o Sampling rate : 365.26/year
#> o Looping on 5 thresholds
#> 
#> o Fit the temporal Poisson process:  non-homogeneous
#> 
#> Number of observations  not used in the estimation process:  0
#> Total number of time observations:  28366
#> Number of events:  772
#> 
#> Convergence code:  0
#> Convergence attained
#> Loglikelihood:  -3537.078
#> 
#> Estimated coefficients: 
#>     b0     b1 
#> -3.757  0.009 
#> Full coefficients: 
#>     b0     b1 
#> -3.757  0.009 
#> attr(,"TypeCoeff")
#> [1] "Fixed: No  fixed parameters"
#> 
#> 
#> o Fit the temporal Poisson process:  non-homogeneous
#> 
#> Number of observations  not used in the estimation process:  0
#> Total number of time observations:  28366
#> Number of events:  671
#> 
#> Convergence code:  0
#> Convergence attained
#> Loglikelihood:  -3165.859
#> 
#> Estimated coefficients: 
#>     b0     b1 
#> -3.912  0.010 
#> Full coefficients: 
#>     b0     b1 
#> -3.912  0.010 
#> attr(,"TypeCoeff")
#> [1] "Fixed: No  fixed parameters"
#> 
#> 
#> o Fit the temporal Poisson process:  non-homogeneous
#> 
#> Number of observations  not used in the estimation process:  0
#> Total number of time observations:  28366
#> Number of events:  565
#> 
#> Convergence code:  0
#> Convergence attained
#> Loglikelihood:  -2762.467
#> 
#> Estimated coefficients: 
#>     b0     b1 
#> -4.087  0.010 
#> Full coefficients: 
#>     b0     b1 
#> -4.087  0.010 
#> attr(,"TypeCoeff")
#> [1] "Fixed: No  fixed parameters"
#> 
#> 
#> o Fit the temporal Poisson process:  non-homogeneous
#> 
#> Number of observations  not used in the estimation process:  0
#> Total number of time observations:  28366
#> Number of events:  450
#> 
#> Convergence code:  0
#> Convergence attained
#> Loglikelihood:  -2299.813
#> 
#> Estimated coefficients: 
#>     b0     b1 
#> -4.337  0.012 
#> Full coefficients: 
#>     b0     b1 
#> -4.337  0.012 
#> attr(,"TypeCoeff")
#> [1] "Fixed: No  fixed parameters"
#> 
#> 
#> o Fit the temporal Poisson process:  non-homogeneous
#> 
#> Number of observations  not used in the estimation process:  0
#> Total number of time observations:  28366
#> Number of events:  323
#> 
#> Convergence code:  0
#> Convergence attained
#> Loglikelihood:  -1749.224
#> 
#> Estimated coefficients: 
#>     b0     b1 
#> -4.754  0.016 
#> Full coefficients: 
#>     b0     b1 
#> -4.754  0.016 
#> attr(,"TypeCoeff")
#> [1] "Fixed: No  fixed parameters"
#> 
Date <- seq(from = as.Date("2022-01-01"), to = as.Date("2043-01-01"), by = "day")
st <- system.time(sim <- simulate(Pgp1, nsim = 8, newdata = Date, trace = 1))
#> tau=0.94 tau=0.95 tau=0.96 tau=0.97 tau=0.98 
#>     0.94     0.95     0.96     0.97     0.98 
#> tau =  0.94
#> Warning: collapsing to unique 'x' values
#> Warning: collapsing to unique 'x' values
#> Warning: collapsing to unique 'x' values
#> Warning: collapsing to unique 'x' values
#> Warning: collapsing to unique 'x' values
#> Warning: collapsing to unique 'x' values
#> Warning: collapsing to unique 'x' values
#> Warning: collapsing to unique 'x' values
#> tau =  0.95
#> Warning: collapsing to unique 'x' values
#> Warning: collapsing to unique 'x' values
#> Warning: collapsing to unique 'x' values
#> Warning: collapsing to unique 'x' values
#> Warning: collapsing to unique 'x' values
#> Warning: collapsing to unique 'x' values
#> Warning: collapsing to unique 'x' values
#> Warning: collapsing to unique 'x' values
#> tau =  0.96
#> Warning: collapsing to unique 'x' values
#> Warning: collapsing to unique 'x' values
#> Warning: collapsing to unique 'x' values
#> Warning: collapsing to unique 'x' values
#> Warning: collapsing to unique 'x' values
#> Warning: collapsing to unique 'x' values
#> Warning: collapsing to unique 'x' values
#> Warning: collapsing to unique 'x' values
#> tau =  0.97
#> Warning: collapsing to unique 'x' values
#> Warning: collapsing to unique 'x' values
#> Warning: collapsing to unique 'x' values
#> Warning: collapsing to unique 'x' values
#> Warning: collapsing to unique 'x' values
#> Warning: collapsing to unique 'x' values
#> Warning: collapsing to unique 'x' values
#> Warning: collapsing to unique 'x' values
#> tau =  0.98
#> Warning: collapsing to unique 'x' values
#> Warning: collapsing to unique 'x' values
#> Warning: collapsing to unique 'x' values
#> Warning: collapsing to unique 'x' values
#> Warning: collapsing to unique 'x' values
#> Warning: collapsing to unique 'x' values
#> Warning: collapsing to unique 'x' values
#> Warning: collapsing to unique 'x' values
autoplot(subset(sim, tau == 0.95))

## more realistic example
if (FALSE) {
    pred <- predict(Pgp1, newdata = Date)
    st <- system.time(sim <- simulate(Pgp1, nsim = 1000, newdata = Date, trace = 1))
    ## compute the maximum 'M' on the period by simulation and by threshold
    M <- with(sim, tapply(TX, list(Sim, tau), max))
    ## quantile over the simulation (by threshold)
    apply(M, 2, quantile, prob = 0.99, na.rm = TRUE)
    ## compare with the computation
    quantile(pred, p = 0.99)
}