#' @importFrom stats dnorm lm optimize pnorm qnorm sd
#' @importFrom graphics abline axis legend lines matplot plot polygon
#' @export murphydiagram murphydiagram_diff expected_score_mean expected_score_quantile fluctuation_test extremal_score apl_score ase_score
# Input checker needed below
input.check <- function(f, y, t, alpha){
if ( (length(t) > 1) | (length(alpha) > 1) | (length(f) != length(y)) | any(c(!is.vector(f), !is.vector(y), !is.vector(t), !is.vector(alpha))) ) stop("invalid input")
}
# Extremal score for quantiles
S.quantile <- function(f, y, t, alpha){
input.check(f, y, t, alpha)
((y < f) - alpha) * ((t < f) - (t < y))
}
# Take the positive part of a vector
pos <- function(x) x*(x >= 0)
# Extremal score for expectiles
S.expectile <- function(f, y, t, alpha){
input.check(f, y, t, alpha)
c1 <- abs((y < f) - alpha)
c2 <- pos(y-t) - pos(f-t) - (y-f)*(t<f)
out <- c1 * c2
out
}
# Newey-West VCV matrix
# u: vector of data
# prewhite: logical, should prewhitening be done?
# k: truncation lag for autocorrelations. If set to NULL, will be chosen automatically.
# meth: either "qs" (Quadratic Spectral, Andrews 1991) or anything else (Bartlett Kernel, Newey/West)
vHAC <- function(u, prewhite = FALSE, k = NULL, meth = "qs"){
if (!is.matrix(u)) u <- as.matrix(u)
n <- nrow(u)
nreg <- ncol(u)
rho <- sigma <- matrix(0, nreg, 1)
# do a VAR(1) prewhitening
if (prewhite == TRUE){
reg.x <- matrix(u[1:(n-1), ], n-1, nreg)
reg.y <- matrix(u[2:n, ], n-1, nreg)
aux <- lm(reg.y~reg.x-1)
beta <- matrix(unname(aux$coefficients), nreg, nreg)
v <- matrix(unname(aux$residuals), n-1, nreg)
} else {
v <- u
beta <- matrix(0, nreg, nreg)
}
nv <- nrow(v)
# choose nr of lags (if not provided)
if (is.null(k)){
for (i in 1:nreg){
aux <- lm(v[2:nv, i]~v[1:(nv-1), i]-1)
rho[i] <- unname(aux$coefficients)
sigma[i] <- sum(unname(aux$residuals)^2) / nv
}
if (meth == "qs"){
# See Eq. (6.4) on page 835 of Andrews (1991) -> Note that his sigma^2 corresponds to our sigma
top <- sum( (4*(rho^2) * (sigma^2)) / ((1-rho)^8) )
bot <- sum( (sigma^2) / ((1-rho)^4) )
k <- ceiling(1.3221*((top/bot)*n)^(0.2))
} else {
top <- sum( (4*(rho^2) * (sigma^2)) / (((1-rho)^6)*((1+rho)^2)) )
bot <- sum( (sigma^2) / ((1-rho)^4) )
k <- min(c(ceiling(1.1447*((top/bot)*n)^(1/3)), round(0.5*n)))
}
}
# compute HAC
vcv <- (t(v) %*% v) / (n-1)
if (k > 0){
if (meth == "qs"){
del <- ((6 * pi)/(5 * k)) * (1:(n-1))
w <- 3 * (sin(del) / del - cos(del)) / (del^2)
if (prewhite == FALSE){
mlag <- n - 1
} else {
mlag <- nv - 1
}
} else {
w <- 1 - (1:k)/(k+1)
mlag <- k
}
for (i in 1:mlag){
cov <- t(v[(i+1):nv, , drop = FALSE]) %*% (v[1:(nv-i), , drop = FALSE]) / (n-1)
vcv <- vcv + w[i]*(cov + t(cov))
}
}
d <- solve(diag(nreg) - t(beta))
hac <- d %*% vcv %*% t(d)
return(list(hac = hac, k = k))
}
get_grid <- function(f1, f2, y, functional, alpha, literal = TRUE, lhdiff = 1e-10){
if (literal == TRUE){
# Choose range of t's (see Corollary 2a, 2b in paper)
if (functional == "quantile"){
tseq <- sort(unique(c(f1, f2, y)))
} else if (functional == "expectile") {
# NOTE: Use y's for Murphy diagrams also in case alpha = 1/2, even if not strictly needed
aux1 <- sort(unique(c(f1, f2, y)))
aux2 <- sort(unique(c(f1, f2))) - lhdiff
tseq <- sort(unique(c(aux1, aux2)))
}
} else {
# Simply choose equally spaced grid of fixed length
aux <- c(f1, f2, y)
tseq <- seq(from = min(aux) - 0.1*sd(aux), to = max(aux) + 0.1*sd(aux), length.out = 100)
}
tseq
}
murphydiagram <- function(f1, f2, y, functional = "expectile", alpha = 0.5, labels = c("Method 1", "Method 2"),
colors = NULL, equally_spaced = FALSE){
cex.gen <- 1.6
# Define function for extremal score
if (functional == "expectile"){
g <- function(f, t) S.expectile(f, y, t, alpha)
} else if (functional == "quantile"){
g <- function(f, t) S.quantile(f, y, t, alpha)
} else {
stop("Please choose either expectile or quantile functional")
}
# Unless specified otherwise: Use colors as in paper
if (is.null(colors)) colors <- c("#D55E00", "#56B4E9", "#000000")
# Grid of theta values
tseq <- get_grid(f1, f2, y, functional, alpha, literal = 1 - equally_spaced)
# Data frame with score entries for all theta values
df <- data.frame(tseq = tseq, s1 = numeric(length(tseq)), s2 = numeric(length(tseq)))
for (j in 1:length(tseq)){
aux1 <- g(f1, tseq[j])
aux2 <- g(f2, tseq[j])
df[j, 2:3] <- c(mean(aux1), mean(aux2))
}
# Plot: Scores for both methods
if (all(y %in% c(0, 1))){
xx <- c(-0.05, 1.05)
} else {
xx <- c(min(tseq) - 0.1, max(tseq) + 0.1)
}
matplot(x = tseq, y = df[,2:3], type = "l", lty = 1, lwd = 4, xlab = expression(paste("Parameter ", theta)), ylab = "", bty = "n", cex.lab = cex.gen,
cex.axis = cex.gen, xlim = xx, ylim = c(0, 1.2*max(df[,2:3])), col = colors)
abline(h = 0, lty = 2)
if (!is.null(labels)) legend("top", labels, col = colors, lwd = 4, bty = "n", horiz = TRUE, cex = 1.2)
}
murphydiagram_diff <- function(f1, f2, y, functional = "expectile", alpha = 0.5, equally_spaced = FALSE, lag_truncate = 0, conf_level = 0.95){
cex.gen <- 1.6
# Some variables
nobs <- length(y)
scl <- abs(qnorm(0.5*(1-conf_level)))
# Define function
if (functional == "expectile"){
g <- function(f, t) S.expectile(f, y, t, alpha)
} else if (functional == "quantile"){
g <- function(f, t) S.quantile(f, y, t, alpha)
} else {
stop("Please choose either expectile or quantile functional")
}
# Choose range of t's
tseq <- get_grid(f1, f2, y, functional, alpha, literal = 1 - equally_spaced)
df <- data.frame(tseq = tseq, s1 = numeric(length(tseq)), s2 = numeric(length(tseq)), lb = numeric(length(tseq)), ub = numeric(length(tseq)))
for (j in 1:length(tseq)){
aux1 <- g(f1, tseq[j])
aux2 <- g(f2, tseq[j])
df[j, 2:3] <- c(mean(aux1), mean(aux2))
aux.v <- try(vHAC(aux1 - aux2, k = lag_truncate, meth = "bartlett")$hac/nobs)
# HAC estimator won't invert for some t in the tails -> Set variance to zero in these cases
if (inherits(aux.v, "try-error")){
aux.v <- 0
}
aux.m <- mean(aux1-aux2)
df[j, 4:5] <- c(aux.m - scl*sqrt(aux.v), aux.m + scl*sqrt(aux.v))
}
# Plot: Score difference + confidence interval
if (all(y %in% c(0, 1))){
xx <- c(-0.05, 1.05)
} else {
xx <- c(min(tseq) - 0.1, max(tseq) + 0.1)
}
matplot(x = df$tseq, y = df[, 4:5], type = "n", ylab = "", xlab = expression(paste("Parameter ", theta)),
bty = "n", cex.axis = cex.gen, cex.lab = cex.gen, xlim = xx, col = 1)
polygon(c(df$t, rev(df$t)), c(df$ub, rev(df$lb)), col = "grey", border = NA)
lines(x = df$t, y = (df$s1-df$s2), type = "l", col = 1, lwd = 2.5)
abline(h = 0, lty = 2)
}
##############################################################################
#
# Helper functions and analytical expressions for simulation study
#
##############################################################################
# cdf and quantile function for unfocused forecaster (forecaster 3)
p3p <- function(q, tau) 0.5*(pnorm(q) + pnorm(q - tau))
p3m <- function(q, tau) 0.5*(pnorm(q) + pnorm(q + tau))
q3p <- function(p, tau) {
if (p == 0) {
return(-Inf)
} else if (p == 1) {
return(Inf)
} else {
optimize(f = function(x) abs(p - p3p(x, tau)), interval = c(0, tau) + qnorm(p))$minimum
}
}
q3m <- function(p, tau) {
if (p == 0) {
return(-Inf)
} else if (p == 1) {
return(Inf)
} else {
optimize(f = function(x) abs(p - p3m(x, tau)), interval = c(0, -tau) + qnorm(p))$minimum
}
}
# analytic expected scores for mean functional
c_theta <- function(t) sqrt(2) * dnorm(t/sqrt(2)) + t * pnorm(t/sqrt(2))
aux_mean <- list()
aux_mean[[1]] <- function(t) 0.5*(c_theta(t) - t * pnorm(t) - dnorm(t)) # perfect (forecaster 1)
aux_mean[[2]] <- function(t) 0.5*(c_theta(t) - t * (t >= 0)) # climatological (forecaster 2)
aux_mean[[3]] <- function(t) { # unfocused (forecaster 3)
tau <- 2
c1 <- c_theta(t)
c2 <- t*(pnorm(t-0.5*tau) + pnorm(t+0.5*tau)) + dnorm(t-0.5*tau) + dnorm(t+0.5*tau)
return(0.5*(c1 - 0.5 * c2))
}
aux_mean[[4]] <- function(t) 0.5*(c_theta(t) - t * pnorm(t) + dnorm(t)) # sign-reversed (forecaster 4)
expected_score_mean <- function(theta, forecaster = "P"){
forecaster_names <- c("P", "C", "U", "SR")
if (! forecaster %in% forecaster_names) stop("Forecaster must be either P, C, U, or SR")
ind <- which(forecaster_names == forecaster)
aux_mean[[ind]](t = theta)
}
# analytic expected scores for quantile / prob functional
aux_quant <- list()
aux_quant[[1]] <- function(t, a) { # perfect (forecaster 1)
bound <- t - qnorm(a)
if (is.finite(bound)) {
c1 <- -a * pnorm(t / sqrt(2))
c2 <- a * pnorm(bound)
c3 <- integrate(function(x) pnorm(t-x) * dnorm(x), bound, Inf)$value
return(c1 + c2 + c3)
} else {
return(0)
}
}
aux_quant[[2]] <- function(t, a) { # climatological (forecaster 2)
x <- pnorm(t/sqrt(2))
if (x < a) {
return(x * (1 - a))
} else {
return(a * (1 - x))
}
}
aux_quant[[3]] <- function(t, a) { # unfocused (forecaster 3)
tau <- 2
bound1 <- t - q3p(a, tau)
bound2 <- t - q3m(a, tau)
if (is.finite(bound1)) {
c1 <- -a * pnorm(t / sqrt(2))
c2 <- a * (pnorm(bound1) + pnorm(bound2))
c31 <- integrate(function(x) pnorm(t-x)*dnorm(x), bound1, Inf)$value
c32 <- integrate(function(x) pnorm(t-x)*dnorm(x), bound2, Inf)$value
return(c1 + 0.5 * (c2 + c31 + c32))
} else {
return(0)
}
}
aux_quant[[4]] <- function(t, a) { # sign-reversed (forecaster 4)
bound <- qnorm(a) - t
if (is.finite(bound)) {
c1 <- -a * pnorm(t / sqrt(2))
c2 <- a * pnorm(-bound)
c3 <- integrate(function(x) pnorm(t-x)*dnorm(x), -Inf, bound)$value
return(c1 + c2 + c3)
} else {
return(0)
}
}
expected_score_quantile <- function(theta, alpha, forecaster = "P") {
forecaster_names <- c("P", "C", "U", "SR")
if (! forecaster %in% forecaster_names) stop("Forecaster must be either P, C, U, or SR")
if (any(alpha < 0 | alpha > 1)) stop("Please provide quantile levels between zero and one")
ind <- which(forecaster_names == forecaster)
E <- function(theta, alpha) aux_quant[[ind]](t = theta, a = alpha)
if (length(alpha) >= 1 & length(theta) == 1) {
out <- sapply(alpha, E, t = theta)
} else if (length(theta) > 1 & length(alpha) == 1) {
out <- sapply(theta, E, a = alpha)
} else {
out <- NA
}
return(out)
}
fluctuation_test <- function(loss1, loss2, mu = 0.5, dmv_fullsample = TRUE,
lag_truncate = 0, time_labels = NULL,
conf_level = 0.05){
# Some input checks
if (length(loss1) != length(loss2)){
stop("Vectors of losses must have the same length")
}
if (all(abs(seq(from = 0.1, to = 0.9, by = 0.1) - mu) > 1e-12)){
stop("mu must be in {0.1, 0.2, ..., 0.9}")
}
if (!lag_truncate %in% 0:5){
stop("lag_truncate must be in {0, 1, .., 5}")
}
if (!is.null(time_labels) & length(time_labels) != length(loss1)){
warning("Specified time labels are inconsistent - simple integers used instead")
time_labels <- NULL
}
if (!conf_level %in% c(0.05, 0.1)){
stop("significance_level must be either 0.05 or 0.1")
}
# Critical values (two-sided test, from Table 1 of Giacomini/Rossi 2010)
# Second/third column correspond to 5/10 percent confidence level
CV <- cbind(seq(from = 0.1, to = 0.9, by = 0.1),
c(3.393, 3.179, 3.012, 2.890, 2.779, 2.634, 2.560, 2.433, 2.248),
c(3.170, 2.948, 2.766, 2.626, 2.500, 2.356, 2.252, 2.130, 1.950))
# Wrapper for HAC variance
vHAC2 <- function(ld, lag_truncate){
vHAC(ld, k = lag_truncate, meth = "Bartlett")$hac
}
# Size for fonts etc
cex_gen <- 1.6
# Window length m
P <- length(loss1) # Notation as in GR2010
m <- round(mu*P)
# Loss differences
ld <- loss1 - loss2
# Vector of statistics
dm_num <- dm_den <- rep(0, P - m + 1)
# Loop over time
for (jj in m:P){
ind <- which(m:P == jj)
ld_tmp <- ld[(jj-m+1):jj]
dm_num[ind] <- mean(ld_tmp)
dm_den[ind] <- sqrt(vHAC2(ld_tmp, lag_truncate)/m)
}
# Construct statistics
# Variant 1: Use full sample for HAC variance
s2hat <- vHAC2(ld, lag_truncate)
dm1 <- sqrt(m)*dm_num/sqrt(s2hat)
# Variant 2: Use rolling samples
dm2 <- dm_num/dm_den
# Choose statistic
if (dmv_fullsample){
dm_final <- dm1
} else {
dm_final <- dm2
}
# Select critical value
if (conf_level == 0.05){
CVs <- CV[abs(CV[, 1] - mu) < 1e-12, 2] * c(-1, 1)
} else if (conf_level == 0.1){
CVs <- CV[abs(CV[, 1] - mu) < 1e-12, 3] * c(-1, 1)
}
# Plot
plot(x = m:length(loss1), y = dm_final, ylim = c(-7, 7), bty = "n", ylab = "",
xlab = "Time (End of Rolling Window)",
type = "l", col = "cornflowerblue", lwd = 2.5, axes = FALSE, cex.lab = cex_gen)
axis(2, cex.axis = cex_gen)
abline(h = CVs, lwd = 3.5)
abline(h = 0, lwd = 1.8, lty = 2)
if (is.null(time_labels)){
time_labels <- 1:length(loss1)
}
inds <- floor(seq(from = m, to = length(time_labels), length.out = 5))
axis(1, at = inds, labels = time_labels[inds], cex.axis = cex_gen)
# Return two-element list which summarizes the results
list(df = data.frame(time = time_labels[m:length(loss1)], dmstat = dm_final), CV = CVs)
}
# Wrapper function for extremal scores
extremal_score <- function(x, y, theta, functional = "expectile", alpha = 0.5){
if (functional == "expectile"){
S.expectile(f = x, y = y, t = theta, alpha = alpha)
} else if (functional == "quantile"){
S.quantile(f = x, y = y, t = theta, alpha = alpha)
} else {
stop("invalid input for functional")
}
}
input.check2 <- function(x, y){
if (length(y) != length(x) | !is.vector(x) | !is.vector(y)) stop("x and y must be vectors of the same length")
}
apl_score <- function(x, y, alpha = 0.5){
input.check2(x, y)
out <- (x - y)
ind <- (out > 0)
out[ind] <- (1-alpha) * out[ind]
out[!ind] <- -alpha * out[!ind]
out
}
ase_score <- function(x, y, alpha = 0.5){
input.check2(x, y)
out <- (x - y)
ind <- (out > 0)
out[ind] <- (1-alpha) * (out[ind]^2)
out[!ind] <- alpha * (out[!ind]^2)
out
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.