# Install FARS from CRAN
#install.packages("FARS")

# Install FARS from GitHub
#install.packages("devtools")
#devtools::install_github("GPEBellocca/FARS")

library(FARS)

library(zoo)

library(readxl)
# Input dep variable
data_input_path <- system.file("extdata", "Data_IMF.xlsx", package = "FARS")
data_input <- read_excel(data_input_path, sheet = "data")
data_ts <- ts(data_input[, -1], frequency = 4)
data_diff <- diff(log(data_ts)) * 400
dep_variable <- data_diff[, 'United States', drop = FALSE] 
dep_variable <- dep_variable[2:60]

# Input data
data_df_path <- system.file("extdata", "DataBase.xlsx", package = "FARS")
data_df <- openxlsx::read.xlsx(data_df_path, sheet = "fulldata", cols = 2:625)
data_df <- data_df[,1:519]
data <- as.matrix(data_df)
dimnames(data) <- NULL

# Generate dates
quarters <- as.yearqtr(seq(from = as.yearqtr("2005 Q2"), by = 0.25, length.out = 59))
dates <- as.Date(quarters)

# MULTI-LEVEL DYNAMIC FACTOR MODEL
# 3 blocks
n_blocks <- 3 # 63 248 208
block_ind <- c(63,311,519)
r <- c(1,0,1,0,1,1,1)

mldfm_result <- mldfm(data, 
                      r=r, 
                      blocks = n_blocks, 
                      block_ind = block_ind , 
                      tol = 1e-6, 
                      max_iter = 1000,
                      method = 0) 

# Plot factors
# plot(mldfm_result, dates = dates) 

# Subsampling
n_samples <- 100
sample_size <- 0.9
mldfm_subsampling_result <- mldfm_subsampling(data, 
                                         r=r, 
                                         blocks = n_blocks, 
                                         block_ind = block_ind , 
                                         tol = 1e-6, 
                                         max_iter = 1000, 
                                         method = 0,
                                         n_samples = n_samples,
                                         sample_size = sample_size,
                                         seed = 123) 

# Create stressed scenario
scenario <- create_scenario(model = mldfm_result,
                                  subsample = mldfm_subsampling_result,
                                  data = data, 
                                  block_ind = block_ind, 
                                  alpha=0.95)

# Compute Quantiles
fars_result <- compute_fars(dep_variable, 
                              mldfm_result$Factors, 
                              scenario = scenario, 
                              h = 1,   
                              edge = 0.05, 
                              min = TRUE) 

# Plot quantiles
plot(fars_result,dates=dates)

# Density 
density <- nl_density(fars_result$Quantiles,  
                             levels = fars_result$Levels,  
                             est_points = 512, 
                             random_samples = 100000,
                             seed = 123)

# Plot density
# plot(density, time_index = dates)

#GaR
GaR0.05 <- quantile_risk(density, QTAU = 0.05)
GaR0.50 <- quantile_risk(density, QTAU = 0.50)
GaR0.95 <- quantile_risk(density, QTAU = 0.95)

# Scenario Density 
scenario_density <- nl_density(fars_result$Scenario_Quantiles,  
                             levels = fars_result$Levels,  
                             est_points = 512, 
                             random_samples = 100000,
                             seed = 123)

#GiS
GiS0.05 <- quantile_risk(scenario_density, QTAU = 0.05)
GiS0.25 <- quantile_risk(scenario_density, QTAU = 0.25)
GiS0.50 <- quantile_risk(scenario_density, QTAU = 0.50)
GiS0.75 <- quantile_risk(scenario_density, QTAU = 0.75)
GiS0.95 <- quantile_risk(scenario_density, QTAU = 0.95)

# Final GaR and GiS plot
library(ggplot2)

time <- dates[-1]

MLGaRGiS <- data.frame(
  time = time,
  dep_variable = as.vector(dep_variable[2:59]),
  GaR.05 = as.vector(GaR0.05[1:58]),
  GaR.95 = as.vector(GaR0.95[1:58]),
  GiS.05 = as.vector(GiS0.05[1:58]), 
  GiS.25 = as.vector(GiS0.25[1:58]),
  GiS.75 = as.vector(GiS0.75[1:58]),
  GiS.95 = as.vector(GiS0.95[1:58])
)


p <- ggplot(MLGaRGiS,aes(x=time,y=dep_variable)) + 
  geom_line() + theme_bw()+
  geom_ribbon(aes(ymin=GiS.05,ymax=GiS.95), alpha=0.1, fill="red") +
  geom_ribbon(aes(ymin=GiS.25,ymax=GiS.75), alpha=0.1, fill="grey10") +
  geom_line(aes(x=time, GaR.05), linewidth=0.1, colour="black", linetype="dashed") +
  geom_line(aes(x=time, GaR.95), linewidth=0.1, colour="black", linetype="dashed") +
  scale_y_continuous("Growth") +
  scale_x_date(date_labels = "%Y", date_breaks = "2 years") +
  theme(axis.text.x = element_text(angle = 90))
fig <- plotly::ggplotly(p)
fig