| Title: | Estimation of Return Curves |
|---|---|
| Description: | Estimates the p-probability return curve proposed by Murphy-Barltrop et al. (2023) <doi:10.1002/env.2797>. Implements pointwise and smooth estimation of the angular dependence function introduced by Wadsworth and Tawn (2013) <doi:10.3150/12-BEJ471>. |
| Authors: | Lídia André [aut, cre], Callum Murphy-Barltrop [aut] |
| Maintainer: | Lídia André <[email protected]> |
| License: | GPL (>= 3) |
| Version: | 1.0.1 |
| Built: | 2025-03-07 06:22:00 UTC |
| Source: | https://github.com/lidiamandre/returncurves |
Implements the estimation of the \(p\)-probability return curve (Murphy-Barltrop et al. 2023), as well as a pointwise and smooth estimation of the angular dependence function (Wadsworth and Tawn 2013).
adf_est: Estimation of the Angular Dependence Function (ADF)
adf_gof: Goodness of fit of the Angular Dependence Function estimates
airdata: Air pollution data
marggpd: Assessing the Marginal Tail Fits
margtransf: Marginal Transformation
rc_est: Return Curve estimation
rc_gof: Goodness of fit of the Return Curve estimates
rc_unc: Uncertainty of the Return Curve estimates
runShiny: Complementary Shiny app for the ReturnCurves package
Maintainer: Lídia André [email protected]
Authors:
Callum Murphy-Barltrop [email protected]
Murphy-Barltrop CJR, Wadsworth JL, Eastoe EF (2023).
“New estimation methods for extremal bivariate return curves.”
Environmetrics, 34(5).
ISSN 1099095X, doi:10.1002/env.2797.
Wadsworth JL, Tawn JA (2013).
“A new representation for multivariate tail probabilities.”
Bernoulli, 19(5B), 2689-2714.
ISSN 13507265, doi:10.3150/12-BEJ471.
library(ReturnCurves) data(airdata) n <- dim(airdata)[1] # Marginal Transformation margdata <- margtransf(airdata) head(margdata@dataexp) # Return Curves estimation prob <- 1/n retcurve <- rc_est(margdata = margdata, p = prob, method = "hill") head(retcurve@rc) # ADF estimation lambda <- adf_est(margdata = margdata, method = "hill") head(lambda@adf)library(ReturnCurves) data(airdata) n <- dim(airdata)[1] # Marginal Transformation margdata <- margtransf(airdata) head(margdata@dataexp) # Return Curves estimation prob <- 1/n retcurve <- rc_est(margdata = margdata, p = prob, method = "hill") head(retcurve@rc) # ADF estimation lambda <- adf_est(margdata = margdata, method = "hill") head(lambda@adf)
Estimation of the angular dependence function \(\lambda(\omega)\) introduced by Wadsworth and Tawn (2013).
adf_est( margdata, w = NULL, method = c("hill", "cl"), q = 0.95, qalphas = rep(0.95, 2), k = 7, constrained = FALSE, tol = 1e-04, par_init = rep(0, k - 1) )adf_est( margdata, w = NULL, method = c("hill", "cl"), q = 0.95, qalphas = rep(0.95, 2), k = 7, constrained = FALSE, tol = 1e-04, par_init = rep(0, k - 1) )
margdata |
An S4 object of class |
w |
Sequence of rays between |
method |
String that indicates which method is used for the estimation of the angular dependence function. Must either be |
q |
Marginal quantile used to define the threshold \(u_\omega\) of the min-projection variable \(T^1\) at ray \(\omega\) \(\left(t^1_\omega = t_\omega - u_\omega | t_\omega > u_\omega\right)\), and/or Hill estimator (Hill 1975). Default is |
qalphas |
A vector containing the marginal quantile used for the Heffernan and Tawn conditional extremes model (Heffernan and Tawn 2004) for each variable, if |
k |
Polynomial degree for the Bernstein-Bezier polynomials used for the estimation of the angular dependence function with the composite likelihood method (Murphy-Barltrop et al. 2024). Default is |
constrained |
Logical. If |
tol |
Convergence tolerance for the composite maximum likelihood procedure. Success is declared when the difference of log-likelihood values between iterations does not exceed this value. Default is |
par_init |
Initial values for the parameters \(\beta\) of the Bernstein-Bezier polynomials used for estimation of the angular dependence function with the composite likelihood method (Murphy-Barltrop et al. 2024). Default is |
The angular dependence function \(\lambda(\omega)\) can be estimated through a pointwise estimator, obtained with the Hill estimator, or via a smoother approach, obtained using Bernstein-Bezier polynomials and estimated via composite likelihood methods.
Knowledge of the conditional extremes framework introduced by Heffernan and Tawn (2004) can be incorporated by setting constrained = TRUE.
Let \(\alpha^1_{x\mid y}=\alpha_{x\mid y} / (1+\alpha_{x\mid y})\) and \(\alpha^1_{y\mid x}=1 /(1+\alpha_{y\mid x})\) with \(\alpha_{x\mid y}\) and \(\alpha_{y\mid x}\)
being the conditional extremes parameters. After obtaining \(\hat{\alpha}_{x\mid y}\) and \(\hat{\alpha}_{y\mid x}\) via maximum likelihood estimation,
\(\lambda(\omega)=\max\lbrace \omega, 1-\omega\rbrace\) for \(\omega \in [0, \hat{\alpha}^1_{x\mid y})\cup (\hat{\alpha}^1_{y\mid x}, 1]\) and
is estimated as before for \(\omega \in [\hat{\alpha}^1_{x\mid y},\hat{\alpha}^1_{y\mid x}]\). For more details see Murphy-Barltrop et al. (2024).
An object of S4 class adf_est.class. This object returns the arguments of the function and two extra slots:
interval: |
A vector containing the maximum likelihood estimates from the conditional extremes model, \(\hat{\alpha}^1_{x\mid y}\) and \(\hat{\alpha}^1_{y\mid x}\), if |
adf: |
A vector containing the estimates of the angular dependence function. |
Heffernan JE, Tawn JA (2004).
“A conditional approach for multivariate extreme values (with discussion).”
Journal of the Royal Statistical Society: Series B (Statistical Methodology), 66(3), 497-546.
doi:10.1111/j.1467-9868.2004.02050.x, https://rss.onlinelibrary.wiley.com/doi/pdf/10.1111/j.1467-9868.2004.02050.x.
Hill BM (1975).
“A Simple General Approach to Inference About the Tail of a Distribution.”
The Annals of Statistics, 3(5), 1163 – 1174.
doi:10.1214/aos/1176343247.
Murphy-Barltrop CJR, Wadsworth JL, Eastoe EF (2024).
“Improving estimation for asymptotically independent bivariate extremes via global estimators for the angular dependence function.”
2303.13237.
Wadsworth JL, Tawn JA (2013).
“A new representation for multivariate tail probabilities.”
Bernoulli, 19(5B), 2689-2714.
ISSN 13507265, doi:10.3150/12-BEJ471.
library(ReturnCurves) data(airdata) n <- dim(airdata)[1] margdata <- margtransf(airdata) lambda <- adf_est(margdata = margdata, method = "hill") plot(lambda) # To see the the S4 object's slots str(lambda) # To access the estimates of the ADF lambda@adf # If constrained = T, the MLE estimates for the conditional extremes model # can be accessed as lambda@intervallibrary(ReturnCurves) data(airdata) n <- dim(airdata)[1] margdata <- margtransf(airdata) lambda <- adf_est(margdata = margdata, method = "hill") plot(lambda) # To see the the S4 object's slots str(lambda) # To access the estimates of the ADF lambda@adf # If constrained = T, the MLE estimates for the conditional extremes model # can be accessed as lambda@interval
Assessment of the goodness of fit of the angular dependence function estimates \(\lambda(\omega)\) following the procedure of Murphy-Barltrop et al. (2024).
adf_gof(adf, ray, blocksize = 1, nboot = 250, alpha = 0.05)adf_gof(adf, ray, blocksize = 1, nboot = 250, alpha = 0.05)
adf |
An S4 object of class |
ray |
Ray \(\omega\) to be considered on the goodness of fit assessment. |
blocksize |
Size of the blocks for the block bootstrap procedure. If |
nboot |
Number of bootstrap samples to be taken. Default is |
alpha |
Significance level to compute the \((1-\alpha)\)% tolerance intervals. Default is |
Define the min-projection variable as \(t^1_\omega = t_\omega - u_\omega | t_\omega > u_\omega\), then variable \(\lambda(\omega)T^1_\omega \sim Exp(1)\) as \(u_\omega \to \infty\) for all \(\omega \in [0,1]\).
Let \(F^{-1}_E\) denote the inverse of the cumulative distribution function of a standard exponential variable and \(T^1_{(i)}\) denote the \(i\)-th ordered increasing statistic, \(i = 1, \ldots, n\).
Function plot shows a QQ plot between the model and empirical exponential quantiles, i.e. points \(\left(F^{-1}_E\left(\frac{i}{n+1}\right), T^1_{(i)}\right)\),
along with the line \(y=x\). Uncertainty is obtained via a (block) bootstrap procedure and shown by the grey region on the plot.
A good fit is shown by agreement of model and empirical quantiles, i.e. points should lie close to the line \(y=x\).
In addition, line \(y = x\) should mainly lie within the \((1-\alpha)\)% tolerance intervals.
We note that, if the grid for \(\omega\) used to estimate the Angular Dependence Function (ADF) does not contain ray, then the closest \(\omega\) in the grid is used to assess the goodness-of-fit of the ADF.
An object of S4 class adf_gof.class. This object returns the arguments of the function and an extra slot gof which is a list containing:
model |
A vector containing the model quantiles. |
empirical |
A vector containing the empirical quantiles. |
lower |
A vector containing the lower bound of the tolerance interval. |
upper |
A vector containing the upper bound of the tolerance interval. |
It is recommended to assess the goodness-of-fit of \(\lambda(\omega)\) for a few values of \(\omega\).
Murphy-Barltrop CJR, Wadsworth JL, Eastoe EF (2024). “Improving estimation for asymptotically independent bivariate extremes via global estimators for the angular dependence function.” 2303.13237.
library(ReturnCurves) data(airdata) n <- dim(airdata)[1] margdata <- margtransf(airdata) lambda <- adf_est(margdata = margdata, method = "hill") # blocksize to account for temporal dependence gof <- adf_gof(adf = lambda, ray = 0.4, blocksize = 10) plot(gof) # To see the the S4 object's slots str(gof) # To access the list of vectors gof@goflibrary(ReturnCurves) data(airdata) n <- dim(airdata)[1] margdata <- margtransf(airdata) lambda <- adf_est(margdata = margdata, method = "hill") # blocksize to account for temporal dependence gof <- adf_gof(adf = lambda, ray = 0.4, blocksize = 10) plot(gof) # To see the the S4 object's slots str(gof) # To access the list of vectors gof@gof
Air pollution data from Marylebone, London (UK), between December and February
of 1998 to 2005. It contains daily measurements of air pollutant concentrations of NOx and
PM10. The dataset is a subset of the data from the openair package and
consists of 1427 observations.
Data frame with 1427 observations (rows) and 2 variables (columns):
Vector containing daily measurements of NOx, in parts per billion.
Vector containing daily measurements of PM10, in ug/m3.
airdata is a subset of the data provided in the openair R package.
Assessment of the marginal tail fits for each margin following the marginal transformation procedure margtransf.
marggpd(margdata, blocksize = 1, nboot = 250, alpha = 0.05)marggpd(margdata, blocksize = 1, nboot = 250, alpha = 0.05)
margdata |
An S4 object of class |
blocksize |
Size of the blocks for the block bootstrap procedure. If |
nboot |
Number of bootstrap samples to be taken. Default is |
alpha |
Significance level to compute the \((1-\alpha)\)% tolerance intervals. Default is |
Let \(X^{GPD}_{(i)}\) denote the \(i\)-th ordered increasing statistic
\((i = 1, \ldots, n)\) of the exceedances, i.e., \(X^{GPD}= (X-u \mid X >u),\)
\(n_{exc}\) denote the sample size of these exceedances, and \(F_{GPD}^{-1}\) denote the
inverse of the cumulative distribution function of a generalised Pareto distribution (GPD).
Function plot shows QQ plots between the model and empirical GPD quantiles for both variables, i.e, for
the first variable points \(\left(F^{-1}_{GPD}\left(\frac{i}{n_{exc}+1}\right) + u, X^{GPD}_{(i)} + u\right)\),
along with the line \(y=x\).
Uncertainty on the empirical quantiles is obtained via a (block) bootstrap procedure and shown by the grey region on the plot. A good fit is shown by agreement of model and empirical quantiles, i.e. points should lie close to the line \(y=x\). In addition, line \(y = x\) should mainly lie within the \((1-\alpha)\)% tolerance intervals.
An object of S4 class marggpd.class. This object returns the arguments of the function and an extra slot marggpd which is a list containing:
model |
A list containing the model quantiles for each variable. |
empirical |
A list containing the empirical quantiles for each variable. |
lower |
A list containing the lower bounds of the tolerance intervals for each variable. |
upper |
A list containing the upper bounds of the tolerance intervals for each variable. |
library(ReturnCurves) data(airdata) n <- dim(airdata)[1] margdata <- margtransf(airdata) # blocksize to account for temporal dependence marggpd <- marggpd(margdata = margdata, blocksize = 10) plot(marggpd) # To see the the S4 object's slots str(marggpd) # To access the list of lists marggpd@marggpdlibrary(ReturnCurves) data(airdata) n <- dim(airdata)[1] margdata <- margtransf(airdata) # blocksize to account for temporal dependence marggpd <- marggpd(margdata = margdata, blocksize = 10) plot(marggpd) # To see the the S4 object's slots str(marggpd) # To access the list of lists marggpd@marggpd
Marginal transformation of a bivariate random vector to standard exponential margins following Coles and Tawn (1991). Variables within each margin are assumed identically distributed.
margtransf(data, qmarg = rep(0.95, 2), constrainedshape = TRUE)margtransf(data, qmarg = rep(0.95, 2), constrainedshape = TRUE)
data |
A matrix containing the data on the original margins. |
qmarg |
A vector containing the marginal quantile used to fit the Generalised Pareto Distribution (GPD) for each variable. Default is |
constrainedshape |
Logical. If |
Given a threshold value \(u\), each stationary random vector is transformed by using the empirical cumulative distribution function (cdf) below \(u\), and a Generalise Pareto Distribution (GPD) fit above \(u\).
The option to constrain \(\xi > -1\) is included as \(\xi \leq -1\) implies that the fitted upper endpoint of the distribution's support is the maximum data point. This situation is rarely encountered in practice.
An object of S4 class margtransf.class. This object returns the arguments of the function, a slot parameters containing a matrix with the shape and scale parameters of the Generalised Pareto Distribution (GPD) for each variable, a slot thresh containing a vector with the threshold \(u\) above which the GPD is fitted, and a slot dataexp containing a matrix with the data on standard exponential margins.
The plot function takes an object of S4 class margtransf.class, and a which argument specifying the type of plot desired (see Examples):
"hist" |
Plots the marginal distributions of the two variables on original and standard exponential margins. |
"ts" |
Plots the time series of the two variables on original and standard exponential margins. |
"joint" |
Plots the joint distribution of the two variables on original and standard exponential margins. |
"all" |
Plots all the above mentioned plots (default). |
Coles SG, Tawn JA (1991). “Modelling Extreme Multivariate Events.” Journal of the Royal Statistical Society. Series B (Methodological), 53(2), 377–392. ISSN 00359246, doi:10.1111/j.2517-6161.1991.tb01830.x.
library(ReturnCurves) data(airdata) n <- dim(airdata)[1] margdata <- margtransf(airdata) # Plots the marginal distributions of X and Y on original vs standard exponential margins plot(margdata, which = "hist") # Plots the time series of X and Y on original vs standard exponential margins plot(margdata, which = "ts") # Plots the joint distribution of X and Y on original vs standard exponential margins plot(margdata, which = "joint") # Plots all the available plots plot(margdata, which = "all") # To see the the S4 object's slots str(margdata) # To access the matrix with the data on standard exponential margins margdata@dataexplibrary(ReturnCurves) data(airdata) n <- dim(airdata)[1] margdata <- margtransf(airdata) # Plots the marginal distributions of X and Y on original vs standard exponential margins plot(margdata, which = "hist") # Plots the time series of X and Y on original vs standard exponential margins plot(margdata, which = "ts") # Plots the joint distribution of X and Y on original vs standard exponential margins plot(margdata, which = "joint") # Plots all the available plots plot(margdata, which = "all") # To see the the S4 object's slots str(margdata) # To access the matrix with the data on standard exponential margins margdata@dataexp
Estimation of the \(p\)-probability return curve following Murphy-Barltrop et al. (2023).
rc_est( margdata, w = NULL, p, method = c("hill", "cl"), q = 0.95, qalphas = rep(0.95, 2), k = 7, constrained = FALSE, tol = 0.001, par_init = rep(0, k - 1) )rc_est( margdata, w = NULL, p, method = c("hill", "cl"), q = 0.95, qalphas = rep(0.95, 2), k = 7, constrained = FALSE, tol = 0.001, par_init = rep(0, k - 1) )
margdata |
An S4 object of class |
w |
Sequence of rays between |
p |
Curve survival probability. Must be \(p < 1-q\) and \(p < 1-q_\alpha\). |
method |
String that indicates which method is used for the estimation of the angular dependence function. Must either be |
q |
Marginal quantile used to define the threshold \(u_\omega\) of the min-projection variable \(T^1\) at ray \(\omega\) \(\left(t^1_\omega = t_\omega - u_\omega | t_\omega > u_\omega\right)\), and/or Hill estimator (Hill 1975). Default is |
qalphas |
A vector containing the marginal quantile used for the Heffernan and Tawn conditional extremes model (Heffernan and Tawn 2004) for each variable, if |
k |
Polynomial degree for the Bernstein-Bezier polynomials used for the estimation of the angular dependence function with the composite likelihood method (Murphy-Barltrop et al. 2024). Default is |
constrained |
Logical. If |
tol |
Convergence tolerance for the composite maximum likelihood procedure. Success is declared when the difference of log-likelihood values between iterations does not exceed this value. Default is |
par_init |
Initial values for the parameters \(\beta\) of the Bernstein-Bezier polynomials used for estimation of the angular dependence function with the composite likelihood method (Murphy-Barltrop et al. 2024). Default is |
Given a probability \(p\) and a joint survival function \(Pr(X>x, Y>y)\), the \(p\)-probability return curve is defined as \[RC(p):=\left\lbrace(x, y) \in R^2: Pr(X>x, Y>y)=p\right\rbrace.\]
This method focuses on estimation of \(RC(p)\) for small \(p\) near \(0\), so that \((X,Y)\) are in the tail of the distribution.
\(Pr(X>x, Y>y)\) is estimated using the angular dependence function \(\lambda(\omega)\) introduced by Wadsworth and Tawn (2013). More details on how to estimate \(\lambda(\omega)\) can be found in adf_est.
The return curve estimation \(\hat{RC}(p)\) is done on standard exponential margins and then back transformed onto the original margins.
An object of S4 class rc_est.class. This object returns the arguments of the function and extra slot rc
interval: |
A vector containing the maximum likelihood estimates from the conditional extremes model, \(\hat{\alpha}^1_{x\mid y}\) and \(\hat{\alpha}^1_{y\mid x}\), if |
rc: |
A matrix with the estimates of the Return Curve. |
Heffernan JE, Tawn JA (2004).
“A conditional approach for multivariate extreme values (with discussion).”
Journal of the Royal Statistical Society: Series B (Statistical Methodology), 66(3), 497-546.
doi:10.1111/j.1467-9868.2004.02050.x, https://rss.onlinelibrary.wiley.com/doi/pdf/10.1111/j.1467-9868.2004.02050.x.
Hill BM (1975).
“A Simple General Approach to Inference About the Tail of a Distribution.”
The Annals of Statistics, 3(5), 1163 – 1174.
doi:10.1214/aos/1176343247.
Murphy-Barltrop CJR, Wadsworth JL, Eastoe EF (2023).
“New estimation methods for extremal bivariate return curves.”
Environmetrics, 34(5).
ISSN 1099095X, doi:10.1002/env.2797.
Murphy-Barltrop CJR, Wadsworth JL, Eastoe EF (2024).
“Improving estimation for asymptotically independent bivariate extremes via global estimators for the angular dependence function.”
2303.13237.
Wadsworth JL, Tawn JA (2013).
“A new representation for multivariate tail probabilities.”
Bernoulli, 19(5B), 2689-2714.
ISSN 13507265, doi:10.3150/12-BEJ471.
library(ReturnCurves) data(airdata) n <- dim(airdata)[1] prob <- 10/n margdata <- margtransf(airdata) retcurve <- rc_est(margdata = margdata, p = prob, method = "hill") plot(retcurve) # To see the the S4 object's slots str(retcurve) # To access the return curve estimation retcurve@rc # If constrained = T, the MLE estimates for the conditional extremes model # can be accessed as retcurve@intervallibrary(ReturnCurves) data(airdata) n <- dim(airdata)[1] prob <- 10/n margdata <- margtransf(airdata) retcurve <- rc_est(margdata = margdata, p = prob, method = "hill") plot(retcurve) # To see the the S4 object's slots str(retcurve) # To access the return curve estimation retcurve@rc # If constrained = T, the MLE estimates for the conditional extremes model # can be accessed as retcurve@interval
Assessment of the goodness-of-fit of the return curve estimates following the approach of Murphy-Barltrop et al. (2023).
rc_gof(retcurve, blocksize = 1, nboot = 250, nangles = 150, alpha = 0.05)rc_gof(retcurve, blocksize = 1, nboot = 250, nangles = 150, alpha = 0.05)
retcurve |
An S4 object of class |
blocksize |
Size of the blocks for the block bootstrap procedure. If |
nboot |
Number of bootstrap samples to be taken. Default is |
nangles |
Number of angles \(m\) in the interval \((0, \pi/2)\) (Murphy-Barltrop et al. 2023). Default is |
alpha |
Significance level to compute the \((1-\alpha)\)% confidence intervals. Default is |
Given a return curve RC(\(p\)), the probability of lying in a survival region is \(p\). Let \[\boldsymbol{\Theta}:= \left\lbrace \frac{\pi(m+1-j)}{2(m+1)} \mid 1\leq j\leq m\right\rbrace\] be a set of angles decreasing from near \(\pi/2\) to \(0\). For each angle \(\theta_j\in \boldsymbol{\Theta,}\) and corresponding point in the estimated return curve \(\lbrace (\hat{x}_{\theta_j}, \hat{y}_{\theta_j}) \rbrace\), the empirical probability \(\hat{p}_j\) of lying in the survival region is given by the proportion of points in the region \((\hat{x}_{\theta_j}, \infty) \times (\hat{y}_{\theta_j}, \infty)\).
Thus, for each angle \(\theta_j\in \boldsymbol{\Theta,}\) a (block) bootstrap procedure to the original data set is applied, and
the empirical probabilities \(\hat{p}_j\) estimated. Then, the median and \((1-\alpha)\)% pointwise confidence intervals are obtained for each \(\theta_j\).
Function plot shows the median of \(\hat{p}_j\), the confidence intervals and the true probability \(p\); ideally, this value should be contained in the confidence region.
We note that due to the use of empirical probabilities, the value of \(p\) should be within the range of the data and not too extreme.
An object of S4 class rc_gof.class. This object returns the arguments of the function and an extra slot gof which is a list containing:
median |
A vector containing the median of the empirical probability of lying in a survival region. |
lower |
A vector containing the lower bound of the confidence interval. |
upper |
A vector containing the upper bound of the confidence interval. |
Murphy-Barltrop CJR, Wadsworth JL, Eastoe EF (2023). “New estimation methods for extremal bivariate return curves.” Environmetrics, 34(5). ISSN 1099095X, doi:10.1002/env.2797.
library(ReturnCurves) data(airdata) n <- dim(airdata)[1] prob <- 10/n margdata <- margtransf(airdata) rc_orig <- rc_est(margdata = margdata, p = prob, method = "hill") # blocksize to account for temporal dependence gof <- rc_gof(retcurve = rc_orig, blocksize = 10) plot(gof) # To see the the S4 object's slots str(gof) # To access the list of vectors gof@goflibrary(ReturnCurves) data(airdata) n <- dim(airdata)[1] prob <- 10/n margdata <- margtransf(airdata) rc_orig <- rc_est(margdata = margdata, p = prob, method = "hill") # blocksize to account for temporal dependence gof <- rc_gof(retcurve = rc_orig, blocksize = 10) plot(gof) # To see the the S4 object's slots str(gof) # To access the list of vectors gof@gof
Uncertainty assessment of the return curve estimates following the procedure of Murphy-Barltrop et al. (2023).
rc_unc(retcurve, blocksize = 1, nboot = 250, nangles = 150, alpha = 0.05)rc_unc(retcurve, blocksize = 1, nboot = 250, nangles = 150, alpha = 0.05)
retcurve |
An S4 object of class |
blocksize |
Size of the blocks for the block bootstrap procedure. If |
nboot |
Number of bootstrap samples to be taken. Default is |
nangles |
Number of angles \(m\) in the interval \((0, \pi/2)\) (Murphy-Barltrop et al. 2023). Default is |
alpha |
Significance level to compute the \((1-\alpha)\)% confidence intervals. Default is |
Define a set of angles \[\boldsymbol{\Theta}:= \left\lbrace \frac{\pi(m+1-j)}{2(m+1)} \mid 1\leq j\leq m\right\rbrace\] decreasing from near \(\pi/2\) to \(0\), and let \(L_\theta:=\left\lbrace(x,y)\in R^2_+ | \tan(\theta)=y/x\right\rbrace\) denote the line segment intersecting the origin with gradient \(\tan(\theta) > 0.\) For each \(\theta\in \boldsymbol{\Theta},\) \(L_\theta\) intersects the estimated \(\hat{RC}(p)\) exactly once, i.e. \(\lbrace(\hat{x}_\theta, \hat{y}_\theta)\rbrace:= \hat{RC}(p)\cap L_\theta.\) Uncertainty of the return curve is then quantified by the distribution of \(\hat{d}_\theta:=(\hat{x}^2_\theta + \hat{y}^2_\theta)^{1/2}\) via a (block) bootstrap procedure.
This procedure is as follows; for \(k = 1, \ldots, \) nboot:
1. (Block) bootstrap the original data set;
2. For each \(\theta\in \boldsymbol{\Theta},\) obtain \(\hat{d}_{\theta,k}\) for the corresponding return curve point estimate.
Full details can be found in Murphy-Barltrop et al. (2023)
An object of S4 class rc_unc.class. This object returns the arguments of the function and an extra slot unc which is a list containing:
median |
A vector containing the median estimates of the return curve. |
mean |
A vector containing the mean estimates of the return curve. |
lower |
A vector containing the lower bound of the confidence interval. |
upper |
A vector containing the upper bound of the confidence interval. |
The plot function takes an object of S4 class rc_unc.class, and a which argument specifying the type of plot desired (see Examples):
"rc" |
Plots the estimated Return Curve and its uncertainty (default). |
"median" |
Plots the median estimates of the Return Curve and its uncertainty. |
"mean" |
Plots the mean estimates of the Return Curve and its uncertainty. |
"all" |
Plots the estimated Return Curve, the median and mean estimates of the Return Curve together, and the associated uncertainty. |
Murphy-Barltrop CJR, Wadsworth JL, Eastoe EF (2023). “New estimation methods for extremal bivariate return curves.” Environmetrics, 34(5). ISSN 1099095X, doi:10.1002/env.2797.
library(ReturnCurves) data(airdata) n <- dim(airdata)[1] prob <- 10/n margdata <- margtransf(airdata) rc_orig <- rc_est(margdata = margdata, p = prob, method = "hill") # Set nboot = 50 for an illustrative example # blocksize to account for temporal dependence unc <- rc_unc(rc_orig, blocksize = 10) # Plots the estimated Return Curve plot(unc, which = "rc") # Plots the median estimates of the Return Curve plot(unc, which = "median") # Plots the mean estimates of the Return Curve plot(unc, which = "mean") # Plots the estimated Return Curve and its the median and mean estimates plot(unc, which = "all") # To see the the S4 object's slots str(unc) # To access the list of vectors unc@unclibrary(ReturnCurves) data(airdata) n <- dim(airdata)[1] prob <- 10/n margdata <- margtransf(airdata) rc_orig <- rc_est(margdata = margdata, p = prob, method = "hill") # Set nboot = 50 for an illustrative example # blocksize to account for temporal dependence unc <- rc_unc(rc_orig, blocksize = 10) # Plots the estimated Return Curve plot(unc, which = "rc") # Plots the median estimates of the Return Curve plot(unc, which = "median") # Plots the mean estimates of the Return Curve plot(unc, which = "mean") # Plots the estimated Return Curve and its the median and mean estimates plot(unc, which = "all") # To see the the S4 object's slots str(unc) # To access the list of vectors unc@unc
Launches the R Shiny app complementary to the ReturnCurves-package.
runShiny()runShiny()