An Introduction to the RecordTest Package"

In RecordTest: Inference Tools in Time Series Based on Record Statistics

library(RecordTest)
library(ggpubr)
data("TX_Zaragoza")
data("ZaragozaSeries")
data("Olympic_records_200m")
library(knitr)
knitr::opts_chunk$set(
  cache = TRUE,
  comment = "#>",
  collapse = TRUE,
  digits = 5,
  tidy = FALSE,
  background = "#FFFF00",
  fig.align = "center",
  warning = FALSE,
  message = FALSE
  )
RNGkind(sample.kind = "Rounding")
options(width = 80, digits = 5)
theme_set(theme_bw())

Introduction

The RecordTest package (Inference Tools in Time Series Based on Record Statistics) contains functions to visualize the behaviour of the record occurrence, functions to calculate a wide variety of distribution-free tests for trend in location, variation and non-stationarity, for change-point detection, and tools to prepare a time series in order to study its records.

Install RecordTest using

install.packages("RecordTest")

The introductory theory and summary for the package is at

help("RecordTest-package")

Here, the main purpose of the package is developed as well as an outline of the functions available in the package.

RecordTest has several functions that attempt to test the classical record model which assumes randomness in its variables, that is, they are continuous independent and identically distributed, by means of hypothesis tests and graphical tools.

Records in the 200-meter Olympic race

Sports Data

To begin with, RecordTest has a dataset Olympic_records_200m containing the record times time and record values value of the Olympic 200-meter, from 1900 to 2020. In this case, only the lower records are available.

Data preparation

The RecordTest functions need a complete series of observations to calculate its records. In order to apply these tools to the series of Olympic records, the function series_record is applied, which generates a series with the same records.

library(RecordTest)
library(ggpubr) # To join plots
data(Olympic_records_200m, package = "RecordTest")

or200m <- series_record(L_lower = Olympic_records_200m$time, 
                        R_lower = Olympic_records_200m$value,
                        Trows = 27)

Some tests and graphical tools

As a preview, the Olympic records series is drawn highlighting its lower records.

records(or200m, type = "points", alpha = c(1,0,1,0)) + ggplot2::ylab("seconds")

The graph below shows the number of accumulated lower records together with its confidence intervals under the null hypothesis, from which we see that the observed sample departs significantly since time $t = 13$, corresponding to the 1960 Olympics.

N.plot(or200m, record = c(0,1,0,0))

Due to the sample size is not very large, an exact one-sided test can be implemented based on the Poisson binomial distribution. The result is highly significant. The number of observed records is 12 while the expected under the null hypothesis is close to 4.

N.test(or200m, record = "lower", distribution = "poisson-binomial")

Records in temperatures and global warming

Temperature Data

RecordTest has a benchmark temperature dataset TX_Zaragoza containing the time series TX of daily maximum temperature at Zaragoza (Spain), from 01/01/1951 to 31/12/2020 measured in tenths of a degree Celsius. In this case, the whole series is available.

Data preparation

As a preview, the temperature series TX_Zaragoza$TX is drawn highlighting its upper and lower records.

data(TX_Zaragoza, package = "RecordTest")
records(TX_Zaragoza$TX, alpha = c(1, 1, 1, 0.1))

A large number of upper records are observed in the first observations and very few lower records. The initial behaviour in the records appear due to the seasonal increment of temperature between January and July, because this series has a strong seasonal component and serial correlation.

To take these characteristics into consideration, splitting the observed series into $M$ uncorrelated subseries is especially useful in the presence of serial correlation and seasonality. series_split splits TX_Zaragoza$TX into 365 subseries, each corresponding to each day of the year. series_uncor selects the larger number of columns or subseries that are not correlated with their adjacent columns.

TxZ365 <- series_split(TX_Zaragoza$TX, Mcols = 365)
TxZ <- series_uncor(TxZ365)
dim(TxZ)

Since the observed series is measured and rounded to tenths of a degree Celsius, ties can occur, which may cause a lower number of records to be identified than actually occurred. In particular, we observe that around $4\%$ of the records are weak records.

series_ties(TxZ365)

We could untie the possible records by adding a random value from a Uniform and independent distribution for each observation.

set.seed(23)
TxZ <- series_untie(TxZ)

Tests and graphics to detect trends

The following plot shows the upper and lower record times in the forward and backward series. Many more points are observed in the plots of the diagonal, giving evidence that there is a positive trend in the series.

L.plot(TxZ365)

The following plots show the mean number of (weighted) upper and lower records in the forward and backward series. Without weights the trend in the forward series is not significant, but backward it is highly significant. With weights, the trend in both directions is significant. This fact is produced because the evolution of temperature in Zaragoza is showing a faster increase since 1980.

ggpubr::ggarrange(N.plot(TxZ), N.plot(TxZ, weights = function(t) t-1),
        ncol = 2, nrow = 1, common.legend = TRUE, legend = "bottom")

A plot that gathers the information about the number of occurrences in the four types of record is obtained as follows.

foster.plot(TxZ) + ggplot2::ylim(-2.5, 2.5)

If we choose incremental weights $\omega_t = t-1$ to the records within the statistic, the trend becomes more significant earlier, then this graphical tool could be useful to identify the first time when non-stationary evidence appears. Several versions of plots of this style can be implemented (see help(foster.plot)).

foster.plot(TxZ, weights = function(t) t-1) + 
  ggplot2::ylim(-85, 85) +
  ggplot2::geom_vline(xintercept = 44, linetype = "dashed")

We can apply the associated test to detect non-stationary behaviour in the records of the series. The result is highly significant.

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

It is possible to use all the series if we compute the p-value with permutations, say 10,000:

set.seed(23)
foster.test(TxZ365, distribution = "normal", weights = function(t) t-1, 
  permutation.test = TRUE, B = 10000)

Under the null hypothesis of randomness, the record probability meets $t p_t = 1$. An exploratory tool can be proposed where $E(t \hat p_t) = \alpha t + \beta$ and also a regression test for the hypothesis [ H_0:\,\alpha=0,\,\beta=1 \qquad \text{and} \qquad H_1:\,\alpha\neq0\,\text{or}\,\beta\neq1. ] For Zaragoza data, plots related to this test detect a clear positive trend that is expressed with more upper records and less lower records in the forward series and the opposite in the backward series.

ggpubr::ggarrange(
  p.plot(TxZ, record = c(1,1,0,0)) + ggplot2::ylim(0, 5),
  p.plot(TxZ, record = c(0,0,1,1)) + ggplot2::ylim(0, 5),
  ncol = 2, nrow = 1, common.legend = TRUE, legend = "bottom")

Those tests can be implemented as follows, where in addition the estimation of the parameters of the line is displayed.

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

Other alternative based on a Monte Carlo approach is to join the information of all previous regression tests. Of the 1000 simulations considered under the null hypothesis, none has a statistic with a value greater than that of the observed series, making the test highly significant.

set.seed(23)
global.test(TxZ, FUN = p.regression.test, B = 1000)

Other powerful tests for trend detection can be implemented as follows:

brown.method(TxZ, weights = function(t) t-1)
N.test(TxZ, weights = function(t) t-1)

or

set.seed(23)
p.chisq.test(TxZ, simulate.p.value = TRUE)
lr.test(TxZ, simulate.p.value = TRUE, B = 10000)
score.test(TxZ)

Other plots:

ggpubr::ggarrange(
  p.plot(TxZ, plot = 1, record = c(1,1,0,0), 
         smooth.method = stats::loess, span = 0.25),
  p.plot(TxZ, plot = 1, record = c(1,1,0,0), 
         smooth.formula = y ~ I(x-1) - 1 + offset(rep(1, length(x)))),
  p.plot(TxZ, plot = 2, record = c(1,1,0,0)),
  p.plot(TxZ, plot = 3, record = c(1,1,0,0)),
  ncol = 2, nrow = 2, common.legend = TRUE, legend = "bottom")

Tests to detect change-points

Signs of climate change have not been present since the beginning of the series. We can implement a test to detect change-points in the series on a daily scale as follows.

change.point(ZaragozaSeries)
change.point(ZaragozaSeries, weights = function(t) sqrt(t), 
  record = "d", simulate.p.value = TRUE, B = 10000)

The change-point is found at time 36 (1986). We can see analogous results for the annual mean temperature series, a change-point is significantly detected with time estimate 38 (1988).

test.result <- change.point(rowMeans(TxZ365, na.rm = TRUE))
test.result

records(rowMeans(TxZ365, na.rm = TRUE)) + 
  ggplot2::geom_vline(xintercept = test.result$estimate, colour = "red")

There are still more tools! Try them yourself.

RecordTest documentation built on Aug. 8, 2023, 1:09 a.m.