#################################################################################################

##---------------------------------------------------------------------------------------------
# Functions for Classes by Prof. Jin Man Lee
# DePaul University
#
# To use any functions here, you need to run the following code to read into your program
# source("http://bigblue.depaul.edu/jlee141/econdata/R/func_tslib.R")
#

#################################################################################################
##---------------------------------------------------------------------------------------------##

adf_fulltest <- function(indata,maxlag){
if (!require("fUnitRoots")) install.packages("fUnitRoots")
maxlag <- maxlag
A <- matrix(nrow=maxlag+1,ncol=7)
for (k in 0:maxlag) {
  k1 <- k + 1
  adf0 <- adfTest(indata, lags = k, type = "nc")    
       A[k1,2] <- round(adf0@test$statistic[1],digit=2) ;  A[k1,3] <- adf0@test$p.value ;  A[k1,1] <- k
  adf0 <- adfTest(indata, lags = k, type = "c")     
       A[k1,4] <- round(adf0@test$statistic[1],digit=2) ;  A[k1,5] <- adf0@test$p.value    
  adf0 <- adfTest(indata, lags = k, type = "ct")   
       A[k1,6] <- round(adf0@test$statistic[1],digit=2) ;  A[k1,7] <- round(adf0@test$p.value,digit=2)     
   }  
colnames(A) <- c("LAG","ADF","P(ADF)","ADF_C","P(ADF_C)","ADF_CT","P(ADF_CT)")
return(A)
}
##---------------------------------------------------------------------------------------------##

#################################################################################################
# Local Projection Model Estimation
##---------------------------------------------------------------------------------------------##

localproj <- function(endog_data,exog_data,lagend,lagexo,maxlag){
results_lin <- lp_lin(endog_data      , 
                      lags_endog_lin = lagend,     # Number of lags for endogenous data
                      exog_data=exog_data ,
                      lags_exog      = lagexo  ,
                      trend          = 0,     # 0 = no trend, 1 = trend, 2 = trend & trend^2    
                      shock_type     = 1,     # 0 = standard deviation shock, 1 = unit shock
                      confint        = 1.96,  # Width of confidence bands: # 1 = 68%, 1.67 = 90%, 1.96 = 95%
                      hor            = maxlag)    # Number of cores to use. When NULL, the number of cores # is chosen automatically    
g1 <- summary(results_lin)

linear_plots <- plot_lin(results_lin)
linear_plots[[1]]
g2 <- linear_plots[[1]]

lin_plots_all <- sapply(linear_plots, ggplotGrob)
g3 <- marrangeGrob(lin_plots_all, nrow = ncol(endog_data), ncol = ncol(endog_data), top = NULL)
gout <- list(g1,g2,g3)
return(gout)
}

################################################################################################
# Automatic Unit Root Tests using tseries
##--------------------------------------------------------------------------------------------##

unitroot_tests <- function(series) {
       if (!require("tseries")) install.packages("tseries")
        adf_test <- adf.test(series, alternative = "stationary")
        kpss_test <- kpss.test(series, null = "Level", lshort = TRUE)
        pp_test <- pp.test(series)
	
        data.frame(
        Test = c("ADF Test", "PP Test", "KPSS Test"),
        Lag_order = c(adf_test$parameter,
                      pp_test$parameter,
                      kpss_test$parameter),
        Statistic = c(adf_test$statistic,
                      pp_test$statistic,
                      kpss_test$statistic),
        Stationary_P_Value = c(adf_test$p.value, 
                               pp_test$p.value, 
                               kpss_test$p.value))
}

# =============================================================================
# A single function to run ADF, PP, KPSS, DF-GLS, and Zivot-Andrews
# with a user-specified lag k (or automatic selection).
# =============================================================================

# --- Dependencies ------------------------------------------------------------
# install.packages(c("tseries", "urca", "forecast"))
library(tseries)
library(urca)
library(forecast)


# =============================================================================
#  FUNCTION: unit_root_tests()
#
#  Arguments
#  ---------
#  x        : numeric vector or ts object — the series to test
#  k        : integer lag length used by ADF and DF-GLS.
#             If NULL (default), lag is chosen automatically:
#               ADF/ADF-GLS → BIC via ur.df / ur.ers
#               (PP and KPSS use kernel bandwidths, not discrete lags)
#  series_name : character label printed in output headers
#  model    : deterministic component assumed by ADF, PP, DF-GLS, ZA
#               "none"  → no constant, no trend
#               "drift" → constant only          (default)
#               "trend" → constant + linear trend
#  za_model : Zivot-Andrews model variant
#               "intercept" → break in intercept (Model A)
#               "trend"     → break in trend     (Model B)
#               "both"      → break in both      (Model C, default)
#  print_details : logical — if TRUE, print full urca summaries;
#                  if FALSE (default), print only the tidy summary table
#
#  Returns
#  -------
#  A named list with elements: adf, pp, kpss, dfgls, za, summary_df
# =============================================================================

unit_root_tests <- function(
    x,
    k            = NULL,
    series_name  = "Series",
    model        = c("drift", "trend", "none"),
    za_model     = c("both", "intercept", "trend"),
    print_details = FALSE
) {
  
  model    <- match.arg(model)
  za_model <- match.arg(za_model)
  
  # Map 'model' to urca and tseries conventions --------------------------------
  urca_type  <- switch(model, none = "nc", drift = "drift", trend = "trend")
  kpss_type  <- switch(model, none = "mu", drift = "mu",    trend = "tau")
  pp_model   <- switch(model, none = "constant", drift = "constant", trend = "trend")
  ers_model  <- switch(model, none = "constant", drift = "constant", trend = "trend")
  
  # Header ---------------------------------------------------------------------
  cat("\n", strrep("=", 68), "\n", sep = "")
  cat("  Unit Root Test Battery  |  Series:", series_name, "\n")
  cat("  Deterministic model    :", model, "\n")
  if (!is.null(k)) cat("  Lag (k)                :", k, "\n")
  else             cat("  Lag (k)                : auto (BIC for ADF/DF-GLS)\n")
  cat(strrep("=", 68), "\n\n", sep = "")
  
  results <- list()
  
  # ============================================================
  # 1. ADF — Augmented Dickey-Fuller
  #    H0: unit root | Ha: stationary
  # ============================================================
  cat("--- [1] Augmented Dickey-Fuller (ADF) ---\n")
  
  if (is.null(k)) {
    adf_fit <- ur.df(x, type = urca_type, selectlags = "BIC")
  } else {
    adf_fit <- ur.df(x, type = urca_type, lags = k)
  }
  
  adf_slot  <- switch(urca_type, nc = "tau1", drift = "tau2", trend = "tau3")
  adf_stat  <- adf_fit@teststat[1, adf_slot]
  adf_cv    <- adf_fit@cval[adf_slot, ]        # 1%, 5%, 10%
  adf_lag   <- length(adf_fit@testreg$coefficients[, 1]) -
    switch(urca_type, nc = 1L, drift = 2L, trend = 3L)
  adf_pval  <- tryCatch(
    suppressWarnings(adf.test(x, k = max(adf_lag, 1))$p.value),
    error = function(e) NA_real_
  )
  
  if (print_details) print(summary(adf_fit))
  
  cat(sprintf("  Statistic (tau): %7.4f    Lag used: %d\n", adf_stat, adf_lag))
  cat(sprintf("  Critical values:  1%% = %.4f  |  5%% = %.4f  |  10%% = %.4f\n",
              adf_cv["1pct"], adf_cv["5pct"], adf_cv["10pct"]))
  cat(sprintf("  Approx. p-value: %.4f\n", adf_pval))
  cat(sprintf("  Decision (5%%):   %s\n\n",
              ifelse(adf_stat < adf_cv["5pct"], "Reject H0 → Stationary",
                     "Fail to Reject H0 → Unit Root")))
  
  results$adf <- list(statistic = adf_stat, cv = adf_cv,
                      lag = adf_lag, pvalue = adf_pval, fit = adf_fit)
  
  # ============================================================
  # 2. PP — Phillips-Perron
  #    H0: unit root | Ha: stationary
  # ============================================================
  cat("--- [2] Phillips-Perron (PP) ---\n")
  
  pp_fit   <- ur.pp(x, type = "Z-tau", model = pp_model, lags = "short")
  pp_stat  <- pp_fit@teststat[1]
  pp_cv    <- pp_fit@cval[1, ]
  
  pp_pval  <- tryCatch(
    suppressWarnings(pp.test(x)$p.value),
    error = function(e) NA_real_
  )
  
  if (print_details) print(summary(pp_fit))
  
  cat(sprintf("  Statistic (Z-tau): %7.4f\n", pp_stat))
  cat(sprintf("  Critical values:   1%% = %.4f  |  5%% = %.4f  |  10%% = %.4f\n",
              pp_cv["1pct"], pp_cv["5pct"], pp_cv["10pct"]))
  cat(sprintf("  Approx. p-value:  %.4f\n", pp_pval))
  cat(sprintf("  Decision (5%%):    %s\n\n",
              ifelse(pp_stat < pp_cv["5pct"], "Reject H0 → Stationary",
                     "Fail to Reject H0 → Unit Root")))
  
  results$pp <- list(statistic = pp_stat, cv = pp_cv,
                     pvalue = pp_pval, fit = pp_fit)
  
  # ============================================================
  # 3. KPSS — Kwiatkowski-Phillips-Schmidt-Shin
  #    H0: STATIONARY | Ha: unit root  ← reversed null!
  # ============================================================
  cat("--- [3] KPSS  (H0 = Stationary) ---\n")
  
  kpss_fit  <- ur.kpss(x, type = kpss_type, lags = "short")
  kpss_stat <- kpss_fit@teststat[1]
  kpss_cv   <- kpss_fit@cval[1, ]
  
  kpss_pval <- tryCatch(
    suppressWarnings(
      kpss.test(x, null = ifelse(kpss_type == "tau",
                                 "Trend", "Level"))$p.value
    ),
    error = function(e) NA_real_
  )
  
  if (print_details) print(summary(kpss_fit))
  
  cat(sprintf("  Statistic: %7.4f\n", kpss_stat))
  cat(sprintf("  Critical values:  1%% = %.4f  |  5%% = %.4f  |  10%% = %.4f\n",
              kpss_cv["1pct"], kpss_cv["5pct"], kpss_cv["10pct"]))
  cat(sprintf("  Approx. p-value: %.4f\n", kpss_pval))
  cat(sprintf("  Decision (5%%):   %s\n\n",
              ifelse(kpss_stat > kpss_cv["5pct"],
                     "Reject H0 → Unit Root",
                     "Fail to Reject H0 → Stationary")))
  
  results$kpss <- list(statistic = kpss_stat, cv = kpss_cv,
                       pvalue = kpss_pval, fit = kpss_fit)
  
  # ============================================================
  # 4. DF-GLS — Elliott, Rothenberg & Stock (1996)
  #    GLS-detrended ADF; near-optimal power
  #    H0: unit root | Ha: stationary
  # ============================================================
  cat("--- [4] DF-GLS  (Elliott, Rothenberg & Stock 1996) ---\n")
  
  lag_max_gls <- if (is.null(k)) floor(12 * (length(x) / 100)^(1/4)) else k
  
  dfgls_fit  <- ur.ers(x, type = "DF-GLS", model = ers_model,
                       lag.max = lag_max_gls)
  dfgls_stat <- dfgls_fit@teststat[1]
  dfgls_cv   <- dfgls_fit@cval[1, ]
  dfgls_lag  <- which.min(abs(dfgls_fit@testreg$coefficients[, 1])) - 1L
  
  if (print_details) print(summary(dfgls_fit))
  
  cat(sprintf("  Statistic: %7.4f    Lag max: %d\n", dfgls_stat, lag_max_gls))
  cat(sprintf("  Critical values:  1%% = %.4f  |  5%% = %.4f  |  10%% = %.4f\n",
              dfgls_cv["1pct"], dfgls_cv["5pct"], dfgls_cv["10pct"]))
  cat(sprintf("  Decision (5%%):   %s\n\n",
              ifelse(dfgls_stat < dfgls_cv["5pct"],
                     "Reject H0 → Stationary",
                     "Fail to Reject H0 → Unit Root")))
  
  results$dfgls <- list(statistic = dfgls_stat, cv = dfgls_cv,
                        lag_max = lag_max_gls, fit = dfgls_fit)
  
  # ============================================================
  # 5. Zivot-Andrews — Unit root with endogenous structural break
  #    H0: unit root (no break) | Ha: trend-stationary with break
  # ============================================================
  cat(sprintf("--- [5] Zivot-Andrews  (structural break | model = '%s') ---\n",
              za_model))
  
  za_lag <- if (is.null(k)) 4L else k
  
  za_fit   <- ur.za(x, model = za_model, lag = za_lag)
  za_stat  <- za_fit@teststat[1]
  za_cv    <- za_fit@cval
  za_break <- za_fit@bpoint
  
  # Recover break date if x is a ts object
  if (inherits(x, "ts")) {
    break_date <- time(x)[za_break]
    break_label <- sprintf("%.2f  (obs %d)", break_date, za_break)
  } else {
    break_label <- sprintf("obs %d", za_break)
  }
  
  if (print_details) print(summary(za_fit))
  
  cat(sprintf("  Statistic: %7.4f    Lag: %d\n", za_stat, za_lag))
  cat(sprintf("  Critical values:  1%% = %.4f  |  5%% = %.4f  |  10%% = %.4f\n",
              za_cv[1], za_cv[2], za_cv[3]))
  cat(sprintf("  Estimated break point: %s\n", break_label))
  cat(sprintf("  Decision (5%%):   %s\n\n",
              ifelse(za_stat < za_cv[2],
                     "Reject H0 → Stationary with structural break",
                     "Fail to Reject H0 → Unit Root")))
  
  results$za <- list(statistic = za_stat, cv = za_cv,
                     break_point = za_break, fit = za_fit)
  
  # ============================================================
  # 6. Tidy Summary Table
  # ============================================================
  cat(strrep("-", 68), "\n")
  cat("  SUMMARY TABLE\n")
  cat(strrep("-", 68), "\n")
  cat(sprintf("  %-10s  %-6s  %-9s  %-8s  %-30s\n",
              "Test", "H0", "Statistic", "p-value", "Decision (5%)"))
  cat(strrep("-", 68), "\n")
  
  fmt_row <- function(name, h0, stat, pval, decision) {
    pval_str <- if (is.na(pval)) "  N/A  " else sprintf("%.4f", pval)
    cat(sprintf("  %-10s  %-6s  %9.4f  %-8s  %s\n",
                name, h0, stat, pval_str, decision))
  }
  
  fmt_row("ADF",    "I(1)",
          adf_stat,  adf_pval,
          ifelse(adf_stat  < adf_cv["5pct"],
                 "Reject → Stationary", "Fail to Reject → I(1)"))
  
  fmt_row("PP",     "I(1)",
          pp_stat,   pp_pval,
          ifelse(pp_stat   < pp_cv["5pct"],
                 "Reject → Stationary", "Fail to Reject → I(1)"))
  
  fmt_row("KPSS",   "I(0)",
          kpss_stat, kpss_pval,
          ifelse(kpss_stat > kpss_cv["5pct"],
                 "Reject → I(1)",       "Fail to Reject → Stationary"))
  
  fmt_row("DF-GLS", "I(1)",
          dfgls_stat, NA_real_,
          ifelse(dfgls_stat < dfgls_cv["5pct"],
                 "Reject → Stationary", "Fail to Reject → I(1)"))
  
  fmt_row("ZA",     "I(1)",
          za_stat,   NA_real_,
          ifelse(za_stat   < za_cv[2],
                 "Reject → TS w/ break", "Fail to Reject → I(1)"))
  
  cat(strrep("-", 68), "\n\n")
  
  # Build tidy data frame
  summary_df <- data.frame(
    Test      = c("ADF", "PP", "KPSS", "DF-GLS", "ZA"),
    H0        = c("I(1)", "I(1)", "I(0)", "I(1)", "I(1)"),
    Statistic = round(c(adf_stat, pp_stat, kpss_stat, dfgls_stat, za_stat), 4),
    p_value   = round(c(adf_pval, pp_pval, kpss_pval, NA, NA), 4),
    CV_5pct   = round(c(adf_cv["5pct"], pp_cv["5pct"], kpss_cv["5pct"],
                        dfgls_cv["5pct"], za_cv[2]), 4),
    stringsAsFactors = FALSE
  )
  
  results$summary_df <- summary_df
  invisible(results)
}

#-----------------------------------------------------------------------------#
# Function for the 
# calculate_price_index(coeff_df, "200801","coeff","yearmm")
# Make sure to have the rest options should have the quotations 
calculate_price_index <- function(data, base_ym, var_input, yearmm_column) {
  # Ensure column names are provided correctly
  if (!yearmm_column %in% names(data) || !var_input %in% names(data)) {
    stop("Column names provided do not exist in the dataframe.")
  }
  
  # Filter to get the base year-month row
  base_val <- data %>%
    filter(get(yearmm_column) == base_ym) %>%
    pull(var_input) %>%
    unlist()
  
  # Check if base_val is empty or NA
  if (length(base_val) == 0 || is.na(base_val)) {
    stop("Base value for specified year-month not found or NA.")
  }
  
  # Calculate the price index using the base value
  index_column <- 100 * data[[var_input]] / base_val[1]
  
  # Return the updated dataframe with the price index
  return(index_column)
}

###########################################################################

##################################################################################################
#  Preparation of Random Forest Data 
prepare_rf_data <- function(data, target_var, lags = lags, date_index = NULL) {
  # Ensure input is a data.frame
  df <- as.data.frame(data)

  # Create lagged features
  for (var in colnames(df)) {
    for (k in 1:lags) {
      df[[paste0(var, "_L", k)]] <- dplyr::lag(df[[var]], k)
    }
  }

  df <- na.omit(df)

  # Create target and remove original target_var from predictors
  df$Target <- df[[target_var]]
  df[[target_var]] <- NULL

  # Attach date index if provided
  if (!is.null(date_index)) {
    df$Date <- tail(date_index, nrow(df))
  }

  return(df)
}
##--------------------------------------------------------------------------------------------##