https://github.com/cran/quantreg
Tip revision: 92269b50ae23297ededa6a58d77823e1302a6b8e authored by Roger Koenker on 06 June 2021, 16:10:02 UTC
version 5.86
version 5.86
Tip revision: 92269b5
ParetoTest.R
ParetoTest <- function(formula, tau = 0.1, data = NULL, flavor = "Hill",
m = 2, cicov = 0.9, ...) {
if(flavor == "Hill")
z <- Hill(formula, tau, data, ...)
else if(flavor == "Pickands")
z <- Pickands(formula, tau, data, m = m, ...)
else
stop(paste("flavor", flavor, "not (yet) implemented"))
summary(z, ci = cicov, ...)
}
# Pickands Estimation and Inference for the Pareto Tail Exponent
Pickands <- function(formula, tau = 0.1, data = NULL, m = 2, ...){
y <- model.response(model.frame(formula, data))
X <- model.matrix(formula, data)
Pickands.fit(X, y, tau, m, ...)
}
Pickands.fit <- function(X, y, tau, m, ...){
if(tau > 0.5){ # Right tail
X <- -X
y <- -y
tau = 1-tau
}
taus = c(tau, m * tau, 2 * m * tau)
F <- rq(y ~ X - 1, tau = taus, ...)
Xbar = colMeans(X)
num = c(crossprod(Xbar, (F$coef[,3] - F$coef[,2])))
denom = c(crossprod(Xbar, (F$coef[,2] - F$coef[,1])))
xi <- (-1/log(m)) * log(num/denom)
z <- list(xi = xi, tau = tau, m = m, F = F)
class(z) <- "Pickands"
z
}
summary.Pickands <- function(object, se = "boot", B = 200, ci = 0.9, ...){
tau <- object$tau
xi <- object$xi
m <- object$m
X <- object$F$x
y <- object$F$y
n <- length(y)
b1 <- object$F$coef[,1]
bm <- object$F$coef[,2]
Xbar <- colMeans(X)
Xg <- X %*% (bm - b1)/c(crossprod(Xbar, bm - b1)) #aargh humbug!
R <- rep(NA, B)
for(i in 1:B){ # Parametric Bootstrap Default Method
ys <- -(rexp(n)^(-xi) - 1)/xi * Xg
R[i] <- Pickands.fit(X, ys, tau, m = m, ...)$xi
}
z <- rep(NA,6)
z[1] <- xi
z[2] <- xi - quantile(R - xi, .5, na.rm=TRUE)
z[3] <- xi - quantile(R - xi, (1+ci)/2, na.rm=TRUE)
z[4] <- xi - quantile(R - xi, (1-ci)/2, na.rm=TRUE)
z[5] <- sqrt(var(R))
z[6] <- B - sum(is.na(R))
names(z) <- c("Estimate", "Bias-Corrected", paste("Lower ", 50*(1-ci),"%",sep=""),
paste("Upper ", 50*(1+ci),"%",sep=""), "StdErr", "Sample Size")
z <- list(z = z, tau = tau)
class(z) = "summary.Pickands"
z
}
print.Pickands <- function(x, ...){
print(x$xi)
}
print.summary.Pickands <- function(x, ...){
cat(paste("Pickands Estimation and Inference: tau = ", x$tau,
"Effective Bootstrap Sample Size = ",x$z[6],"\n\n"))
print(x$z[1:5])
}
if(FALSE){
# Test problem
require(quantreg)
n = 500
x = rnorm(n)
y = x + rt(n,2)
P = Pickands(y ~ x, .9)
S = summary(P)
}
# Hill Estimation and Inference for the Pareto Tail Exponent
Hill <- function(formula, tau = 0.1, data = NULL, ...){
y <- model.response(model.frame(formula, data))
X <- model.matrix(formula, data)
Hill.fit(X, y, tau)
}
Hill.fit <- function(X, y, tau, ...){
if(tau > 0.5){ # Right tail
X <- -X
y <- -y
tau = 1-tau
}
F <- rq.fit(X, y, tau = tau, ...)
yhat <- F$fitted.values
z <- log(ifelse(y <= yhat, abs(y/yhat), 1))
z <- z[is.finite(z)]
z <- list(xi = sum(z)/(length(z)*tau), tau = tau, F = F)
class(z) <- "Hill"
z
}
summary.Hill <- function(object, se = "boot", B = 200, m = 2, ci = 0.9, ...){
tau <- object$tau
xi <- object$xi
X <- object$F$x
y <- object$F$y
n <- length(y)
b0 <- object$F$coef
bm <- rq.fit(X, y, m * tau, ...)$coef
Xbar <- colMeans(X)
Xg <- X %*% (bm - b0)/c(crossprod(Xbar, bm - b0)) #aargh humbug!
R <- rep(NA, B)
for(i in 1:B){ # Parametric Bootstrap Default Method
ys <- -(rexp(n)^(-xi) - 1)/xi * Xg
R[i] <- Hill.fit(X, ys, tau, ...)$xi
}
z <- rep(NA,6)
z[1] <- xi
z[2] <- xi - quantile(R - xi, .5, na.rm=TRUE)
z[3] <- xi - quantile(R - xi, (1+ci)/2, na.rm=TRUE)
z[4] <- xi - quantile(R - xi, (1-ci)/2, na.rm=TRUE)
z[5] <- sqrt(var(R))
z[6] <- B - sum(is.na(R))
names(z) <- c("Estimate", "Bias-Corrected", paste("Lower ", 50*(1-ci),"%",sep=""),
paste("Upper ", 50*(1+ci),"%",sep=""), "BS_StdErr", "Sample Size")
z <- list(z = z, tau = tau)
class(z) = "summary.Hill"
z
}
print.Hill <- function(x, ...){
print(x$xi)
}
print.summary.Hill <- function(x, ...){
cat(paste("Hill Estimation and Inference: tau = ", x$tau,
"Effective Bootstrap Sample Size = ",x$z[6],"\n\n"))
print(x$z[1:5])
}