R/nsp_tvreg.R

In nsp: Inference for Multiple Change-Points in Linear Models

#' Narrowest Significance Pursuit algorithm with general covariates
#'
#' This function runs the Narrowest Significance Pursuit (NSP) algorithm on data sequence \code{y} and design matrix \code{x}
#' to return localised regions (intervals) of the domain in which the parameters of the linear regression model y_t = beta(t) x_t + z_t significantly
#' depart from constancy (e.g. by containing change-points), at the global significance level \code{alpha}. For any interval considered by the algorithm, 
#' significant departure from parameter constancy is achieved if the multiscale
#' deviation measure (see Details for the literature reference) exceeds a threshold, which is either provided as input
#' or determined from the data (as a function of \code{alpha}).
#' The function works best when the errors z_t in the linear regression formulation y_t = beta(t) x_t + z_t are independent and
#' identically distributed Gaussians.
#' 
#' The NSP algorithm is described in P. Fryzlewicz (2021) "Narrowest Significance Pursuit: inference for multiple change-points in linear
#' models", preprint.
#' 
#' @param y A vector containing the data sequence being the response in the linear model y_t = beta(t) x_t + z_t.
#' @param x The design matrix in the regression model above, with the regressors as columns.
#' @param M The minimum number of intervals considered at each recursive stage, unless the number of all intervals is smaller, in which case all intervals
#' are used.
#' @param thresh.val Numerical value of the significance threshold (lambda in the paper); or \code{NULL} if the threshold is to be determined from
#' the data (see \code{thresh.type}).
#' @param sigma The standard deviation of the errors z_t; if \code{NULL} then will be estimated from the data via the MOLS estimator described in the paper.
#' @param alpha Desired maximum probability of obtaining an interval that does not contain a change-point (the significance threshold will be 
#' determined as a function of this parameter).
#' @param power A parameter for the MOLS estimator of sigma; the span of the moving window in the MOLS estimator is \code{min(n, max(round(n^power), min.size))},
#' where \code{n} is the length of \code{y}.
#' @param min.size (See immediately above.)
#' @param overlap If \code{FALSE}, then on discovering a significant interval, the search continues recursively to the left and to the right of that 
#' interval. If \code{TRUE}, then the search continues to the left and to the right of the midpoint of that interval.
#' @return A list with the following components:
#' \item{intervals}{A data frame containing the estimated intervals of significance: \code{starts} and \code{ends} is where the intervals start and end, 
#'	respectively; \code{values} are the values of the deviation measure on each given interval; \code{midpoints} are their midpoints.}
#' \item{threshold.used}{The threshold value.}
#' @author Piotr Fryzlewicz, \email{p.fryzlewicz@@lse.ac.uk}
#' @seealso \code{\link{nsp}}, \code{\link{nsp_poly}}, \code{\link{nsp_poly_ar}}, \code{\link{nsp_selfnorm}}, \code{\link{nsp_poly_selfnorm}}
#' @examples
#' set.seed(1)
#' f <- c(1:100, 100:1, 1:100)
#' y <- f + stats::rnorm(300) * 15
#' x <- matrix(0, 300, 2)
#' x[,1] <- 1
#' x[,2] <- seq(from = 0, to = 1, length = 300)
#' nsp_tvreg(y, x, 100)
#' @importFrom stats rnorm
#' @export




nsp_tvreg <- function(y, x, M = 1000, thresh.val = NULL, sigma = NULL, alpha = 0.1, power = 1/2, min.size = 20, overlap = FALSE) {
			
	n <- length(y)
	
	if (is.null(thresh.val)) {
		
		if (is.null(sigma)) sigma <- est_var(y, x, power, min.size)
		if (is.null(alpha)) alpha <- 0.1
		
		thresh.val <- sigma * thresh_kab(n, alpha)
	}
	
	nsp(y, x, M, thresh.val, overlap)
	
}