The Non-Stationary Return Level (NSRL) of Parey et al is implemented
for Time-Varying GEV models in the predictUncond
function
of the NSGEV package. The package has a function named
check_predictUncond
that can be used to check the results
of NSGEV::predictUncond
. For now, the check does not check
the computation of the confidence limits, but only the value of the RL
as given by the Quant
column in the data frame returned by
predictUncond
. Since the NSRL is given by an equation with
no explicit solution, it is worth checking that the returned value is
consistent with the definition of the NSRL.
The check relies on the simulate
method for the
"TVGEV"
class. Firstly, an unconditional prediction is
obtained for each period in a set of periods of interest. All these
periods share the same beginning date, say \(t_0\). The period with number \(j\) contains \(B_j\) blocks where \(B_j\) is the maximal horizon chosen via
period
(or set by default if period
is not
given). Remind that predictUncond
generally computes the
results for several periods with the same origin \(t_0\) as required here. Secondly, a large
number, say \(m\), of paths \(Y^{[i]}(t)\) is generated from the model.
For each predicted period \(j\) and
each path \(i\), we compute the number
\(N^{[i]}(u_j)\) of exceedances over
the quantile \(u_j\) corresponding to
the period. We eventually take the average of the \(m\) values \(N^{[i]}(u_j)\) for \(i= 1\) to \(m\) and compare it to the expectation. If
\(m\) is large, then the average should
be close to the expectation.
The example run with example(TVGEV)
defines a
TVGEV
model for the annual maxima of daily maximum
temperature (TX) at Dijon (Fr) in an object named res2
.
Consider as period of interest the period starting at year 2020. The following code computes the NSRL for several periods of time.
pu <- predictUncond(res2, newdateFrom = "2020-01-01", confintMethod = "proflik")
kable(pu, digits = 2)
Date | Period | Level | Quant | L | U |
---|---|---|---|---|---|
2020-01-01 | 5 | 95% | 37.15 | 36.15 | 38.22 |
2020-01-01 | 10 | 95% | 38.23 | 37.16 | 39.49 |
2020-01-01 | 20 | 95% | 39.31 | 38.07 | 40.99 |
2020-01-01 | 30 | 95% | 40.04 | 38.48 | 42.00 |
2020-01-01 | 40 | 95% | 40.65 | 38.76 | 42.87 |
2020-01-01 | 50 | 95% | 41.22 | 39.07 | 43.69 |
2020-01-01 | 70 | 95% | 42.30 | 39.63 | 45.35 |
2020-01-01 | 100 | 95% | 43.90 | 40.30 | 47.85 |
plot(pu)
## Warning in plot.predict.TVGEV(pu): Inasmuch the plot is actually a 'ggplot', it
## is better to use the 'autoplot' method for consistency
For instance, for the period of \(B_j =
40\) years beginning at the date "2020-01-01"
, the
level \(u_j =40.65\) Celsius is
expected to be exceeded once (in average) across a large number of
paths.
set.seed(12345)
check_predictUncond(res2, newdateFrom = "2020-01-01")
## 5 10 20 30 40 50 70 100
## 1.0175000 1.0113333 0.9903333 0.9911667 1.0005000 0.9916667 1.0031667 1.0260000
We use the default number of paths \(m =
5000\), but only the 100
first paths are displayed
on the plot. The results are in good accordance with the definition.