cmfs_R_package

# CMFS Package: Capillarity Model of Financial Stability # Reference: Hammad (2025), https://zenodo.org/records/17943080 # Main Wrapper: Run Full Analysis with Prompts run_cmfs_analysis <- function() { cat("Welcome to CMFS Package: Capillarity Model Analysis\n") cat("Reference: Rayan S Hammad (2025), https://zenodo.org/records/17943080\n\n") inputs <- cmfs_input_dates() data <- cmfs_get_data(inputs) prep_data <- cmfs_prepare_data(data) sti <- cmfs_compute_sti(prep_data) external_var <- data$volume[-1] # Example; prompt for custom cat("Use default external variable (volume)? (y/n): ") if (tolower(readline()) == "n") { cat("Enter custom external variable vector: ") external_var <- eval(parse(text = readline())) # Advanced; assume vector } setar_model <- cmfs_fit_setar(sti, external_var) cmfs_output(setar_model, sti) cmfs_robustness(setar_model, ts(sti)) cat("\nAnalysis complete! Check outputs and files.\n") } # Block 1: Date/Frequency Input (with Validation) cmfs_input_dates <- function() { tryCatch({ cat("Enter start date (YYYY-MM-DD): ") start_date <- as.Date(readline()) if (is.na(start_date)) stop("Invalid start date.") cat("Enter end date (YYYY-MM-DD): ") end_date <- as.Date(readline()) if (is.na(end_date)) stop("Invalid end date.") if (start_date > end_date) stop("Start before end.") cat("Enter frequency (daily, weekly, monthly, quarterly, annual): ") freq <- tolower(readline()) if (!freq %in% c("daily", "weekly", "monthly", "quarterly", "annual")) stop("Invalid frequency.") return(list(start = start_date, end = end_date, freq = freq)) }, error = function(e) { cat("Error:", e$message, "\nRetry.\n") cmfs_input_dates() # Recursive retry }) } # Block 2: Data Source (with File Checks) cmfs_get_data <- function(inputs) { tryCatch({ cat("Enter source (yahoo, fred, csv, excel): ") source_type <- tolower(readline()) if (source_type %in% c("yahoo", "fred")) { cat("Enter symbol (e.g., SPY): ") symbol <- readline() data <- getSymbols(symbol, from = inputs$start, to = inputs$end, auto.assign = FALSE) prices <- Cl(data) volume <- Vo(data) } else if (source_type == "csv") { cat("Enter CSV path: ") file_path <- readline() if (!file.exists(file_path)) stop("File not found.") data <- read.csv(file_path) prices <- data$Close volume <- data$Volume } else if (source_type == "excel") { cat("Enter Excel path: ") file_path <- readline() if (!file.exists(file_path)) stop("File not found.") cat("Enter sheet: ") sheet <- readline() data <- read_excel(file_path, sheet = sheet) prices <- data$Close volume <- data$Volume } else stop("Invalid source.") if (inputs$freq != "daily") { data_xts <- xts(cbind(prices, volume), order.by = as.Date(index(data))) data_agg <- to.period(data_xts, period = inputs$freq) prices <- Cl(data_agg) volume <- Vo(data_agg) } return(list(prices = prices, volume = volume)) }, error = function(e) { cat("Error:", e$message, "\nRetry.\n") cmfs_get_data(inputs) }) } # Block 3: Prepare Data cmfs_prepare_data <- function(data) { tryCatch({ if (length(data$prices) < 2) stop("Insufficient data.") traded_value <- data$prices * data$volume return(list(prices = data$prices, traded_value = traded_value)) }, error = function(e) { cat("Error in data prep:", e$message, "\n") return(NULL) }) } # Block 4: Compute STI (with GARCH Variants) cmfs_compute_sti <- function(prep_data) { tryCatch({ returns <- diff(log(prep_data$prices))[-1] amihud <- abs(returns) / prep_data$traded_value[-1] cat("Choose GARCH variant (sGARCH, eGARCH, gjrGARCH; default sGARCH): ") variant <- tolower(readline()) if (!variant %in% c("sgarch", "egarch", "gjrgarch")) variant <- "sgarch" garch_spec <- ugarchspec(variance.model = list(model = toupper(variant), garchOrder = c(1,1)), mean.model = list(armaOrder = c(0,0))) garch_fit <- ugarchfit(garch_spec, data = returns) if (garch_fit@fit$convergence != 0) warning("GARCH issue; approximate results.") sigma_sq <- fitted(garch_fit)^2 sti <- prep_data$traded_value[-1] / sigma_sq return(sti) }, error = function(e) { cat("Error in STI:", e$message, "\nFallback to variance.\n") sigma_sq <- var(returns) sti <- prep_data$traded_value[-1] / sigma_sq return(sti) }) } # Block 5: Fit SETAR cmfs_fit_setar <- function(sti, external_var) { tryCatch({ cat("Enter SETAR lags (m, default 1): ") m <- as.integer(readline()) if (is.na(m) || m < 1) m <- 1 sti_ts <- ts(sti) setar_model <- setar(sti_ts, m = m, thDelay = 1, thVar = external_var) print(summary(setar_model)) return(setar_model) }, error = function(e) { cat("Error in SETAR:", e$message, "\nAdjust params.\n") return(NULL) }) } # Block 6: Output Results cmfs_output <- function(setar_model, sti) { if (is.null(setar_model)) stop("No model.") cat("SETAR Summary:\n") print(summary(setar_model)) threshold <- quantile(sti, 0.05) ggplot(data.frame(Date = seq_along(sti), STI = sti), aes(x = Date, y = STI)) + geom_bar(stat = "identity") + geom_hline(yintercept = threshold, color = "red") + labs(title = "Systemic Tenacity Index (STI)", subtitle = "Red = 5th Percentile Threshold") ggsave("sti_plot.png") cat("Plot saved.\n") cat("Save to CSV? (y/n): ") save <- tolower(readline()) if (save == "y") { write.csv(data.frame(STI = sti), "sti_results.csv") cat("Saved.\n") } } # Block 7: Robustness Checks cmfs_robustness <- function(setar_model, ts_data) { if (is.null(setar_model)) stop("No model.") cat("Running checks...\n") # Linearity lin_test <- linearityTest(ts_data, m = setar_model$m) print(lin_test) resids <- residuals(setar_model) # Normality jb <- jarque.bera.test(resids) print(jb) # Autocorrelation lb <- Box.test(resids, lag = 20, type = "Ljung-Box") print(lb) # Heteroskedasticity bp <- bptest(resids ~ fitted(setar_model)) print(bp) # Threshold sensitivity th_est <- setar_model$th sens_low <- setar(ts_data, m = setar_model$m, th = th_est * 0.9) sens_high <- setar(ts_data, m = setar_model$m, th = th_est * 1.1) print(c(Original_AIC = AIC(setar_model), Low = AIC(sens_low), High = AIC(sens_high))) # Out-of-sample vs. AR/GARCH train_end <- floor(0.8 * length(ts_data)) train <- ts_data[1:train_end] test <- ts_data[(train_end+1):length(ts_data)] train_setar <- setar(train, m = setar_model$m, thDelay = setar_model$thDelay) fc_setar <- predict(train_setar, n.ahead = length(test)) mse_setar <- mean((fc_setar - test)^2, na.rm = TRUE) ar_bench <- auto.arima(train) fc_ar <- forecast(ar_bench, h = length(test))$mean mse_ar <- mean((fc_ar - test)^2, na.rm = TRUE) garch_spec <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1,1))) garch_fit <- ugarchfit(garch_spec, data = train) fc_garch <- ugarchforecast(garch_fit, n.ahead = length(test))@forecast$sigmaFor^2 mse_garch <- mean((fc_garch - test)^2, na.rm = TRUE) print(c(SETAR_MSE = mse_setar, AR_MSE = mse_ar, GARCH_MSE = mse_garch)) }