Package 'qrmtools'

Computing allocations

Description

Computing (capital) allocations.

Usage

## For elliptical distributions under certain assumptions
alloc_ellip(total, loc, scale)

## Nonparametrically
conditioning(x, level, risk.measure = "VaR_np", ...)
alloc_np(x, level, risk.measure = "VaR_np", include.conditional = FALSE, ...)

Arguments

total	total to be allocated (typically the risk measure of the sum of the underlying loss random variables).
loc	location vector of the elliptical distribution of the loss random vector.
scale	scale (covariance) matrix of the elliptical distribution of the loss random vector.
x	$(n, d)$-matrix containing $n$ iid $d$-dimensional losses.
level	either one or two confidence level(s) for risk.measure; in the former case the upper bound on the conditioning region is determined by confidence level 1.
risk.measure	character string or function specifying the risk measure to be computed on the row sums of x based on the given level(s) in order to determine the conditioning region.
include.conditional	logical indicating whether the computed sub-sample of x is to be returned, too.
...	additional arguments passed to risk.measure.

Details

The result of alloc_ellip() for loc = 0 can be found in McNeil et al. (2015, Corollary 8.43). Otherwise, McNeil et al. (2015, Theorem 8.28 (1)) can be used to derive the result.

Value

$d$-vector of allocated amounts (the allocation) according to the Euler principle under the assumption that the underlying loss random vector follows a $d$-dimensional elliptical distribution with location vector loc ($\bm{mu}$ in the reference) and scale matrix scale ($\Sigma$ in the reference, a covariance matrix) and that the risk measure is law-invariant, positive-homogeneous and translation invariant.

Author(s)

Marius Hofert

References

McNeil, A. J., Frey, R. and Embrechts, P. (2015). Quantitative Risk Management: Concepts, Techniques, Tools. Princeton University Press.

Examples

### Ellipitical case ###########################################################

## Construct a covariance matrix
sig <- 1:3 # standard deviations
library(copula) # for p2P() here
P <- p2P(c(-0.5, 0.3, 0.5)) # (3, 3) correlation matrix
Sigma <- P * sig %*% t(sig) # corresponding covariance matrix
stopifnot(all.equal(cov2cor(Sigma), P)) # sanity check

## Compute the allocation of 1.2 for a joint loss L ~ E_3(0, Sigma, psi)
AC <- alloc_ellip(1.2, loc = 0, scale = Sigma) # allocated amounts
stopifnot(all.equal(sum(AC), 1.2)) # sanity check
## Be careful to check whether the aforementioned assumptions hold.


### Nonparametrically ##########################################################

## Generate data
set.seed(271)
X <- qt(rCopula(1e5, copula = gumbelCopula(2, dim = 5)), df = 3.5)

## Estimate an allocation via MC based on a sub-sample whose row sums have a
## nonparametric VaR with confidence level in ...
alloc_np(X, level = 0.9) # ... (0.9, 1]
CA  <- alloc_np(X, level = c(0.9, 0.95)) # ... in (0.9, 0.95]
CA. <- alloc_np(X, level = c(0.9, 0.95), risk.measure = VaR_np) # providing a function
stopifnot(identical(CA, CA.))

---

Black–Scholes formula and the Greeks

Description

Compute the Black–Scholes formula and the Greeks.

Usage

Black_Scholes(t, S, r, sigma, K, T, type = c("call", "put"))
Black_Scholes_Greeks(t, S, r, sigma, K, T, type = c("call", "put"))

Arguments

t	initial or current time $t$ (in years).
S	stock price at time $t$.
r	risk-free annual interest rate.
sigma	annual volatility (standard deviation).
K	strike.
T	maturity (in years).
type	character string indicating whether a call (the default) or a put option is considered.

Details

Note again that t is time in years. In the context of McNeil et al. (2015, Chapter 9), this is $\tau_t = t/250$.

Value

Black_Scholes() returns the value of a European-style call or put option (depending on the chosen type) on a non-dividend paying stock.

Black_Scholes_Greeks() returns the first-order derivatives delta, theta, rho, vega and the second-order derivatives gamma, vanna and vomma (depending on the chosen type) in this order.

Author(s)

Marius Hofert

References

McNeil, A. J., Frey, R., and Embrechts, P. (2015). Quantitative Risk Management: Concepts, Techniques, Tools. Princeton University Press.

---

Brownian and Related Motions

Description

Simulate paths of dependent Brownian motions, geometric Brownian motions and Brownian bridges based on given increment copula samples. And extract copula increments from paths of dependent Brownian motions and geometric Brownian motions.

Usage

rBrownian(N, t, d = 1, U = matrix(runif(N * n * d), ncol = d),
          drift = 0, vola = 1, type = c("BM", "GBM", "BB"), init = 1)
deBrowning(x, t, drift = 0, vola = 1, type = c("BM", "GBM"))

Arguments

N	number $N$ of paths to simulate (positive integer).
x	$n+1$-vector containing one path of the specified stochastic process or $(n+1, d)$-matrix containing one path of the specified $d$ stochastic processes or $(N, n+1, d)$-array containing $N$ paths of the specified $d$ stochastic processes.
t	$n+1$-vector of the form $(t_0,\dots,t_n)$ with $0 = t_0 < \dots < t_n$ containing the time points where the stochastic processes are considered.
d	number $d$ of stochastic processes to simulate (positive integer).
U	$(N\cdot n, d)$-matrix of copula realizations to be converted to the joint increments of the stochastic processes.
drift	$d$-vector or number (then recycled to a $d$-vector) of drifts (typically denoted by $\mu$). Note that risk-neutral drifts are $r - \sigma^2/2$, where $r$ is the risk-free interest rate and $\sigma$ the volatility.
vola	$d$-vector or number (then recycled to a $d$-vector) of volatilities (typically denoted by $\sigma$).
type	character string indicating whether a Brownian motion ("BM"), geometric Brownian motion ("GBM") or Brownian bridge ("BB") is to be considered.
init	$d$-vector or number (then recycled to a $d$-vector) of initial values (typically stock prices at time 0) for type = "GBM".

Value

rBrownian() returns an $(N, n+1, d)$-array containing the $N$ paths of the specified $d$ stochastic processes simulated at the $n+1$ time points ($t_0 = 0, t_1,\dots, t_n$).

deBrowning() returns an $(N, n, d)$-array containing the $N$ paths of the copula increments of the $d$ stochastic processes over the $n+1$ time points ($t_0 = 0, t_1,\dots, t_n$).

Author(s)

Marius Hofert

Examples

## Setup
d <- 3 # dimension
library(copula)
tcop <- tCopula(iTau(tCopula(), tau = 0.5), dim = d, df = 4) # t_4 copula
vola <- seq(0.05, 0.20, length.out = d) # volatilities sigma
r <- 0.01 # risk-free interest rate
drift <- r - vola^2/2 # marginal drifts
init <- seq(10, 100, length.out = d) # initial stock prices
N <- 100 # number of replications
n <- 25 # number of time intervals
t <- 0:n/n # time points 0 = t_0 < ... < t_n

## Simulate N paths of a cross-sectionally dependent d-dimensional
## (geometric) Brownian motion ((G)BM) over n time steps
set.seed(271)
U <- rCopula(N * n, copula = tcop) # for dependent increments
X <- rBrownian(N, t = t, d = d, U = U, drift = drift, vola = vola) # BM
S <- rBrownian(N, t = t, d = d, U = U, drift = drift, vola = vola,
               type = "GBM", init = init) # GBM
stopifnot(dim(X) == c(N, n+1, d), dim(S) == c(N, n+1, d))

## DeBrowning
Z.X <- deBrowning(X, t = t, drift = drift, vola = vola) # BM
Z.S <- deBrowning(S, t = t, drift = drift, vola = vola, type = "GBM") # GBM
stopifnot(dim(Z.X) == c(N, n, d), dim(Z.S) == c(N, n, d))
## Note that for BMs, one loses one observation as X_{t_0} = 0 (or some other
## fixed value, so there is no random increment there that can be deBrowned.


    ## If we map the increments back to their copula sample, do we indeed
    ## see the copula samples again?
    U.Z.X <- pnorm(Z.X) # map to copula sample
    U.Z.S <- pnorm(Z.S) # map to copula sample
    stopifnot(all.equal(U.Z.X, U.Z.S)) # sanity check
    ## Visual check
    pairs(U.Z.X[,1,], gap = 0) # check at the first time point of the BM
    pairs(U.Z.X[,n,], gap = 0) # check at the last time point of the BM
    pairs(U.Z.S[,1,], gap = 0) # check at the first time point of the GBM
    pairs(U.Z.S[,n,], gap = 0) # check at the last time point of the GBM
    ## Numerical check
    ## First convert the (N * n, d)-matrix U to an (N, n, d)-array but in
    ## the right way (array(U, dim = c(N, n, d)) would use the U's in the
    ## wrong order)
    U. <- aperm(array(U, dim = c(n, N, d)), perm = c(2,1,3))
    ## Now compare
    stopifnot(all.equal(U.Z.X, U., check.attributes = FALSE))
    stopifnot(all.equal(U.Z.S, U., check.attributes = FALSE))



    ## Generate dependent GBM sample paths with quasi-random numbers
    library(qrng)
    set.seed(271)
    U.. <- cCopula(to_array(sobol(N, d = d * n, randomize = "digital.shift"), f = n),
                   copula = tcop, inverse = TRUE)
    S. <- rBrownian(N, t = t, d = d, U = U.., drift = drift, vola = vola,
                    type = "GBM", init = init)
    pairs(S [,2,], gap = 0) # pseudo-samples at t_1
    pairs(S.[,2,], gap = 0) # quasi-samples at t_1
    pairs(S [,n+1,], gap = 0) # pseudo-samples at t_n
    pairs(S.[,n+1,], gap = 0) # quasi-samples at t_n



    ## Generate paths from a Brownian bridge
    B <- rBrownian(N, t = t, type = "BB")
    plot(NA, xlim = 0:1, ylim = range(B),
         xlab = "Time t", ylab = expression("Brownian bridge path"~(B[t])))
    for(i in 1:N)
        lines(t, B[i,,], col = adjustcolor("black", alpha.f = 25/N))

---

Catching Results, Warnings and Errors Simultaneously

Description

Catches results, warnings and errors.

Usage

catch(expr)

Arguments

expr	expression to be evaluated, typically a function call.

Details

This function is particularly useful for large(r) simulation studies to not fail until finished.

Value

list with components:

value	value of expr or NULL in case of an error.
warning	warning message (see simpleWarning or warning()) or NULL in case of no warning.
error	error message (see simpleError or stop()) or NULL in case of no error.

Author(s)

Marius Hofert (based on doCallWE() and tryCatch.W.E() in the R package simsalapar).

Examples

catch(log(2))
catch(log(-1))
catch(log("a"))

---

Fitting ARMA-GARCH Processes

Description

Fail-safe componentwise fitting of univariate ARMA-GARCH processes.

Usage

fit_ARMA_GARCH(x, ugarchspec.list = ugarchspec(), solver = "hybrid",
               verbose = FALSE, ...)

Arguments

x	matrix-like data structure, possibly an xts object.
ugarchspec.list	object of class uGARCHspec (as returned by ugarchspec()) or a list of such. In case of a list, its length has to be equal to the number of columns of x. ugarchspec.list provides the ARMA-GARCH specifications for each of the time series (columns of x).
solver	string indicating the solver used; see ?ugarchfit.
verbose	logical indicating whether verbose output is given.
...	additional arguments passed to the underlying ugarchfit().

Value

If x consists of one column only (e.g. a vector), ARMA_GARCH() returns the fitted object; otherwise it returns a list of such.

Author(s)

Marius Hofert

See Also

fit_GARCH_11() for fast(er) and numerically more robust fitting of GARCH(1,1) processes.

Examples

library(rugarch)
library(copula)

## Read the data, build -log-returns
data(SMI.12) # Swiss Market Index data
stocks <- c("CSGN", "BAER", "UBSN", "SREN", "ZURN") # components we work with
x <- SMI.12[, stocks]
X <- -returns(x)
n <- nrow(X)
d <- ncol(X)

## Fit ARMA-GARCH models to the -log-returns
## Note: - Our choice here is purely for demonstration purposes.
##         The models are not necessarily adequate
##       - The sample size n is *too* small here for properly capturing GARCH effects.
##         Again, this is only for demonstration purposes here.
uspec <- c(rep(list(ugarchspec(distribution.model = "std")), d-2), # ARMA(1,1)-GARCH(1,1)
           list(ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(2,2)),
                           distribution.model = "std")),
           list(ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(2,1)),
                           mean.model = list(armaOrder = c(1,2), include.mean = TRUE),
                           distribution.model = "std")))
system.time(fitAG <- fit_ARMA_GARCH(X, ugarchspec.list = uspec))
str(fitAG, max.level = 1) # list with components fit, warning, error
## Now access the list to check

## Not run: 
## Pick out the standardized residuals, plot them and fit a t copula to them
## Note: ugarchsim() needs the residuals to be standardized; working with
##       standardize = FALSE still requires to simulate them from the
##       respective standardized marginal distribution functions.
Z <- sapply(fitAG$fit, residuals, standardize = TRUE)
U <- pobs(Z)
pairs(U, gap = 0)
system.time(fitC <- fitCopula(tCopula(dim = d, dispstr = "un"), data = U,
                              method = "mpl"))

## Simulate (standardized) Z
set.seed(271)
U. <- rCopula(n, fitC@copula) # simulate from the fitted copula
nu <- sapply(1:d, function(j) fitAG$fit[[j]]@fit$coef["shape"]) # extract (fitted) d.o.f. nu
Z. <- sapply(1:d, function(j) sqrt((nu[j]-2)/nu[j]) * qt(U.[,j], df = nu[j])) # Z

## Simulate from fitted model
X. <- sapply(1:d, function(j)
    fitted(ugarchsim(fitAG$fit[[j]], n.sim = n, m.sim = 1, startMethod = "sample",
                     rseed = 271, custom.dist = list(name = "sample",
                                                     distfit = Z.[,j, drop = FALSE]))))

## Plots original vs simulated -log-returns
opar <- par(no.readonly = TRUE)
layout(matrix(1:(2*d), ncol = d)) # layout
ran <- range(X, X.)
for(j in 1:d) {
    plot(X[,j],  type = "l", ylim = ran, ylab = paste(stocks[j], "-log-returns"))
    plot(X.[,j], type = "l", ylim = ran, ylab = "Simulated -log-returns")
}
    par(opar)

## End(Not run)

---

Fast(er) and Numerically More Robust Fitting of GARCH(1,1) Processes

Description

Fast(er) and numerically more robust fitting of GARCH(1,1) processes according to Zumbach (2000).

Usage

fit_GARCH_11(x, init = NULL, sig2 = mean(x^2), delta = 1,
             distr = c("norm", "st"), control = list(), ...)
tail_index_GARCH_11(innovations, alpha1, beta1,
                    interval = c(1e-6, 1e2), ...)

Arguments

x	vector of length $n$ containing the data (typically log-returns) to be fitted a GARCH(1,1) to.
init	vector of length 2 giving the initial values for the likelihood fitting. Note that these are initial values for $z_{corr}$ and $z_{ema}$ as in Zumbach (2000).
sig2	annualized variance (third parameter of the reparameterization according to Zumbach (2000)).
delta	unit of time (defaults to 1 meaning daily data; for yearly data, use 250).
distr	character string specifying the innovation distribution ("norm" for N(0,1) or "st" for a standardized $t$ distribution).
control	see ?optim().
innovations	random variates from the innovation distribution; for example, obtained via rnorm() or rt(, df = nu) * sqrt((nu-2)/nu) where nu are the d.o.f. of the $t$ distribution.
alpha1	nonnegative GARCH(1,1) coefficient $alpha[1]$ satisfying $alpha[1] + beta[1] < 1$.
beta1	nonnegative GARCH(1,1) coefficient $beta[1]$ satisfying $alpha[1] + beta[1] < 1$.
interval	initial interval for computing the tail index; passed to the underlying uniroot().
...	fit_GARCH_11(): additional arguments passed to the underlying optim(). tail_index_GARCH_11(): additional arguments passed to the underlying uniroot().

Value

fit_GARCH_11():

coef:

estimated coefficients $\alpha_0$, $\alpha_1$, $\beta_1$ and, if distr = "st" the estimated degrees of freedom.

logLik:

maximized log-likelihood.

counts:

number of calls to the objective function; see ?optim.

convergence:

convergence code ('0' indicates successful completion); see ?optim.

message:

see ?optim.

sig.t:

vector of length $n$ giving the conditional volatility.

Z.t:

vector of length $n$ giving the standardized residuals.

tail_index_GARCH_11():

The tail index $alpha$ estimated by Monte Carlo via McNeil et al. (2015, p. 576), so the $alpha$ which solves

$E({(\alpha_1Z^2 + \beta_1)}^{\alpha/2}) = 1$

, where $Z$ are the innovations. If no solution is found (e.g. if the objective function does not have different sign at the endpoints of interval), NA is returned.

Author(s)

Marius Hofert

References

Zumbach, G. (2000). The pitfalls in fitting GARCH (1,1) processes. Advances in Quantitative Asset Management 1, 179–200.

McNeil, A. J., Frey, R. and Embrechts, P. (2015). Quantitative Risk Management: Concepts, Techniques, Tools. Princeton University Press.

See Also

fit_ARMA_GARCH() based on rugarch.

Examples

### Example 1: N(0,1) innovations ##############################################

## Generate data from a GARCH(1,1) with N(0,1) innovations
library(rugarch)
uspec <- ugarchspec(variance.model = list(model = "sGARCH",
                                          garchOrder = c(1, 1)),
                    distribution.model = "norm",
                    mean.model = list(armaOrder = c(0, 0)),
                    fixed.pars = list(mu = 0,
                                      omega = 0.1, # alpha_0
                                      alpha1 = 0.2, # alpha_1
                                      beta1 = 0.3)) # beta_1
X <- ugarchpath(uspec, n.sim = 1e4, rseed = 271) # sample (set.seed() fails!)
X.t <- as.numeric(X@path$seriesSim) # actual path (X_t)

## Fitting via ugarchfit()
uspec. <- ugarchspec(variance.model = list(model = "sGARCH",
                                           garchOrder = c(1, 1)),
                     distribution.model = "norm",
                     mean.model = list(armaOrder = c(0, 0)))
fit <- ugarchfit(uspec., data = X.t)
coef(fit) # fitted mu, alpha_0, alpha_1, beta_1
Z <- fit@fit$z # standardized residuals
stopifnot(all.equal(mean(Z), 0, tol = 1e-2),
          all.equal(var(Z),  1, tol = 1e-3))

## Fitting via fit_GARCH_11()
fit. <- fit_GARCH_11(X.t)
fit.$coef # fitted alpha_0, alpha_1, beta_1
Z. <- fit.$Z.t # standardized residuals
stopifnot(all.equal(mean(Z.), 0, tol = 5e-3),
          all.equal(var(Z.),  1, tol = 1e-3))

## Compare
stopifnot(all.equal(fit.$coef, coef(fit)[c("omega", "alpha1", "beta1")],
                    tol = 5e-3, check.attributes = FALSE)) # fitted coefficients
summary(Z. - Z) # standardized residuals


### Example 2: t_nu(0, (nu-2)/nu) innovations ##################################

## Generate data from a GARCH(1,1) with t_nu(0, (nu-2)/nu) innovations
uspec <- ugarchspec(variance.model = list(model = "sGARCH",
                                          garchOrder = c(1, 1)),
                    distribution.model = "std",
                    mean.model = list(armaOrder = c(0, 0)),
                    fixed.pars = list(mu = 0,
                                      omega = 0.1, # alpha_0
                                      alpha1 = 0.2, # alpha_1
                                      beta1 = 0.3, # beta_1
                                      shape = 4)) # nu
X <- ugarchpath(uspec, n.sim = 1e4, rseed = 271) # sample (set.seed() fails!)
X.t <- as.numeric(X@path$seriesSim) # actual path (X_t)

## Fitting via ugarchfit()
uspec. <- ugarchspec(variance.model = list(model = "sGARCH",
                                           garchOrder = c(1, 1)),
                     distribution.model = "std",
                     mean.model = list(armaOrder = c(0, 0)))
fit <- ugarchfit(uspec., data = X.t)
coef(fit) # fitted mu, alpha_0, alpha_1, beta_1, nu
Z <- fit@fit$z # standardized residuals
stopifnot(all.equal(mean(Z), 0, tol = 1e-2),
          all.equal(var(Z),  1, tol = 5e-2))

## Fitting via fit_GARCH_11()
fit. <- fit_GARCH_11(X.t, distr = "st")
c(fit.$coef, fit.$df) # fitted alpha_0, alpha_1, beta_1, nu
Z. <- fit.$Z.t # standardized residuals
stopifnot(all.equal(mean(Z.), 0, tol = 2e-2),
          all.equal(var(Z.),  1, tol = 2e-2))

## Compare
fit.coef <- coef(fit)[c("omega", "alpha1", "beta1", "shape")]
fit..coef <- c(fit.$coef, fit.$df)
stopifnot(all.equal(fit.coef, fit..coef, tol = 7e-2, check.attributes = FALSE))
summary(Z. - Z) # standardized residuals

---

Parameter Estimators of the Generalized Extreme Value Distribution

Description

Quantile matching estimator, probability weighted moments estimator, log-likelihood and maximum-likelihood estimator for the parameters of the generalized extreme value distribution (GEV).

Usage

fit_GEV_quantile(x, p = c(0.25, 0.5, 0.75), cutoff = 3)
fit_GEV_PWM(x)

logLik_GEV(param, x)
fit_GEV_MLE(x, init = c("shape0", "PWM", "quantile"),
            estimate.cov = TRUE, control = list(), ...)

Arguments

x	numeric vector of data. In the block maxima method, these are the block maxima.
p	numeric(3) specifying the probabilities whose quantiles are matched.
cutoff	positive $z$ after which $\exp(-z)$ is truncated to 0.
param	numeric(3) containing the value of the shape $\xi$ (a real), location $\mu$ (a real) and scale $\sigma$ (positive real) parameters of the GEV distribution in this order.
init	character string specifying the method for computing initial values. Can also be numeric(3) for directly providing $\xi$, $\mu$, $\sigma$.
estimate.cov	logical indicating whether the asymptotic covariance matrix of the parameter estimators is to be estimated (inverse of observed Fisher information (negative Hessian of log-likelihood evaluated at MLE)) and standard errors for the estimators of $\xi$, $\mu$, $\sigma$ returned, too.
control	list; passed to the underlying optim().
...	additional arguments passed to the underlying optim().

Details

fit_GEV_quantile() matches the empirical p-quantiles.

fit_GEV_PWM() computes the probability weighted moments (PWM) estimator of Hosking et al. (1985); see also Landwehr and Wallis (1979).

fit_GEV_MLE() uses, as default, the case $\xi = 0$ for computing initial values; this is actually a small positive value since Nelder–Mead could fail otherwise. For the other available methods for computing initial values, $\sigma$ (obtained from the case $\xi = 0$) is doubled in order to guarantee a finite log-likelihood at the initial values. After several experiments (see the source code), one can safely say that finding initial values for fitting GEVs via MLE is non-trivial; see also the block maxima method script about the Black Monday event on https://qrmtutorial.org.

Caution: See Coles (2001, p. 55) for how to interpret $\xi\le -0.5$; in particular, the standard asymptotic properties of the MLE do not apply.

Value

fit_GEV_quantile() and fit_GEV_PWM() return a numeric(3) giving the parameter estimates for the GEV distribution.

logLik_GEV() computes the log-likelihood of the GEV distribution (-Inf if not admissible).

fit_GEV_MLE() returns the return object of optim() (by default, the return value value is the log-likelihood) and, appended, the estimated asymptotic covariance matrix and standard errors of the parameter estimators, if estimate.cov.

Author(s)

Marius Hofert

References

McNeil, A. J., Frey, R. and Embrechts, P. (2015). Quantitative Risk Management: Concepts, Techniques, Tools. Princeton University Press.

Hosking, J. R. M., Wallis, J. R. and Wood, E. F. (1985). Estimation of the Generalized Extreme-Value Distribution by the Method of Probability-Weighted Moments. Technometrics 27(3), 251–261.

Landwehr, J. M. and Wallis, J. R. (1979). Probability Weighted Moments Compared With Some Traditional Techniques in Estimating Gumbel Parameters and Quantiles. Water Resourches Research 15(5), 1055–1064.

Coles, S. (2001). An Introduction to Statistical Modeling of Extreme Values. Springer-Verlag.

Examples

## Simulate some data
xi <- 0.5
mu <- -2
sig <- 3
n <- 1000
set.seed(271)
X <- rGEV(n, shape = xi, loc = mu, scale = sig)

## Fitting via matching quantiles
(fit.q <- fit_GEV_quantile(X))
stopifnot(all.equal(fit.q[["shape"]], xi,  tol = 0.12),
          all.equal(fit.q[["loc"]],   mu,  tol = 0.12),
          all.equal(fit.q[["scale"]], sig, tol = 0.005))

## Fitting via PWMs
(fit.PWM <- fit_GEV_PWM(X))
stopifnot(all.equal(fit.PWM[["shape"]], xi,  tol = 0.16),
          all.equal(fit.PWM[["loc"]],   mu,  tol = 0.15),
          all.equal(fit.PWM[["scale"]], sig, tol = 0.08))

## Fitting via MLE
(fit.MLE <- fit_GEV_MLE(X))
(est <- fit.MLE$par) # estimates of xi, mu, sigma
stopifnot(all.equal(est[["shape"]], xi,  tol = 0.07),
          all.equal(est[["loc"]],   mu,  tol = 0.12),
          all.equal(est[["scale"]], sig, tol = 0.06))
fit.MLE$SE # estimated asymp. variances of MLEs = std. errors of MLEs

## Plot the log-likelihood in the shape parameter xi for fixed
## location mu and scale sigma (fixed as generated)
xi. <- seq(-0.1, 0.8, length.out = 65)
logLik <- sapply(xi., function(xi..) logLik_GEV(c(xi.., mu, sig), x = X))
plot(xi., logLik, type = "l", xlab = expression(xi),
     ylab = expression("GEV distribution log-likelihood for fixed"~mu~"and"~sigma))
## => Numerically quite challenging (for this seed!)

## Plot the profile likelihood for these xi's
## Note: As initial values for the nuisance parameters mu, sigma, we
##       use their values in the case xi = 0 (for all fixed xi = xi.,
##       in particular those xi != 0). Furthermore, for the given data X
##       and xi = xi., we make sure the initial value for sigma is so large
##       that the density is not 0 and thus the log-likelihood is finite.
pLL <- sapply(xi., function(xi..) {
    scale.init <- sqrt(6 * var(X)) / pi
    loc.init <- mean(X) - scale.init * 0.5772157
    while(!is.finite(logLik_GEV(c(xi.., loc.init, scale.init), x = X)) &&
          is.finite(scale.init)) scale.init <- scale.init * 2
    optim(c(loc.init, scale.init), fn = function(nuis)
                logLik_GEV(c(xi.., nuis), x = X),
    		        control = list(fnscale = -1))$value
})
plot(xi., pLL, type = "l", xlab = expression(xi),
     ylab = "GEV distribution profile log-likelihood")

---

Parameter Estimators of the Generalized Pareto Distribution

Description

Method-of-moments estimator, probability weighted moments estimator, log-likelihood and maximum-likelihood estimator for the parameters of the generalized Pareto distribution (GPD).

Usage

fit_GPD_MOM(x)
fit_GPD_PWM(x)

logLik_GPD(param, x)
fit_GPD_MLE(x, init = c("PWM", "MOM", "shape0"),
            estimate.cov = TRUE, control = list(), ...)

Arguments

x	numeric vector of data. In the peaks-over-threshold method, these are the excesses (exceedances minus threshold).
param	numeric(2) containing the value of the shape $\xi$ (a real) and scale $\beta$ (positive real) parameters of the GPD in this order.
init	character string specifying the method for computing initial values. Can also be numeric(2) for directly providing $\xi$ and $\beta$.
estimate.cov	logical indicating whether the asymptotic covariance matrix of the parameter estimators is to be estimated (inverse of observed Fisher information (negative Hessian of log-likelihood evaluated at MLE)) and standard errors for the estimators of $\xi$ and $\beta$ returned, too.
control	list; passed to the underlying optim().
...	additional arguments passed to the underlying optim().

Details

fit_GPD_MOM() computes the method-of-moments (MOM) estimator.

fit_GPD_PWM() computes the probability weighted moments (PWM) estimator of Hosking and Wallis (1987); see also Landwehr et al. (1979).

fit_GPD_MLE() uses, as default, fit_GPD_PWM() for computing initial values. The former requires the data x to be non-negative and adjusts $\beta$ if $\xi$ is negative, so that the log-likelihood at the initial value should be finite.

Value

fit_GEV_MOM() and fit_GEV_PWM() return a numeric(3) giving the parameter estimates for the GPD.

logLik_GPD() computes the log-likelihood of the GPD (-Inf if not admissible).

fit_GPD_MLE() returns the return object of optim() and, appended, the estimated asymptotic covariance matrix and standard errors of the parameter estimators, if estimate.cov.

Author(s)

Marius Hofert

References

McNeil, A. J., Frey, R. and Embrechts, P. (2015). Quantitative Risk Management: Concepts, Techniques, Tools. Princeton University Press.

Hosking, J. R. M. and Wallis, J. R. (1987). Parameter and Quantile Estimation for the Generalized Pareto Distribution. Technometrics 29(3), 339–349.

Landwehr, J. M., Matalas, N. C. and Wallis, J. R. (1979). Estimation of Parameters and Quantiles of Wakeby Distributions. Water Resourches Research 15(6), 1361–1379.

Examples

## Simulate some data
xi <- 0.5
beta <- 3
n <- 1000
set.seed(271)
X <- rGPD(n, shape = xi, scale = beta)

## Fitting via matching moments
(fit.MOM <- fit_GPD_MOM(X))
stopifnot(all.equal(fit.MOM[["shape"]], xi,   tol = 0.52),
          all.equal(fit.MOM[["scale"]], beta, tol = 0.24))

## Fitting via PWMs
(fit.PWM <- fit_GPD_PWM(X))
stopifnot(all.equal(fit.PWM[["shape"]], xi,   tol = 0.2),
          all.equal(fit.PWM[["scale"]], beta, tol = 0.12))

## Fitting via MLE
(fit.MLE <- fit_GPD_MLE(X))
(est <- fit.MLE$par) # estimates of xi, mu, sigma
stopifnot(all.equal(est[["shape"]], xi,   tol = 0.12),
          all.equal(est[["scale"]], beta, tol = 0.11))
fit.MLE$SE # estimated asymp. variances of MLEs = std. errors of MLEs

## Plot the log-likelihood in the shape parameter xi for fixed
## scale beta (fixed as generated)
xi. <- seq(-0.1, 0.8, length.out = 65)
logLik <- sapply(xi., function(xi..) logLik_GPD(c(xi.., beta), x = X))
plot(xi., logLik, type = "l", xlab = expression(xi),
     ylab = expression("GPD log-likelihood for fixed"~beta))

## Plot the profile likelihood for these xi's
## (with an initial interval for the nuisance parameter beta such that
##  logLik_GPD() is finite)
pLL <- sapply(xi., function(xi..) {
    ## Choose beta interval for optimize()
    int <- if(xi.. >= 0) {
               ## Method-of-Moment estimator
               mu.hat <- mean(X)
               sig2.hat <- var(X)
               shape.hat <- (1-mu.hat^2/sig2.hat)/2
               scale.hat <- mu.hat*(1-shape.hat)
               ## log-likelihood always fine for xi.. >= 0 for all beta
               c(1e-8, 2 * scale.hat)
           } else { # xi.. < 0
               ## Make sure logLik_GPD() is finite at endpoints of int
               mx <- max(X)
               -xi.. * mx * c(1.01, 100) # -xi * max(X) * scaling
               ## Note: for shapes xi.. closer to 0, the upper scaling factor
               ##       needs to be chosen sufficiently large in order
               ##       for optimize() to find an optimum (not just the
               ##       upper end point). Try it with '2' instead of '100'.
           }
    ## Optimization
    optimize(function(nuis) logLik_GPD(c(xi.., nuis), x = X),
             interval = int, maximum = TRUE)$maximum
})
plot(xi., pLL, type = "l", xlab = expression(xi),
     ylab = "GPD profile log-likelihood")

---

Tools for Getting and Working with Data

Description

Download (and possibly) merge data from freely available databases.

Usage

get_data(x, from = NULL, to = NULL,
         src = c("yahoo", "quandl", "oanda", "FRED", "google"),
         FUN = NULL, verbose = TRUE, warn = TRUE, ...)

Arguments

x	vector of ticker symbols (e.g. "^GSPC" if src = "yahoo" or "EUR/USD" if src = "oanda").
from	start date as a Date object or character string (in international date format "yyyy-mm-dd"); if NULL, the earliest date with available data is picked.
to	end date as a Date object or character string (in international date format "yyyy-mm-dd"); if NULL, the last date with available data is picked.
src	character string specifying the data source (e.g. "yahoo" for stocks or "oanda" for FX data); see getSymbols() and Quandl().
FUN	function to be applied to the data before being returned. This can be the identity: if the data could not be retrieved (and is thus replaced by NA); the given FUN: if FUN has been provided; a useful default: if FUN = NULL; the default uses the adjusted close price Ad() if src = "yahoo", the close price Cl() if src = "google" and the identity otherwise.
verbose	logical indicating whether progress monitoring should be done.
warn	logical indicating whether a warning is given showing the error message when fetching x fails.
...	additional arguments passed to the underlying getSymbols() from quantmod or Quandl() from Quandl (if src = "quandl").

Details

FUN is typically one of quantmod's Op, Hi, Lo, Cl, Vo, Ad or one of the combined functions OpCl, ClCl, HiCl, LoCl, LoHi, OpHi, OpLo, OpOp.

Value

xts object containing the data with column name(s) adjusted to be the ticker symbol (in case lengths match; otherwise the column names are not adjusted); NA if data is not available.

Author(s)

Marius Hofert

Examples

## Not run: 
## Note: This needs a working internet connection
## Get stock and volatility data (for all available trading days)
dat <- get_data(c("^GSPC", "^VIX")) # note: this needs a working internet connection
## Plot them (Alternative: plot.xts() from xtsExtra)
library(zoo)
plot.zoo(dat, screens = 1, main = "", xlab = "Trading day", ylab = "Value")

## End(Not run)

---

Generalized Extreme Value Distribution

Description

Density, distribution function, quantile function and random variate generation for the generalized extreme value distribution (GEV).

Usage

dGEV(x, shape, loc = 0, scale = 1, log = FALSE)
pGEV(q, shape, loc = 0, scale = 1, lower.tail = TRUE, log.p = FALSE)
qGEV(p, shape, loc = 0, scale = 1, lower.tail = TRUE, log.p = FALSE)
rGEV(n, shape, loc = 0, scale = 1)

Arguments

x, q	vector of quantiles.
p	vector of probabilities.
n	number of observations.
shape	GEV shape parameter $\xi$, a real.
loc	GEV location parameter $\mu$, a real.
scale	GEV scale parameter $\sigma$, a positive real.
lower.tail	logical; if TRUE (default) probabilities are $P(X \le x)$ otherwise, $P(X > x)$.
log, log.p	logical; if TRUE, probabilities p are given as log(p).

Details

The distribution function of the generalized extreme value distribution is given by

$F(x) = \left\{ \begin{array}{ll} \exp(-(1-\xi(x-\mu)/\sigma)^{-1/\xi}), & \xi\neq 0,\ 1+\xi(x-\mu)/\sigma>0,\\ \exp(-e^{-(x-\mu)/\sigma}), & \xi = 0, \end{array}\right.$

where $\sigma>0$.

Value

dGEV() computes the density, pGEV() the distribution function, qGEV() the quantile function and rGEV() random variates of the generalized extreme value distribution.

Author(s)

Marius Hofert

References

McNeil, A. J., Frey, R., and Embrechts, P. (2015). Quantitative Risk Management: Concepts, Techniques, Tools. Princeton University Press.

Examples

## Basic sanity checks
plot(pGEV(rGEV(1000, shape = 0.5), shape = 0.5)) # should be U[0,1]
curve(dGEV(x, shape = 0.5), from = -3, to = 5)

---

Fitted GEV Shape as a Function of the Threshold

Description

Fit GEVs to block maxima and plot the fitted GPD shape as a function of the block size.

Usage

GEV_shape_plot(x, blocksize = tail(pretty(seq_len(length(x)/20), n = 64), -1),
               estimate.cov = TRUE, conf.level = 0.95,
               CI.col = adjustcolor(1, alpha.f = 0.2),
               lines.args = list(), xlim = NULL, ylim = NULL,
               xlab = "Block size",  ylab = NULL,
               xlab2 = "Number of blocks", plot = TRUE, ...)

Arguments

x	vector of numeric data.
blocksize	numeric vector of block sizes for which to fit a GEV to the block maxima.
estimate.cov	logical indicating whether confidence intervals are to be computed.
conf.level	confidence level of the confidence intervals if estimate.cov.
CI.col	color of the pointwise asymptotic confidence intervals (CIs); if NA, no CIs are shown.
lines.args	list of arguments passed to the underlying lines() for drawing the shape parameter as a function of the block size.
xlim, ylim, xlab, ylab	see plot().
xlab2	label of the secondary x-axis.
plot	logical indicating whether a plot is produced.
...	additional arguments passed to the underlying plot().

Details

Such plots can be used in the block maxima method for determining the optimal block size (as the smallest after which the plot is (roughly) stable).

Value

Invisibly returns a list containing the block sizes considered, the corresponding block maxima and the fitted GEV distribution objects as returned by the underlying fit_GEV_MLE().

Author(s)

Marius Hofert

Examples

set.seed(271)
X <- rPar(5e4, shape = 4)
GEV_shape_plot(X)
abline(h = 1/4, lty = 3) # theoretical xi = 1/shape for Pareto

---

(Generalized) Pareto Distribution

Description

Density, distribution function, quantile function and random variate generation for the (generalized) Pareto distribution (GPD).

Usage

dGPD(x, shape, scale, log = FALSE)
pGPD(q, shape, scale, lower.tail = TRUE, log.p = FALSE)
qGPD(p, shape, scale, lower.tail = TRUE, log.p = FALSE)
rGPD(n, shape, scale)

dPar(x, shape, scale = 1, log = FALSE)
pPar(q, shape, scale = 1, lower.tail = TRUE, log.p = FALSE)
qPar(p, shape, scale = 1, lower.tail = TRUE, log.p = FALSE)
rPar(n, shape, scale = 1)

Arguments

x, q	vector of quantiles.
p	vector of probabilities.
n	number of observations.
shape	GPD shape parameter $\xi$ (a real number) and Pareto shape parameter $\theta$ (a positive number).
scale	GPD scale parameter $\beta$ (a positive number) and Pareto scale parameter $\kappa$ (a positive number).
lower.tail	logical; if TRUE (default) probabilities are $P(X \le x)$ otherwise, $P(X > x)$.
log, log.p	logical; if TRUE, probabilities p are given as log(p).

Details

The distribution function of the generalized Pareto distribution is given by

$F(x) = \left\{ \begin{array}{ll} 1-(1+\xi x/\beta)^{-1/\xi}, & \xi \neq 0,\\ 1-\exp(-x/\beta), & \xi = 0, \end{array}\right.$

where $\beta>0$ and $x\ge0$ if $\xi\ge 0$ and $x\in[0,-\beta/\xi]$ if $\xi<0$.

The distribution function of the Pareto distribution is given by

$F(x) = 1-(1+x/\kappa)^{-\theta},\ x\ge 0,$

where $\theta > 0$, $\kappa > 0$.

In contrast to dGPD(), pGPD(), qGPD() and rGPD(), the functions dPar(), pPar(), qPar() and rPar() are vectorized in their main argument and the parameters.

Value

dGPD() computes the density, pGPD() the distribution function, qGPD() the quantile function and rGPD() random variates of the generalized Pareto distribution.

Similary for dPar(), pPar(), qPar() and rPar() for the Pareto distribution.

Author(s)

Marius Hofert

References

McNeil, A. J., Frey, R., and Embrechts, P. (2015). Quantitative Risk Management: Concepts, Techniques, Tools. Princeton University Press.

Examples

## Basic sanity checks
curve(dGPD(x, shape = 0.5, scale = 3), from = -1, to = 5)
plot(pGPD(rGPD(1000, shape = 0.5, scale = 3), shape = 0.5, scale = 3)) # should be U[0,1]

---

Fitted GPD Shape as a Function of the Threshold

Description

Fit GPDs to various thresholds and plot the fitted GPD shape as a function of the threshold.

Usage

GPD_shape_plot(x, thresholds = seq(quantile(x, 0.5), quantile(x, 0.99),
                                   length.out = 65),
               estimate.cov = TRUE, conf.level = 0.95,
               CI.col = adjustcolor(1, alpha.f = 0.2),
               lines.args = list(), xlim = NULL, ylim = NULL,
               xlab = "Threshold", ylab = NULL,
               xlab2 = "Excesses", plot = TRUE, ...)

Arguments

x	vector of numeric data.
thresholds	numeric vector of thresholds for which to fit a GPD to the excesses.
estimate.cov	logical indicating whether confidence intervals are to be computed.
conf.level	confidence level of the confidence intervals if estimate.cov.
CI.col	color of the pointwise asymptotic confidence intervals (CIs); if NA, no CIs are shown.
lines.args	list of arguments passed to the underlying lines() for drawing the shape parameter as a function of the threshold.
xlim, ylim, xlab, ylab	see plot().
xlab2	label of the secondary x-axis.
plot	logical indicating whether a plot is produced.
...	additional arguments passed to the underlying plot().

Details

Such plots can be used in the peaks-over-threshold method for determining the optimal threshold (as the smallest after which the plot is (roughly) stable).

Value

Invisibly returns a list containing the thresholds considered, the corresponding excesses and the fitted GPD objects as returned by the underlying fit_GPD_MLE().

Author(s)

Marius Hofert

Examples

set.seed(271)
X <- rt(1000, df = 3.5)
GPD_shape_plot(X)

---

GPD-Based Tail Distribution (POT method)

Description

Density, distribution function, quantile function and random variate generation for the GPD-based tail distribution in the POT method.

Usage

dGPDtail(x, threshold, p.exceed, shape, scale, log = FALSE)
pGPDtail(q, threshold, p.exceed, shape, scale, lower.tail = TRUE, log.p = FALSE)
qGPDtail(p, threshold, p.exceed, shape, scale, lower.tail = TRUE, log.p = FALSE)
rGPDtail(n, threshold, p.exceed, shape, scale)

Arguments

x, q	vector of quantiles.
p	vector of probabilities.
n	number of observations.
threshold	threshold $u$ in the POT method.
p.exceed	probability of exceeding the threshold u; for the Smith estimator, this is mean(x > threshold) for x being the data.
shape	GPD shape parameter $\xi$ (a real number).
scale	GPD scale parameter $\beta$ (a positive number).
lower.tail	logical; if TRUE (default) probabilities are $P(X \le x)$ otherwise, $P(X > x)$.
log, log.p	logical; if TRUE, probabilities p are given as log(p).

Details

Let $u$ denote the threshold (threshold), $p_u$ the exceedance probability (p.exceed) and $F_{GPD}$ the GPD distribution function. Then the distribution function of the GPD-based tail distribution is given by

$F(q) = 1-p_u(1-F_{GPD}(q - u))$

. The quantile function is

$F^{-1}(p) = u + F_GPD^{-1}(1-(1-p)/p_u)$

and the density is

$f(x) = p_u f_{GPD}(x - u)$

, where $f_{GPD}$ denotes the GPD density.

Note that the distribution function has a jumpt of height $P(X \le u)$ (1-p.exceed) at $u$.

Value

dGPDtail() computes the density, pGPDtail() the distribution function, qGPDtail() the quantile function and rGPDtail() random variates of the GPD-based tail distribution in the POT method.

Author(s)

Marius Hofert

References

McNeil, A. J., Frey, R., and Embrechts, P. (2015). Quantitative Risk Management: Concepts, Techniques, Tools. Princeton University Press.

Examples

## Generate data to work with
set.seed(271)
X <- rt(1000, df = 3.5) # in MDA(H_{1/df}); see MFE (2015, Section 16.1.1)

## Determine thresholds for POT method
mean_excess_plot(X[X > 0])
abline(v = 1.5)
u <- 1.5 # threshold

## Fit GPD to the excesses (per margin)
fit <- fit_GPD_MLE(X[X > u] - u)
fit$par
1/fit$par["shape"] # => close to df

## Estimate threshold exceedance probabilities
p.exceed <- mean(X > u)

## Define corresponding densities, distribution function and RNG
dF <- function(x) dGPDtail(x, threshold = u, p.exceed = p.exceed,
                           shape = fit$par["shape"], scale = fit$par["scale"])
pF <- function(q) pGPDtail(q, threshold = u, p.exceed = p.exceed,
                           shape = fit$par["shape"], scale = fit$par["scale"])
rF <- function(n) rGPDtail(n, threshold = u, p.exceed = p.exceed,
                           shape = fit$par["shape"], scale = fit$par["scale"])

## Basic check of dF()
curve(dF, from = u - 1, to = u + 5)

## Basic check of pF()
curve(pF, from = u, to = u + 5, ylim = 0:1) # quite flat here
abline(v = u, h = 1-p.exceed, lty = 2) # mass at u is 1-p.exceed (see 'Details')

## Basic check of rF()
set.seed(271)
X. <- rF(1000)
plot(X., ylab = "Losses generated from the fitted GPD-based tail distribution")
stopifnot(all.equal(mean(X. == u), 1-p.exceed, tol = 7e-3)) # confirms the above
## Pick out 'continuous part'
X.. <- X.[X. > u]
plot(pF(X..), ylab = "Probability-transformed tail losses") # should be U[1-p.exceed, 1]

---

Construction of Hierarchical Matrices

Description

Constructing hierarchical matrices, used, for example, for hierarchical dependence models, clustering, etc.

Usage

hierarchical_matrix(x, diagonal = rep(1, d))

Arguments

x	list of length 2 or 3 containing the homogeneous numeric entry of the current block of the hierarchical matrix, the integer components belongning to the current block (or NULL) and, possibly, another (nested) list of the same type.
diagonal	diagonal elements of the hierarchical matrix.

Details

See the examples for how to use.

Value

A hierarchical matrix of the structure as specified in x with off-diagonal entries as specified in x and diagonal entries as specified in diagonal.

Author(s)

Marius Hofert

Examples

rho <- c(0.2, 0.3, 0.5, 0.8) # some entries (e.g., correlations)

## Test homogeneous case
x <- list(rho[1], 1:6)
hierarchical_matrix(x)

## Two-level case with one block of size 2
x <- list(rho[1], 1, list(rho[2], 2:3))
hierarchical_matrix(x)

## Two-level case with one block of size 2 and a larger homogeneous block
x <- list(rho[1], 1:3, list(rho[2], 4:5))
hierarchical_matrix(x)

## Test two-level case with three blocks of size 2
x <- list(rho[1], NULL, list(list(rho[2], 1:2),
                             list(rho[3], 3:4),
                             list(rho[4], 5:6)))
hierarchical_matrix(x)

## Test three-level case
x <- list(rho[1], 1:3, list(rho[2], NULL, list(list(rho[3], 4:5),
                                               list(rho[4], 6:8))))
hierarchical_matrix(x)

## Test another three-level case
x <- list(rho[1], c(3, 6, 1), list(rho[2], c(9, 2, 7, 5),
                                   list(rho[3], c(8, 4))))
hierarchical_matrix(x)

---

Hill Estimator and Plot

Description

Compute the Hill estimator and Hill plot.

Usage

Hill_estimator(x, k = c(10, length(x)), conf.level = 0.95)
Hill_plot(x, k = c(10, length(x)), conf.level = 0.95, Hill.estimator = NULL,
          log = "x", xlim = NULL, ylim = NULL,
          xlab = "Order statistics", ylab = "Tail index",
          CI.col = adjustcolor(1, alpha.f = 0.2), lines.args = list(),
          xaxis2 = TRUE, xlab2 = "Empirical probability", ...)

Arguments

x	vector of numeric data.
k	vector of length 2, determining the smallest and largest number of order statistics of x to compute the Hill estimator for (the smallest needs to be >= 2). If k is of length 1, k is expanded by length(x).
conf.level	confidence level of the confidence intervals.
Hill.estimator	object as returned by Hill_estimator().
log, xlim, ylim, xlab, ylab	see plot().
CI.col	color of the pointwise asymptotic confidence intervals (CIs); if NA, no CIs are shown.
lines.args	list of additional arguments for the underlying lines() call to draw the Hill estimator.
xaxis2	logical indicating whether a third axis is drawn that shows the empirical probabilities 1-(k-1)/length(x) corresponding to k, so the value of the column k.prob as returned by Hill_estimator().
xlab2	label of the secondary x-axis.
...	additional arguments passed to the underlying plot().

Details

See McNeil et al. (2015, Section 5.2.4, (5.23))

Value

Hill_estimator():

A five-column matrix containing the indices k, their corresponding empirical probabilities k.prob, the estimated tail indices tail.index, and the lower and upper CI endpoints CI.low and CI.up.

Hill_plot():

Hill plot by side-effect.

Author(s)

Marius Hofert

References

McNeil, A. J., Frey, R. and Embrechts, P. (2015). Quantitative Risk Management: Concepts, Techniques, Tools. Princeton University Press.

Examples

set.seed(271)
X <- rt(1000, df = 3.5)
Y <- X[X > 0]
Hill_plot(Y)
Hill_plot(Y, log = "", CI.col = NA)

---

Density Plot of the Values from a Lower Triangular Matrix

Description

Density plot of all values in the lower triangular part of a matrix.

Usage

matrix_density_plot(x, xlab = "Entries in the lower triangular matrix",
                    main = "", text = NULL, side = 4, line = 1, adj = 0, ...)

Arguments

x	matrix-like object.
xlab	x-axis label.
main	title.
text	see mtext(). The text = "", it is omitted.
side	see mtext().
line	see mtext().
adj	see mtext().
...	additional arguments passed to the underlying plot().

Details

matrix_density_plot() is typically used for symmetric matrices (like correlation matrices, matrices of pairwise Kendall's tau or tail dependence parameters) to check the distribution of their off-diagonal entries.

Value

invisible().

Author(s)

Marius Hofert

Examples

## Generate a random correlation matrix
d <- 50
L <- diag(1:d)
set.seed(271)
L[lower.tri(L)] <- runif(choose(d,2))
Sigma <- L 
P <- cor(Sigma)
## Density of its lower triangular entries
matrix_density_plot(P)

---

Graphical Tool for Visualizing Matrices

Description

Plot of a matrix.

Usage

matrix_plot(x, ran = range(x, na.rm = TRUE), ylim = rev(c(0.5, nrow(x) + 0.5)),
            xlab = "Column", ylab = "Row",
            scales = list(alternating = c(1,1), tck = c(1,0),
                          x = list(at = pretty(1:ncol(x)), rot = 90),
                          y = list(at = pretty(1:nrow(x)))),
            at = NULL, colorkey = NULL, col = c("royalblue3", "white", "maroon3"),
            col.regions = NULL, ...)

Arguments

x	matrix-like object.
ran	range (can be used to enforce (-1,1), for example).
ylim	y-axis limits in reverse order (for the rows to appear 'top down').
xlab	x-axis label.
ylab	y-axis label.
scales	see levelplot(); if NULL, labels and ticks are omitted.
at	see levelplot(). If NULL, a useful default is computed based on the given values in x.
colorkey	see levelplot(). If NULL, a useful default is computed based on at.
col	vector of length two (if all values of x are non-positive or all are non-negative; note that also a vector of length three is allowed in this case) or three (if x contains negative and positive values) providing the color key's default colors.
col.regions	see levelplot(). If NULL, a useful default is computed based on at.
...	additional arguments passed to the underlying levelplot().

Details

Plot of a matrix.

Value

The plot, a Trellis object.

Author(s)

Marius Hofert

Examples

## Generate a random correlation matrix
d <- 50
L <- diag(1:d)
set.seed(271)
L[lower.tri(L)] <- runif(choose(d,2)) # random Cholesky factor
Sigma <- L 
P <- cor(Sigma)

## Default
matrix_plot(P)
matrix_plot(P, ran = c(-1, 1)) # within (-1, 1)
matrix_plot(abs(P)) # if nonnegative
L. <- L
diag(L.) <- NA
matrix_plot(L.) # Cholesky factor without diagonal

## Default if nonpositive
matrix_plot(-abs(P))

## Changing colors
matrix_plot(P, ran = c(-1, 1),
            col.regions = grey(c(seq(0, 1, length.out = 100),
                                 seq(1, 0, length.out = 100))))

## An example with overlaid lines
library(lattice)
my_panel <- function(...) {
    panel.levelplot(...)
    panel.abline(h = c(10, 20), v = c(10, 20), lty = 2)
}
matrix_plot(P, panel = my_panel)

---

Mean Excess

Description

Sample mean excess function, mean excess function of a GPD and sample mean excess plot.

Usage

mean_excess_np(x, omit = 3)
mean_excess_plot(x, omit = 3,
                 xlab = "Threshold", ylab = "Mean excess over threshold", ...)
mean_excess_GPD(x, shape, scale)

Arguments

x	mean_excess_GPD(): numeric vector of evaluation points of the mean excess function of the GPD. otherwise: numeric vector of data.
omit	number $\ge 1$ of unique last observations to be omitted from the sorted data (as mean excess plot becomes unreliable for these observations as thresholds).
xlab	x-axis label.
ylab	y-axis label.
...	additional arguments passed to the underlying plot().
shape	GPD shape parameter $\xi$.
scale	GPD scale parameter $\beta$.

Details

Mean excess plots can be used in the peaks-over-threshold method for choosing a threshold. To this end, one chooses the smallest threshold above which the mean excess plot is roughly linear.

Value

mean_excess_np() returns a two-column matrix giving the sorted data without the omit-largest unique values (first column) and the corresponding values of the sample mean excess function (second column). It is mainly used in mean_excess_plot().

mean_excess_plot() returns invisible().

mean_excess_GPD() returns the mean excess function of a generalized Pareto distribution evaluated at x.

Author(s)

Marius Hofert

Examples

## Generate losses to work with
set.seed(271)
X <- rt(1000, df = 3.5) # in MDA(H_{1/df}); see MFE (2015, Section 16.1.1)

## (Sample) mean excess plot and threshold choice
mean_excess_plot(X[X > 0]) # we only use positive values here to see 'more'
## => Any value in [0.8, 2] seems reasonable as threshold at first sight
##    but 0.8 to 1 turns out to be too small for the degrees of
##    freedom implied by the GPD estimator to be close to the true value 3.5.
## => We go with threshold 1.5 here.
u <- 1.5 # thresholds

## An alternative way
ME <- mean_excess_np(X[X > 0])
plot(ME, xlab = "Threshold", ylab = "Mean excess over threshold")

## Mean excess plot with mean excess function of the fitted GPD
fit <- fit_GPD_MLE(X[X > u] - u)
q <- seq(u, ME[nrow(ME),"x"], length.out = 129)
MEF.GPD <- mean_excess_GPD(q-u, shape = fit$par[["shape"]], scale = fit$par[["scale"]])
mean_excess_plot(X[X > 0]) # mean excess plot for positive losses...
lines(q, MEF.GPD, col = "royalblue", lwd = 1.4) # ... with mean excess function of the fitted GPD

---

Graphical Tool for Visualizing NAs in a Data Set

Description

Plot NAs in a data set.

Usage

NA_plot(x, col = c("black", "white"), xlab = "Time", ylab = "Component",
        text = "Black: NA; White: Available data",
        side = 4, line = 1, adj = 0, ...)

Arguments

x	matrix (ideally an xts object).
col	bivariate vector containing the colors for missing and available data, respectively.
xlab	x-axis label.
ylab	y-axis label.
text	see mtext(). The text = "", it is omitted.
side	see mtext().
line	see mtext().
adj	see mtext().
...	additional arguments passed to the underlying image().

Details

Indicate NAs in a data set.

Value

invisible().

Author(s)

Marius Hofert

Examples

## Generate data
n <- 1000 # sample size
d <- 100 # dimension
set.seed(271) # set seed
x <- matrix(runif(n*d), ncol = d) # generate data

## Assign missing data
k <- ceiling(d/4) # fraction of columns with some NAs
j <- sample(1:d, size = k) # columns j with NAs
i <- sample(1:n, size = k) # 1:i will be NA in each column j
X <- x
for(k. in seq_len(k)) X[1:i[k.], j[k.]] <- NA # put in NAs

## Plot NAs
NA_plot(X) # indicate NAs

---

P-P and Q-Q Plots

Description

Probability-probability plots and quantile-quantile plots.

Usage

pp_plot(x, FUN, pch = 20, xlab = "Theoretical probabilities",
        ylab = "Sample probabilities", ...)
qq_plot(x, FUN = qnorm, method = c("theoretical", "empirical"),
        pch = 20, do.qqline = TRUE, qqline.args = NULL,
        xlab = "Theoretical quantiles", ylab = "Sample quantiles",
        ...)

Arguments

x	data vector.
FUN	function. For pp_plot(): The distribution function (vectorized). qq_plot(): The quantile function (vectorized).
pch	plot symbol.
xlab	x-axis label.
ylab	y-axis label.
do.qqline	logical indicating whether a Q-Q line is plotted.
method	method used to construct the Q-Q line. If "theoretical", the theoretically true line with intercept 0 and slope 1 is displayed; if "empirical", the intercept and slope are determined with qqline(). The former helps deciding whether x comes from the distribution specified by FUN exactly, the latter whether x comes from a location-scale transformed distribution specified by FUN.
qqline.args	list containing additional arguments passed to the underlying abline() functions. Defaults to list(a = 0, b = 1) if method = "theoretical" and list() if method = "empirical".
...	additional arguments passed to the underlying plot().

Details

Note that Q-Q plots are more widely used than P-P plots (as they highlight deviations in the tails more clearly).

Value

invisible().

Author(s)

Marius Hofert

Examples

## Generate data
n <- 1000
mu <- 1
sig <- 3
nu <- 3.5
set.seed(271) # set seed
x <- mu + sig * sqrt((nu-2)/nu) * rt(n, df = nu) # sample from t_nu(mu, sig^2)

## P-P plot
pF <- function(q) pt((q - mu) / (sig * sqrt((nu-2)/nu)), df = nu)
pp_plot(x, FUN = pF)

## Q-Q plot
qF <- function(p) mu + sig * sqrt((nu-2)/nu) * qt(p, df = nu)
qq_plot(x, FUN = qF)

## A comparison with R's qqplot() and qqline()
qqplot(qF(ppoints(length(x))), x) # the same (except labels)
qqline(x, distribution = qF) # slightly different (since *estimated*)

## Difference of the two methods
set.seed(271)
z <- rnorm(1000)
## Standardized data
qq_plot(z, FUN = qnorm) # fine
qq_plot(z, FUN = qnorm, method = "empirical") # fine
## Location-scale transformed data
mu <- 3
sig <- 2
z. <- mu+sig*z
qq_plot(z., FUN = qnorm) # not fine (z. comes from N(mu, sig^2), not N(0,1))
qq_plot(z., FUN = qnorm, method = "empirical") # fine (as intercept and slope are estimated)

---

Computing Returns and Inverse Transformation

Description

Compute log-returns, simple returns and basic differences (or the inverse operations) from given data.

Usage

returns(x, method = c("logarithmic", "simple", "diff"), inverse = FALSE,
        start, start.date)
returns_qrmtools(x, method = c("logarithmic", "simple", "diff"),
                 inverse = FALSE, start, start.date)

Arguments

x	matrix or vector (possibly a xts object) to be turned into returns (if inverse = FALSE) or returns to be turned into the original data (if inverse = TRUE).
method	character string indicating the method to be used (log-returns (logarithmic changes), simple returns (relative changes), or basic differences). Note that this can also be a vector of such methods of length equal to the number of columns of x.
inverse	logical indicating whether the inverse transformation (data from given returns) shall be computed (if TRUE, this requires start to be specified).
start	if inverse = TRUE, the last available value of the time series to be constructed from the given returns x.
start.date	character or Date object to be used as the date corresponding to the value start; currently only used for xts objects.

Details

If inverse = FALSE and x is an xts object, the returned object is an xts, too.

Note that the R package timeSeries also contains a function returns() (and hence the order in which timeSeries and qrmtools are loaded matters to get the right returns()). For this reason, returns_qrmtools() is an alias for returns() from qrmtools.

Value

vector or matrix with the same number of columns as x just one row less if inverse = FALSE or one row more if inverse = TRUE.

Author(s)

Marius Hofert

Examples

## Generate two paths of a geometric Brownian motion
S0 <- 10 # current stock price S_0
r <- 0.01 # risk-free annual interest rate
sig <- 0.2 # (constant) annual volatility
T <- 2 # maturity in years
N <- 250 # business days per year
t <- 1:(N*T) # time points to be sampled
npath <- 2 # number of paths
set.seed(271) # for reproducibility
S <- replicate(npath, S0 * exp(cumsum(rnorm(N*T, # sample paths of S_t
                                            mean = (r-sig^2/2)/N,
                                            sd = sqrt((sig^2)/N))))) # (N*T, npath)

## Turn into xts objects
library(xts)
sdate <- as.Date("2000-05-02") # start date
S.  <- as.xts(S, order.by = seq(sdate, length.out = N*T, by = "1 week"))
plot(S.[,1], main = "Stock 1")
plot(S.[,2], main = "Stock 2")


### Log-returns ################################################################

## Based on S[,1]
X <- returns(S[,1]) # build log-returns (one element less than S)
Y <- returns(X, inverse = TRUE, start = S[1,1]) # transform back
stopifnot(all.equal(Y, S[,1]))

## Based on S
X <- returns(S) # build log-returns (one element less than S)
Y <- returns(X, inverse = TRUE, start = S[1,]) # transform back
stopifnot(all.equal(Y, S))

## Based on S.[,1]
X <- returns(S.[,1])
Y <- returns(X, inverse = TRUE, start = S.[1,1], start.date = sdate)
stopifnot(all.equal(Y, S.[,1], check.attributes = FALSE))

## Based on S.
X <- returns(S.)
Y <- returns(X, inverse = TRUE, start = S.[1], start.date = sdate)
stopifnot(all.equal(Y, S., check.attributes = FALSE))

## Sign-adjusted (negative) log-returns
X <- -returns(S) # build -log-returns
Y <- returns(-X, inverse = TRUE, start = S[1,]) # transform back
stopifnot(all.equal(Y, S))


### Simple returns #############################################################

## Simple returns based on S
X <- returns(S, method = "simple")
Y <- returns(X, method = "simple", inverse = TRUE, start = S[1,])
stopifnot(all.equal(Y, S))

## Simple returns based on S.
X <- returns(S., method = "simple")
Y <- returns(X, method = "simple", inverse = TRUE, start = S.[1,],
             start.date = sdate)
stopifnot(all.equal(Y, S., check.attributes = FALSE))

## Sign-adjusted (negative) simple returns
X <- -returns(S, method = "simple")
Y <- returns(-X, method = "simple", inverse = TRUE, start = S[1,])
stopifnot(all.equal(Y, S))


### Basic differences ##########################################################

## Basic differences based on S
X <- returns(S, method = "diff")
Y <- returns(X, method = "diff", inverse = TRUE, start = S[1,])
stopifnot(all.equal(Y, S))

## Basic differences based on S.
X <- returns(S., method = "diff")
Y <- returns(X, method = "diff", inverse = TRUE, start = S.[1,],
             start.date = sdate)
stopifnot(all.equal(Y, S., check.attributes = FALSE))

## Sign-adjusted (negative) basic differences
X <- -returns(S, method = "diff")
Y <- returns(-X, method = "diff", inverse = TRUE, start = S[1,])
stopifnot(all.equal(Y, S))


### Vector-case of 'method' ####################################################

X <- returns(S., method = c("logarithmic", "diff"))
Y <- returns(X, method = c("logarithmic", "diff"), inverse = TRUE, start = S.[1,],
             start.date = sdate)
stopifnot(all.equal(Y, S., check.attributes = FALSE))

---

Risk Measures

Description

Computing risk measures.

Usage

## Value-at-risk
VaR_np(x, level, names = FALSE, type = 1, ...)
VaR_t(level, loc = 0, scale = 1, df = Inf)
VaR_t01(level, df = Inf)
VaR_GPD(level, shape, scale)
VaR_Par(level, shape, scale = 1)
VaR_GPDtail(level, threshold, p.exceed, shape, scale)

## Expected shortfall
ES_np(x, level, method = c(">", ">="), verbose = FALSE, ...)
ES_t(level, loc = 0, scale = 1, df = Inf)
ES_t01(level, df = Inf)
ES_GPD(level, shape, scale)
ES_Par(level, shape, scale = 1)
ES_GPDtail(level, threshold, p.exceed, shape, scale)

## Range value-at-risk
RVaR_np(x, level, ...)

## Multivariate geometric value-at-risk and expectiles
gVaR(x, level, start = colMeans(x),
     method = if(length(level) == 1) "Brent" else "Nelder-Mead", ...)
gEX(x, level, start = colMeans(x),
    method = if(length(level) == 1) "Brent" else "Nelder-Mead", ...)

Arguments

x	gVaR(), gEX(): matrix of (rowwise) multivariate losses. VaR_np(), ES_np(), RVaR_np(): if x is a matrix then rowSums() is applied first (so value-at-risk and expected shortfall of the sum is computed). otherwise: vector of losses.
level	RVaR_np(): vector of length 1 or 2 giving the lower and upper confidence level; if of length 1, it is interpreted as the lower confidence level and the upper one is taken to be 1. gVaR(), gEX(): vector or matrix of (rowwise) confidence levels $\alpha$ (all in $[0,1]$). otherwise: confidence level $\alpha\in[0,1]$.
names	see ?quantile.
type	see ?quantile.
loc	location parameter $\mu$.
shape	VaR_GPD(), ES_GPD(): GPD shape parameter $\xi$, a real number. VaR_Par(), ES_Par(): Pareto shape parameter $\theta$, a positive number.
scale	VaR_t(), ES_t(): $t$ scale parameter $\sigma$, a positive number. VaR_GPD(), ES_GPD(): GPD scale parameter $\beta$, a positive number. VaR_Par(), ES_Par(): Pareto scale parameter $\kappa$, a positive number.
df	degrees of freedom, a positive number; choose df = Inf for the normal distribution. For the standardized $t$ distributions, df has to be greater than 2.
threshold	threhold $u$ (used to estimate the exceedance probability based on the data x).
p.exceed	exceedance probability; typically mean(x > threshold) for x being the data modeled with the peaks-over-threshold (POT) method.
start	vector of initial values for the underlying optim().
method	ES_np(): character string indicating the method for computing expected shortfall. gVaR(), gEX(): the optimization method passed to the underlying optim().
verbose	logical indicating whether verbose output is given (in case the mean is computed over (too) few observations).
...	VaR_np(): additional arguments passed to the underlying quantile(). ES_np(), RVaR_np(): additional arguments passed to the underlying VaR_np(). gVaR(), gEX(): additional arguments passed to the underlying optim().

Details

The distribution function of the Pareto distribution is given by

$F(x) = 1-(\kappa/(\kappa+x))^{\theta},\ x\ge 0,$

where $\theta > 0$, $\kappa > 0$.

Value

VaR_np(), ES_np(), RVaR_np() estimate value-at-risk, expected shortfall and range value-at-risk non-parametrically. For expected shortfall, if method = ">=" (method = ">", the default), losses greater than or equal to (strictly greater than) the nonparametric value-at-risk estimate are averaged; in the former case, there might be no such loss, in which case NaN is returned. For range value-at-risk, losses greater than the nonparametric VaR estimate at level level[1] and less than or equal to the nonparametric VaR estimate at level level[2] are averaged.

VaR_t(), ES_t() compute value-at-risk and expected shortfall for the $t$ (or normal) distribution. VaR_t01(), ES_t01() compute value-at-risk and expected shortfall for the standardized $t$ (or normal) distribution, so scaled $t$ distributions to have mean 0 and variance 1; note that they require a degrees of freedom parameter greater than 2.

VaR_GPD(), ES_GPD() compute value-at-risk and expected shortfall for the generalized Pareto distribution (GPD).

VaR_Par(), ES_Par() compute value-at-risk and expected shortfall for the Pareto distribution.

gVaR(), gEX() compute the multivariate geometric value-at-risk and expectiles suggested by Chaudhuri (1996) and Herrmann et al. (2018), respectively.

Author(s)

Marius Hofert

References

McNeil, A. J., Frey, R. and Embrechts, P. (2015). Quantitative Risk Management: Concepts, Techniques, Tools. Princeton University Press.

Chaudhuri, P. (1996). On a geometric notion of quantiles for multivariate data. Journal of the American Statistical Assosiation 91(434), 862–872.

Herrmann, K., Hofert, M. and Mailhot, M. (2018). Multivariate geometric expectiles. Scandinavian Actuarial Journal, 2018(7), 629–659.

Examples

### 1 Univariate measures ######################################################

## Generate some losses and (non-parametrically) estimate VaR_alpha and ES_alpha
set.seed(271)
L <- rlnorm(1000, meanlog = -1, sdlog = 2) # L ~ LN(mu, sig^2)
## Note: - meanlog = mean(log(L)) = mu, sdlog = sd(log(L)) = sig
##       - E(L) = exp(mu + (sig^2)/2), var(L) = (exp(sig^2)-1)*exp(2*mu + sig^2)
##         To obtain a sample with E(L) = a and var(L) = b, use:
##         mu = log(a)-log(1+b/a^2)/2 and sig = sqrt(log(1+b/a^2))
VaR_np(L, level = 0.99)
ES_np(L,  level = 0.99)

## Example 2.16 in McNeil, Frey, Embrechts (2015)
V <- 10000 # value of the portfolio today
sig <- 0.2/sqrt(250) # daily volatility (annualized volatility of 20%)
nu <- 4 # degrees of freedom for the t distribution
alpha <- seq(0.001, 0.999, length.out = 256) # confidence levels
VaRnorm <- VaR_t(alpha, scale = V*sig, df = Inf)
VaRt4 <- VaR_t(alpha, scale = V*sig*sqrt((nu-2)/nu), df = nu)
ESnorm <- ES_t(alpha, scale = V*sig, df = Inf)
ESt4 <- ES_t(alpha, scale = V*sig*sqrt((nu-2)/nu), df = nu)
ran <- range(VaRnorm, VaRt4, ESnorm, ESt4)
plot(alpha, VaRnorm, type = "l", ylim = ran, xlab = expression(alpha), ylab = "")
lines(alpha, VaRt4, col = "royalblue3")
lines(alpha, ESnorm, col = "darkorange2")
lines(alpha, ESt4, col = "maroon3")
legend("bottomright", bty = "n", lty = rep(1,4), col = c("black",
       "royalblue3", "darkorange3", "maroon3"),
       legend = c(expression(VaR[alpha]~~"for normal model"),
                  expression(VaR[alpha]~~"for "*t[4]*" model"),
                  expression(ES[alpha]~~"for normal model"),
                  expression(ES[alpha]~~"for "*t[4]*" model")))


### 2 Multivariate measures ####################################################

## Setup
library(copula)
n <- 1e4 # MC sample size
nu <- 3 # degrees of freedom
th <- iTau(tCopula(df = nu), tau = 0.5) # correlation parameter
cop <- tCopula(param = th, df = nu) # t copula
set.seed(271) # for reproducibility
U <- rCopula(n, cop = cop) # copula sample
theta <- c(2.5, 4) # marginal Pareto parameters
stopifnot(theta > 2) # need finite 2nd moments
X <- sapply(1:2, function(j) qPar(U[,j], shape = theta[j])) # generate X
N <- 17 # number of angles (rather small here because of run time)
phi <- seq(0, 2*pi, length.out = N) # angles
r <- 0.98 # radius
alpha <- r * cbind(alpha1 = cos(phi), alpha2 = sin(phi)) # vector of confidence levels

## Compute geometric value-at-risk
system.time(res <- gVaR(X, level = alpha))
gvar <- t(sapply(seq_len(nrow(alpha)), function(i) {
    x <- res[[i]]
    if(x[["convergence"]] != 0) # 0 = 'converged'
        warning("No convergence for alpha = (", alpha[i,1], ", ", alpha[i,2],
                ") (row ", i, ")")
    x[["par"]]
})) # (N, 2)-matrix

## Compute geometric expectiles
system.time(res <- gEX(X, level = alpha))
gex <- t(sapply(seq_len(nrow(alpha)), function(i) {
    x <- res[[i]]
    if(x[["convergence"]] != 0) # 0 = 'converged'
        warning("No convergence for alpha = (", alpha[i,1], ", ", alpha[i,2],
                ") (row ", i, ")")
    x[["par"]]
})) # (N, 2)-matrix

## Plot geometric VaR and geometric expectiles
plot(gvar, type = "b", xlab = "Component 1 of geometric VaRs and expectiles",
     ylab = "Component 2 of geometric VaRs and expectiles",
     main = "Multivariate geometric VaRs and expectiles")
lines(gex, type = "b", col = "royalblue3")
legend("bottomleft", lty = 1, bty = "n", col = c("black", "royalblue3"),
       legend = c("geom. VaR", "geom. expectile"))
lab <- substitute("MC sample size n ="~n.*","~t[nu.]~"copula with Par("*th1*
                  ") and Par("*th2*") margins",
                  list(n. = n, nu. = nu, th1 = theta[1], th2 = theta[2]))
mtext(lab, side = 4, line = 1, adj = 0)

---