#' The aira main class.
#'
#' @field bootstrap_iterations the number of bootstrap iterations to do for determining the significance of the effects
#' @field horizon the number of steps to look in the future
#' @field var_model the var model to perform the calculations on
#' @field orthogonalize use orthogonalized IRF
#' @importFrom methods setRefClass
#' @export Aira
#' @exportClass Aira
Aira <- setRefClass('Aira',
fields = c(
"bootstrap_iterations",
"horizon",
"var_model",
"orthogonalize",
"irf_cache",
"vars_functions",
"name"
),
methods = list(
initialize = function(bootstrap_iterations, horizon, var_model, orthogonalize, name="", reverse_order=FALSE) {
#rm(list=ls(pos='.GlobalEnv',all=TRUE),pos='.GlobalEnv')
irf_cache <<- list()
vars_functions <<- VarsFunctions$new(bootstrap_iterations = bootstrap_iterations,
horizon = horizon,
var_model = var_model,
orthogonalize = orthogonalize,
reverse_order = reverse_order
)
callSuper(bootstrap_iterations= bootstrap_iterations, horizon = horizon,
var_model = var_model, orthogonalize = orthogonalize, name = name)
},
determine_best_node_from_all = function(negative_variables = c()) {
"Returns the total effect a variable has on all other variables in the network.
If bootstrap iterations provided to aira is 0, we will not run any bootstrapping.
If bootstrap iterations >0, we will only consider the significant effects in the response. If negative_variables
are provided, we will convert those variables to positive ones (i.e., depression will become -1 * depression)
@param negative_variables the variables to invert to positibe variables"
total <- list()
for (variable in 1:var_model$K) {
variable_name <- .get_variable_name(variable)
total[[variable_name]] <- 0
for (response_name in 1:var_model$K) {
if (variable == response_name) next
response_name <- .get_variable_name(response_name)
# If the current response variable is a negative variable, we'd like to invert the result
multiplier <- ifelse(response_name %in% negative_variables, -1, 1)
irf <- multiplier * .calculate_irf(variable_name=variable_name, response=response_name)
total[[variable_name]] <- total[[variable_name]] + irf
}
}
total
},
determine_effect_network = function(include_autoregressive_effects = FALSE) {
"Returns the summed effect each node has on the other nodes node
@param include_autoregressive_effects if enabled, autoregressive effects are used (default FALSE). Not yet supported!"
K <- var_model$K
effect_matrix <- matrix(0, K, K)
colnames(effect_matrix) <- get_all_variable_names()
rownames(effect_matrix) <- get_all_variable_names()
for (variable in 1:var_model$K) {
variable_name <- .get_variable_name(variable)
for (response_name in 1:var_model$K) {
if (variable == response_name & !include_autoregressive_effects) next # Don't consider autocorrelation
response_name <- .get_variable_name(response_name)
result <- .calculate_irf(variable_name = variable_name, response = response_name)
effect_matrix[variable_name, response_name] <- result
}
}
effect_matrix
},
determine_percentage_effect = function(variable_to_improve, percentage) {
"Returns the percentage for each variable in the network (other then the provided variable)
to be changed in order to change the variable_to_improve with the given percentage.
@param variable_to_improve the name of the variable in the network which we'd like to improve
@param percentage the percentage with which we'd like to improve the variable to improve"
total <- list()
for (variable in 1:var_model$K) {
variable_name <- .get_variable_name(variable)
if (variable_name == variable_to_improve)
next
effect <- .calculate_irf(variable_name, variable_to_improve)
if (abs(effect) < 0.0001) {
total[[variable_name]] <- list(percentage = Inf, length_of_effect = Inf)
next
}
length_of_effect <- determine_length_of_effect(variable_name, variable_to_improve, 1, first_effect_only=FALSE, plot_results=FALSE)
length_of_effect <- ceiling(length_of_effect$effective_horizon)
needed_difference <- mean(var_model$y[,variable_to_improve]) * length_of_effect * (percentage / 100) * sd(var_model$y[,variable_name])
print(paste('Needed difference:', needed_difference, ', effect: ', effect, ' SD of var to use:', sd(var_model$y[,variable_name])))
denominator <- mean(var_model$y[,variable_name]) * effect * sd(var_model$y[,variable_to_improve])
print(paste('Numerator', needed_difference))
print(paste('Denominator', denominator))
needed_difference <- needed_difference / denominator
total[[variable_name]] <- list(percentage = needed_difference * 100, length_of_effect = length_of_effect)
}
total
},
determine_length_of_effect = function(variable_name, response, measurement_interval, first_effect_only=FALSE, plot_results=FALSE) {
"Returns the time in minues a variable is estimated to have an effect on another variable.
@param variable_name the name of the variable to receive the shock
@param response the name of the variable to respond to the shock
@param measurement interval the time in minutes between two measurements"
# TODO: The effect is currently not cached.
if(bootstrap_iterations <= 0) {
result <- vars_functions$irf(from=variable_name, to=response)
lower <- result$irf[[variable_name]]
upper <- lower
} else {
print('Determining with bootstrap')
result <- vars_functions$bootstrapped_irf(from=variable_name, to=response)
lower <- result$Lower[[variable_name]]
upper <- result$Upper[[variable_name]]
}
if(plot_results) plot(result)
# In order to be able to support normal irf (not only bootstrapped irf) we should incorporate rounding errors.
threshold <- 1e-4
low <- (lower > threshold)
high <- (upper < -threshold)
begin <- 1
end <- begin
effect_started <- FALSE
effective_horizon <- 0
exact_length <- 0
for (i in 1:horizon) {
if (low[i] | high[i]) {
# If the beginning of the effect is not on the first measurement, we should interpolate.
if (i > 1 & !effect_started) {
if (low[i]){
begin = i - ((lower[i] - threshold) / (lower[i] - lower[i - 1]))
} else if(high[i]) { # high[i] == true
begin = i - ((upper[i] + threshold) / (upper[i] - upper[i - 1]))
}
}
effect_started = TRUE
} else {
if (effect_started) {
if (low[i-1]){
end = i - 1 + (1 - ((lower[i] - threshold) / (lower[i] - lower[i - 1])))
} else if(high[i-1]) { # high[i-1] == true
end = i - 1 + (1 - ((upper[i] + threshold) / (upper[i] - upper[i - 1])))
}
effective_horizon <- end
effect_started <- FALSE
exact_length <- exact_length + (end - begin)
if(first_effect_only) break
}
}
}
# The effect has never stopped?
if(effect_started) {
exact_length <- exact_length + ((horizon - 1) - begin)
effective_horizon <- horizon
}
data.frame(
length_in_minutes = exact_length * measurement_interval,
length_of_effect = exact_length,
effective_horizon = effective_horizon
)
},
get_all_variable_names = function() {
"returns all variables in the var model"
dimnames(var_model$y)[[2]]
},
.get_variable_name = function(id) {
"Returns a variable name based on its id"
get_all_variable_names()[[id]]
},
.calculate_irf = function(variable_name, response = NULL, plot_results = FALSE){
"Calculates IRF and returns the total effect"
resulting_score <- 0
result <- ''
key <- paste(variable_name, response, sep="|")
# If we have processed this call before, return it from the cache
if((key %in% names(irf_cache)) & !plot_results) return(irf_cache[[key]])
if (bootstrap_iterations > 0) {
## TERRIBLE;
set_exo(var_model)
result <- vars_functions$bootstrapped_irf(from=variable_name, to=response)
low <- (result$Lower[[variable_name]] * (result$Lower[[variable_name]] > 0))
high <- (result$Upper[[variable_name]] * (result$Upper[[variable_name]] < 0))
sign_effects <- (low + high)[, !dimnames(result$Lower[[variable_name]])[[2]] %in% variable_name]
resulting_score <- sum(sign_effects)
} else {
result <- vars_functions$irf(from=variable_name, to=response)
resulting_score <- result$irf[[variable_name]][, !dimnames(result$irf[[variable_name]])[[2]] %in% variable_name, drop=FALSE]
resulting_score <- as.numeric(colSums(resulting_score))
}
if (plot_results) plot(result)
irf_cache[[key]] <<- resulting_score
resulting_score
}
)
)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.