#################################################### ### 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)