Introduction to SNSeg and Examples
SNSeg supports change-points estimation for both univariate and
multivariate time series (including high-dimensional time series with
dimension greater than 10.) using Self-Normalization (SN) based
framework. Please read Zhao, Jiang and Shao (2022)
<doi.org/10.1111/rssb.12552> for details of the SN-based
algorithms.
The package contain three functions for change-points estimation:
- The function
SNSeg_Uni works for a univariate time
series.
- The function
SNSeg_Multi works for multivariate time
series with dimension up to ten.
- The function
SNSeg_HD works for high-dimensional time
series with dimension above ten.
All functions contain two important input arguments:
grid_size_scale and grid_size.
grid_size_scale: the trimming parameter \(\epsilon\), a parameter to control the
local window size grid_size. The range of
grid_size_scale should be between 0.05 and 0.5. Any input
grid_size_scale smaller than 0.05 will be automatically
changed to 0.05, and any input grid_size_scale greater than
0.5 will be automatically changed to 0.5.grid_size: the local window size \(h\) for local scanning for each time point.
According to Zhao et al. (2022), grid_size is a product of
grid_size_scale and the length of time.
Users can set their own grid_size or leave it as
NULL in the input arguments. If grid_size is
NULL, these functions will compute the SN test statistic
using the value of grid_size_scale. If
grid_size is set by users, the functions will first
calculate grid_size_scale by diving the length of time, and
then compute the SN test statistic.
For the other input arguments: * ts: Users should enter
a time series for this argument. For SNSeg_Uni, the
dimension of ts should be exactly one in most cases (or two
when paras_to_test = 'bivcor'); for
SNSeg_Multi, the dimension must be at least two but no more
than ten; for SNSeg_HD, the dimension must be greater than
ten. * paras_to_test: This argument in functions
SNSeg_Uni and SNSeg_Multi allows users to
enter the parameter(s) they would like to test the change in. *
confidence: Users should choose a confidence level among
0.9, 0.95, 0.99, 0.995 and 0.995. A smaller confidence level is easier
to reject the null hypothesis and detect a change-point.
The package also offers an option to plot the time series (this only
works for the univariate time series cases!) Users need to set
plot_SN = TRUE to visualize the time series, and
est_cp_loc = TRUE to add the estimated change-point
locations in the plot.
SN Test Statistic Plot: max_SNsweep
To visualize the computed SN-based test statistic at each time point,
users can apply the function max_SNsweep by plugging in the
output object from one of the functions SNSeg_Uni,
SNSeg_Multi and SNSeg_HD. The options
est_cp_loc = TRUE and critical_loc = TRUE are
provided to draw the estimated change-point locations and the critical
value threshold inside the test statistic plot.
max_SNsweep also returns the SN-based test statistic for
each time point. A large number of test statistics in the output can be
messy for users who only seek to generate a plot. To hide the test
statistics output, users can create an arbitrary variable name and set
it to the max_SNsweep function, e.g.,
SN_stat <- max_SNsweep(...).
Parameter Estimates of Each Segment Separated by the Detected
Change-Points: SNSeg_estimate
The function SNSeg_estimate allows users to compute the
parameter estimates (e.g., mean, variance, acf, quantile, etc.) of each
of the segments separated by the estimated change-points. To use this
function, users should use the output of the functions
SNSeg_Uni, SNSeg_Multi and
SNSeg_HD as the input of SNSeg_estimate.
S3 methods: summary, print and plot
The typical S3 methods summary, print and
plot are available to SNSeg_Uni,
SNSeg_Multi and SNSeg_HD objects. The
summary method displays the parameter to be tested, the
estimated change-point amount and locations, the grid_size,
confidence level as well as the critical value of the SN-based test. The
print method shows the change-point locations. The
plot method plots the time series, and similar to the
argument plot_SN = TRUE, the plot method
allows users to generate time series segmentation plot(s) with the
estimated change-point locations. It also provides ts_index
option to allow users to plot any individual time series they want.
Users can apply their preferred color of the change-point(s) within the
plot(s) by setting cpts.col to any color.
We then provide the examples of the functions SNSeg_Uni,
SNSeg_Multi and SNSeg_HD.
Examples of SNSeg_Uni:
The function SNSeg_Uni detect change-points for a
univariate time series based on the change in a single or
parameters.
- For a univariate time series, the input argument
paras_to_test allows one or a combination of multiple
parameters from “mean”, “variance”, “acf” and a numeric quantile value
between 0 and 1. For instance, to test the change in autocorrelation and
80th quantile, users should set paras_to_test to c(‘acf’,
0.8).
- To test the change in correlation between two univariate time
series, the input time series must be two-dimensional and
paras_to_test should be set to bivcor.
- Users can also set
paras_to_test using their own
parameters. In this case, the input of paras_to_test needs
to be a function which returns a numeric value.
We provide examples for different cases:
# Please run the following function before running examples:
mix_GauGPD <- function(u,p,trunc_r,gpd_scale,gpd_shape){
indicator <- u<p
rv <- rep(0, length(u))
rv[indicator>0] <- qtruncnorm(u[indicator>0]/p,a=-Inf,b=trunc_r)
rv[indicator<=0] <- qgpd((u[indicator<=0]-p)/(1-p), loc=trunc_r, scale=gpd_scale,shape=gpd_shape)
return(rv)
}
Test in a single parameter
Segmentation for Mean
set.seed(7)
n <- 2000
reptime <- 2
cp_sets <- round(n*c(0,cumsum(c(0.5,0.25)),1))
mean_shift <- c(0.4,0,0.4)
rho <- -0.7
ts <- MAR(n, reptime, rho)
no_seg <- length(cp_sets)-1
for(index in 1:no_seg){ # Mean shift
tau1 <- cp_sets[index]+1
tau2 <- cp_sets[index+1]
ts[tau1:tau2,] <- ts[tau1:tau2,] + mean_shift[index]
}
ts <- ts[,2]
# grid_size undefined
result <- SNSeg_Uni(ts, paras_to_test = "mean", confidence = 0.9,
grid_size_scale = 0.05, grid_size = NULL,
plot_SN = FALSE, est_cp_loc = FALSE)
# grid_size defined & generate time series segmentation plot
result <- SNSeg_Uni(ts, paras_to_test = "mean", confidence = 0.9,
grid_size_scale = 0.05, grid_size = 116,
plot_SN = TRUE, est_cp_loc = TRUE)
# Estimated change-point locations
result$est_cp
#> [1] 1018 1487
# Parameter estimates (mean) of each segment
SNSeg_estimate(result)
#> $mean
#> [1] 0.39444849 0.02107613 0.38362925
#>
#> attr(,"class")
#> [1] "SNSeg_estimate"
# plot the SN-based test statistic
SN_stat <- max_SNsweep(result, plot_SN = TRUE, est_cp_loc = TRUE,
critical_loc = TRUE)
We then show how to use the S3 methods summary,
print and plot.
summary(result)
#> There are 2 change-point(s) detected at 90th confidence level based on the change in the single mean parameter.
#>
#> The critical value of SN-based test is 135.39283594
#>
#> The detected change-point location(s) are 1018,1487 with a grid_size of 116
print(result)
#> The detected change-point location(s) are 1018,1487
plot(result, cpts.col = 'red')
We can see both the plot_SN = TRUE option and the
plot.SN function generates the same plot. Users can select
any choice based on their preferences. The class of the argument
result is SNSeg_Uni.
Segmentation for Variance
set.seed(7)
ts <- MAR_Variance(2, "V1")
ts <- ts[,2]
# grid_size defined
result <- SNSeg_Uni(ts, paras_to_test = "variance", confidence = 0.9,
grid_size_scale = 0.05, grid_size = NULL,
plot_SN = FALSE, est_cp_loc = TRUE)
# Estimated change-point locations
result$est_cp
# Parameter estimates (variance) of each segment
SNSeg_estimate(result)
# plot the SN-based test statistic
SN_stat <- max_SNsweep(result, plot_SN = TRUE, est_cp_loc = TRUE,
critical_loc = TRUE)
# built-in functions
# summary statistic
summary(result)
# print
print(result)
# plot
plot(result, cpts.col = 'red')
Segmentation for Autocorrelation
set.seed(7)
ts <- MAR_Variance(2, "A3")
ts <- ts[,2]
# grid_size defined
result <- SNSeg_Uni(ts, paras_to_test = "acf", confidence = 0.9,
grid_size_scale = 0.05, grid_size = 92, plot_SN = FALSE,
est_cp_loc = TRUE)
# Estimated change-point locations
result$est_cp
# Parameter estimates (acf) of each segment
SNSeg_estimate(result)
# plot the SN-based test statistic
SN_stat <- max_SNsweep(result, plot_SN = TRUE, est_cp_loc = TRUE,
critical_loc = TRUE)
# built-in functions
# summary statistic
summary(result)
# print
print(result)
# plot
plot(result, cpts.col = 'red')
Segmentation for bivariate correlation
library(mvtnorm)
set.seed(7)
n <- 1000
sigma_cross <- list(4*matrix(c(1,0.8,0.8,1), nrow=2),
matrix(c(1,0.2,0.2,1), nrow=2),
matrix(c(1,0.8,0.8,1), nrow=2))
cp_sets <- round(c(0,n/3,2*n/3,n))
noCP <- length(cp_sets)-2
rho_sets <- rep(0.5, noCP+1)
ts <- MAR_MTS_Covariance(n, 2, rho_sets, cp_sets, sigma_cross)
ts <- ts[1][[1]]
# grid_size defined
result <- SNSeg_Uni(ts, paras_to_test = "bivcor", confidence = 0.9,
grid_size_scale = 0.05, grid_size = 77,
plot_SN = FALSE, est_cp_loc = TRUE)
# Estimated change-point locations
result$est_cp
# Parameter estimates (bivariate correlation) of each segment
SNSeg_estimate(result)
# plot the SN-based test statistic
SN_stat <- max_SNsweep(result, plot_SN = TRUE, est_cp_loc = TRUE,
critical_loc = TRUE)
# built-in functions
# summary statistic
summary(result)
# print
print(result)
# plot
plot(result, cpts.col = 'red')
Segmentation for quantile
Please download the packages truncnorm and
evd before running the following code.
library(truncnorm)
library(evd)
set.seed(7)
n <- 1000
cp_sets <- c(0,n/2,n)
noCP <- length(cp_sets)-2
reptime <- 2
rho <- 0.2
# AR time series with no change-point (mean, var)=(0,1)
ts <- MAR(n, reptime, rho)*sqrt(1-rho^2)
trunc_r <- 0
p <- pnorm(trunc_r)
gpd_scale <- 2
gpd_shape <- 0.125
for(ts_index in 1:reptime){
ts[(cp_sets[2]+1):n, ts_index] <- mix_GauGPD(pnorm(ts[(cp_sets[2]+1):n, ts_index]),
p,trunc_r,gpd_scale,gpd_shape)
}
ts <- ts[,2]
# grid_size undefined
# test in 90% quantile
result <- SNSeg_Uni(ts, paras_to_test = 0.9, confidence = 0.9,
grid_size_scale = 0.066, grid_size = NULL,
plot_SN = FALSE, est_cp_loc = TRUE)
# Estimated change-point locations
result$est_cp
# Parameter estimates (90th quantile) of each segment
SNSeg_estimate(result)
# plot the SN-based test statistic
SN_stat <- max_SNsweep(result, plot_SN = TRUE, est_cp_loc = TRUE,
critical_loc = TRUE)
# built-in functions
# summary statistic
summary(result)
# print
print(result)
# plot
plot(result, cpts.col = 'red')
Test in a general functional
To perform the SN-based test using any general functional, users
should define their own function as the input of
paras_to_test. The input of this function should consist of
the followings:
ts: A univariate time series provided by the user.
The user-defined function should return a numeric value as the
output.
set.seed(7)
n <- 500
reptime <- 2
cp_sets <- round(n*c(0,cumsum(c(0.5,0.25)),1))
mean_shift <- c(0.4,0,0.4)
rho <- -0.7
ts <- MAR(n, reptime, rho)
no_seg <- length(cp_sets)-1
for(index in 1:no_seg){ # Mean shift
tau1 <- cp_sets[index]+1
tau2 <- cp_sets[index+1]
ts[tau1:tau2,] <- ts[tau1:tau2,] + mean_shift[index]
}
ts <- ts[,2]
# Define a general functional for the input 'paras_to_test'
paras_to_test = function(ts){
mean(ts)
}
result.SNCP.general <- SNSeg_Uni(ts, paras_to_test = paras_to_test,
confidence = 0.9, grid_size_scale = 0.05,
grid_size = NULL, plot_SN = FALSE,
est_cp_loc = TRUE)
# Estimated change-point locations
result.SNCP.general$est_cp
# Parameter estimates (general functional) of each segment
SNSeg_estimate(result.SNCP.general)
# plot the SN-based test statistic
SN_stat <- max_SNsweep(result.SNCP.general, plot_SN = TRUE, est_cp_loc = TRUE,
critical_loc = TRUE)
# built-in functions
# summary statistic
summary(result)
# print
print(result)
# plot
plot(result, cpts.col = 'red')
Test in multiple parameters
set.seed(7)
n <- 1000
cp_sets <- c(0,333,667,1000)
no_seg <- length(cp_sets)-1
rho <- 0
# AR time series with no change-point (mean, var)=(0,1)
ts <- MAR(n, 2, rho)*sqrt(1-rho^2)
no_seg <- length(cp_sets)-1
sd_shift <- c(1,1.6,1)
for(index in 1:no_seg){ # Mean shift
tau1 <- cp_sets[index]+1
tau2 <- cp_sets[index+1]
ts[tau1:tau2,] <- ts[tau1:tau2,]*sd_shift[index]
}
d <- 2
ts <- ts[,2]
# Test in 90th and 95th quantile with grid_size undefined
result <- SNSeg_Uni(ts, paras_to_test = c(0.9, 0.95), confidence = 0.9,
grid_size_scale = 0.05, grid_size = NULL,
plot_SN = FALSE, est_cp_loc = FALSE)
# Test in 90th quantile and the variance with grid_size undefined
result <- SNSeg_Uni(ts, paras_to_test = c(0.9, 'variance'),
confidence = 0.95, grid_size_scale = 0.078,
grid_size = NULL, plot_SN = FALSE,
est_cp_loc = FALSE)
# Test in 90th quantile, variance and acf with grid_size undefined
result <- SNSeg_Uni(ts, paras_to_test = c(0.9,'variance', "acf"),
confidence = 0.9, grid_size_scale = 0.064,
grid_size = NULL, plot_SN = TRUE,
est_cp_loc = TRUE)
# Estimated change-point locations
result$est_cp
# Test in 60th quantile, mean, variance and acf with grid_size defined
result.last <- SNSeg_Uni(ts, paras_to_test = c(0.6, 'mean', "variance",
"acf"), confidence = 0.9, grid_size_scale = 0.05,
grid_size = 83, plot_SN = FALSE, est_cp_loc = TRUE)
# Estimated change-point locations
result.last$est_cp
# Parameter estimates (of the last example) of each segment
SNSeg_estimate(result.last)
# SN-based test statistic segmentation plot
SN_stat <- max_SNsweep(result.last, plot_SN = TRUE, est_cp_loc = TRUE,
critical_loc = TRUE)
# built-in functions
# summary statistic
summary(result)
# print
print(result)
# plot
plot(result, cpts.col = 'red')
Examples: SNSeg_Multi
The function SNSeg_Multi detects change-points for
multivariate time series based on the change in multivariate means or
covariances. Users can set paras_to_test to
mean or covariance for change-points
estimation.
- The dimension of parameters to be tested must be at most ten. If the
parameter is
mean, the dimension is equivalent to the
number of time series. If the parameter is covariance, the
dimension is equivalent to the number of unknown parameters within the
covariance matrix.
For examples of multivariate means and covariances, please check the
following commands:
# Please run this function before running examples:
exchange_cor_matrix <- function(d, rho){
tmp <- matrix(rho, d, d)
diag(tmp) <- 1
return(tmp)
}
Segmentation for Multivariate Mean
library(mvtnorm)
set.seed(10)
d <- 5
n <- 1000
cp_sets <- round(n*c(0,cumsum(c(0.075,0.3,0.05,0.1,0.05)),1))
mean_shift <- c(-3,0,3,0,-3,0)/sqrt(d)
mean_shift <- sign(mean_shift)*ceiling(abs(mean_shift)*10)/10
rho_sets <- 0.5
sigma_cross <- list(exchange_cor_matrix(d,0))
ts <- MAR_MTS_Covariance(n, 2, rho_sets, cp_sets=c(0,n), sigma_cross)
noCP <- length(cp_sets)-2
no_seg <- length(cp_sets)-1
for(rep_index in 1:2){
for(index in 1:no_seg){ # Mean shift
tau1 <- cp_sets[index]+1
tau2 <- cp_sets[index+1]
ts[[rep_index]][,tau1:tau2] <- ts[[rep_index]][,tau1:tau2] + mean_shift[index]
}
}
ts <- ts[1][[1]]
# grid_size undefined
result <- SNSeg_Multi(ts, paras_to_test = "mean", confidence = 0.95,
grid_size_scale = 0.079, grid_size = NULL,
plot_SN = FALSE, est_cp_loc = TRUE)
# grid_size defined
result <- SNSeg_Multi(ts, paras_to_test = "mean", confidence = 0.99,
grid_size_scale = 0.05, grid_size = 65,
plot_SN = FALSE, est_cp_loc = TRUE)
# Estimated change-point locations
result$est_cp
# Parameter estimates (multivariate mean) of each segment
SNSeg_estimate(result)
# SN-based test statistic segmentation plot
SN_stat <- max_SNsweep(result, plot_SN = TRUE, est_cp_loc = TRUE,
critical_loc = TRUE)
# built-in functions
# summary statistic
summary(result)
# print
print(result)
# plot
par(mfrow=c(2,3))
plot(result, cpts.col = 'red')
Segmentation for Multivariate Covariance
library(mvtnorm)
set.seed(10)
reptime <- 2
d <- 4
n <- 1000
sigma_cross <- list(exchange_cor_matrix(d,0.2),
2*exchange_cor_matrix(d,0.5),
4*exchange_cor_matrix(d,0.5))
rho_sets <- c(0.3,0.3,0.3)
mean_shift <- c(0,0,0) # with mean change
cp_sets <- round(c(0,n/3,2*n/3,n))
ts <- MAR_MTS_Covariance(n, reptime, rho_sets, cp_sets, sigma_cross)
noCP <- length(cp_sets)-2
no_seg <- length(cp_sets)-1
for(rep_index in 1:reptime){
for(index in 1:no_seg){ # Mean shift
tau1 <- cp_sets[index]+1
tau2 <- cp_sets[index+1]
ts[[rep_index]][,tau1:tau2] <- ts[[rep_index]][,tau1:tau2] +
mean_shift[index]
}
}
ts <- ts[[1]]
# grid_size undefined
result <- SNSeg_Multi(ts, paras_to_test = "covariance",
confidence = 0.9, grid_size_scale = 0.05,
grid_size = NULL, plot_SN = FALSE, est_cp_loc = FALSE)
# grid_size defined
result <- SNSeg_Multi(ts, paras_to_test = "covariance",
confidence = 0.9, grid_size_scale = 0.05,
grid_size = 81, plot_SN = FALSE,
est_cp_loc = TRUE)
# Estimated change-point locations
result$est_cp
# Parameter estimates (covariance estimate) of each segment
SNSeg_estimate(result)
# SN-based test statistic segmentation plot
SN_stat <- max_SNsweep(result, plot_SN = TRUE, est_cp_loc = TRUE,
critical_loc = TRUE)
# built-in functions
# summary statistic
summary(result)
# print
print(result)
# plot
par(mfrow=c(1,1))
plot(result, cpts.col = 'red')
Examples: SNSeg_HD
The function SNSeg_HD performs change-points estimation
for high-dimensional time series (dimension greater than 10) using the
change in high-dimensional means. For usage examples of
SNSeg_HD, we provide the followig code:
n <- 600
p <- 100
nocp <- 5
cp_sets <- round(seq(0,nocp+1,1)/(nocp+1)*n)
num_entry <- 5
kappa <- sqrt(4/5) # Wang et al(2020)
mean_shift <- rep(c(0,kappa),100)[1:(length(cp_sets)-1)]
set.seed(1)
ts <- matrix(rnorm(n*p,0,1),n,p)
no_seg <- length(cp_sets)-1
for(index in 1:no_seg){ # Mean shift
tau1 <- cp_sets[index]+1
tau2 <- cp_sets[index+1]
ts[tau1:tau2,1:num_entry] <- ts[tau1:tau2,1:num_entry] +
mean_shift[index] # sparse change
}
# SN segmentation plot
# grid_size undefined (plot the first 3 time series)
result <- SNSeg_HD(ts, confidence = 0.9, grid_size_scale = 0.05,
grid_size = NULL, plot_SN = FALSE, est_cp_loc = TRUE,
ts_index = c(1:3))
# grid_size defined (plot the 1st, 3rd and 5th time series)
result <- SNSeg_HD(ts, confidence = 0.9, grid_size_scale = 0.05,
grid_size = 52, plot_SN = FALSE, est_cp_loc = TRUE,
ts_index = c(1,3,5))
# Estimated change-point locations
result$est_cp
# Parameter estimates (high-dimensional means) of each segment
summary.stat <- SNSeg_estimate(result)
# SN-based test statistic segmentation plot
SN_stat <- max_SNsweep(result, plot_SN = TRUE, est_cp_loc = TRUE,
critical_loc = TRUE)
# built-in functions
# summary statistic
summary(result)
# print
print(result)
# plot
plot(result, cpts.col = 'red', ts_index = c(1:3))
We note that only the selected time series were plotted by setting
values for ts_index.