Nothing
setClass("SSM", representation(Info = "list", Bootstrap = "list",
show = "matrix"))
setMethod("show", "SSM", function(object) {
cat(paste("\n------------------------------------------"))
cat(paste("\n- Superior Set of Models -"))
cat(paste("\n------------------------------------------\n"))
show(object@show)
cat(paste("\nDetails"))
cat(paste("\n------------------------------------------\n"))
cat(paste("\nNumber of eliminated models\t:\t"))
cat(object@Info$n_elim)
cat(paste("\nStatistic\t:\t"))
cat(object@Info$statistic)
cat(paste("\nElapsed Time\t:\t"))
print(object@Info$elapsed.time)
})
boot.block <- function(x, v, n, k) {
startIndexes = sample(1:(n - k), v + 1)
blocks = do.call(c, lapply(startIndexes, function(p, k) p +
seq_len(k), k = k))
return(blocks[1:n])
}
d_b_i_mean.fun <- function(x, M, model.names) {
export = sapply(1:M, function(l) sum(x[paste(model.names[l],
model.names[-l], sep = ".")])/(M - 1))
return(export)
}
LossVaR <- function(realized, evaluated, which = "asymmetricLoss",
type = "normal", delta = 25, tau) {
if (delta <= 0) {
stop("delta must be greater than zero")
}
if (which != "asymmetricLoss") {
stop("Not valid choice of \"which\". Valid choices are \"asymmetricLoss\"")
}
if (which == "asymmetricLoss" & all(type != c("normal", "differentiable"))) {
stop("Not valid choice of \"type\". Valid choices are \"normal\" and \"differentiable\" ")
}
if (tau < 0 || tau > 1) {
stop("tau must be in (0,1)")
}
if (is.matrix(evaluated)) {
m = ncol(evaluated)
} else {
m = 1
}
if (!is.matrix(evaluated)) {
evaluated = as.matrix(evaluated)
}
if (!is.matrix(realized)) {
realized = as.matrix(realized)
}
if (ncol(realized) > 1) {
stop("realized must of one column")
}
if (nrow(realized) != nrow(evaluated)) {
stop("the number of realized and evaluated observation must be the same")
}
if (which == "asymmetricLoss") {
if (type == "normal") {
dummies = matrix(0, nrow = length(realized), ncol = m)
for (j in 1:m) {
dummies[realized < evaluated[, j], j] = 1
}
loss = (tau - dummies) * (realized - evaluated)
}
if (type == "differentiable") {
loss = (tau - (1 + exp(delta * (realized - evaluated)))^(-1)) *
(realized - evaluated)
}
}
return(loss)
}
LossVol <- function(realized, evaluated, which = "SE1") {
if (!which %in% c("SE1", "SE2", "QLIKE", "R2LOG", "AE1",
"AE2")) {
stop("Not valid choice of \"which\". Valid choices are \"SE1\",\"SE2\",\"QLIKE\",\"R2LOG\",\"AE1\",\"AE2\" ")
}
if (is.matrix(evaluated)) {
m = ncol(evaluated)
} else {
m = 1
}
if (!is.matrix(evaluated)) {
evaluated = as.matrix(evaluated)
}
if (!is.matrix(realized)) {
realized = as.matrix(realized)
}
if (ncol(realized) > 1) {
stop("realized must of one column")
}
if (nrow(realized) != nrow(evaluated)) {
stop("the number of realized and evaluated observation must be the same")
}
if (which == "SE1") {
loss = (evaluated - realized)^2
}
if (which == "SE2") {
loss = (evaluated^2 - realized^2)^2
}
if (which == "QLIKE") {
loss = log(evaluated^2) + realized^2 * evaluated^(-2)
}
if (which == "R2LOG") {
loss = (log(realized^2 * evaluated^-2))^2
}
if (which == "AE1") {
loss = abs(evaluated - realized)
}
if (which == "AE2") {
loss = abs(evaluated^2 - realized^2)
}
return(loss)
}
LossLevel <- function(realized, evaluated, which = "SE") {
if (all(which != c("SE", "AE")))
stop("Not valid choice of \"which\". Valid choices are \"SE\",\"AE\" ")
if (is.matrix(evaluated)) {
m = ncol(evaluated)
} else {
m = 1
}
if (!is.matrix(evaluated)) {
evaluated = as.matrix(evaluated)
}
if (!is.matrix(realized)) {
realized = as.matrix(realized)
}
if (ncol(realized) > 1)
stop("realized must of one column")
if (nrow(realized) != nrow(evaluated))
stop("the number of realized and evaluated observation must be the same")
if (which == "SE") {
loss = (evaluated - realized)^2
}
if (which == "AE") {
loss = abs(evaluated - realized)
}
return(loss)
}
MCSprocedure <- function(Loss, alpha = 0.15, B = 5000, cl = NULL,
ram.allocation = TRUE, statistic = "Tmax", k = NULL, min.k = 3,
verbose = TRUE) {
time.start = Sys.time()
Loss = as.matrix(Loss)
M_start = ncol(Loss)
# no points in colnames(Loss)
colnames(Loss) = gsub(".", "_", colnames(Loss), fixed = T)
if (is.null(colnames(Loss))) {
colnames(Loss) = paste("model", 1:ncol(Loss), sep = "_")
}
if (is.null(cl)) {
cl_number = 0
}
if (!is.null(cl)) {
max.cores = length(cl)
}
B = round(B)
if (B < 1000) {
cat(paste("Warning: B is small"))
}
if (any(is.na(Loss))) {
stop("NAs in Loss are not allowed")
}
if (any(abs(Loss) == Inf)) {
stop("Inf in Loss are not allowed")
}
N = nrow(Loss)
repeat {
M = ncol(Loss)
model.names = colnames(Loss)
col.names.d = do.call(c,
lapply(1:M, function(x) {
paste(model.names[x], model.names[-x], sep = ".")
}))
d = do.call(cbind,
lapply(1:M, function(x) {
Loss[, x] - Loss[, -x]
}))
colnames(d) = col.names.d
d_ij_mean = colMeans(d)
d_i_mean = sapply(1:M, function(x) {
sum(d_ij_mean[paste(model.names[x], model.names[-x], sep = ".")])/(M - 1)
})
names(d_i_mean) = model.names
block.length = NULL
foo = expand.grid(1:M, 1:M)[, 2:1]
foo = foo[foo[, 1] != foo[, 2], ]
index = col.names.d[foo[, 1] < foo[, 2]]
if (ncol(d) > 2) {
d_cut = as.list(as.data.frame(d[, index]))
if (!is.null(cl)) {
k = max(na.omit(as.numeric(parSapply(cl, d_cut, function(x) {
try(ar(x)$order, silent = TRUE)
}))))
}
if (is.null(cl)) {
k = max(na.omit(as.numeric(sapply(d_cut, function(x) {
try(ar(x)$order, silent = TRUE)
}))))
}
} else {
k = ar(d[, 1])$order
}
if (k < min.k) {
k = min.k
}
v = ceiling(nrow(d)/k)
n = nrow(d)
if (!is.null(cl)) {
indexes_b = parLapply(cl, 1:B, boot.block, v = v, n = n, k = k)
}
if (is.null(cl)) {
indexes_b = lapply(1:B, boot.block, v = v, n = n, k = k)
}
if (!is.null(cl)) {
if (ram.allocation) {
weigth.in.cl = (as.numeric(object.size(d))/1e+06) *
3.2
used.memory = sum(sapply(ls(), function(x) object.size(get(x))))/1e+06
cl_number = min(as.numeric(ceiling((memory.limit() * 0.8 - used.memory)/weigth.in.cl) - 1), max.cores)
} else {
cl_number = max.cores
}
} else {
cl_number = 1
}
d = as.data.frame(d)
if (cl_number > 1) {
d_ij_avg_resampled = parLapply(cl[1:cl_number], indexes_b,
function(x, d, N) {
colSums((d[x, ,drop = FALSE]))/N
}, d = d, N = N)
} else {
d_ij_avg_resampled = lapply(indexes_b, function(x, d, N) {
colSums((d[x, ,drop = FALSE]))/N
}, d = d, N = N)
}
if (!is.null(cl)) {
d_b_i_mean = do.call(rbind, parLapply(cl, d_ij_avg_resampled, d_b_i_mean.fun, M = M, model.names = model.names))
}
if (is.null(cl)) {
d_b_i_mean = do.call(rbind, lapply(d_ij_avg_resampled, d_b_i_mean.fun, M = M, model.names = model.names))
}
d_b_i_var = colSums(t(apply(d_b_i_mean, 1, function(z) (z - d_i_mean)^2)))/B
names(d_b_i_var) = model.names
d_ij_avg_resampled = do.call(rbind, d_ij_avg_resampled)
d_var_ij = matrix(colSums(t(apply(d_ij_avg_resampled, 1, function(z) (z - d_ij_mean)^2)))/B, nrow = 1)
colnames(d_var_ij) = names(d_ij_mean)
TR = max((abs(d_ij_mean))/(d_var_ij^0.5))
TM = max(d_i_mean/(d_b_i_var^0.5))
Tb_R = sapply(1:B, function(i) {
max(abs(d_ij_avg_resampled[i, ] - d_ij_mean)/(d_var_ij^0.5))
})
Tb_M = sapply(1:B, function(i) {
max((d_b_i_mean[i, ] - d_i_mean)/(d_b_i_var^0.5))
})
Pr = length(which(TR < Tb_R))/B
Pm = length(which(TM < Tb_M))/B
# elimination rules
v_i_M = d_i_mean/(d_b_i_var^0.5)
v_i_M_order = order(v_i_M, decreasing = TRUE)
TR_all = ((d_ij_mean))/(d_var_ij^0.5)
v_i_R = sapply(1:M, function(x) {
max(TR_all[1, paste(model.names[x], model.names[-x], sep = ".")])
})
names(v_i_R) = model.names
v_i_R_order = order(v_i_R, decreasing = TRUE)
Pm_h0_mk = NULL
for (i in 1:M) {
if (i == 1) {
model_temp = model.names
} else {
model_temp = model.names[-v_i_M_order[1:(i - 1)]]
}
TM_temp = max(d_i_mean[model_temp]/(d_b_i_var[model_temp]^0.5))
Pm_h0_mk = c(Pm_h0_mk, length(which(TM_temp < Tb_M))/B)
}
mcs_Pm = sapply(1:M, function(x) max(Pm_h0_mk[1:x]))
# For TR statistic
combine.names = names(d_ij_mean)
Pr_h0_mk = NULL
for (i in 1:M) {
remove = NULL
if (i == 1) {
model_temp = combine.names
} else {
remove = do.call(c, lapply(1:(i - 1), function(x) {
which(gsub(model.names[v_i_R_order[x]], "", combine.names) != combine.names)
}))
model_temp = combine.names[-remove]
}
if (i < M) {
TR_temp = max((abs(d_ij_mean[model_temp]))/(d_var_ij[1, model_temp]^0.5))
Pr_h0_mk = c(Pr_h0_mk, length(which(TR_temp < Tb_R))/B)
} else {
Pr_h0_mk = c(Pr_h0_mk, 1)
}
}
mcs_Pr = sapply(1:M, function(x) max(Pr_h0_mk[1:x]))
matrix_show = matrix(NA, M, 7, dimnames = list(model.names,
c("Rank_M", "v_M", "MCS_M", "Rank_R", "v_R", "MCS_R",
"Loss")))
matrix_show[, "v_M"] = v_i_M
matrix_show[names(sort(v_i_M)), "Rank_M"] = 1:M
matrix_show[model.names[v_i_M_order], "MCS_M"] = mcs_Pm
matrix_show[, "v_R"] = v_i_R
matrix_show[names(sort(v_i_R)), "Rank_R"] = 1:M
matrix_show[model.names[v_i_R_order], "MCS_R"] = mcs_Pr
matrix_show[, "Loss"] = colMeans(Loss)
rm(list = "indexes_b")
if (!is.null(cl))
clusterEvalQ(cl, {
rm(list = ls())
gc()
})
gc()
if (statistic == "Tmax") {
p2test = Pm
}
if (statistic == "TR") {
p2test = Pr
}
if (p2test > alpha | all(d_var_ij == 0)) {
if (verbose) {
cat(paste("\n###########################################################################################################################\n"))
cat(paste("Superior Set Model created\t:\n"))
show(matrix_show)
cat(paste("p-value\t:\n"))
print(p2test)
cat(paste("\n###########################################################################################################################"))
}
break
} else {
if (statistic == "Tmax")
eliminate = which(v_i_M == max(v_i_M))
if (statistic == "TR")
eliminate = which(v_i_R == max(v_i_R))
if (verbose)
cat(paste("\nModel", model.names[eliminate],
"eliminated", Sys.time()))
Loss = as.matrix(Loss[, -eliminate])
colnames(Loss) = model.names[-eliminate]
if (ncol(Loss) == 1) {
if (verbose) {
cat(paste("\n###########################################################################################################################\n"))
cat(paste("Superior Set Model created\t:\n"))
matrix_show = matrix(matrix_show[-eliminate, ], nrow = 1, dimnames = list(colnames(Loss),
colnames(matrix_show)))
show(matrix_show)
cat(paste("p-value\t:\n"))
show(p2test)
cat(paste("\n###########################################################################################################################"))
}
break
}
}
}
elapsed.time = Sys.time() - time.start
n_elim = M_start - nrow(matrix_show)
out = new("SSM", show = matrix_show, Info = list(model.names = rownames(matrix_show),
elapsed.time = elapsed.time, statistic = statistic, n_elim = n_elim,
mcs_pvalue = p2test, alpha = alpha, B = B, k = k), Bootstrap = list(TR = list(Stat = TR,
BootDist = Tb_R), Tmax = list(Stat = TM, BootDist = Tb_M)))
return(out)
}
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.