Package 'FatTailsR'

Package FatTailsR

Description

This package includes Kiener distributions K1, K2, K3, K4 and K7 and two estimation algorithms to characterize with a high precision symmetric or asymmetric distributions with left and right fat tails that appear in market finance, neuroscience and many other disciplines. The estimation of the distribution parameters, quantiles, value-at-risk and expected shortfall is usually very accurate. Two datasets are provided, as well as power hyperbolas and power hyperbolic functions which are simplified versions of symmetric distribution K1.

Download the pdf cited in the references to get an overview of the theoretical part and several examples on stocks and indices.

A commercial package, FatTailsRplot, with advanced plotting functions and calculation of matrix of stocks over rolling windows is also developped by the author.

Details

With so many functions, this package could look fat. But it's not! It's rather agile and easy to use! The various functions included in this package can be assigned to the following groups:

1. Two datasets presented in different formats: list, data.frame, matrix, timeSeries, xts, zoo:

   • getDSdata.

   • extractData, dfData, mData, tData, xData, zData.

2. Functions to check the dimensions of vector, matrix, array, list:

   • dimdim, dimdim1, dimdimc.

3. Functions to calculate (positive, negative) prices to returns on vector, matrix, array, list, data.frame, timeSeries, xts, zoo:

   • elevate.

   • fatreturns, logreturns.

4. Several predefined vectors of probability. One function to check them. A conversion function from probabilities to characters

   • pprobs0, pprobs1, pprobs2, ..., pprobs9.

   • checkquantiles.

   • getnamesk.

5. Miscellaneous functions related to the logistic function:

   • logit, invlogit, ltmlogis, rtmlogis, eslogis.

6. Power hyperbolas, power hyperbolic functions and their reciprocal functions:

   • exphp, coshp, sinhp, tanhp, sechp, cosechp, cotanhp.

   • loghp, acoshp, asinhp, atanhp, asechp, acosechp, acotanhp.

   • kashp, dkashp_dx, ashp.

7. Logishp function, kogit and invkogit = logistic function + power hyperbolas:

   • d, p, q, r, dp, dq, l, dl, ql logishp.

   • kogit, invkogit.

8. Conversion functions between parameters related to Kiener distributions K1, K2, K3, K4:

   • aw2k, aw2d, aw2e, ad2e, ad2k, ad2w, ae2d, ae2k, ae2w, ak2e, ak2w, de2a, de2k, de2w, dk2a, dk2e, dw2a, dw2e, dw2k, ek2a, ak2d, ek2w, aw2a, aw2d, ew2a, aw2d, ew2k, kd2a, kd2e, kd2w, ke2a, ke2d, ke2w, kw2a, kw2d, kw2e.

   • pk2pk.

9. Kiener distributions K1, K2, K3, K4 and the new K7 (introduced in v1.7-0):

   • d, p, q, r, dp, dq, l, dl, ql, var, ltm, rtm, dtmq, es kiener1,

   • d, p, q, r, dp, dq, l, dl, ql, var, ltm, rtm, dtmq, es kiener2,

   • d, p, q, r, dp, dq, l, dl, ql, var, ltm, rtm, dtmq, es kiener3,

   • d, p, q, r, dp, dq, l, dl, ql, var, ltm, rtm, dtmq, es kiener4,

   • d, p, q, r, dp, dq, l, dl, ql, var, ltm, rtm, dtmq, es kiener7.

10. Quantile (VaR) corrective function (as a multiplier of the logistic function). Expected shortfall corrective function (as a multiplier of the expected shortfall of the logistic distribution):

    • ckiener1, ckiener2, ckiener3, ckiener4, ckiener7.

    • hkiener1, hkiener2, hkiener3, hkiener4, hkiener7.

11. Moments of the distribution estimated from the dataset and from the regression parameters:

    • xmoments.

    • kmoments, kmoment, kcmoment, kmean, kstandev, kvariance, kskewness, kkurtosis, kekurtosis.

12. Regression and estimation functions to estimate Kiener distribution parameters on a given dataset. *fit* and *param* are wrappers of algorithms reg and estim. reg uses an unweighted nonlinear regression function. estim uses a fast estimation based on quantiles:

    • regkienerLX, laplacegaussnorm.

    • fitkienerX.

    • paramkienerX, paramkienerX5, paramkienerX7.

13. Functions related to paramkienerX:

    • elevenprobs, sevenprobs, fiveprobs.

    • estimkiener11, estimkiener7, estimkiener5.

    • roundcoefk.

    • checkcoefk.

14. Predefined subsets of parameters to extract them from the long vector fitk obtained after regression/estimation regkienerLX, fitkienerX :

    • exfit0, ..., exfit7.

For a quick start, jump to the functions regkienerLX, fitkienerX and run the examples. Then, download and read the documents in pdf format cited in the references to get an overview on the major functions. Finally, explore the other examples.

Author(s)

Maintainer: Patrice Kiener [email protected] (ORCID)

References

P. Kiener, Explicit models for bilateral fat-tailed distributions and applications in finance with the package FatTailsR, 8th R/Rmetrics Workshop and Summer School, Paris, 27 June 2014. Download it from: https://www.inmodelia.com/exemples/2014-0627-Rmetrics-Kiener-en.pdf

P. Kiener, Fat tail analysis and package FatTailsR, 9th R/Rmetrics Workshop and Summer School, Zurich, 27 June 2015. Download it from: https://www.inmodelia.com/exemples/2015-0627-Rmetrics-Kiener-en.pdf

See Also

Useful links:

• https://www.inmodelia.com/fattailsr-en.html

Examples

require(graphics)
require(minpack.lm)
require(timeSeries)

### Load the datasets and select one number (1-16)
DS     <- getDSdata()
j      <- 5

### and run this block
X      <- DS[[j]]
nameX  <- names(DS)[j]
reg    <- regkienerLX(X)
lgn    <- laplacegaussnorm(X)
lleg   <- c("logit(0.999) = 6.9", "logit(0.99)   = 4.6", 
           "logit(0.95)   = 2.9", "logit(0.50)   = 0", 
           "logit(0.05)   = -2.9", "logit(0.01)   = -4.6", 
           "logit(0.001) = -6.9  ")
pleg   <- c( paste("m =",  reg$coefr4[1]), paste("g  =", reg$coefr4[2]), 
             paste("k  =", reg$coefr4[3]), paste("e  =", reg$coefr4[4]) )

## Main plot
op     <- par(mfrow = c(1,1), mgp = c(1.5,0.8,0), mar = c(3,3,2,1))
plot(reg$dfrXP, main = nameX)
legend("top", legend = pleg, cex = 0.9, inset = 0.02 )
lines(reg$dfrEP, col = 2, lwd = 2)
points(reg$dfrQkPk, pch = 3, col = 2, lwd = 2, cex = 1.5)
lines(lgn$dfrXPn, col = 7, lwd = 2)

## Plot F(X) > 0,97
front = c(0.06, 0.39, 0.50, 0.95)
par(fig = front, new = TRUE, mgp = c(1.5, 0.6, 0), las = 0)
plot( reg$dfrXP[which(reg$dfrXP$P > 0.97),] , pch = 1, xlab = "", ylab = "", main = "F(X) > 0,97" )
lines(reg$dfrEP[which(reg$dfrEP$P > 0.97),], type="l", col = 2, lwd = 3 )
lines(lgn$dfrXPn[which(lgn$dfrXPn$Pn > 0.97),], type = "l", col = 7, lwd= 2 )
points(reg$dfrQkPk, pch = 3, col = 2, lwd = 2, cex = 1.5)
points(lgn$dfrQnPn, pch = 3, col = 7, lwd = 2, cex = 1)

## Plot F(X) < 0,03
front = c(0.58, 0.98, 0.06, 0.61)
par(fig = front, new = TRUE, mgp = c(0.5, 0.6, 0), las = 0 )
plot( reg$dfrXP[which(reg$dfrXP$P < 0.03),] , pch = 1, xlab = "", ylab = "", main = "F(X) < 0,03")
lines(reg$dfrEP[which(reg$dfrEP$P < 0.03),], type = "l", col = 2, lwd = 3 )
lines(lgn$dfrXPn[which(lgn$dfrXPn$Pn < 0.03),], type = "l", col= 7, lwd= 2 )
points(reg$dfrQkPk, pch = 3, col = 2, lwd = 2, cex = 1.5)
points(lgn$dfrQnPn, pch = 3, col = 7, lwd = 2, cex = 1)

## Moments from the parameters (k) and from the Dataset (X)
round(cbind("k" = kmoments(reg$coefk, lengthx = nrow(reg$dfrXL)), "X" = xmoments(X)), 2)
attributes(reg)
### End block

---

Local Conversion Functions Between Kiener Distribution Parameters

Description

Conversion functions between parameters a, k, w, d, e used in Kiener distributions K2, K3 and K4.

Usage

aw2k(a, w)

aw2d(a, w)

aw2e(a, w)

ad2e(a, d)

ad2k(a, d)

ad2w(a, d)

ae2d(a, e)

ae2k(a, e)

ae2w(a, e)

ak2d(a, k)

ak2e(a, k)

ak2w(a, k)

de2a(d, e)

de2k(d, e)

de2w(d, e)

dk2a(d, k)

dk2e(d, k)

dk2w(d, k)

dw2a(d, w)

dw2e(d, w)

dw2k(d, w)

ek2a(e, k)

ek2d(e, k)

ek2w(e, k)

ew2a(e, w)

ew2d(e, w)

ew2k(e, w)

kd2a(k, d)

kd2e(k, d)

kd2w(k, d)

ke2a(k, e)

ke2d(k, e)

ke2w(k, e)

kw2a(k, w)

kw2d(k, w)

kw2e(k, w)

Arguments

a	a numeric value.
w	a numeric value.
d	a numeric value.
e	a numeric value.
k	a numeric value.

Details

a (alpha) is the left tail parameter, w (omega) is the right tail parameter, d (delta) is the distortion parameter, e (epsilon) is the eccentricity parameter. k (kappa) is the harmonic mean of a and w and describes a global tail parameter. They are defined by:

$aw2k(a, w) = k = 2 / (1/a + 1/w) = \frac{2}{\frac{1}{a} +\frac{1}{w}}$

$aw2d(a, w) = d = (-1/a + 1/w) / 2 = \frac{-\frac{1}{a} +\frac{1}{w}}{2}$

$aw2e(a, w) = e = (a - w) / (a + w) = \frac{a-w}{a+w}$

$kd2a(k, d) = a = 1 / ( 1/k - d) = \frac{1}{\frac{1}{k} - d}$

$kd2w(k, d) = w = 1 / ( 1/k + d) = \frac{1}{\frac{1}{k} + d}$

$ke2a(k, e) = a = k / (1 - e) = \frac{k}{1-e}$

$ke2w(k, e) = w = k / (1 + e) = \frac{k}{1+e}$

$ke2d(k, e) = d = e / k = \frac{e}{k}$

$kd2e(k, d) = e = k * d$

$de2k(k, e) = k = e / d = \frac{e}{d}$

See Also

The asymmetric Kiener distributions K2, K3, K4: kiener2, kiener3, kiener4

Examples

aw2k(4, 6); aw2d(4, 6); aw2e(4, 6)
outer(1:6, 1:6, aw2k)

---

Check Coefk

Description

Check that coefk is either a vector of length 7 or a matrix with 7 columns or an array with length of last dimension equal to 7.

Usage

checkcoefk(coefk, dim = c(1, 2), STOP = TRUE)

Arguments

coefk	numeric, matrix or data.frame representing parameters c(m,g,a,k,w,d,e).
dim	numeric. Accepted dimension(s) for coefk: 1 for vector, 2 for matrix, 3 for array. List is not accepted. Default is c(1, 2).
STOP	boolean. If an error is encountered, TRUE stops the function and returns an error message. FALSE just returns FALSE.

Examples

(coefk <- paramkienerX(getDSdata()))
checkcoefk(coefk)
checkcoefk(t(coefk), STOP = FALSE)

---

Check Quantiles and Probabilities

Description

Check that quantiles (or probabilities) are all different from each other and correctly ordered. If proba = TRUE, check that values are in range (0, 1).

Usage

checkquantiles(x, proba = FALSE, acceptNA = FALSE, STOP = TRUE)

Arguments

x	vector of quantiles.
proba	boolean. If TRUE, check range (0,1).
acceptNA	boolean. If FALSE, NA value are not accepted.
STOP	boolean. If an error is encountered, TRUE stops the function and returns an error message. FALSE just returns FALSE.

Examples

lst <- list(
  0.8,
  c(0.1, 0.5, 0.8),
  c(0.1, 0.5, 0.8, 0.2),
  c(2, 3, 1),
  c(2, 3),
  -0.01,
  NA,
  c(NA, NA),
  c(0.1, NA),
  c(0.1, NA, 0.5, 0.8),
  c(0.1, NA, 0.8, NA, 0.5),
  c(12, NA)
)

## Evaluate
for (i in seq_along(lst)) { 
  cat(i, lst[[i]], " : ",
      checkquantiles(lst[[i]], proba = FALSE, STOP = FALSE), 
      checkquantiles(lst[[i]], proba = TRUE, STOP = FALSE), 
      checkquantiles(lst[[i]], proba = FALSE, acceptNA = TRUE, STOP = FALSE), 
      checkquantiles(lst[[i]], proba = TRUE,  acceptNA = TRUE, STOP = FALSE), 
	     "\n") 
}

sapply(lst, checkquantiles, proba = TRUE, acceptNA = TRUE, STOP = FALSE)

## Not run: 
checkquantiles(matrix((1:12)/16, ncol=3), proba = TRUE, STOP = FALSE)
## End(Not run)

---

Quantile (VaR) and Expected Shortfall Corrective Functions

Description

Quantile functions (or VaR) and Expected Shortfall of Kiener distributions K1, K2, K3 and K4, usually calculated at pprobs2 = c(0.01, 0.025, 0.05, 0.95, 0.975, 0.99), can be expressed as:

1. Quantile of the logit function multiplied by a fat tail (c)orrective function ckiener1234;

2. Expected s(h)ortfall of the logistic function multiplied by a corrective function hkiener1234.

Both functions ckiener1234 and hkiener1234 are independant from the scale parameter g and are indirect measures of the tail curvature. A value close to 1 indicates a model similar to the logistic function with almost no curvature and probably parameter k > 8. When k (or a,w) decreases, the values of c and h increase and indicate some more pronounced symmetric or asymmetric curvature, depending on values of d,e. Note that if (min(a,k,w) <= 1), ckiener1234 still exists but the expected shortfall and hkiener1234 become undefined (NA).

Some financial applications use threshold values on ckiener1234 or hkiener1234 to select or discard stocks over time as they become less or more risky.

Usage

hkiener1(p, m = 0, g = 1, k = 3.2, lower.tail = TRUE, log.p = FALSE)

hkiener2(p, m = 0, g = 1, a = 3.2, w = 3.2, lower.tail = TRUE,
  log.p = FALSE)

hkiener3(p, m = 0, g = 1, k = 3.2, d = 0, lower.tail = TRUE,
  log.p = FALSE)

hkiener4(p, m = 0, g = 1, k = 3.2, e = 0, lower.tail = TRUE,
  log.p = FALSE)

hkiener7(p, coefk = c(0, 1, 3.2, 3.2, 3.2, 0, 0), lower.tail = TRUE,
  log.p = FALSE)

ckiener1(p, k = 3.2, lower.tail = TRUE, log.p = FALSE)

ckiener2(p, a = 3.2, w = 3.2, lower.tail = TRUE, log.p = FALSE)

ckiener3(p, k = 3.2, d = 0, lower.tail = TRUE, log.p = FALSE)

ckiener4(p, k = 3.2, e = 0, lower.tail = TRUE, log.p = FALSE)

ckiener7(p, coefk = c(0, 1, 3.2, 3.2, 3.2, 0, 0), lower.tail = TRUE,
  log.p = FALSE)

Arguments

p	numeric or vector of probabilities.
m	numeric. parameter m used in model K1, K2, K3 and K4.
g	numeric. parameter g used in model K1, K2, K3 and K4.
k	numeric. parameter k used in model K1, K3 and K4.
lower.tail	logical. If TRUE, use p. If FALSE, use 1-p.
log.p	logical. If TRUE, probabilities p are given as log(p).
a	numeric. parameter a used in model K2.
w	numeric. parameter w used in model K2.
d	numeric. parameter d used in model K3.
e	numeric. parameter e used in model K4.
coefk	vector or 7 columns-matrix representing parameters c(m,g,a,k,w,d,e) obtained from paramkienerX.

See Also

logit, qkiener1, qkiener2, qkiener3, qkiener4, fitkienerX.

---

Datasets dfData, mData, tData, xData, zData, extractData : dfData

Description

A list of datasets in data.frame, matrix, timeSeries, xts and zoo formats. This is the data.frame format. Visit extractData for more information.

---

Length and Dimensions of Vector, Matrix, Array, Data.Frame, List

Description

Dimensions and length of vector, matrix, array, data.frame and list. A friendly version of dim that returns the true dimension rather than the sometimes unexpected NULL value. The number of dimensions appears first, then the length in each dimension. A special case is list: the list's length (number of items) is turned into a negative integer and the dimension/length of each item is either positive if the item is a vector, matrix, array or data.frame or negative if the item is itself a list. Only the first level of the list is explored.

dimdim1 and dimdimc return the first item of dimdim, thus the true dimension, either as an integer or as a character and, in this latest case, always "-1" for lists.

Notes: From version 1.6.2 (April 2016), dimdim(NULL) = c(0, 0). (before c(1, 0)). Hence, dimdim1(NULL) = 0 and dimdimc(NULL) = "0". Some problems may occur with S4 objects like dimdim(qualityTools::fracDesign(k = 3, gen = "C = AB")).

Usage

dimdim(x)

dimdim1(x)

dimdimc(x)

Arguments

x	vector, matrix, array, data.frame, list.

Examples

require(timeSeries)

dimdim(NULL) 
dimdim(NA); dimdim(NaN); dimdim(Inf); dimdim(TRUE); dimdim(FALSE)
dimdim(11:39)
dimdim(LETTERS[1:8])
dimdim(matrix(1:60, ncol=5))
dimdim(extractData())
dimdim(as.data.frame(extractData()))
dimdim(data.frame(X=1:2, Y=1:4, Z=LETTERS[1:8]))
dimdim(array(1:240, c(8,6,5)))
dimdim(array(1:240, c(4,2,6,5)))
dimdim(getDSdata())
dimdim(zData)
dimdim(xData)
dimdim(tData)

dimdim1(matrix(1:60, ncol=5))
dimdimc(matrix(1:60, ncol=5))
dimdim1(tData)
dimdimc(tData)

---

Elevate

Description

A transformation to turn negative prices into positive prices and maintain at the same time the hierachy between all prices.

Usage

elevate(X, e = NULL)

Arguments

X	The prices.
e	numeric. The focal point of the hyperbola.

Details

Negative prices in financial markets, like interest rates in Europe, are a nightmare as the rough calculation of the returns generates non-sense values. elevate uses an hyperbola and implements the following formula:

$elevate(x, e) = (x + sqrt(x*x + e*e)) / 2$

There is currently no rule of thumb to calculate e. When $e = NULL$, there is no change and the output is identical to the input. When $e = 0$, all negative values are turned to 0.

Examples

require(graphics)

X <- (-50:100)/5
plot( X, elevate(X, e = 5), type = "l", ylim = c(0, 20) )
lines(X, elevate(X, e = 2),   col = 2)
lines(X, elevate(X, e = 1),   col = 3)
lines(X, elevate(X, e = 0.5), col = 4)
lines(X, elevate(X, e = 0),   col = 1)

---

Eleven, Seven, Five Probabilities

Description

Extract from a dataset X a vector of 11, 7 or 5 probabilities:

• c(p1, p2, p3, 0.25, 0.35, 0.50, 0.65, 0.75, 1-p3, 1-p2, 1-p1)

• c(p1, p2, 0.25, 0.50, 0.75, 1-p2, 1-p1)

• c(p1, 0.25, 0.50, 0.75, 1-p1)

where p1, p2 and p3 are the most extreme probabilities with values finishing by ..01, ..025 or ..05 that can be extracted from the dataset X. Parameters names are displayed if parnames = TRUE.

From version 1.8-0, p1 and 1-p1 can be associated to the i-th and (N-i)-th element.

Usage

elevenprobs(X, parnames = FALSE)

sevenprobs(X, parnames = FALSE)

fiveprobs(X, i = 4, parnames = FALSE)

Arguments

X	numeric. Vector of quantiles.
parnames	boolean. Output parameter vector with or without names.
i	integer. The i-th and (N-i)-th elements for which the probabilities p1 and 1-p1 are calculated. If (i == 0), the method used before version 1.8-0 : the extreme finishing by ..01, ..025 or ..05.

See Also

fitkienerX, estimkiener11.

Examples

require(timeSeries)

## DS
DS  <- getDSdata()
for (j in 1:16) { print(round(elevenprobs(DS[[j]]), 6)) }
z   <- cbind(t(sapply(DS, elevenprobs)), sapply(DS, length))
colnames(z) <- c("p1","p2","p3","p.25","p.35","p.50","p.65","p.75","1-p3","1-p2","1-p1","length")
z

## Choose j in 1:16
j   <- 1
X   <- sort(DS[[j]])
leX <- logit(eX <- elevenprobs(X))
lpX <- logit(ppoints(length(X), a = 0))
plot(X, lpX)
abline(h = leX, lty = 3)
mtext(eX, side = 4, at = leX, las = 1, line = -3.3)

---

Estimation Functions with 5, 7 or 11 Quantiles

Description

Several functions to estimate the parameters of asymmetric Kiener distributions with just 5, 7 or 11 quantiles.

Usage

estimkiener11(x11, p11, ord = 7, maxk = 10)

estimkiener7(x7, p7, maxk = 10)

estimkiener5(x5, p5, maxk = 20, maxe = 0.9)

Arguments

ord	integer. Option for probability selection and treatment.
maxk	numeric. Maximum value for k (kappa).
x5, x7, x11	vector of 5, 7 or 11 quantiles.
p5, p7, p11	vector of 5, 7 or 11 probabilities.
maxe	numeric. Maximum value for abs(e) (epsilon). Maximum is maxe = 1.

Details

These functions, called by paramkienerX5, paramkienerX7, paramkienerX, use 5, 7 or 11 probabilites and quantiles to estimate the parameters of Kiener distributions.

p5, x5 are obtained with functions fiveprobs(X) and quantile(p5).

p7, x7 are obtained with functions sevenprobs(X) and quantile(p7).

p11, x11 are obtained with functions elevenprobs(X) and quantile(p11).

The extraction of the 11 probabilities is controlled with the option ord which can take 12 integer values, ord = 7 being the default. Small dataset should consider ord = 5 and large dataset can consider ord = 12:

1. c(p1, 0.35, 0.50, 0.65, 1-p1)

2. c(p2, 0.35, 0.50, 0.65, 1-p2)

3. c(p1, p2, 0.35, 0.50, 0.65, 1-p2, 1-p1)

4. c(p1, p2, p3, 0.35, 0.50, 0.65, 1-p3, 1-p2, 1-p1)

5. c(p1, 0.25, 0.50, 0.75, 1-p1)

6. c(p2, 0.25, 0.50, 0.75, 1-p2)

7. c(p1, p2, 0.25, 0.50, 0.75, 1-p2, 1-p1)

8. c(p1, p2, p3, 0.25, 0.50, 0.75, 1-p3, 1-p2, 1-p1)

9. c(p1, 0.25, 0.35, 0.50, 0.65, 0.75, 1-p1)

10. c(p2, 0.25, 0.35, 0.50, 0.65, 0.75, 1-p2)

11. c(p1, p2, 0.25, 0.35, 0.50, 0.65, 0.75, 1-p2, 1-p1)

12. c(p1, p2, p3, 0.25, 0.35, 0.50, 0.65, 0.75, 1-p3, 1-p2, 1-p1)

p5 = fiveprobs(X) corresponds to c(p1, 0.25, 0.50, 0.75, 1-p1).

p7 = sevenprobs(X) corresponds to c(p1, p2, 0.25, 0.50, 0.75, 1-p2, 1-p1).

The above probabilities are then transfered to the quantile function whose parameter type can change significantly the extracted quantiles. Our experience is that type = 6 is appropriate when k > 1.9 and type = 5 is appropriate when k < 1.9. Other types type = 8 and type = 9 can be considered as well. The other types should be ignored. (Note: when k < 1.5, algorithm algo = "reg" returns better results).

Parameter maxk controls the maximum allowed value for estimated parameter k. Reasonnable values are maxk = 10, 15, 20. Default is maxk = 10 to be consistent with regkienerLX.

See Also

elevenprobs, paramkienerX, quantile, roundcoefk.

Examples

require(timeSeries)

## Choose j in 1:16. Choose ord in 1:12 (7 is default)
j    <- 5
ord  <- 5
DS   <- getDSdata()
p11  <- elevenprobs(DS[[j]])
x11  <- quantile(DS[[j]], probs = p11, na.rm = TRUE, names = TRUE, type = 6) 
round(estimkiener11(x11, p11, ord), 3)

## Compare the results obtained with the 12 different values of ord on stock j
compare <- function(ord, x11, p11) {estimkiener11(x11, p11, ord)}
coefk   <- t(sapply(1:12, compare, x11, p11)) 
rownames(coefk) <- 1:12
mcoefk  <- apply(coefk, 2, mean) # the mean of the 12 results above
roundcoefk(rbind(coefk, mcoefk), 13)

---

Parameter Subsets

Description

Some vectors of parameter names to be used with parameter exfitk in functions regkienerLX(.., exfitk = ...) and fitkienerX(.., exfitk = ...) or to subset the vector (or matrix) fitk obtained after regression fitk <- regkienerLX(..)$fitk or estimation fitk <- fitkienerX(..). Visit fitkienerX for details on each parameter.

exfit0 <- c("lh", "ret")

exfit1 <- c("m", "g", "a", "k", "w", "d", "e")

exfit2 <- c("m1", "sd", "sk", "ke", "m1x", "sdx", "skx", "kex")

exfit3 <- c("q.01", "q.05", "q.95", "q.99", "ltm.025", "rtm.975")

exfit4 <- c("VaR.01", "VaR.05", "VaR.95", "VaR.99", "ES.025", "ES.975")

exfit5 <- c("c.01", "c.05", "c.95", "c.99", "h.025", "h.975")

exfit6 <- c(exfit1, exfit2, exfit3, exfit4, exfit5)

exfit7 <- c(exfit0, exfit1, exfit2, exfit3, exfit4, exfit5)

Usage

exfit0

exfit1

exfit2

exfit3

exfit4

exfit5

exfit6

exfit7

Format

An object of class character of length 2.

An object of class character of length 7.

An object of class character of length 8.

An object of class character of length 6.

An object of class character of length 6.

An object of class character of length 6.

An object of class character of length 33.

An object of class character of length 35.

Examples

require(minpack.lm)
require(timeSeries)

### Load the datasets and select one number j in 1:16
j      <- 5
DS     <- getDSdata()
(fitk  <- regkienerLX(DS[[j]])$fitk)
fitk[exfit3]
fitkienerX(DS[[j]], exfitk = exfit3)

---

Power Hyperbolas and Power Hyperbolic Functions

Description

These functions define the power hyperbola exphp and the associated power hyperbolic cosine, sine, tangent, secant, cosecant, cotangent. They are similar to the traditional hyperbolic functions with term x receiving a nonlinear transformation via the function kashp.

Usage

exphp(x, k = 1)

coshp(x, k = 1)

sinhp(x, k = 1)

tanhp(x, k = 1)

sechp(x, k = 1)

cosechp(x, k = 1)

cotanhp(x, k = 1)

Arguments

x	a numeric value, vector or matrix.
k	a numeric value, preferably strictly positive.

Details

exphp function is defined for x in (-Inf, +Inf) by:

$exphp(x, k) = exp( kashp(x, k) ) = exp( k * asinh(x / 2 / k) )$

coshp function is defined for x in (-Inf, +Inf) by:

$coshp(x, k) = cosh( kashp(x, k) )$

sinhp function is defined for x in (-Inf, +Inf) by:

$sinhp(x, k) = sinh( kashp(x, k) )$

tanhp function is defined for x in (-Inf, +Inf) by:

$tanhp(x, k) = tanh( kashp(x, k) )$

sechp function is defined for x in (-Inf, +Inf) by:

$sechp(x, k) = 1 / coshp(x, k)$

cosechp function is defined for x in (-Inf, 0) U (0, +Inf) by:

$cosechp(x, k) = 1 / sinhp(x, k)$

cotanhp function is defined for x in (-Inf, 0) U (0, +Inf) by:

$cotanhp(x, k) = 1 / tanhp(x, k)$

The undesired case k = 0 returns 0 for sinhp and tanhp, 1 for exphp, coshp and sechp, Inf for cosechp and cotanhp.

If k is a vector of length > 1, then the use of the function outer is recommanded.

See Also

The nonlinear transformation kashp, the inverse power hyperbolas and the inverse power hyperbolic functions loghp.

Examples

### Example 1
x  <- (-3:3)*3 
exphp(x, k = 4)
coshp(x, k = 4)
sinhp(x, k = 4) 
tanhp(x, k = 4)

### Example 2 outer + plot(exphp, coshp, sinhp, tanhp)
xmin  <- -10
xd    <- 0.5
x     <- seq(xmin, -xmin, xd) ; names(x) <- x
k     <- c(0.6, 1, 1.5, 2, 3.2, 10) ; names(k) <- k
olty  <- c(2, 1, 2, 1, 2, 1, 1)
olwd  <- c(1, 1, 2, 2, 3, 4, 2)
ocol  <- c(2, 2, 4, 4, 3, 3, 1)
op    <- par(mfrow = c(2,2), mgp = c(1.5,0.8,0), mar = c(3,3,2,1))

## exphp(x, k)
Texphp <- ts(cbind(outer(-x, k, exphp), "exp(-x/2)" = exp(-x/2)), 
             start = xmin, deltat = xd)
plot(Texphp, plot.type = "single", ylim = c(0,20), 
       lty = olty, lwd = olwd, col = ocol, xaxs = "i", yaxs = "i", xlab = "", 
       ylab = "", main = "exphp(-x, k)" )
legend("topright", title = expression(kappa), legend = colnames(Texphp), 
       inset = 0.02, lty = olty, lwd = olwd, col = ocol, cex = 0.7 )

## coshp(x, k)
Tcoshp <- ts(cbind(outer(x, k, coshp), "cosh(x/2)" = cosh(x/2)), 
             start = xmin, deltat = xd)
plot(Tcoshp, plot.type = "single", ylim = c(0,20), 
       lty = olty, lwd = olwd, col = ocol, xaxs = "i", yaxs = "i", 
       xlab = "", ylab = "", main = "coshp(x, k)" )
legend("top", title = expression(kappa), legend = colnames(Tcoshp), 
       inset = 0.02, lty = olty, lwd = olwd, col = ocol, cex = 0.7 )

## sinhp(x, k)
Tsinhp <- ts(cbind(outer(x, k, sinhp), "sinh(x/2)" = sinh(x/2)), 
             start = xmin, deltat=xd)
plot(Tsinhp, plot.type = "single", ylim = c(-10,10), 
       lty = olty, lwd = olwd, col = ocol, xaxs = "i", yaxs = "i", 
       xlab = "", ylab = "", main = "sinhp(x, k)" )
legend("topleft", title = expression(kappa), legend = colnames(Tsinhp), 
       inset = 0.02, lty = olty, lwd=  olwd, col = ocol, cex = 0.7 )

## tanhp(x, k)
Ttanhp <- ts(cbind(outer(x, k, tanhp), "tanh(x/2)" = tanh(x/2)), 
             start = xmin, deltat = xd)
plot(Ttanhp, plot.type = "single", ylim = c(-1,1), 
       lty = olty, lwd = olwd, col = ocol, xaxs = "i", yaxs = "i", xlab = "", 
       ylab = "", main = "tanhp(x, k)" )
legend("topleft", title = expression(kappa), legend = colnames(Ttanhp), 
       inset = 0.02, lty = olty, lwd = olwd, col = ocol, cex = 0.7 )
### End Example 3

---

Datasets dfData, mData, tData, xData, zData, extractData : extractData

Description

dfData, mData, tData, xData, zData are datasets made of lists of data.frame, matrix, timeSeries, xts and zoo components. They describe prices and returns of 10 financial series used in the documents and demos presented at 8th and 9th R/Rmetrics conferences (2014, 2015). See the references. The last serie (CHF, interest rates in Switzerland) exhibits negative prices. All distributions of logreturns exhibit fat tails. Function extractData converts subsets of mData, tData, xData, zData.

Usage

extractData(pr = "p", ft = "tss", start = "2007-01-01",
  end = "2013-12-31")

Arguments

pr	character. Extract prices or returns: c("p","r","prices","returns").
ft	character. Output format among c("tss","xts","zoo","dfr","bfr","mat").
start	character. Start date.
end	character. End date.

Examples

library(zoo) 
library(xts) 
library(timeSeries) 

### dfData, tData, xData, zData : prices only
attributes(dfData); attributes(tData); attributes(xData); attributes(zData) 
lapply(dfData, head, 3)
lapply( mData, head, 3)
lapply( tData, head, 3)
lapply( xData, head, 3)
lapply( zData, head, 3)

### extractData : prices and logreturns
head(ptD <- extractData("p", "tss", "2009-01-01", "2012-12-31")) ; tail(ptD)
head(rtD <- extractData("r", "tss")) 
head(pxD <- extractData("p", "xts")) 
head(rxD <- extractData("r", "xts")) 
head(pzD <- extractData("p", "zoo")) 
head(rzD <- extractData("r", "zoo")) 
head(pbD <- extractData("p", "bfr")) 
head(rbD <- extractData("r", "bfr")) 
head(pmD <- extractData("p", "mat")) 
head(rmD <- extractData("r", "mat")) 

### Remove item CHF (negative prices) from dfData, tData, xData, zData
Z <- dfData[names(dfData)[1:9]]; attributes(Z)
Z <- tData[names(tData)[1:9]]; attributes(Z)
Z <- xData[names(xData)[1:9]]; attributes(Z)
Z <- zData[names(zData)[1:9]]; attributes(Z)

---

Simple and Elaborated Prices to Returns

Description

fatreturns is an elaborated function to compute prices to returns. It includes a pre-treatment for negative prices. It computes either log-returns (default) or percentage-returns. It handles properly NA values in the input vector, replacing them by 0 in the output vector. Doing so, it warrants that the sum of the log-returns (when selected) is equal to the difference of the log-prices. It works with vector, matrix, data.frame, timeSeries, xts, zoo, list, list of lists and even list of vector, data.frame, timeSeries, xts, zoo mixed together. The returned object is of same dimension and same class than the input object with the first line filled with 0. The results may be as per one, per cent (default), per thousand and per ten thousand.

logreturns is an improved version of function 100*diff(log(x)) to handle vector, matrix, data.frame and list. It handles properly the first line and the NA values. It does not control time, rownames and colnames but may return them.

Usage

fatreturns(x, log = TRUE, per = "cent", e = NULL, dfrcol = 1,
  na.rm = TRUE)

logreturns(x)

replaceNA(x)

Arguments

x	The prices (vector, data.frame, matrix, timeSeries, xts, zoo, list).
log	boolean. log returns or percentage returns.
per	character. Either "one", "cent, "thousand", "tenthousand" or "o", "c", "th", "te". Multiply the result by 1, 100, 1000, 10000.
e	NULL or positive numeric. NULL is for no change f(x)=x. A positive numeric designates the focal point of the hyperbola to turn negative prices into positive prices, keeping the hierarchy: f(x)=(x+sqrt(x*x+e*e))/2. There is currently no rules of thumb for the optimal value of e.
dfrcol	integer. For data.frame only, designates the column that handles the time and must be processed separately. Use dfrcol = 0 if all columns must be processed and there is no time (or turn the data.frame to a matrix).
na.rm	boolean. Replace x[t]=NA with the previous non-NA value available in the price serie such that (x[t-1], x[t]=x[t-1], x[t+1]) and calculate the returns accordingly. Force 0 in the first line of the returns if x[1]=NA.

Details

Negative prices in financial markets, like interest rates in Europe, are a nightmare as the rough calculation of the returns generates non-sense values. elevate uses an hyperbola and implements the following formula:

$elevate(x, e) = (x + sqrt(x*x + e*e)) / 2$

There is currently no rule of thumb to calculate e. When $e = NULL$, there is no change and the output is identical to the input. When $e = 0$, all negative values are turned to 0.

Examples

fatreturns(extractData())
logreturns(getDSdata())

---

Estimation and Regression Functions for Kiener Distributions

Description

Several functions to estimate the parameters of asymmetric Kiener distributions and display the results in a numeric vector or in a matrix. Algorithm "reg" (the default) uses a nonlinear regression and handle difficult cases. Algorithm "estim" has been completely rewritten in version 1.8-0 and is now very accurate, even for k<1. Adjustement on extreme quantiles can be controlled very precisely.

Usage

fitkienerX(X, algo = c("r", "reg", "e", "estim"), ord = 7, maxk = 10,
  mink = 1.53, maxe = 0.5, probak = pprobs2, dgts = NULL,
  exfitk = NULL, dimnames = FALSE, ncores = 1)

paramkienerX(X, algo = c("r", "reg", "e", "estim"), ord = 7, maxk = 10,
  mink = 1.53, maxe = 0.5, dgts = 3, parnames = TRUE,
  dimnames = FALSE, ncores = 1)

paramkienerX7(X, dgts = 3, n = 10, maxk = 20, maxe = 0.9,
  parnames = TRUE, dimnames = FALSE, ncores = 1)

paramkienerX5(X, dgts = 3, i = 4, maxk = 20, maxe = 0.9,
  parnames = TRUE, dimnames = FALSE, ncores = 1)

Arguments

X	numeric. Vector, matrix, array or list of quantiles.
algo	character. The algorithm used: "r" or "reg" for regression (default) and "e" or "estim" for quantile estimation.
ord	integer. Option for probability selection and treatment.
maxk	numeric. The maximum value of tail parameter k.
mink	numeric. The minimum value of tail parameter k.
maxe	numeric. The maximum value of absolute tail parameter |e|.
probak	numeric. Ordered vector of probabilities.
dgts	integer. The rounding of output parameters.
exfitk	character. A vector of parameter names to subset the output.
dimnames	boolean. Display dimnames.
ncores	integer. The number of cores for parallel processing of arrays.
parnames	boolean. Display parameter names.
n	integer. The 1:n and (N+i-n):N elements of X used to calculate synthetic quantiles at probability levels p1 and 1-p1.
i	integer. The i-th and (N-i)-th elements of X used to extract probabilities p1 and 1-p1 and quantiles x(p) and x(1-p).

Details

FatTailsR package currently uses two different algorithms to estimate the parameters of Kiener distributions K1, K2, K3 and K4.

• Functions fitkienerX(algo = "reg"), paramkienerX(algo = "reg") and regkienerLX use an unweighted nonlinear regression from logit(p) to X over the whole dataset. Depending the size of the dataset, calculation can be slow but is usually accurate and describes very well the last 1-10 points in the tails (except if there is a huge outlier).

• Functions fitkienerX(algo = "estim"), paramkienerX(algo = "estim"), paramkienerX5 and paramkienerX7 estimate the parameters with just 5 to 11 quantiles, 5 being the minimum. For averaging purpose, 11 quantiles are proposed (see below). Computation is almost instantaneous and reasonnably accurate. This is the recommanded method for intensive computation.

A typical input is a numeric vector or a matrix that describes the returns of a stock. A matrix must be in the format DS with DATES as rownames, STOCKS as colnames and (log-)returns as the content of the matrix. An array must be in the format DSL with DATES as rownames, STOCKS as colnames LAGS in the third dimension and (log-)returns as the content of the array. A list can be a list of numeric but neither a list of matrix, a list of data.frame or a list of arrays.

Conversion from a (possible) time series format to a sorted numeric vector is done automatically and without any check of the initial format. Empirical probabilities of each point in the sorted dataset is calculated with the function ppoints whose parameter a has been set to a = 0 as large datasets are very common in finance. The lowest acceptable size of a dataset is not clear at this moment. A minimum of 11 points has been set in "reg" algorithm and a minimum of 15 points has been set in "estim" algorithm. It might change in the future. If possible, use at least 21 points.

Parameter algo controls the algorithm used. Default is "reg".

When algo = "reg" (or algo = "r"), a nonlinear regression is performed with nlsLM from the logit of the empirical probabilities logit(p) over the quantiles X with the function qlkiener4. The maximum value of the tail parameter k is controlled by maxk. An upper value maxk = 10 is appropriate for datasets of low and medium size, less than 20.000 or 50.000 points. For larger datasets, the upper limit can be extended up to maxk = 20. When this limit is reached, the shape of the distribution is very similar to the logistic distribution (at least when e = 0) and the use of this distribution should be considered. Remember that value k < 2 describes a distribution with no stable variance and k < 1 describes a distribution with no stable mean.

When algo = "estim" (or algo = "e"), 5 to 11 quantiles are used to estimate the parameters. The minimum is 5 quantiles : the median x.50, two quantiles at medium distance to the median, usually x.25 and x.75 and two quantiles located close to the extremes of the dataset, for instance x.01 and x.99 if the dataset X has more than 100 points, x.0001 and x.9999 if the dataset X has more than 10.000 points and so on if the dataset is larger. These quantiles are extracted with function fiveprobs. Small datasets must contain at least 15 different points.

With the idea of averaging the results (but without any guarantee of better estimates), calculation has been extended to 11 probabilities extracted from X with the function elevenprobs where p1, p2 and p3 are the most extreme probabilities of the dataset X with values finishing either by .x01 or .x025 or .x05:

• p11 = c(p1, p2, p3, 0.25, 0.35, 0.50, 0.65, 0.75, 1-p3, 1-p2, 1-p1)

Selection of subsets among these 11 probabilities is controlled with the option ord which can take 12 different values. For instance, the default ord = 7 computes the parameters at probabilities c(p1, 0.25, 0.50, 0.75, 1-p1) and c(p2, 0.25, 0.50, 0.75, 1-p2). Parameters d and k are averaged first and the results of these averages are used to compute the other parameters g, a, w, e. Small dataset should consider ord = 5 and large dataset can consider ord = 12. The 12 possible values of ord are:

1. c(p1, 0.35, 0.50, 0.65, 1-p1)

2. c(p2, 0.35, 0.50, 0.65, 1-p2)

3. c(p1, p2, 0.35, 0.50, 0.65, 1-p2, 1-p1)

4. c(p1, p2, p3, 0.35, 0.50, 0.65, 1-p3, 1-p2, 1-p1)

5. c(p1, 0.25, 0.50, 0.75, 1-p1)

6. c(p2, 0.25, 0.50, 0.75, 1-p2)

7. c(p1, p2, 0.25, 0.50, 0.75, 1-p2, 1-p1)

8. c(p1, p2, p3, 0.25, 0.50, 0.75, 1-p3, 1-p2, 1-p1)

9. c(p1, 0.25, 0.35, 0.50, 0.65, 0.75, 1-p1)

10. c(p2, 0.25, 0.35, 0.50, 0.65, 0.75, 1-p2)

11. c(p1, p2, 0.25, 0.35, 0.50, 0.65, 0.75, 1-p2, 1-p1)

12. c(p1, p2, p3, 0.25, 0.35, 0.50, 0.65, 0.75, 1-p3, 1-p2, 1-p1)

paramkienerX5 is a simplified version of paramkienerX with predefined values algo = "estim", ord = 5, maxk = 10 and direct access to internal subfunctions. It uses the following probabilities:

• p5 = c(p1, 0.25, 0.50, 0.75, 1-p1)

paramkienerX7 is a simplified version of paramkienerX with predefined values algo = "estim", ord = 7, maxk = 10 and direct access to internal subfunctions. It uses the following probabilities:

• p7 = c(p1, p2, 0.25, 0.50, 0.75, 1-p2, 1-p1)

The quantiles corresponding to the above probabilities are then extracted with the function quantile whose parameter type has been set to type = 6 as it returns the closest values to the true quantiles (according to our experience) for all k > 1.9. (Note: when k < 1.5, algorithm algo = "reg" returns better results). Both probabilities and quantiles are then transfered to estimkiener11 for calculation.

probak controls the probabilities at which the model is tested with the parameter estimates. fitkienerX and regkienerLX share the same subroutines. The default for fitkienerX and regkienerLX is pprobs2 = c(0.01, 0.025, 0.05, 0.95, 0.975, 0.99) as those values are usual in finance. Other sets of values are provided at pprobs0.

Rounding the results is useful to display nice results, especially in a matrix or in a data.frame. dgts = 13 is recommanded as a, k, w are usually significant at 1 digit.

• dgts = NULL does not perform any rounding.

• dgts = 0 to 9 rounds all parameters at the same level.

• dgts = 10 to 27 rounds the parameters at various levels for nice display. See roundcoefk for the details. (Note: the rounding 10 to 27 currently works with paramkienerX, paramkienerX5, paramkienerX7 but not yet with fitkienerX).

Extracting the most useful parameters from the (quite long) vector/matrix fitk is controlled by parameter exfitk that calls user-defined or predefined parameter subsets like exfit0, ..., exfit7. IMPORTANT: never subset fitk by rank number as new items may be added in the future and rank may vary.

Calculation of vectors, matrices and lists is not parallelized. Parallelization of code for arrays was introduced in version 1.5-0 and improved in version 1.5-1. ncores controls the number of cores allowed to the process (through parApply which runs on Unices and Windows and requires about 2 seconds to start). ncores = 1 means no parallelization. ncores = 0 is the recommanded option. It uses the maximum number of cores available on the computer, as detected by detectCores, minus 1 core, which gives the best performance in most cases. Although appealing, this automatic selection may be sometimes dangerous. For instance, the instruction f(X, ncores_max) - f(X, ncores_max), a nice way to compute an array of 0, will call 2 ncores_max and crash R. ncores = 2,..,99 sets manually the number of cores. If the requested value is larger than the maximum number of cores, this value is automatically reduced (with a warning) to this maximum. Hence, this latest option provides one core more than option ncores = 0.

NOTE: fitkienerLX, regkienerX, estimkiener(X,5,7) were introduced in v1.2-0 and replaced in version v1.4-1 by fitkienerX and paramkiener(X,5,7) to accomodate vector, matrix, arrays and lists. We apologize to early users who need to rewrite their codes.

Value

paramkienerX: a vector (or a matrix) of parameter estimates c(m, g, a, k, w, d, e).

fitkienerX: a vector (or a matrix) made of several parts:

• ret : the return over the period calculated with sum(x). Thus, assume log-returns.

• m, g, a, k, w, d, e : the parameter estimates.

• m1, sd, sk, ke : the mean, standard deviation, skewness and excess of kurtosis computed from the parameter estimates.

• m1x, sdx, skx, kex : The mean, standard deviation, skewness and excess of kurtosis computed from the dataset.

• lh : the length of the dataset over the period.

• q. : quantile estimated with the parameter estimates.

• VaR. : Value-at-Risk, positive in most cases.

• c. : corrective tail coefficient = (q - m) / (q_logistic_function - m).

• ltm. : left tail mean (signed ES on the left tail, usually negative).

• rtm. : right tail mean (signed ES on the right tail, usually positive).

• dtmq. : (p<=0.5 left, p>0.5 right) tail mean minus quantile.

• ES. : expected shortfall, positive in most cases.

• h. : corrective ES = (ES - m) / (ES_logistic_function - m).

• desv. : ES - VaR, usually positive.

• l. : quantile estimated by the tangent logistic function.

• dl. : quantile - quantile_logistic_function.

• g. : quantile estimated by the Laplace-Gauss function.

• dg. : quantile - quantile_Laplace_Gauss_function.

IMPORTANT : if you need to subset fitk, always subset it by parameter names and never subset it by rank number as new items may be added in the future and rank may vary. Use for instance exfit0, ..., exfit7.

References

P. Kiener, Fat tail analysis and package FatTailsR, 9th R/Rmetrics Workshop and Summer School, Zurich, 27 June 2015. https://www.inmodelia.com/exemples/2015-0627-Rmetrics-Kiener-en.pdf

See Also

regkienerLX, estimkiener11, roundcoefk, exfit6.

Examples

require(minpack.lm)
require(timeSeries)

### Load the datasets and choose j in 1:16
DS     <- getDSdata()
j      <- 5

### and run this block
probak <- c(0.01, 0.05, 0.95, 0.99)
X      <- DS[[j]] ; names(DS)[j]
elevenprobs(X)
fitkienerX(X, algo = "reg", dgts = 3, probak = probak)
fitkienerX(X, algo = "estim", ord = 5, probak = probak, dgts = 3)
paramkienerX(X)
paramkienerX5(X)

### Compare the 12 values of paramkienerX(ord/row = 1:12) and paramkienerX (row 13)
compare <- function(ord, X) { paramkienerX(X, ord, algo = "estim", dgts = 13) }
rbind(t(sapply( 1:12, compare, X)), paramkienerX(X, algo = "reg", dgts = 13))

### Analyze DS in one step
t(sapply(DS, paramkienerX, algo = "reg", dgts = 13))
t(sapply(DS, paramkienerX, algo = "estim", dgts = 13))
paramkienerX(DS, algo = "reg", dgts = 13)
paramkienerX(DS, algo = "estim", dgts = 13)
system.time(fitk_rDS <- fitkienerX(DS, algo = "r", probak = pprobs2, dgts = 3))
system.time(fitk_eDS <- fitkienerX(DS, algo = "e", probak = pprobs2, dgts = 3))
fitk_rDS
fitk_eDS

### Subset rDS and eDS with exfit0,..,exfit7
fitk_rDS[,exfit4]
fitk_eDS[,exfit7]
fitkienerX(DS, algo = "e", probak = pprobs2, dgts = 3, exfitk = exfit7)

### Array (new example introduced in v1.5-1)
### Increase the number of cores and crash R.
## Not run:
arr <- array(rkiener1(3000), c(4,3,250))
paramkienerX7(arr, ncores = 2)
## paramkienerX7(arr, ncores = 2) - paramkienerX(arr, ncores = 2)
## End(Not run)

### End

---

Get DS Dataset

Description

A function to extract the log-returns of 16 financial series and time series provided by the packages datasets (EuStockMarkets, sunspot.year) and timeSeries (USDCHF, MSFT, LPP2005REC). The 16 datasets are converted to a list of numeric without any reference to the original dates. This list is usually called DS, hence the name.

Usage

getDSdata()

Details

The dataset is usually created by the instruction DS <- getDSdata(). Then, it is used with a call to DS[[j]] with j in 1:16.

1. "USDCHF" (USDCHF, timeSeries)

2. "MSFT" (MSFT, timeSeries)

3. "DAX" (EuStockMarkets, datasets)

4. "SMI" (EuStockMarkets, datasets)

5. "CAC" (EuStockMarkets, datasets)

6. "FTSE" (EuStockMarkets, datasets)

7. "SBI" (LPP2005REC, timeSeries)

8. "SPI" (LPP2005REC, timeSeries)

9. "SII" (LPP2005REC, timeSeries)

10. "LMI" (LPP2005REC, timeSeries)

11. "MPI" (LPP2005REC, timeSeries)

12. "ALT" (LPP2005REC, timeSeries)

13. "LPP25" (LPP2005REC, timeSeries)

14. "LPP40" (LPP2005REC, timeSeries)

15. "LPP60" (LPP2005REC, timeSeries)

16. "sunspot" (sunspot.year, datasets)

Note that sunspot.year is regularly updated with each new version of R. The generated dataset is logreturn(sunspot.year + 1000).

See Also

EuStockMarkets, sunspot.year, TimeSeriesData, regkienerLX, fitkienerX

Examples

require(timeSeries) 

getDSdata
DS  <- getDSdata()
attributes(DS)
sapply(DS, length)
sapply(DS, head)

---

Generate a list of vectors of characters from a vector of probabilities

Description

Generate vector of characters from a vector of probabilities, replacing 0. by letters:

• p. : probability.

• q. : quantile.

• VaR. : Value-at-Risk, positive in most cases.

• c. : corrective tail coefficient = (q - m) / (q_logistic_function - m).

• ltm. : left tail mean (signed ES on the left tail, usually negative).

• rtm. : right tail mean (signed ES on the right tail, usually positive).

• dtmq. : (p<=0.5 left, p>0.5 right) tail mean minus quantile.

• ES. : expected shortfall, positive in most cases.

• h. : corrective ES = (ES - m) / (ES_logistic_function - m).

• desv. : ES - VaR, usually positive.

• l. : quantile of the tangent logistic function.

• dl. : quantile - quantile_logistic_function.

• g. : quantile of the Laplace-Gauss function.

• dg. : quantile - quantile_Laplace_Gauss_function.

, q., VaR., c., ltm., rtm., ES., h., l., dl., g., dg.. The result is a list of vectors.

Usage

getnamesk(probak = pprobs2, check = TRUE)

getnprobak(probak = pprobs2, check = TRUE)

Arguments

probak	a vector of ordered probabilities with 0 and 1 excluded.
check	boolean. Apply checkquantiles function.

See Also

Probabilities: pprobs0

Examples

getnamesk(pprobs1)
getnamesk(pprobs8)

---

Kashp Function

Description

kashp, which stands for kappa times arc-sine-hyperbola-power is the nonlinear transformation of x at the heart of power hyperbolas, power hyperbolic functions and symmetric Kiener distributions. dkashp_dx is its derivative with respect to x. ashp is provided for convenience.

Usage

kashp(x, k = 1)

dkashp_dx(x, k = 1)

ashp(x, k = 1)

Arguments

x	a numeric value, vector or matrix.
k	a numeric value or vector, preferably strictly positive.

Details

ashp function is defined for x in (-Inf, +Inf) by:

$ashp(x, k) = asinh(x / 2 / k)$

kashp function is defined for x in (-Inf, +Inf) by:

$kashp(x, k) = k * asinh(x / 2 / k)$

dkashp_dx function is defined for x in (-Inf, +Inf) by:

$dkashp_dx(x, k) = 1 / sqrt( x*x/k/k + 4 ) = 1 / 2 / cosh( ashp(x, k) )$

If k is a vector, then the use of the function outer is recommanded.

The undesired case k=0 returns 0 for kashp and dkashp_dx, 1 for exphp, -Inf, NaN, +Inf for ashp.

See Also

The power hyperbolas and the power hyperbolic functions exphp.

Examples

require(graphics)

### Example 1
x    <- (-3:3)*3 ; names(x) <- x
kashp(x, k=2)
k    <- c(-2, 0, 1, 2, 3, 5, 10) ; names(k) <- k
outer(x, k, kashp)
outer(x, k, exphp)

### Example 2
xmin   <- -10
xd     <- 0.5
x      <- seq(xmin, -xmin, xd) ; names(x) <- x
k      <- c(0.6, 1, 1.5, 2, 3.2, 10) ; names(k) <- k
olty   <- c(2, 1, 2, 1, 2, 1, 1)
olwd   <- c(1, 1, 2, 2, 3, 4, 2)
ocol   <- c(2, 2, 4, 4, 3, 3, 1)
op     <- par(mfrow = c(2,2), mgp = c(1.5,0.8,0), mar = c(3,3,2,1))

Tkashp <- ts(cbind(outer(x, k, kashp), "x/2" = x/2), start = xmin, deltat = xd)
plot(Tkashp, plot.type = "single", ylim = c(-5, +5), 
       lty = olty, lwd = olwd, col = ocol, xaxs = "i", yaxs = "i", xlab = "", 
       ylab = "", main = "kashp(x, k)" )
legend("topleft", title = expression(kappa), legend = colnames(Tkashp), 
       inset = 0.02, lty = olty, lwd = olwd, col = ocol, cex = 0.7 )

Tdkashp <- ts(cbind(outer(x, k, dkashp_dx)), start = xmin, deltat = xd)
plot(Tdkashp, plot.type = "single", ylim = c(0, 0.8), 
       lty = olty, lwd = olwd, col = ocol, xaxs = "i", yaxs = "i", xlab = "", 
       ylab = "", main="dkashp_dx(x, k)" )
legend("topleft", title = expression(kappa), legend = colnames(Tdkashp), 
       inset = 0.02, lty = olty, lwd = olwd, col = ocol, cex = 0.7 )

Tashp <- ts(cbind(outer(x, k, ashp), "x/2" = x/2), start = xmin, deltat = xd)
plot(Tashp, plot.type = "single", ylim = c(-5, +5), 
       lty = olty, lwd = olwd, col = ocol, xaxs = "i", yaxs = "i", xlab = "", 
       ylab = "", main = "ashp(x, k)" )
legend("topleft", title = expression(kappa), legend = colnames(Tashp), 
       inset = 0.02, lty = olty, lwd = olwd, col = ocol, cex = 0.7 )
### End example 2

---

Symmetric Kiener Distribution K1

Description

Density, distribution function, quantile function, random generation, value-at-risk, expected shortfall (+ signed left/right tail mean) and additional formulae for symmetric Kiener distribution K1. This distribution is similar to the power hyperbola logistic distribution but with additional parameters for location (m) and scale (g).

Usage

dkiener1(x, m = 0, g = 1, k = 3.2, log = FALSE)

pkiener1(q, m = 0, g = 1, k = 3.2, lower.tail = TRUE, log.p = FALSE)

qkiener1(p, m = 0, g = 1, k = 3.2, lower.tail = TRUE, log.p = FALSE)

rkiener1(n, m = 0, g = 1, k = 3.2)

dpkiener1(p, m = 0, g = 1, k = 3.2, log = FALSE)

dqkiener1(p, m = 0, g = 1, k = 3.2, log = FALSE)

lkiener1(x, m = 0, g = 1, k = 3.2)

dlkiener1(lp, m = 0, g = 1, k = 3.2, log = FALSE)

qlkiener1(lp, m = 0, g = 1, k = 3.2, lower.tail = TRUE)

varkiener1(p, m = 0, g = 1, k = 3.2, lower.tail = TRUE,
  log.p = FALSE)

ltmkiener1(p, m = 0, g = 1, k = 3.2, lower.tail = TRUE,
  log.p = FALSE)

rtmkiener1(p, m = 0, g = 1, k = 3.2, lower.tail = TRUE,
  log.p = FALSE)

dtmqkiener1(p, m = 0, g = 1, k = 3.2, lower.tail = TRUE,
  log.p = FALSE)

eskiener1(p, m = 0, g = 1, k = 3.2, lower.tail = TRUE, log.p = FALSE,
  signedES = FALSE)

Arguments

x	vector of quantiles.
m	numeric. The median.
g	numeric. The scale parameter, preferably strictly positive.
k	numeric. The tail parameter, preferably strictly positive.
log	logical. If TRUE, densities are given in log scale.
q	vector of quantiles.
lower.tail	logical. If TRUE, use p. If FALSE, use 1-p.
log.p	logical. If TRUE, probabilities p are given as log(p).
p	vector of probabilities.
n	number of observations. If length(n) > 1, the length is taken to be the number required.
lp	vector of logit of probabilities.
signedES	logical. FALSE (default) returns positive numbers for left and right tails. TRUE returns negative number (= ltmkiener1) for left tail and positive number (= rtmkiener1) for right tail.

Details

Kiener distributions use the following parameters, some of them being redundant. See aw2k and pk2pk for the formulas and the conversion between parameters:

• m (mu) is the median of the distribution,.

• g (gamma) is the scale parameter.

• a (alpha) is the left tail parameter.

• k (kappa) is the harmonic mean of a and w and describes a global tail parameter.

• w (omega) is the right tail parameter.

• d (delta) is the distortion parameter.

• e (epsilon) is the eccentricity parameter.

Kiener distributions K1(m, g, k, ...) describe distributions with symmetric left and right fat tails with tail parameter k. This parameter is the power exponent mentionned in Pareto formula and Karamata theorems.

m is the median of the distribution. g is the scale parameter and the inverse of the density at the median: $g = 1 / 8 / f(m)$. As a first estimate, it is approximatively one fourth of the standard deviation $g \approx \sigma / 4$ but is independant from it.

dkiener1 function is defined for x in (-Inf, +Inf) by:

$dkiener1(x, m, g, k) = 1 / 4 / g / cosh( ashp((x - m)/g, k) ) / (1 + cosh( kashp((x - m)/g, k)))$

pkiener1 function is defined for q in (-Inf, +Inf) by:

$pkiener1(q, m, g, k) = 1/(1 + exp(- kashp((q - m)/g, k)))$

qkiener1 function is defined for p in (0, 1) by:

$qkiener1(p, m, g, k) = m + 2 * g * k * sinh( logit(p)/k )$

rkiener1 generates n random quantiles.

In addition to the classical d, p, q, r functions, the prefixes dp, dq, l, dl, ql are also provided.

dpkiener1 is the density function calculated from the probability p. It is defined for p in (0, 1) by:

$dpkiener1(p, m, g, k) = p * (1 - p) / 2 / g / cosh( logit(p)/k )$

dqkiener1 is the derivate of the quantile function calculated from the probability p. It is defined for p in (0, 1) by:

$dqkiener1(p, m, g, k) = 2 * g / p / (1 - p) * cosh( logit(p)/k )$

lkiener1 function is equivalent to kashp function but with additional parameters m and g. Being computed from the x (or q) vector, it can be compared to the logit of the empirical probability logit(p) through a nonlinear regression with ordinary or weighted least squares to estimate the distribution parameters. It is defined for x in (-Inf, +Inf) by:

$lkiener1(x, m, g, k) = kashp((x - m)/g, k)$

dlkiener1 is the density function calculated from the logit of the probability lp = logit(p). It is defined for lp in (-Inf, +Inf) by:

$dlkiener1(lp, m, g, k) = p * (1 - p) / 2 / g / cosh( lp/k )$

qlkiener1 is the quantile function calculated from the logit of the probability lp = logit(p). It is defined for lp in (-Inf, +Inf) by:

$qlkiener1(lp, m, g, k) = m + g * k * 2 * sinh( lp/k )$

varkiener1 designates the Value a-risk and turns negative numbers into positive numbers with the following rule:

$varkiener1 <- if(p <= 0.5) (- qkiener1) else (qkiener1)$

Usual values in finance are p = 0.01, p = 0.05, p = 0.95 and p = 0.99. lower.tail = FALSE uses 1-p rather than p.

ltmkiener1, rtmkiener1 and eskiener1 are respectively the left tail mean, the right tail mean and the expected shortfall of the distribution (sometimes called average VaR, conditional VaR or tail VaR). Left tail mean is the integrale from -Inf to p of the quantile function qkiener1 divided by p. Right tail mean is the integrale from p to +Inf of the quantile function qkiener1 divided by 1-p. Expected shortfall turns negative numbers into positive numbers with the following rule:

$eskiener1 <- if(p <= 0.5) (- ltmkiener1) else (rtmkiener1)$

Usual values in finance are p = 0.01, p = 0.025, p = 0.975 and p = 0.99. lower.tail = FALSE uses 1-p rather than p.

dtmqkiener1 is the difference between the left tail mean and the quantile when (p <= 0.5) and the difference between the right tail mean and the quantile when (p > 0.5). It is in quantile unit and is an indirect measure of the tail curvature.

References

P. Kiener, Explicit models for bilateral fat-tailed distributions and applications in finance with the package FatTailsR, 8th R/Rmetrics Workshop and Summer School, Paris, 27 June 2014. Download it from: https://www.inmodelia.com/exemples/2014-0627-Rmetrics-Kiener-en.pdf

P. Kiener, Fat tail analysis and package FatTailsR, 9th R/Rmetrics Workshop and Summer School, Zurich, 27 June 2015. Download it from: https://www.inmodelia.com/exemples/2015-0627-Rmetrics-Kiener-en.pdf

C. Acerbi, D. Tasche, Expected shortfall: a natural coherent alternative to Value at Risk, 9 May 2001. Download it from: https://www.bis.org/bcbs/ca/acertasc.pdf

See Also

Power hyperbola logistic distribution logishp, asymmetric Kiener distributions K2, K3, K4 and K7 kiener2, kiener3, kiener4, kiener7, regression function regkienerLX.

Examples

require(graphics)

### Example 1
pp <- c(ppoints(11, a = 1), NA, NaN) ; pp
qkiener1(p = pp, k = 4)


### Example 2: Try various value of k = 1.5, 3, 5, 10
k       <- 5  # 1.5, 3, 5, 10
set.seed(2014)
mainTC  <- paste("qkiener1(p, m = 0, g = 1, k = ", k, ")")
mainsum <- paste("cumulated qkiener1(p, m = 0, g = 1, k = ", k, ")")
T       <- 500
C       <- 4
TC      <- qkiener1(p = runif(T*C), m = 0, g = 1, k = k)
matTC   <- matrix(TC, nrow = T, ncol = C, dimnames = list(1:T, letters[1:C]))
head(matTC)
plot.ts(matTC, main = mainTC)
#
matsum  <- apply(matTC, MARGIN=2, cumsum) 
head(matsum)
plot.ts(matsum, plot.type = "single", main = mainsum)
### End example 2


### Example 3 (four plots: probability, density, logit, logdensity)
x  <- q  <- seq(-15, 15, length.out=101)
k     <- c(0.6, 1, 1.5, 2, 3.2, 10) ; names(k) <- k ; k
olty  <- c(2, 1, 2, 1, 2, 1, 1)
olwd  <- c(1, 1, 2, 2, 3, 3, 2)
ocol  <- c(2, 2, 4, 4, 3, 3, 1)
lleg  <- c("logit(0.999) = 6.9", "logit(0.99)   = 4.6", "logit(0.95)   = 2.9", 
           "logit(0.50)   = 0", "logit(0.05)   = -2.9", "logit(0.01)   = -4.6", 
           "logit(0.001) = -6.9  ")
op    <- par(mfrow=c(2,2), mgp=c(1.5,0.8,0), mar=c(3,3,2,1))

plot(x, plogis(x, scale = 2), type = "b", lwd = 2, ylim = c(0, 1),
     xaxs = "i", yaxs = "i", xlab = "", ylab = "", 
     main = "pkiener1(q, m, g, k)")
for (i in 1:length(k)) lines(x, pkiener1(x, k = k[i]), 
       lty = olty[i], lwd = olwd[i], col = ocol[i] )
legend("topleft", title = expression(kappa), legend = c(k, "logistic"), 
       cex = 0.7, inset = 0.02, lty = olty, lwd = olwd, col = ocol )

plot(x, dlogis(x, scale = 2), type = "b", lwd = 2, ylim = c(0, 0.14), 
     xaxs = "i", yaxs = "i", xlab = "", ylab = "", main = "dkiener1(x, m, g, k)" )
for (i in 1:length(k)) lines(x, dkiener1(x, k = k[i]), 
       lty = olty[i], lwd = olwd[i], col = ocol[i] )
legend("topright", title = expression(kappa), legend = c(k, "logistic"), 
       cex = 0.7, inset = 0.02, lty = olty, lwd = olwd, col = ocol )

plot(x, x/2, type = "b", lwd = 2, ylim = c(-7.5, 7.5), yaxt="n", xaxs = "i", 
     yaxs = "i", xlab = "", ylab = "", main = "logit(pkiener1(q, m, g, k))")
axis(2, las=1, at=c(-6.9, -4.6, -2.9, 0, 2.9, 4.6, 6.9) )
for (i in 1:length(k)) lines(x, lkiener1(x, k = k[i]),  
       lty = olty[i], lwd = olwd[i], col = ocol[i] )
lines(x, logit(pnorm(x, 0, 3.192)), type="l", lty=1, lwd=3, col=7) # erfx
legend("topleft", legend = lleg, cex = 0.7, inset = 0.02 )
legend("bottomright", title = expression(kappa), 
       legend = c(k, "logistic", "Gauss"), cex = 0.7, inset = 0.02, 
       lty = c(olty, 1), lwd = c(olwd, 3), col = c(ocol , 7) )

plot(x, log(dlogis(x, scale = 2)), lwd = 2, type = "b", ylim = c(-8, -1.5), 
     xaxs = "i", yaxs = "i", xlab = "", ylab = "", main = "log(dkiener1(x, m, g, k))") 
for (i in 1:length(k)) lines(x, log(dkiener1(x, k = k[i])),  
       lty = olty[i], lwd = olwd[i], col = ocol[i] )
lines(x, dnorm(x, 0, 3.192, log = TRUE), type = "l", lty = 1, lwd = 3, col = 7)
legend("bottom", title = expression(kappa), legend = c(k, "logistic", "Gauss"), 
       cex = 0.7, inset = 0.02, lty = c(olty, 1), lwd = c(olwd, 3), col = c(ocol , 7) )
### End example 3


### Example 4 (four plots: quantile, derivate, density and quantiles from p)
p   <- ppoints(199, a=0)
k   <- c(0.6, 1, 1.5, 2, 3.2, 10) ; names(k) <- k ; k
op  <- par(mfrow=c(2,2), mgp=c(1.5,0.8,0), mar=c(3,3,2,1))
plot(p, qlogis(p, scale = 2), type = "o", lwd = 2, ylim = c(-15, 15),
     xaxs = "i", yaxs = "i", xlab = "", ylab = "", 
     main = "qkiener1(p, m, g, k)")
for (i in 1:length(k)) lines(p, qkiener1(p, k = k[i]), 
          lty = olty[i], lwd = olwd[i], col = ocol[i] )
legend("topleft", title = expression(kappa), legend = c(k, "qlogis(x/2)"), 
          inset = 0.02, lty = olty, lwd = olwd, col = ocol, cex = 0.7 )

plot(p, 2/p/(1-p), type = "o", lwd = 2, xlim = c(0, 1), ylim = c(0, 100),
     xaxs = "i", yaxs = "i", xlab = "", ylab = "", 
     main = "dqkiener1(p, m, g, k)")
for (i in 1:length(k)) lines(p, dqkiener1(p, k = k[i]), 
          lty = olty[i], lwd = olwd[i], col = ocol[i] )
legend("top", title = expression(kappa), legend = c(k, "p*(1-p)/2"), 
          inset = 0.02, lty = olty, lwd = olwd, col = ocol, cex = 0.7 )

plot(qlogis(p, scale = 2), p*(1-p)/2, type = "o", lwd = 2, xlim = c(-15, 15), 
     ylim = c(0, 0.14), xaxs = "i", yaxs = "i", xlab = "", ylab = "", 
     main = "qkiener1, dpkiener1(p, m, g, k)")
for (i in 1:length(k)) lines(qkiener1(p, k = k[i]), dpkiener1(p, k = k[i]),
          lty = olty[i], lwd = olwd[i], col = ocol[i] )
legend("topleft", title = expression(kappa), legend = c(k, "p*(1-p)/2"), 
          inset = 0.02, lty = olty, lwd = olwd, col = ocol, cex = 0.7 )
### End example 4


### Example 5 (q and VaR, ltm, rtm, and ES)
pp <- c(0.001, 0.0025, 0.005, 0.01, 0.025, 0.05, 
        0.10, 0.20, 0.35, 0.5, 0.65, 0.80, 0.90,
        0.95, 0.975, 0.99, 0.995, 0.9975, 0.999)
m <- -10 ; g <- 1 ; k <- 4
round(c(m = m, g = g, a = k, k = k, w = k, d = 0, e = 0), 2) 
plot(qkiener1(pp, m, g, k), pp, type = "b")
round(cbind(p = pp, "1-p" = 1-pp, 
	q   =   qkiener1(pp, m, g, k), 
	ltm = ltmkiener1(pp, m, g, k), 
	rtm = rtmkiener1(pp, m, g, k), 
	es  =  eskiener1(pp, m, g, k), 
	VaR = varkiener1(pp, m, g, k)), 4)
round(kmean(c(m, g, k), model = "K1"), 4) # limit value of ltm, rtm
round(cbind(p = pp, "1-p" = 1-pp, 
	q   =   qkiener1(pp, m, g, k, lower.tail = FALSE), 
	ltm = ltmkiener1(pp, m, g, k, lower.tail = FALSE), 
	rtm = rtmkiener1(pp, m, g, k, lower.tail = FALSE), 
	es  =  eskiener1(pp, m, g, k, lower.tail = FALSE), 
	VaR = varkiener1(pp, m, g, k, lower.tail = FALSE)), 4)
### End example 5

---

Asymmetric Kiener Distribution K2

Description

Density, distribution function, quantile function, random generation, value-at-risk, expected shortfall (+ signed left/right tail mean) and additional formulae for asymmetric Kiener distribution K2.

Usage

dkiener2(x, m = 0, g = 1, a = 3.2, w = 3.2, log = FALSE)

pkiener2(q, m = 0, g = 1, a = 3.2, w = 3.2, lower.tail = TRUE,
  log.p = FALSE)

qkiener2(p, m = 0, g = 1, a = 3.2, w = 3.2, lower.tail = TRUE,
  log.p = FALSE)

rkiener2(n, m = 0, g = 1, a = 3.2, w = 3.2)

dpkiener2(p, m = 0, g = 1, a = 3.2, w = 3.2, log = FALSE)

dqkiener2(p, m = 0, g = 1, a = 3.2, w = 3.2, log = FALSE)

lkiener2(x, m = 0, g = 1, a = 3.2, w = 3.2)

dlkiener2(lp, m = 0, g = 1, a = 3.2, w = 3.2, log = FALSE)

qlkiener2(lp, m = 0, g = 1, a = 3.2, w = 3.2, lower.tail = TRUE)

varkiener2(p, m = 0, g = 1, a = 3.2, w = 3.2, lower.tail = TRUE,
  log.p = FALSE)

ltmkiener2(p, m = 0, g = 1, a = 3.2, w = 3.2, lower.tail = TRUE,
  log.p = FALSE)

rtmkiener2(p, m = 0, g = 1, a = 3.2, w = 3.2, lower.tail = TRUE,
  log.p = FALSE)

dtmqkiener2(p, m = 0, g = 1, a = 3.2, w = 3.2, lower.tail = TRUE,
  log.p = FALSE)

eskiener2(p, m = 0, g = 1, a = 3.2, w = 3.2, lower.tail = TRUE,
  log.p = FALSE, signedES = FALSE)

Arguments

x	vector of quantiles.
m	numeric. The median.
g	numeric. The scale parameter, preferably strictly positive.
a	numeric. The left tail parameter, preferably strictly positive.
w	numeric. The right tail parameter, preferably strictly positive.
log	logical. If TRUE, densities are given in log scale.
q	vector of quantiles.
lower.tail	logical. If TRUE, use p. If FALSE, use 1-p.
log.p	logical. If TRUE, probabilities p are given as log(p).
p	vector of probabilities.
n	number of observations. If length(n) > 1, the length is taken to be the number required.
lp	vector of logit of probabilities.
signedES	logical. FALSE (default) returns positive numbers for left and right tails. TRUE returns negative number (= ltmkiener4) for left tail and positive number (= rtmkiener4) for right tail.

Details

Kiener distributions use the following parameters, some of them being redundant. See aw2k and pk2pk for the formulas and the conversion between parameters:

• m (mu) is the median of the distribution,.

• g (gamma) is the scale parameter.

• a (alpha) is the left tail parameter.

• k (kappa) is the harmonic mean of a and w and describes a global tail parameter.

• w (omega) is the right tail parameter.

• d (delta) is the distortion parameter.

• e (epsilon) is the eccentricity parameter.

Kiener distributions K2(m, g, a, w) are distributions with asymmetrical left and right fat tails described by the parameters a (alpha) for the left tail and w (omega) for the right tail. These parameters correspond to the power exponent that appear in Pareto formula and Karamata theorems.

As a and w are highly correlated, the use of Kiener distributions (K3(..., k, d) K4 (K4(..., k, e) is an alternate solution.

m is the median of the distribution. g is the scale parameter and the inverse of the density at the median: $g = 1 / 8 / f(m)$. As a first estimate, it is approximatively one fourth of the standard deviation $g \approx \sigma / 4$ but is independant from it.

The d, p functions have no explicit forms. They are provided here for convenience. They are estimated from a reverse optimization on the quantile function and can be (very) slow, depending the number of points to estimate. We recommand to use the quantile function as far as possible. WARNING: Results may become inconsistent when a or w are smaller than 1. Hopefully, this case seldom happens in finance.

qkiener2 function is defined for p in (0, 1) by:

$qkiener2(p, m, g, a, w) = m + g * k * (- exp(-logit(p)/a) + exp(logit(p)/w) )$

where k is the harmonic mean of the tail parameters a and w calculated by $k = aw2k(a, w)$.

rkiener2 generates n random quantiles.

In addition to the classical d, p, q, r functions, the prefixes dp, dq, l, dl, ql are also provided.

dpkiener2 is the density function calculated from the probability p. It is defined for p in (0, 1) by:

$dpkiener2(p, m, g, a, w) = p * (1 - p) / k / g / ( exp(-logit(p)/a)/a + exp(logit(p)/w)/w$

dqkiener2 is the derivate of the quantile function calculated from the probability p. It is defined for p in (0, 1) by:

$dqkiener2(p, m, g, a, w) = k * g / p / (1 - p) * ( exp(-logit(p)/a)/a + exp(logit(p)/w)/w )$

lkiener2 function is estimated from a reverse optimization and can be (very) slow depending the number of points to estimate. Initialization is done by assuming a symmetric distribution lkiener1 around the harmonic mean k, then optimization is performed to take into account the true values a and w. The result can be then compared to the empirical probability logit(p). WARNING: Results may become inconsistent when a or w are smaller than 1. Hopefully, this case seldom happens in finance.

dlkiener2 is the density function calculated from the logit of the probability lp = logit(p). it is defined for lp in (-Inf, +Inf) by:

$dlkiener2(lp, m, g, a, w) = p * (1 - p) / k / g / ( exp(-lp/a)/a + exp(lp/w)/w )$

qlkiener2 is the quantile function calculated from the logit of the probability. It is defined for lp in (-Inf, +Inf) by:

$qlkiener2(lp, m, g, a, w) = m + g * k * ( - exp(-lp/a) + exp(lp/w) )$

varkiener2 designates the Value a-risk and turns negative numbers into positive numbers with the following rule:

$varkiener2 <- if(p <= 0.5) { - qkiener2 } else { qkiener2 }$

Usual values in finance are p = 0.01, p = 0.05, p = 0.95 and p = 0.99. lower.tail = FALSE uses 1-p rather than p.

ltmkiener2, rtmkiener2 and eskiener2 are respectively the left tail mean, the right tail mean and the expected shortfall of the distribution (sometimes called average VaR, conditional VaR or tail VaR). Left tail mean is the integrale from -Inf to p of the quantile function qkiener2 divided by p. Right tail mean is the integrale from p to +Inf of the quantile function qkiener2 divided by 1-p. Expected shortfall turns negative numbers into positive numbers with the following rule:

$eskiener2 <- if(p <= 0.5) { - ltmkiener2 } else { rtmkiener2 }$

Usual values in finance are p = 0.01, p = 0.025, p = 0.975 and p = 0.99. lower.tail = FALSE uses 1-p rather than p.

dtmqkiener2 is the difference between the left tail mean and the quantile when (p <= 0.5) and the difference between the right tail mean and the quantile when (p > 0.5). It is in quantile unit and is an indirect measure of the tail curvature.

References

P. Kiener, Explicit models for bilateral fat-tailed distributions and applications in finance with the package FatTailsR, 8th R/Rmetrics Workshop and Summer School, Paris, 27 June 2014. Download it from: https://www.inmodelia.com/exemples/2014-0627-Rmetrics-Kiener-en.pdf

P. Kiener, Fat tail analysis and package FatTailsR, 9th R/Rmetrics Workshop and Summer School, Zurich, 27 June 2015. Download it from: https://www.inmodelia.com/exemples/2015-0627-Rmetrics-Kiener-en.pdf

C. Acerbi, D. Tasche, Expected shortfall: a natural coherent alternative to Value at Risk, 9 May 2001. Download it from: https://www.bis.org/bcbs/ca/acertasc.pdf

See Also

Symmetric Kiener distribution K1 kiener1, asymmetric Kiener distributions K3, K4 and K7 kiener3, kiener4, kiener7, conversion functions aw2k, estimation function fitkienerX, regression function regkienerLX.

Examples

require(graphics)

### Example 1
pp <- c(ppoints(11, a = 1), NA, NaN) ; pp
lp <- logit(pp) ; lp
qkiener2(  p = pp, m = 2, g = 1.5, a = 4, w = 6)
qkiener2(  p = pp, m = 2, g = 1.5, a = 4, w = 6)
qlkiener2(lp = lp, m = 2, g = 1.5, a = 4, w = 6)
dpkiener2( p = pp, m = 2, g = 1.5, a = 4, w = 6)
dlkiener2(lp = lp, m = 2, g = 1.5, a = 4, w = 6)
dqkiener2( p = pp, m = 2, g = 1.5, a = 4, w = 6)


### Example 2
a       <- 6
w       <- 4
set.seed(2014)
mainTC  <- paste("qkiener2(p, m = 0, g = 1, a = ", a, ", w = ", w, ")")
mainsum <- paste("cumulated qkiener2(p, m = 0, g = 1, a = ", a, ", w = ", w, ")")
T       <- 500
C       <- 4
TC      <- qkiener2(p = runif(T*C), m = 0, g = 1, a = a, w = w)
matTC   <- matrix(TC, nrow = T, ncol = C, dimnames = list(1:T, letters[1:C]))
head(matTC)
plot.ts(matTC, main = mainTC)
#
matsum  <- apply(matTC, MARGIN=2, cumsum) 
head(matsum)
plot.ts(matsum, plot.type = "single", main = mainsum)
### End example 2


### Example 3 (four plots: probability, density, logit, logdensity)
x     <- q  <- seq(-15, 15, length.out=101)
w     <- c(0.6, 1, 1.5, 2, 3.2, 10) ; names(w) <- w
olty  <- c(2, 1, 2, 1, 2, 1, 1)
olwd  <- c(1, 1, 2, 2, 3, 3, 2)
ocol  <- c(2, 2, 4, 4, 3, 3, 1)
lleg  <- c("logit(0.999) = 6.9", "logit(0.99)   = 4.6", "logit(0.95)   = 2.9", 
           "logit(0.50)   = 0", "logit(0.05)   = -2.9", "logit(0.01)   = -4.6", 
           "logit(0.001) = -6.9  ")
op    <- par(mfrow=c(2,2), mgp=c(1.5,0.8,0), mar=c(3,3,2,1))

plot(x, plogis(x, scale = 2), type = "n", lwd = 2, ylim = c(0, 1),
     xaxs = "i", yaxs = "i", xlab = "", ylab = "", 
     main = "pkiener2(q, m, g, a=2, w=...)")
for (i in 1:length(w)) lines(x, pkiener2(x, a = 2, w = w[i]), 
       lty = olty[i], lwd = olwd[i], col = ocol[i] )
legend("topleft", title = expression(omega), legend = c(w), 
       cex = 0.7, inset = 0.02, lty = olty, lwd = olwd, col = ocol )

plot(x, dlogis(x, scale = 2), type = "n", lwd = 2, ylim = c(0, 0.17),
     xaxs = "i", yaxs = "i", xlab = "", ylab = "", 
     main = "dkiener2(q, m, g, a=2, w=...)")
for (i in 1:length(w)) lines(x, dkiener2(x, a = 2, w = w[i]), 
       lty = olty[i], lwd = olwd[i], col = ocol[i] )
legend("topright", title = expression(omega), legend = c(w), 
       cex = 0.7, inset = 0.02, lty = olty, lwd = olwd, col = ocol )

plot(x, x/2, type = "n", lwd = 1, ylim = c(-7.5, 7.5), yaxt="n", xaxs = "i", 
     yaxs = "i", xlab = "", ylab = "", 
     main = "logit(pkiener2(q, m, g, a=2, w=...))")
axis(2, las=1, at=c(-6.9, -4.6, -2.9, 0, 2.9, 4.6, 6.9) )
for (i in 1:length(w)) lines(x, lkiener2(x, a = 2, w = w[i]),  
       lty = olty[i], lwd = olwd[i], col = ocol[i] )
legend("topleft", legend = lleg, cex = 0.7, inset = 0.02 )
legend("bottomright", title = expression(omega), legend = c(w), 
       cex = 0.7, inset = 0.02, lty = c(olty), lwd = c(olwd), col = c(ocol) )

plot(x, dlogis(x, scale = 2, log=TRUE), type = "n", lwd = 2, ylim = c(-8, -1.5),
     xaxs = "i", yaxs = "i", xlab = "", ylab = "", 
     main = "log(dkiener2(q, m, g, a=2, w=...))")
for (i in 1:length(w)) lines(x, dkiener2(x, a = 2, w = w[i], log=TRUE), 
       lty = olty[i], lwd = olwd[i], col = ocol[i] )
legend("bottom", title = expression(omega), legend = c(w), 
       cex = 0.7, inset = 0.02, lty = olty, lwd = olwd, col = ocol )
### End example 3


### Example 4 (four plots: quantile, derivate, density and quantiles from p)
p     <- ppoints(199, a=0)
w     <- c(0.6, 1, 1.5, 2, 3.2, 10) ; names(w) <- w ; w
op    <- par(mfrow=c(2,2), mgp=c(1.5,0.8,0), mar=c(3,3,2,1))

plot(p, qlogis(p, scale = 2), type = "l", lwd = 2, xlim = c(0, 1), 
     ylim = c(-15, 15), xaxs = "i", yaxs = "i", xlab = "", ylab = "", 
     main = "qkiener2(p, m, g, a=2, w=...)")
for (i in 1:length(w)) lines(p, qkiener2(p, a = 2, w = w[i]), 
          lty = olty[i], lwd = olwd[i], col = ocol[i] )
legend("topleft", title = expression(omega), legend = c(w, "qlogis(x/2)"), 
          inset = 0.02, lty = olty, lwd = olwd, col = ocol, cex = 0.7 )

plot(p, 2/p/(1-p), type = "l", lwd = 2, xlim = c(0, 1), ylim = c(0, 100),
     xaxs = "i", yaxs = "i", xlab = "", ylab = "", 
     main = "dqkiener2(p, m, g, a=2, w=...)")
for (i in 1:length(w)) lines(p, dqkiener2(p, a = 2, w = w[i]), 
          lty = olty[i], lwd = olwd[i], col = ocol[i] )
legend("top", title = expression(omega), legend = c(w, "p*(1-p)/2"), 
          inset = 0.02, lty = olty, lwd = olwd, col = ocol, cex = 0.7 )

plot(qlogis(p, scale = 2), p*(1-p)/2, type = "l", lwd = 2, xlim = c(-15, 15), 
     ylim = c(0, 0.18), xaxs = "i", yaxs = "i", xlab = "", ylab = "", 
     main = "qkiener2, dpkiener2(p, m, g, a=2, w=...)")
for (i in 1:length(w)) { 
     lines(qkiener2(p, a = 2, w = w[i]), dpkiener2(p, a = 2, w = w[i]),
           lty = olty[i], lwd = olwd[i], col = ocol[i] ) }
legend("topleft", title = expression(omega), legend = c(w, "p*(1-p)/2"), 
          inset = 0.02, lty = olty, lwd = olwd, col = ocol, cex = 0.7 )

plot(qlogis(p, scale = 2), p, type = "l", lwd = 2, xlim = c(-15, 15), 
     ylim = c(0, 1), xaxs = "i", yaxs = "i", xlab = "", ylab = "", 
     main = "inverse axis qkiener2(p, m, g, a=2, w=...)")
for (i in 1:length(w)) lines(qkiener2(p, a = 2, w = w[i]), p,
          lty = olty[i], lwd = olwd[i], col = ocol[i] )
legend("topleft", title = expression(omega), legend = c(w, "qlogis(x/2)"), 
          inset = 0.02, lty = olty, lwd = olwd, col = ocol, cex = 0.7 )
### End example 4


### Example 5 (q and VaR, ltm, rtm, and ES)
pp <- c(0.001, 0.0025, 0.005, 0.01, 0.025, 0.05, 
        0.10, 0.20, 0.35, 0.5, 0.65, 0.80, 0.90,
        0.95, 0.975, 0.99, 0.995, 0.9975, 0.999)
m <- -10 ; g <- 1 ; a <- 5 ; w = 3 
k <- aw2k(a, w) ; d <- aw2d(a, w) ; e <- aw2e(a, w)
round(c(m = m, g = g, a = a, k = k, w = w, d = d, e = e), 2) 
plot(qkiener2(pp, m, g, a, w), pp, type = "b")
round(cbind(p = pp, "1-p" = 1-pp,
	q   =   qkiener2(pp, m, g, a, w), 
	ltm = ltmkiener2(pp, m, g, a, w), 
	rtm = rtmkiener2(pp, m, g, a, w), 
	ES  =  eskiener2(pp, m, g, a, w), 
	VaR = varkiener2(pp, m, g, a, w)), 4)
round(kmean(c(m, g, a, w), model = "K2"), 4) # limit value for ltm and rtm
round(cbind(p = pp, "1-p" = 1-pp, 
	q   =   qkiener2(pp, m, g, a, w, lower.tail = FALSE), 
	ltm = ltmkiener2(pp, m, g, a, w, lower.tail = FALSE), 
	rtm = rtmkiener2(pp, m, g, a, w, lower.tail = FALSE), 
	ES  =  eskiener2(pp, m, g, a, w, lower.tail = FALSE), 
	VaR = varkiener2(pp, m, g, a, w, lower.tail = FALSE)), 4)
### End example 5

---