# 4.2. Computational Performance ------------------------------------------

library("here")
library("tidyverse")

sample_size <- seq(100, 1000, 100)

# read intermediate files
r_elapsed <- readr::read_csv("simulation/R/elapsed-r.txt")
eviews_elapsed <- readr::read_csv("simulation/eviews/elapsed-eviews.txt", col_names = FALSE) %>% 
  map_dbl(median)
matlab_elapsed <- readr::read_csv("simulation/matlab/elapsed-matlab.txt", col_names = FALSE) %>% 
  map_dbl(median)
stata_elapsed <- readr::read_csv("simulation/stata/elapsed-stata.txt", col_names = FALSE) %>% 
  map_dbl(median)

# Tabulate simulation results
tbl0 <- r_elapsed %>% 
  mutate(
    id = sample_size, 
    eviews = eviews_elapsed * 1000, 
    matlab = matlab_elapsed * 1000, 
    stata = stata_elapsed * 1000) %>% 
  select(id, everything())


# plot the time performance
pkgs <- c("MultipleBubbles", "psymonitor",  "exuber", "rtadf (Eviews)", "PSY (Matlab)", "radf (Stata)")
pkgs_reorder <- pkgs[c(4, 1, 2, 6, 5, 3)]

tbl0 %>% 
  set_names(c("ID", pkgs)) %>% 
  pivot_longer(cols = all_of(pkgs), values_to = "TimeElapsed", names_to = "Package") %>%
  group_by(ID) %>% 
  mutate(Package = factor(Package, levels = pkgs_reorder)) %>% 
  ggplot(aes(ID, TimeElapsed, col = Package, linetype = Package)) +
  geom_line() +
  geom_point() +
  labs(
    y = expression("Milliseconds (log scale)"), 
    x = "Sample size") +
  scale_x_continuous(breaks = seq(100, 1000, 100))  +
  scale_y_continuous(
    labels = scales::scientific_format(digits = 2), 
    trans = "log", 
    n.breaks = 7
  ) +
  scale_linetype_manual(values = c(6, 2, 1, 3, 4, 7)) +
  theme_bw() +
  theme(
    legend.title = element_blank(), 
    legend.key = element_blank(), 
    legend.position = c(0.12, 0.87), #c(0.8, 0.2), #
    legend.background = element_blank(), 
    panel.grid.minor = element_blank(), 
    panel.grid.major = element_line(linetype = "dashed")
  )

# Relative time comparison

options(pillar.sigfig = 5)
tbl0 %>% filter(id == 1000) %>%
  mutate_at(vars(-id), ~.x/1000)


# 5.1. Simulated Bubble Series --------------------------------------------

library("exuber")

sim <- sim_psy1(100, te = 50, tf = 70, seed = 145)

# Figure 2
figure2 <- autoplot(sim)
figure2

radf_sim <- radf(sim)

summary(radf_sim)

# Figure 3
library("ggplot2")
autoplot_sim_gsadf <- autoplot(radf_sim) + 
  labs(title = NULL)
autoplot_sim_sadf <- autoplot(radf_sim, option = "sadf") + 
  labs(title = NULL)
library("patchwork")
figure3 <- autoplot_sim_sadf + autoplot_sim_gsadf
figure3

datestamp(radf_sim)
datestamp(radf_sim, option = "sadf")

# 5.2. International House Prices -----------------------------------------

library("dplyr")
library("tidyr")

hprices <- ihpdr::ihpd_get(version = "1501") %>% 
  filter(country != "Aggregate") %>%
  select(Date, country, rhpi) %>% 
  pivot_wider(names_from = country, values_from = rhpi)
hprices

# Figure 4
figure4 <- ihpdr::ihpd_get(version = "1501") %>% 
  filter(Date <= "2015-12-31", country != "Aggregate") %>%
  select(Date, country, rhpi) %>% 
  group_by(Date) %>% 
  summarise(
    mean = mean(rhpi), 
    q25 = quantile(rhpi, probs = 0.25), 
    q75 = quantile(rhpi, probs = 0.75)
    ) %>% 
  ungroup() %>% 
  ggplot(aes(Date, mean)) +
  geom_line() +
  geom_ribbon(aes(ymin = q25, ymax = q75), fill = "#174B97", alpha = 0.5) +
  labs(x = NULL, y = NULL) +
  theme_bw()
figure4

min_dur <- psy_ds(hprices)
min_dur

radf_hprices <- radf(hprices, minw = 36, lag = 1)
mc_critical_values <- radf_mc_cv(nrow(hprices), minw = 36, seed = 145)

summary_results <- summary(radf_hprices, mc_critical_values)
summary_results[c("US", "UK")]

datestamp_results <- datestamp(radf_hprices, mc_critical_values, min_duration = min_dur)
datestamp_results[c("US", "UK")]

# Figure 5
figure5 <- autoplot(radf_hprices, mc_critical_values, select_series = c("US", "UK"), min_duration = min_dur)
figure5

tidy(radf_hprices)
tidy(mc_critical_values)

diagnostics(radf_hprices, mc_critical_values)
diagnostics(radf_hprices, mc_critical_values, option = "sadf")

# Figure 6
figure6 <- autoplot(datestamp_results)
figure6

# Application: Panel ------------------------------------------------------

panel_critical_values <- radf_sb_cv(hprices, minw = 36, lag = 1, seed = 145)
summary(radf_hprices, panel_critical_values)

# Figure 7
figure7 <- autoplot(radf_hprices, cv = panel_critical_values, min_duration = min_dur)
figure7

datestamp(radf_hprices, panel_critical_values, min_duration = min_dur)