# R/getEconomicValue.R In matteodefelice/eneaR:

#' @description Compute Economic Value following Richardson definitions. For further details cfr Richardson, Economic Value & Skills, 2003
#! @param pp A vector containing probabilistic forecast
#! @param pp_ref A vector of an observed binary event (0,1)
#! @param COST the vector of the COST/LOSS ratio to consider. The default is seq(0.1, 0.9, 0.05)
#! @param method NA or 'bootstrap' to compute a bootstrap for the significativity
#' @export
#' @author M. De Felice
getEconomicValue = function(pp, pp_ref, COST = c(), method = NA, NREP = 100) {

if (length(COST) == 0) {
SCALE = seq(0.1, 0.9, by = 0.05)
} else {
SCALE = COST
}

N = length(pp)
pp_ref = as.logical(pp_ref) #
V = rep(NA, length(SCALE))
# Discretizing probabilistic forecast
# Assuming that optimal disc. threshold is equal to C/L -> C
for (i in 1:length(SCALE)) {
pp_disc = (pp > SCALE[i])
# Setting L (Loss) as 1 we explore various values of C (cost)

C = SCALE[i]
s = mean(pp_ref * 1)
Eperfect = s*C
Eclimate = min(C, s)

a = sum(1*(pp_disc & pp_ref)) / N
b = sum(1*(pp_disc & !pp_ref)) / N
c = sum(1*(!pp_disc & pp_ref)) / N
d = sum(1*(!pp_disc & !pp_ref)) / N
Eforecast = a*C +b*C+c # Eq. (8.6)

V[i] = (Eclimate - Eforecast)/(Eclimate - Eperfect)
# print(sprintf('C: %2.2f: Eclim: %2.2f Eforec: %2.2f Eperf: %2.2f',
#                   C, Eclimate, Eforecast, Eperfect))
}
# Bootstrap
if (!is.na(method)) {
bsV = matrix(NA, nrow = length(SCALE), ncol = NREP)
for (n in 1:NREP) {
for (i in 1:length(SCALE)) {
bspp = sample(pp, length(pp), replace = T)
pp_disc = (bspp > SCALE[i])
# Setting L (Loss) as 1 we explore various values of C (cost)
C = SCALE[i]
s = mean(pp_ref * 1)
Eperfect = s*C
Eclimate = min(C, s)

a = sum(1*(pp_disc & pp_ref)) / N
b = sum(1*(pp_disc & !pp_ref)) / N
c = sum(1*(!pp_disc & pp_ref)) / N
d = sum(1*(!pp_disc & !pp_ref)) / N
Eforecast = a*C +b*C+c # Eq. (8.6)

bsV[i,n] = (Eclimate - Eforecast)/(Eclimate - Eperfect)
}
}
# Compute significativity
QQs = c(0.01, 0.025, 0.05)
SIGs = rep(NA, length(SCALE))
for (l in 1:length(SCALE)) {
for (q in QQs) {
qq = quantile(bsV[l,], c(q, 1-q))
if (V[i] <= qq[1] || V[i]>=qq[2]) {
SIGs[l] = q
break
}
}
}
# If bootstrap is enable this function exits here
return(list(x = SCALE, v = V, SIG = SIGs))
} else {
return(list(x = SCALE, v = V))
}

}
matteodefelice/eneaR documentation built on May 21, 2019, 12:59 p.m.