library(tidyverse)
library(tidymodels)
library(modeltime)
library(timetk)
library(tidyquant)


#Coinbase Global
df_coin <- 
  tq_get("COIN") %>% 
  select(date, Coinbase = close)

#Bitcoin
df_btc <- 
  tq_get("BTC-USD") %>% 
  select(date, Bitcoin = close)

#Merging the datasets
df_merged <- 
  df_coin %>% 
  left_join(df_btc) %>% 
  drop_na() %>% 
  filter(date >= last(date) - months(36)) %>% 
  pivot_longer(-date,
               names_to = "id",
               values_to = "close") %>% 
  mutate(id = as_factor(id)) 


#Split Data
splits <- 
  time_series_split(
  df_merged,
  assess     = "15 days",
  cumulative = TRUE
)

#Create & Fit Forecasting Models

#Recipe
recipe_ml <- 
  recipe(close ~ ., training(splits)) %>%
  step_timeseries_signature(date) %>%
  step_rm(date) %>%
  step_dummy(all_nominal_predictors(), one_hot = TRUE) %>% 
  step_zv(all_predictors()) %>% 
  step_normalize(all_numeric_predictors())


#Model & Workflow
model_xgb <- 
  boost_tree("regression") %>%
  set_engine("xgboost")


wflw_fit_xgb <- 
  workflow() %>%
  add_model(model_xgb) %>%
  add_recipe(recipe_ml) %>%
  fit(training(splits))


#Adding fitted models to a Model Table
models_tbl <- modeltime_table(
  wflw_fit_xgb
)


#Calibrating the model to a testing set
calibration_tbl <- 
  models_tbl %>%
  modeltime_calibrate(
    new_data = testing(splits), 
    id       = "id"
  )


#Accuracy of the finalized model
calibration_tbl %>%
  modeltime_accuracy(metric_set = metric_set(rmse, rsq, mape),
                     acc_by_id = TRUE) %>% 
  table_modeltime_accuracy() 


#Conformal Split Method
#https://business-science.github.io/modeltime/articles/modeltime-conformal-prediction.html

forecast_tbl <- 
  calibration_tbl %>%
  modeltime_forecast(
    new_data      = testing(splits),
    actual_data   = df_merged %>% filter(date >= as.Date("2025-07-23")),
    conf_interval = 0.95,
    conf_method   = "conformal_split", # Split Conformal Method
    conf_by_id    = TRUE, # TRUE = local CI by ID, FALSE = global CI
    keep_data     = TRUE
  )

#Plotting prediction intervals
forecast_tbl %>%
  group_by(id) %>%
  plot_modeltime_forecast(
    .facet_ncol  = 1, 
    .line_size = 1.5,
    .interactive = FALSE
  ) +
  labs(title = "<span style = 'color:dimgrey;'>Conformal Prediction Intervals</span> of <span style = 'color:red;'>XGBoost</span> Model", 
       y = "", 
       x = "") + 
  scale_y_continuous(labels = scales::label_currency()) +
  scale_x_date(labels = scales::label_date("%b %d"),
               date_breaks = "4 days") +
  theme_tq(base_family = "Roboto Slab", base_size = 16) +
  theme(plot.title = ggtext::element_markdown(face = "bold", 
                                              hjust = 0.5,
                                              size = 18),
        strip.text = element_text(face = "bold"),
        legend.position = "none")