Nothing
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Library General Public
# License as published by the Free Software Foundation; either
# version 2 of the License, or (at your option) any later version.
#
# This library is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Library General Public License for more details.
#
# You should have received a copy of the GNU Library General
# Public License along with this library; if not, write to the
# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston,
# MA 02111-1307 USA
################################################################################
# FUNCTION: DESCRIPTION:
# 'fTHETA' Class representation for extremal index
# show.fTHETA S4: Print Method for extremal index
# thetaSim Simulates a time series with known theta
# FUNCTION: DESCRIPTION:
# blockTheta Computes theta from Block Method
# clusterTheta Computes theta from Reciprocal Cluster Method
# runTheta Computes theta from Run Method
# ferrosegersTheta Computes Theta according to Ferro and Segers
# FUNCTION: DESCRIPTION:
# exindexesPlot Computes and Plot Theta(1,2,3)
# exindexPlot Computes Theta(1,2) and Plot Theta(1)
################################################################################
setClass("fTHETA",
representation(
call = "call",
data = "list",
theta = "data.frame",
title = "character",
description = "character")
)
# ------------------------------------------------------------------------------
setMethod("show", "fTHETA",
function(object)
{ # A function implemented by Diethelm Wuertz
# FUNCTION:
# Unlike print the argument for show is 'object'.
x = object
# Title:
cat("\nTitle:\n ", x@title, "\n", sep = "")
# Call:
cat("\nCall:\n ", deparse(x@call), "\n", sep = "")
# Extremal Index:
cat("\nExtremal Index:\n")
print(object@theta)
# Description:
cat("\nDescription:\n ", x@description, sep = "")
cat("\n\n")
# Return Value:
invisible()
})
# ------------------------------------------------------------------------------
thetaSim =
function(model = c("max", "pair"), n = 1000, theta = 0.5)
{ # A function implemented by Diethelm Wuertz
# Description:
# Simulates a time series with known theta
# Arguments:
# model - a character string denoting the model
# "max" - Max Frechet Series
# "pair" - Paired Exponential Series
# FUNCTION:
# Model Argument:
model = match.arg(model)
# Simulate:
model = model[1]
X = rep(0, n)
if (model == "max") {
# Frechet rvs:
eps = 1/(-log(runif(n)))
X[1] = theta*eps[1]
for ( i in 2:n ) X[i] = max( (1-theta)*X[i-1], theta*eps[i] )
} else if (model == "pair") {
theta = 0.5
eps = rexp(n+1)
for ( i in 1:n ) X[i] = max(eps[i], eps[i+1])
}
# As time series:
X = as.ts(X)
attr(X, "control") = c(model = model, theta = as.character(theta))
# Return Value:
X
}
################################################################################
blockTheta =
function (x, block = 22, quantiles = seq(0.950, 0.995, length = 10),
title = NULL, description = NULL)
{ # A function implemented by Diethelm Wuertz
# Description:
# Calculates theta from Block method, i.e. theta1.
# Example:
# blockTheta(thetaSim(n=10000))
# FUNCTION:
# Check if block is numeric:
stopifnot(is.numeric(block))
# Number of blocks and number of data points:
X = as.vector(x)
ordered = sort(X)
k = floor(length(X)/block)
n = k*block
# Now organize your X:
# 1) truncate the rest of the time series,
# 2) arrange them in matrix form,
# 3) sort them in reverse order, ie. from high (pos) to low (neg)
X = matrix(X[1:(k*block)], ncol = block, byrow = TRUE)
# Threshold values associated to quantiles:
thresholds = ordered[floor(quantiles*length(X))]
# Presettings:
theta1 = rep(0, times = length(quantiles))
# Calculate Extremal Imdex:
run = 0
keepK = keepN = NULL
for ( u in thresholds ) {
run = run + 1
# N # of exceedances | K # of blocks with exceedances:
N = length(X[X > u])
K = floor(sum(sign(apply(X, 1, max) - u) + 1) / 2)
if (K/k < 1) {
theta1[run] = (k/n) * log(1-K/k) / log(1-N/n)
} else {
theta1[run] = NA
}
keepK = c(keepK, K)
keepN = c(keepN, N)
}
# Theta Values:
ans = data.frame(quantiles = quantiles, thresholds = thresholds,
N = keepN, K = keepK, theta = theta1)
# Add title and description:
if (is.null(title)) title = "Extremal Index from Block Method"
if (is.null(description)) description = description()
# Return Value:
new("fTHETA",
call = match.call(),
data = list(x = x, block = block),
theta = ans,
title = title,
description = description)
}
# ------------------------------------------------------------------------------
clusterTheta =
function (x, block = 22, quantiles = seq(0.950, 0.995, length = 10),
title = NULL, description = NULL)
{ # A function implemented by Diethelm Wuertz
# Description:
# Calculates theta from Reciprocal Mean Cluster Size method, i.e. theta2.
# Example:
# clusterTheta(thetaSim(n=10000))
# FUNCTION:
# Check if block is numeric:
stopifnot(is.numeric(block))
# Number of blocks and number of data points:
X = as.vector(x)
ordered = sort(X)
k = floor(length(X)/block)
n = k*block
# Now organize your X:
# 1) truncate the rest of the time series,
# 2) arrange them in matrix form,
# 3) sort them in reverse order, ie. from high (pos) to low (neg)
X = matrix(X[1:(k*block)], ncol = block, byrow = TRUE)
# Threshold values associated to quantiles:
thresholds = ordered[floor(quantiles*length(X))]
# Presettings:
theta2 = rep(0, times = length(quantiles))
# Calculate Extremal Imdex:
run = 0
keepK = keepN = NULL
for ( u in thresholds ) {
run = run + 1
# N # of exceedances | K # of blocks with exceedances:
N = length(X[X > u])
K = floor(sum(sign(apply(X, 1, max) - u) + 1) / 2)
theta2[run] = K/N
keepK = c(keepK, K)
keepN = c(keepN, N)
}
# Theta Values:
ans = data.frame(quantiles = quantiles, thresholds = thresholds,
N = keepN, K = keepK, theta = theta2)
# Add title and description:
if (is.null(title))
title = "Extremal Index from Reciprocal Cluster Method"
if (is.null(description)) description = description()
# Return Value:
new("fTHETA",
call = match.call(),
data = list(x = x, block = block),
theta = ans,
title = title,
description = description)
}
# ------------------------------------------------------------------------------
runTheta =
function (x, block = 22, quantiles = seq(0.950, 0.995, length = 10),
title = NULL, description = NULL)
{ # A function implemented by Diethelm Wuertz
# Description:
# Calculates theta from Run method, i.e. theta3.
# Example:
# runTheta(thetaSim(n=10000))
# FUNCTION:
# Check if block is numeric:
stopifnot(is.numeric(block))
# Number of blocks and number of data points:
X = as.vector(x)
ordered = sort(X)
k = floor(length(X)/block)
n = k*block
Count = 1:n
# Now organize your X:
# 1) truncate the rest of the time series,
# 2) arrange them in matrix form,
# 3) sort them in reverse order, ie. from high (pos) to low (neg)
X = matrix(X[1:(k*block)], ncol = block, byrow = TRUE)
# Threshold values associated to quantiles:
thresholds = ordered[floor(quantiles*length(X))]
# Presettings:
theta3 = rep(0, times = length(quantiles))
# Calculate Extremal Imdex:
run = 0
keepN = NULL
for ( u in thresholds ) {
run = run + 1
# N # of exceedances | K # of blocks with exceedances:
N = length(X[X > u])
Y = diff(Count[X > u])
Y = Y[Y > block]
theta3[run] = length(Y)/N
keepN = c(keepN, N)
}
# Theta Values:
ans = data.frame(quantiles = quantiles, thresholds = thresholds,
N = keepN, theta = theta3)
# Add title and description:
if (is.null(title))
title = "Extremal Index from Run Method"
if (is.null(description)) description = description()
# Return Value:
new("fTHETA",
call = match.call(),
data = list(x = x, block = block),
theta = ans,
title = title,
description = description)
}
# ------------------------------------------------------------------------------
ferrosegersTheta =
function (x, quantiles = seq(0.950, 0.995, length= 10),
title = NULL, description = NULL)
{
# Description:
# Estimates the extremal index based on the intervals estimator
# due to Ferro and Segers (2003).
# Note:
# Adapted from function 'extremalindex' in contributed R-package
# 'extRemes' written and maintained by ...
# FUNCTION:
# Settings:
x = as.vector(x)
n = length(x)
N = floor(quantiles*n)
sorted = sort(x)
U = sorted[N]
ans = NULL
# Extremal Index:
for ( u in U ) {
msg = 0
id = x > u
N = sum(id)
S = (1:n)[id]
TT = diff(S)
if (!any(TT > 2)) {
theta = 2*sum(TT, na.rm = TRUE)^2/((N-1) * sum(TT^2, na.rm = TRUE))
# msg = paste("theta.hat used because no values of T>2.")
msg = msg + 1
if (theta > 1) {
theta = 1
# msg = paste(msg, "Using theta = 1 because theta.hat > 1.",
# sep = "\n")
msg = msg + 10
}
} else {
theta = 2 * sum(TT-1, na.rm = TRUE)^2/((N-1) * sum((TT-1) *
(TT-2), na.rm = TRUE))
# msg = paste("theta.tilde used because a value(s) exists of T>2.")
msg = msg + 100
if (theta > 1) {
theta = 1
# msg = paste(msg, "Using theta = 1 as theta.hat > 1.")
msg = msg + 1000
}
}
K = ifelse(round(theta*N) != theta*N, floor(theta*N) + 1, theta*N)
T.order = order(TT, na.last = TRUE, decreasing = TRUE)
T.ranked = TT[T.order]
T.K = T.ranked[K]
if (sum(TT == T.K, na.rm = TRUE) > 1) {
for (i in 1:K) {
K = K - 1
T.K = T.ranked[K]
if (sum(TT == T.K, na.rm = TRUE) > 1) {
next
} else {
break
}
}
}
ans = rbind(ans, c(T.K, K, msg, theta))
}
# Result:
ans = data.frame(quantiles, U, ans)
colnames(ans) = c("Threshold", "Quantiles",
"RunLength", "Clusters", "messageNo", "theta")
# Add title and description:
if (is.null(title))
title = "Extremal Index from Ferro-Segers Method"
if (is.null(description)) description = description()
# Return Value:
new("fTHETA",
call = match.call(),
data = list(x = x),
theta = ans,
title = title,
description = description)
}
################################################################################
exindexesPlot =
function (x, block = 22, quantiles = seq(0.950, 0.995, length = 10),
doplot = TRUE, labels = TRUE, ...)
{ # A function implemented by Diethelm Wuertz
# Description:
# Calculates and Plots Theta(1,2,3) for numeric block lengths
# Areguments:
# x - an univariate time series, or any other object which can be
# transformed by the function as.vector into a numeric vector.
# block - an integer value which denotes the length of the blocks.
# quantiles - a numeric vector of quantile values.
# doplot - alogical flag. Should a plot be produced?
# Example:
# exindexesPlot(as.timeSeries(data(bmwRet)), 20)
# FUNCTION:
# Settings:
if (!is.numeric(block)) stop("Argument block must be an integer value.")
doprint = FALSE
# Block Size:
blocklength = block # argument renamed
# Note, in finance the x's should be residuals
resid = as.vector(x)
# Extremal Index - Theta_1, Theta_2 and Theta_3
k = floor(length(resid)/blocklength) # Number of blocks
n = k*blocklength # Number of data points
# Now organize your residuels:
# 1) truncate the rest of the time series,
# 2) arrange them in matrix form,
# 3) sort them in reverse order, ie. from high (pos) to low (neg)
resid1 = resid[1:(k*blocklength)]
resid1 = matrix(resid1, ncol = blocklength, byrow = TRUE)
ordered1 = sort(resid1)
# Threshold values associated to quantiles:
z0 = ordered1[floor(quantiles*length(resid1))]
# Presettings:
theta1 = theta2 = theta3 = rep(0, times = length(quantiles))
# Calculate Extremal Imdex:
run = 0
for ( z in z0 ) {
run = run + 1
# N - number of exceedances:
N = length(resid1[resid1 > z])
# K - number of blocks with exceedances:
# DW: floor()
K = floor(sum(sign(apply(resid1, 1, max)-z)+1) / 2)
if (K/k < 1) {
theta1[run] = (k/n) * log(1-K/k) / log(1-N/n)
} else {
theta1[run] = NA
}
theta2[run] = K/N
x = 1:n
xx = diff(x[resid1 > z])
xx = xx[xx > blocklength]
theta3[run] = length(xx)/N
# Printing:
if (doprint) {
print(c(N, K, quantiles[run], z))
print(c(theta1[run], theta2[run], theta3[run]))
}
}
# Plotting:
if (doplot) {
plot(quantiles, theta1,
xlim = c(quantiles[1], quantiles[length(quantiles)]),
ylim = c(0, 1.2), type = "b", pch = 1,
xlab = "", ylab = "", main = "", ...)
points(quantiles, theta2, pch = 2, col = 3)
points(quantiles, theta3, pch = 4, col = 4)
if (labels) {
title(main = "Extremal Index")
title(xlab = "Quantile", ylab = "Theta 1,2,3")
mtext("Threshold", side = 3, line = 3)
grid()
mtext(text = paste("Blocklength: ", as.character(block)),
adj = 0, side = 4, cex = 0.7)
}
}
# Theta Values:
ans = data.frame(quantiles = quantiles, thresholds = z0,
theta1 = theta1, theta2 = theta2, theta3 = theta3)
# Return Value:
ans
}
# ------------------------------------------------------------------------------
exindexPlot =
function(x, block = c("monthly", "quarterly"), start = 5, end = NA,
doplot = TRUE, plottype = c("thresh", "K"), labels = TRUE, ...)
{
# Example:
# exindexPlot(as.timeSeries(data(bmwRet)), 20)
# exindexPlot(as.timeSeries(data(bmwRet)), "monthly")
# exindexPlot(as.vector(as.timeSeries(data(bmwRet))), 20)
# Settings:
plottype = match.arg(plottype)
reverse = FALSE
if (plottype == "K") reverse = TRUE
# Extremal Index - following A. McNeil:
b.maxima = rev(sort(as.vector(blockMaxima(x, block))))
data = as.vector(x)
sorted = rev(sort(data))
n = length(sorted)
if (!is.numeric(block)) block = round(length(data)/length(b.maxima))
k = round(n/block)
un = unique(b.maxima)[-1]
K = match(un, b.maxima) - 1
N = match(un, sorted) - 1
if (is.na(end)) end = k
cond = (K < end) & (K >= start)
un = un[cond]
K = K[cond]
N = N[cond]
theta2 = K/N
theta = logb(1 - K/k)/(block * logb(1 - N/n))
ans = data.frame(N = N, K = K, un = un, theta2 = theta2, theta = theta)
yrange = range(theta)
index = K
if (reverse) index = - K
# Plot:
if (doplot) {
plot(index, theta, ylim = yrange, type = "b", xlab = "", ylab = "",
axes = FALSE, ...)
IDX = round(seq(1, length(index), length = 10))
axis(1, at = index[IDX], labels = paste(K)[IDX])
axis(2)
axis(3, at = index[IDX], labels = paste(format(signif(un, 3)))[IDX])
box()
if (labels) {
ylabel =
paste("theta (", k, " blocks of size ", block, ")", sep = "")
title(xlab = "K", ylab = ylabel)
mtext("Threshold", side = 3, line = 3)
lines(index, theta, col = "steelblue")
grid()
mtext(text = paste("Blocklength: ", as.character(block)),
adj = 0, side = 4, cex = 0.7)
}
}
# Return Value:
ans
}
################################################################################
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.