Non-stationary Copula Objects

Description

The purpose of this function is creating a NScopula object. This object contains information of the copula: "family", "type", "gamma", and "name". This object is needed to use some of the functions in the package.

Usage

NScopula(
  family = c("GUM", "GAL", "HUS", "CLA", "JOE"),
  type = c("0", "1", "2"),
  gamma,
  name = "Copula"
)

Arguments

Details

The purpose of this function is creating a NScopula object. This object contains information of the copula; "family" refers to the copula family, Gumbel, Galambos, Hüsler-Reiss, Clayton, or Joe; "type" refers to the relationship between the copula parameter and the covariate, stationary, linear, or quadratic; "gamma" is a vector containing coefficients of the parameter (\gamma_0, \gamma_1, \gamma_2). To ensure that \gamma is in an adequate range for each copula family, different transformations must be applied to obtain it. Let \nu be the covariate of the data series. The gamma parameter for Gumbel and Joe copula families is given by

\gamma(\nu) = 1 + exp(\gamma_0 + \gamma_1 \nu + \gamma_2 \nu^2),

for the Clayton copula

\gamma(\nu) = exp(\gamma_0 + \gamma_1 \nu + \gamma_2 \nu^2) - 1,

and for Hüsler-Reiss and Galambos

\gamma(\nu) = exp(\gamma_0 + \gamma_1 \nu + \gamma_2 \nu^2)

.

For more information on the estimation see copulaNSestimation.

Value

A NScopula object

Examples

# A Gumbel copula, with a quadratic non-stationary parameter
copula <- NScopula(family = "GUM", type = "2", gamma = c(1, 0.1, 0.03))

Non-stationary Marginal Objects

Description

The purpose of this function is creating a NSmarginal object. This object contains information of the marginal distribution: "family", "type", "mu", "sig", "k", and"name". This object is needed to use other functions of the package.

Usage

NSmarginal(
  family = c("GEV", "LN2", "LN3"),
  type = c("00", "10", "20", "11", "21"),
  mu,
  sig,
  k,
  name = "Marginal"
)

Arguments

Details

The purpose of this function is creating a NSmarginal object. This object contains information of the marginal distribution; "family" refers to the distribution, Generalized Extreme Value, Log-normal or three-parameter Log normal; "type" refers to the relationship between the distribution parameters and the covariate, stationary, linear, or quadratic; "mu" is a vector containing coefficients of the location parameter \mu = \mu_0 + \mu_1\nu + \mu_2\nu^2, where \nu is the covariate; "sig" is vector containing the coefficients for the scale parameter \sigma = \sigma_0 + \sigma_1\nu; "k" is a vector with the coefficients for the shape parameter in the GEV distribution and threshold in the LN3 distributions k = k_0 + k_1\nu; lastly the optional "name" is used to name the marginal distribution. Based on the non zero inputs of "mu", "sig", and "k" parameters the "type" of marginal might be overwritten.

For more information on the estimation see marginalNSestimation.

Value

A NSmarginal object

References

• El Adlouni, S., Ouarda, T. B., Zhang, X., Roy, R., & Bobée, B. (2007). Generalized maximum likelihood estimators for the nonstationary generalized extreme value model. Water Resources Research, 43(3).

• Villarini, G., Smith, J.A., Serinaldi, F., Bales, J., Bates, P.D., Krajewski, W.F., 2009. Flood frequency analysis for nonstationary annual peak records in an urban drainage basin. Advances in Water Resources, 32(8): 1255-1266.

Examples

# A GEV marginal, with quadratic location parameter
marginal <- NSmarginal(family = "GEV", type = "20",
                        mu = c(0.5, 0.3, 0.001),
                        sig = 1, k = 0.5)

Synthetic Data

Description

Synthetic data modeled using copulas with the purpose of exampling the functions of the package, 50 two dimensional observations that follow a Clayton copula with a linear non stationary parameter.

Usage

SampData

Format

SampData

A data frame with 50 rows and 2 columns:

X

Numeric, first variable

Y

Numeric, second variable

Source

Goutali and Chebana (2024) doi:10.1016/j.envsoft.2024.106090

Clayton Copula Fitting

Description

Fitting of a Clayton class copula using Pseudo-maximum Likelihood, the parameters are estimated as stationary, linear, or quadratic, relating to the covariate (\nu)

\gamma = exp(\gamma_0 + \gamma_1 \nu + \gamma_2 \nu^2) - 1

.

It returns a data frame with the coefficients of each estimated model, and its corresponding AIC (or AICc), BIC and negative logarithmic pseudo-likelihood. It also returns a 'NScopula' object with filled in attributes according to the best model, selected on minimizing AIC or AICc.

Usage

claytonEstimation(data, covar = NULL)

Arguments

Details

The bivariate Clayton copula function is given by:

C_{\gamma}(u, v)=( \text{max}\{ u^{-\gamma} + v^{-\gamma} - 1; 0\})^{-1/\gamma}

where u and v are the pseudo-observations of the data series and \gamma \in [-1, \infty) \setminus \{0\} is the copula parameter.

To measure the fit of the models two measures are used the AIC and BIC. In case the series contains less than 25 observations the small-sample corrected AIC (AICc) is to be used.

Note: During optimization, numerical warnings may occur when the likelihood is evaluated at invalid parameter values. These warnings do not affect the final parameter estimates.

Value

• A 'NScopula' object containing the parameters of the best fitted model. NAs are used when a model doesn't take the coefficient.

• Data frame containing the Negative Log-likelihood, AIC, BIC, and the estimated parameter coefficients for each model.

References

Clayton, David G. (1978). "A model for association in bivariate life tables and its application in epidemiological studies of familial tendency in chronic disease incidence". Biometrika. 65 (1): 141–151.

Examples

# Load sample data

data("SampData")

# Estimate parameters for Clayton copula

temp <- claytonEstimation(SampData)

temp$all_fits

Copula Parameter Estimation

Description

This function uses Maximum Likelihood Estimation to fit to a given data series, the parameters, stationary or parametric, of the following copulas. Gumbel, Galambos, Hüsler-Reiss, Joe, and Clayton. It returns a 'NScopula' object with filled in attributes according to the best model selected on minimizing AIC (or AICc).

Also returns a table of fits with the negative log-likelihood, AIC (or small-sample corrected AIC for series with less than 25 observations), BIC, and estimated parameters for each distribution. If declared, saves to an Excel workspace said table.

Usage

copulaNSestimation(data, covar = NULL, file = NULL)

Arguments

Details

Chebana and Ouarda (2021) consider three copula classes due to their applications in the field of hydrology; logistic Gumbel, Galambos, and Hüsler-Reiss copulas. In addition the Clayton and Joe copulas are also implemented as alternatives.

All copulas depend depend on the \gamma parameter, \gamma should be bigger than 1 for the Gumbel and Joe copulas, for the Clayton copula bigger than -1, and bigger than 0 for the rest. To ensure this holds the following transformations are applied when fitting the models, \phi_{\gamma}(\nu) = ln(\gamma(\nu) + 1) (Gumbel and Joe), \phi_{\gamma}(\nu) = ln(\gamma(\nu) - 1) (Clayton), and \phi_{\gamma}(\nu) = ln(\gamma(\nu) ) (Galambos and Hüsler-Reiss). This results in the following expressions of the parameter that allow for variations in stationariety:

\gamma(\nu) = 1 + exp(\gamma_0 + \gamma_1 \nu + \gamma_2 \nu^2),

\gamma(\nu) = exp(\gamma_0 + \gamma_1 \nu + \gamma_2 \nu^2) - 1, \quad \text{and}

\gamma(\nu) = exp(\gamma_0 + \gamma_1 \nu + \gamma_2 \nu^2)

where \gamma_n are the coefficients to be estimated and \nu the covariate.

To estimate the coefficient the Maximum Pseudo-likelihood (MLP) is preferred. It is based on the pseudo-observations u_{ij} rather than the original.

L_n(\gamma) = \prod_{j=1}^{n}c_{\gamma}(u_{1j}, u_{2j}; \gamma(\nu_j))

where c_{\gamma} represents the density function corresponding to the copula C_{\gamma}.

To assess the fit of the models two metrics are to be used, the Akaike Information Criterion and the Bayesian Information Criterion. The one with the lowest AIC is selected as the best.

Note: During optimization, numerical warnings may occur when the likelihood is evaluated at invalid parameter values. These warnings do not affect the final parameter estimates.

Value

• Data frame containing the Negative Log-likelihood, AIC, BIC, and the estimated parameter coefficients for each model (All 'NAs' in the data frame indicate that the model doesn't take said parameter).

• 'NScopula' object containing the parameters of the best fitted model.

• If declared, '.xlsx' file with the evaluation results.

References

Fateh Chebana, Taha B.M.J. Ouarda, 2021, Multivariate non-stationary hydrological frequency analysis, Journal of Hydrology, Volume 593, doi:10.1016/j.jhydrol.2020.125907.

See Also

• gumbelEstimation, galambosEstimation, huslerReissEstimation, claytonEstimation

• copula: a comprehensive package that allows for copula estimation, analysis, and distribution.

Examples

# Load sample data
data("SampData")

# Estimate parameters for different copula classes
temp <- copulaNSestimation(SampData, covar = NULL, file = NULL)

temp$all_fits

# Get 'NScopula' object with the best fit
copula <- temp$best

Galambos Copula Fitting

Description

Fitting of a Galambos class copula using Pseudo-maximum Likelihood, the parameters are estimated as stationary, linear, or quadratic, relating to the covariate (\nu)

\gamma = exp(\gamma_0 + \gamma_1 \nu + \gamma_2 \nu^2)

.

It returns a data frame with the coefficients of each estimated model, and its corresponding AIC (or AICc), BIC and negative log-pseudo-likelihood. It also returns a 'NScopula' object with filled in attributes according to the best model, selected on minimizing AIC or AICc.

Usage

galambosEstimation(data, covar = NULL)

Arguments

Details

The bivariate Galambos copula function is given by:

C_\gamma(u, v) = u \ v \ exp \left( [(-log u)^{-\gamma} + (-log v)^{-\gamma}]^{-1/\gamma} \right)

where u and v are the pseudo-observations of the data series and \gamma \geq 0 is the copula parameter.

To measure the fit of the models two measures are used the AIC and BIC. In case the series contains less than 25 observations the small-sample corrected AIC (AICc) is to be used.

Note: During optimization, numerical warnings may occur when the likelihood is evaluated at invalid parameter values. These warnings do not affect the final parameter estimates.

Value

• A 'NScopula' object containing the parameters of the best fitted model. NAs are used when a model doesn't take the coefficient.

• Data frame containing the Negative Log-likelihood, AIC, BIC, and the estimated parameter coefficients for each model.

References

Galambos, J., 1975. Order statistics of samples from multivariate distributions. J. Am. Stat. Assoc. 70 (351), 674–680.

Examples

# Load sample data

data("SampData")

# Estimate parameters for Galambos copula

temp <- galambosEstimation(SampData)

temp$all_fits

Gumbel Copula Fitting

Description

Fitting of a Gumbel class copula using Pseudo-maximum Likelihood, the parameters are estimated as stationary, linear, or quadratic, relating to the covariate (\nu)

\gamma = 1 + exp(\gamma_0 + \gamma_1 \nu + \gamma_2 \nu^2)

.

It returns a data frame with the coefficients of each estimated model, and its corresponding AIC (or AICc), BIC and negative log-pseudo-likelihood. It also returns a 'NScopula' object with filled in attributes according to the best model, selected on minimizing AIC or AICc.

Usage

gumbelEstimation(data, covar = NULL)

Arguments

Details

The bivariate Gumbel copula function is given by:

C_\gamma(u,v) = exp \biggl( -[(-log u)^\gamma + (-log v)^\gamma]^{1/\gamma}\biggr)

where u and v are the pseudo-observations of the data series and \gamma \geq 1 is the copula parameter.

To measure the fit of the models two measures are used, the AIC and BIC. In case the series contains less than 25 observations the small-sample corrected AIC (AICc) is to be used.

Note: During optimization, numerical warnings may occur when the likelihood is evaluated at invalid parameter values. These warnings do not affect the final parameter estimates.

Value

• A 'NScopula' object containing the parameters of the best fitted model. NAs are used when a model doesn't take the coefficient.

• Data frame containing the Negative Log-likelihood, AIC, BIC, and the estimated parameter coefficients for each model.

Examples

# Load sample data

data("SampData")

# Estimate parameters for Gumbel copula

temp <- gumbelEstimation(SampData)

temp$all_fits

Hüsler-Reiss Copula Fitting

Description

Fitting of a Joe class copula using Pseudo-maximum Likelihood, the parameters are estimated as stationary, linear, or quadratic, relating to the covariate (\nu)

\gamma = 1 + exp(\gamma_0 + \gamma_1 \nu + \gamma_2 \nu^2)

.

It returns a data frame with the coefficients of each estimated model, and its corresponding AIC (or AICc), BIC and negative log-pseudo-likelihood. It also returns a 'NScopula' object with filled in attributes according to the best model, selected on minimizing AIC or AICc.

Usage

huslerReissEstimation(data, covar = NULL)

Arguments

Details

The Hüsler-Reiss copula function is given by:

C_\gamma(u,v) = exp \Biggl[(\text{log} u) \Phi \Bigl(\frac{1}{\gamma} + \frac{\gamma}{2} \text{log}(\frac{\text{log} u}{\text{log} v})\Bigr) + (\text{log} v) \Phi \Bigl(\frac{1}{\gamma} + \frac{\gamma}{2} \text{log}(\frac{\text{log} v}{\text{log} u})\Bigr) \Biggr]

where \Phi is the cumulative function of the standard normal distribution, u and v the pseudo-observations of the data series, and \gamma \geq 0 the copula parameter.

To measure the fit of the models two measures are used the AIC and BIC. In case the series contains less than 25 observations the small-sample corrected AIC (AICc) is to be used.

Note: During optimization, numerical warnings may occur when the likelihood is evaluated at invalid parameter values. These warnings do not affect the final parameter estimates.

Value

• A 'NScopula' object containing the parameters of the best fitted model. NAs are used when a model doesn't take the coefficient.

• Data frame containing the Negative Log-likelihood, AIC, BIC, and the estimated parameter coefficients for each model.

References

Hüsler, J., Reiss, R.-D., 1989. Maxima of normal random vectors: between independence and complete dependence. Statistics Probability Lett. 7 (4), 283–286.

Examples

# Load sample data

data("SampData")

# Estimate parameters for Husler-Reiss copula

temp <- huslerReissEstimation(SampData)

temp$all_fits

Joe Copula Fitting

Description

Fitting of a Joe class copula using Pseudo-maximum Likelihood, the parameters are estimated as stationary, linear, or quadratic, relating to the covariate (\nu)

\gamma = 1 + exp(\gamma_0 + \gamma_1 \nu + \gamma_2 \nu^2)

.

It returns a data frame with the coefficients of each estimated model, and its corresponding AIC (or AICc), BIC and negative log-pseudo-likelihood. It also returns a 'NScopula' object with filled in attributes according to the best model, selected on minimizing AIC or AICc.

Usage

joeEstimation(data, covar = NULL)

Arguments

Details

The Joe copula function is given by:

C_\gamma(u, v) = 1 - \bigl[(1-u)^\gamma + (1-v)^\gamma - (1-u)^\gamma(1-v)^\gamma\bigl]^{1/\gamma}

where u and v are the data series pseudo-observations, and \gamma \geq 1 is the copula parameter.

To measure the fit of the models two measures are used the AIC and BIC. In case the series contains less than 25 observations the small-sample corrected AIC (AICc) is to be used.

Note: During optimization, numerical warnings may occur when the likelihood is evaluated at invalid parameter values. These warnings do not affect the final parameter estimates.

Value

• A 'NScopula' object containing the parameters of the best fitted model. NAs are used when a model doesn't take the coefficient.

• Data frame containing the Negative Log-likelihood, AIC, BIC, and the estimated parameter coefficients for each model.

Examples

# Load sample data

data("SampData")

# Estimate parameters for Joe copula

temp <- joeEstimation(SampData)

temp$all_fits

Marginal Distribution Estimation

Description

This function uses Maximum Likelihood Estimation to fit to a given data series, the parameters, stationary or parametric, of the following marginal distributions: Generalized Extreme Value, Log-normal, and three parameter Log-normal. It prints or saves to an Excel workspace a table of fits with the negative log-likelihood, AIC (or small-sample corrected AIC for series with less than 25 points), BIC, and estimated parameters for each distribution and returns a 'marginal' object with filled in attributes according to the best model selected on minimizing AIC (or AICc).

Usage

marginalNSestimation(data, covar = NULL, file = NULL)

Arguments

Details

The function fits a distribution and parameters, stationary or not, to the provided data. Due to their applications in hydrology, three distributions are implemented Chebana and Ouarda (2021).

Generalized Extreme Value:

The cumulative distribution funciton of the GEV distribution is given by:

F_{GEV}(x; \mu, \sigma, k) = exp \left[ -(1 + k(\frac{x-\mu}{\sigma}))^{\frac{-1}{k}} \right] \quad \text{if } k \neq 0

\text{and} \quad exp \left[ -exp(\frac{x-\mu}{\sigma})\right] \quad \text{if } k = 0

where \mu \in \mathbb{R},\ \sigma > 0 and k \in \mathbb{R} are respectively the location, scale, and shape parameters.

Log-normal distribution:

The density function for the two parameter LN2 distribution is given by:

f_{LN2}(x; \mu, \sigma) = \frac{1}{x \sigma \sqrt{2 \pi}} exp(\frac{-(\text{log}x-\mu)^2}{2\sigma^2}) \quad x>0

where \mu \in \mathbb{R} and \sigma > 0 are respectively the mean and standard deviation of \text{log}x.

Three parameter Log-normal:

The density function for the three parameter LN3 distribution is given by:

f_{LN3}(x;\mu, \sigma, k) = f_{LN2}(x-k;\mu\sigma) \quad x>k

where \mu and \sigma are as in LN2 and k \in \mathbb{R} is the threshold parameter or lower bound.

To incorporate the possible non-stationarity Chebana and Ouarda (2021) propose to link the distribution parameters to the covariate in the shape:

\alpha(\nu) = \alpha_0 + \alpha_1\nu + \alpha_2\nu^2

where \alpha is the parameter where trend is inducted and \nu is is the covariate. Various models; stationary, with linear or quadratic trend, are to be fitted.

To estimate the coefficients of any parameter the maximum likelihood approach is used. To assess the fit of the model the Akaike Information Criterion and Bayesia Information Criterion are computed for each model. The one with the lowest AIC will be selected as best.

Note: During optimization, numerical warnings may occur when the likelihood is evaluated at invalid parameter values. These warnings do not affect the final parameter estimates.

Value

• Data frame containing the Negative Log-likelihood, AIC, BIC and the estimated parameter coefficients for each model (All 'NAs' in the data frame indicate that the model doesn't take said parameter).

• 'NSmarginal' object containing the parameters of the best fitted model.

• If declared, '.xlsx' file with the evaluation results.

References

• Fateh Chebana, Taha B.M.J. Ouarda, 2021, Multivariate non-stationary hydrological frequency analysis, Journal of Hydrology, Volume 593, doi:10.1016/j.jhydrol.2020.125907.

• Aitchison, J., Brown, J.A.C., 1957. The Lognormal Distribution. Cambridge University Press, London, p. 176.

Examples

# Load sample data
data("SampData")
X <- SampData[, 1]
Y <- SampData[, 2]

# Estimate parameters for different marginal distributions
temp <- marginalNSestimation(X, covar = NULL)

temp$all_fits

# Get 'NSmarginal' object with the best fit
marginal <- temp$best

Multivariate Quantile Curve

Description

This file contains the functions for computing the bivariate quantile curve for a given model. It defines a grid on [0, 1] \times [0, 1] and calls contourLines to find the values that yield a quantile at level p.

Usage

mqc(
  seriesX,
  seriesY,
  copula,
  covariate = 1,
  p = 0.95,
  event = c("E.AND.E", "E.OR.E", "nE.AND.nE", "nE.OR.nE"),
  step = 0.01,
  file = NULL,
  c_name = "Time"
)

Arguments

Details

This function calculates the p-level bivariate quantile curve using the MQC approach. It first sets up a grid on [0, 1] \times [0, 1], calculates the (u,v) probabilities using the copula, computes the proper event based on the chosen event type, and then obtains the contour lines to find the (u,v) points corresponding to the desired quantile level p. Finally, the inverse marginal distributions are applied to convert (u,v) points to (X,Y) coordinates in the original data space.

Value

A 'MQC' object containing the X and Y coordinates of the p-level bivariate quantile curve and all the input information.

Examples

# Load sample data
data(SampData)

X <- SampData[, 1]
Y <- SampData[, 2]

# Fit marginal distributions
temp <- marginalNSestimation(X, covar = NULL, file = NULL)
marx <- temp$best

temp <- marginalNSestimation(Y, covar = NULL, file = NULL)
mary <- temp$best

# Fit copula model
temp <- copulaNSestimation(SampData, covar = NULL)
cop <- temp$best

# Compute the multivariate quantile curves

MQC <- mqc(marx, mary, cop,
           covariate = seq(1:50), p = 0.05,
           event = "E.AND.E",
           step = 0.01,
           file = NULL,
           c_name = "Time")

Non-stationary Quantile Curve Plot

Description

This file contains the functions for plotting the bivariate quantile curves for a given model. It needs the results from mqc, and the data used to estimate the marginals and copula.

Usage

quantilePlotNS(
  MQC,
  data = NULL,
  ylim = NULL,
  xlim = NULL,
  color = c("cool", "warm", "dark", "blue"),
  file = NULL
)

Arguments

Value

The quantile curves plot, if declared a png file

Examples

# Load sample data
data(SampData)

X <- SampData[, 1]
Y <- SampData[, 2]

# Fit marginal distributions
temp <- marginalNSestimation(X, covar = NULL, file = NULL)
marx <- temp$best

temp <- marginalNSestimation(Y, covar = NULL, file = NULL)
mary <- temp$best

# Fit copula model
temp <- copulaNSestimation(SampData, covar = NULL)
cop <- temp$best

# Compute the multivariate quantile curves

MQC <- mqc(marx, mary, cop,
           covariate = seq(1:50), p = 0.05,
           event = "E.AND.E",
           step = 0.01,
           file = NULL,
           c_name = "Time")

# Plot the quantile curve
quantilePlotNS(MQC, data = SampData)