####################################################
### R code to reproduce examples of
###
### "RecordTest: An R Package to Analyze
###  Non-Stationarity in the Extremes Based on
###  Record-Breaking Events"
###
####################################################

### SECTION 4. An example ###

## Dataset
library("RecordTest")
Tx <- TX_Zaragoza$TX / 10

## Data preparation and exploratory analysis
records(Tx) + 
  ggplot2::geom_smooth(formula = y ~ x, method = stats::loess, 
    mapping = ggplot2::aes(y = Tx), se = FALSE, col = "red") 

TxZ365 <- series_split(Tx, Mcols = 365)
TxZ    <- series_uncor(TxZ365)
dim(TxZ)

series_uncor(TxZ365, return.value = "indexes")
range(diff(series_uncor(TxZ365, return.value = "indexes")))

lapply(series_ties(TxZ, record = "upper"), round, digits = 2)
lapply(series_ties(TxZ, record = "lower"), round, digits = 2)

records(TxZ[, 30], direction = "both") +
  ggplot2::scale_x_continuous(name = "Year", 
    breaks = c(10, 30, 50, 70), labels = c("1960", "1980", "2000", "2020"))

L.plot(TxZ) 


## Inference tools to study non-stationarity in the extremes

## Tests to detect deviations from stationarity
N.test(TxZ, weights = function(t) t - 1)
set.seed(1234)
N.test(TxZ, weights = function(t) t - 1, simulate.p.value = TRUE)

N.test(TxZ, weights = function(t) t - 1, record = "lower",  
  alternative = "less")

foster.test(TxZ, weights = function(t) t - 1, statistic = "U")

foster.test(TxZ, weights = function(t) t - 1, statistic = "D")

foster.test(TxZ, weights = function(t) t - 1, statistic = "D",
  distribution = "t")

brown.method(TxZ, weights = function(t) t - 1,
  alternative = c("FU" = "greater", "FL" = "less", 
                  "BU" = "less", "BL" = "greater"))


## Graphical tools to detect deviations from stationarity
ggpubr::ggarrange(
  N.plot(TxZ, 
    point.col = c("FU" = "red", "FL" = "red", "BU" = "blue", "BL" = "blue"), 
    point.shape = c("FU" = 19, "FL" = 4, "BU" = 19, "BL" = 4)) + 
    ggplot2::scale_x_continuous(name = "Year (forward)",
      breaks = c(10, 30, 50, 70), 
      labels=c("1960", "1980", "2000", "2020"), 
      sec.axis = ggplot2::sec_axis(~ 2021 - ., name = "Year (backward)",
                                   breaks = 1950 + c(10, 30, 50, 70))) +
    ggplot2::theme(axis.title.x = ggplot2::element_text(colour = "red"), 
      axis.text.x = ggplot2::element_text(colour = "red"),
      axis.title.x.top = ggplot2::element_text(colour = "blue"), 
      axis.text.x.top = ggplot2::element_text(colour = "blue")),
  N.plot(TxZ, weights = function(t) t - 1,
    point.col = c("FU" = "red", "FL" = "red", "BU" = "blue", "BL" = "blue"), 
    point.shape = c("FU" = 19, "FL" = 4, "BU" = 19, "BL" = 4)) + 
    ggplot2::scale_x_continuous(name = "Year (forward)",
      breaks = c(10, 30, 50, 70), 
      labels=c("1960", "1980", "2000", "2020"), 
      sec.axis = ggplot2::sec_axis(~ 2021 - ., name = "Year (backward)",
                                   breaks = 1950 + c(10, 30, 50, 70))) +
    ggplot2::theme(axis.title.x = ggplot2::element_text(colour = "red"), 
      axis.text.x = ggplot2::element_text(colour = "red"),
      axis.title.x.top = ggplot2::element_text(colour = "blue"), 
      axis.text.x.top = ggplot2::element_text(colour = "blue")),
  nrow = 1, common.legend = TRUE, legend = "bottom")

foster.plot(TxZ, weights = function(t) t - 1) + 
  ggplot2::ylim(-75,75) +
  ggplot2::scale_x_continuous(name = "Year", 
    breaks = c(10, 30, 50, 70), labels = c("1960", "1980", "2000", "2020"))

p.plot(TxZ, record = c("FU" = 1, "FL" = 0, "BU" = 0, "BL" = 0),
  smooth.formula = y ~ poly(x, degree = 2))

p.plot(TxZ, record = c("FU" = 0, "FL" = 1, "BU" = 0, "BL" = 1), 
  point.col = c("FU" = NA, "FL" = "blue", "BU" = NA,  "BL" = "red"))

p.regression.test(TxZ, formula = y ~ poly(x, degree = 2))

p.regression.test(series_rev(TxZ), record = "lower")


## Tests for change-point detection
TxZmean <- rowMeans(TxZ365, na.rm = TRUE)

records(TxZmean) + 
  ggplot2::scale_x_continuous(name = "Year", breaks = c(10, 30, 50, 70), 
    labels = c("1960", "1980", "2000", "2020")) +   
  ggplot2::geom_vline(xintercept = change.point(TxZmean)$estimate, 
    color = "red")

change.point(TxZmean)

change.point(TxZ)