We consider a problem of detection, estimation, and segmentation of local nonlinear signals imbedded in a sequence of observations. As a model, we assume for u = 1, . . . , T
where t 1 < t 2 < . . . < t M define the locations and f the "shape," of the local signals.
For asymptotic analysis given below, we assume that the t k are scaled, so t k /T = t 0,k .
Initially we assume the X u are fixed, but for some applications they are random. The u are independent mean 0, normally distributed errors. For the moment we assume their variances are known and equal one, and discuss later how they should be estimated. The nuisance parameters µ(X u ) depend on the variable X u and may be constant or a parameterized regression function. They might, for example, be the mean values for a control group, for which it is assumed that all of the ξ j = 0. For our general theory we assume ρ is an unknown constant, but we estimate ρ differently in different problems.
An important special case that we return to in examples below is X u = u and µ u = α + β[(u − (T + 1)/2)/T ], so with this notation the model becomes
It is occasionally convenient to simplify several basic calculations by considering an alternative continuous time model, for which (2) can be conveniently written
for 0 ≤ u ≤ T , where dW defines white noise residuals.
Specific choices of the function f appear to be appropriate for a number of applications, some of which are discussed below. The special case f (x) = I{x > 0} is a frequently discussed "change-point" model, which has been applied to a variety of problems, usually with ρ assumed equal to 0. See, for example, Olshen et al. (2004) , Robbins, Gallagher, and Lund (2016) and references given there. Fang, Li and Siegmund (2018) and Fryzlewicz (2014) address a version of the problem that involves the possibility of multiple changepoints, and a major goal is segmentation of the observations to identify their number and locations, while controlling the probability of false positive detections. In this paper we extend the methods of Fang, Li and Siegmund (2018) to deal with more general signals and with observations that may have autoregressive dependence. See also Baranowski, Chen, and Fryzlewicz (2019) .
An interesting special case is a broken line regression model, where f (x) = x + , so the model is given by
This model has a long history, again primarily for the case of independent observations and at most one break-point in a regression line. An application to monitor kidney function after transplant has been discussed in a series of papers by A. F. M. Smith and others (e.g. Smith and Cook (1980) ). Davies (1987) and Knowles and Siegmund (1989) give theoretical analyses for independent observations. Toms and Lesperance (2003) Shin et al. (2013) , Schwartzman et al. (2013) ). In some applications it is often appropriate to add a scale parameter τ , so at u the signal located at t is of the form ξf [(u − t)/(τ T )] (cf. Siegmund and Worsley (1995) ). The model of "paired change-points"
of Olshen et al. (2004) can be described similarly, with f (x) = I{0 < x ≤ 1}.
In the special case that µ is a constant, X u = Y u−1 , and f (x) is the indicator that
x ≤ 0, (1) is the simplest example of a threshold autoregression, as introduced by Tong and studied by Chan and Tong and others (e.g., ) in numerous projects. In a still different context, X u may be a regressor, say a biomarker, and the local signal can be used to study whether a subset of individuals, defined by the value of that biomarker, differ from others in some respect, e.g., response to a treatment.
In some applications X u may be multidimensional or the noise distribution may not be normal, with the Poisson distribution representing a particularly interesting alternative.
Although the focus or our paper is on segmentation of multiple signals, we also develop confidence regions for t and jointly for t and ξ. After beginning with some basic calculations and a discussion of testing for at most one change in Section 2, segmentation is discussed in Section 3 and illustrated in 3.1. New approximations to control the false positive error rate given in Section 3, which combine smooth and not smooth stochastic processes, play 4 an important role. Sections 4 and 5 briefly consider a number of other special cases of our general framework. Additional theory and examples are contained in an online Supplement.
We have used simple models to minimize issues of over-fitting and data re-use that could arguably complicate interpretation of our segmentations. Once segmented, the models become simple linear models, so we also consider a second stage analysis to estimate more accurately the nuisance parameters and the ξ j , which then allows us to re-evaluate the accuracy of the segmentation.
Although large sample theory suggests use of maximum likelihood estimators to estimate the nuisance parameters ρ and σ 2 , we find that in some cases the signals themselves introduce large bias into those estimators. We discuss some ad hoc alternatives below.
Change.
To introduce our notation and provide some basic results, we begin with the model given by (2), with M = 1 or 0, and we consider a test of the hypothesis ξ = 0. The log likelihood is given by (ξ, t, θ)
where θ = (α, β, ρ) , and we consider a test based on the maximum with respect to t of the standardized score statistic
where ξ = d /dξ andθ is the maximum likelihood estimator of the nuisance parameters under the hypothesis ξ = 0; and σ 2 (t) is the asymptotic variance of the numerator.
We begin with the standard large sample expansion of the numerator in (5) given by
where I ·,· denotes elements of the Fisher information matrix and by the law of large numbers all quantities on the right hand side are evaluated at ξ = 0 and true values of the other parameters. This expansion is valid in large samples up to terms that are o(T ) in probability. To emphasize the structure of this approximation we find it convenient to use the notation E 0 to denote expectation under the hypothesis ξ = 0 and write (6) in the form
is not random although it depends on t; θ does not depend on t, and under the hypothesis ξ = 0 it has mean value 0 and covariance matrix A −1 = I θ,θ .
Remark. It may be shown that the decomposition given in (6) does not depend on the normality of the j , although we can no longer use the terminology of likelihood, efficient score, etc., and must rely on the central limit theorem to justify the asymptotic normality of probability calculations given below. The special case of Poisson observaions is discussed briefly in Section 5.
Calculations yield a number of simple propositions.
Proposition 1.
In particular the variance of V t is σ 2 (t) = Σ(t, t).
Remark. In the case the local signal contains a scale parameter, so it takes the form
, the covariance is similar in its general formulation, but slightly more complicated to compute in examples.
Writing E t,ξ to denote expectation under the alternative, where t and ξ can be vectors, we have
, and consider the test statistic max t |Z t |. The false positive probability for the test that detects a signal when max t |Z t | is large is given by
To evaluate this probability, we distinguish essentially two cases, one where f is continuous and the second where it is discontinuous. The case where f is the indicator of a half line or an interval (of unspecified length) is discussed in Fang, Li and Siegmund (2018) .
Here we assume that f is continuous and piecewise differentiable. An approximation based on Rice's formula is given by
whereŻ t denotes the derivative of Z at t. For the model where the local signal at u has the form ξf [(u − t)/τ T ], and we maximize over both t and a range of values of τ , a first order approximation becomes
whereŻ = ∇Z t,τ . This approximation can be improved by adding terms involving edge effects and corrections for curvature. The most important is the boundary correction at the minimum value of τ , which equals (8π) −1/2 ϕ(b)
. It is often convenient when T is large to ignore edge effects, which simplifies integration over 7 [T 0 , T 1 ]. See Siegmund and Worsley (1995) or Adler and Taylor (2007) for more detailed approximations.
To evaluate (8) and (9), the following result is useful.
Finally observe that by Propositions 1 and 2, we have the "matched filter" conditions, and consequently
Proposition 5. Given t (or t, τ ), the variable V t (or V t,τ ) is sufficient for ξ.
For the specific case of (4), it simplifies calculations somewhat to consider the continuous time version given in (3), which amounts to replacing certain sums by integrals, leading to
As expected, this result does not depend on the nuisance parameters α, β. It is perhaps surprising that it also does not depend on ρ. See the online Supplement for details.
It is possible to obtain a confidence region for the parameters t or a joint region for t and ξ. Calculation shows that the expected value of V s equals ξ times the covariance of V s and V t (cf. Proposition 2). An easy consequence is that Z t = V t /σ(t) is sufficient for ξ, so the conditional distribution of Z s given Z t is the same as it would be under the null hypothesis, ξ = 0. Now we must consider two cases, depending on whether the process has a jump at t.
The case of a jump was studied by Fang, Li and Siegmund (2018) , so here we consider 8 the continuous case. For a similar approach based on a different probability evaluation see Knowles, Siegmund and Zhang (1991) . We use the Kac-Slepian model process, which for a mean 0 unit variance stationary, differentiable, Gaussian process, say U s , gives a parabolic approximation for the process in a neighborhood of s 0 conditional on U s 0 assuming a large value. As a consequence of those arguments, in a neighborhood of s 0
with errors that converge to 0 in probability. It follows that the conditional distribution of max s (Z 2 s − Z 2 t ) is asymptotically χ 2 with one degree of freedom, so a 1 − α confidence region is given by the set of all s that satisfy Z 2 s ≥ Z 2 t − χ 2 1−α , where χ 2 1−α is the 1 − α quantile of the χ 2 distribution with one degree of freedom. This is in effect the same result one would obtain by inverting the log likelihood ratio statistic if Z were the log likelihood of a parameter satisfying standard regularity conditions. A joint confidence region for ξ and t is also easily obtained. To the condition that Z 2 t must be within a given distance of max s Z 2 s , we also require that |Z t − ξσ(t)| is sufficiently small. The pairs ξ, t that satisfy max δ [Z 2 t+δ − Z 2 t ] + (Z t − ξσ(t)) 2 ≤ c 2 provide a joint confidence region having the confidence coefficient
which equals the distribution of a χ 2 random variable with two degrees of freedom, again as if standard regularity conditions had been satisfied.
Looking ahead to the possibility of finding joint confidence regions for two (or more) change-points t 1 < t 2 , we face a more subtle argument. First note that with t = (t 1 , t 2 ), and ξ = (ξ 1 , ξ 2 ), a straightforward calculation extending the result given in Proposition
. This produces by a second calculation the conclusion that the log
. Hence the log likelihood ratio statistic, given t, for testing that ξ is the 0 vector equals
whereZ t 1 ,t 2 = Σ(t 1 , t 2 ) −1/2 (V t 1 , V t 2 ) has a standard bivariate normal distribution when ξ = 0, and || · || denotes the Euclidean norm. It follows from sufficiency and the Kac-Slepian argument used above that
has a χ 2 distribution with two degrees of freedom and can be used as above to obtain a joint confidence region for (t 1 , t 2 ).
A joint confidence region for t, ξ follows by an argment similar to that given above.
Remark. The confidence regions for the change-points are easy to visualize and seem reasonable. The estimates of ξ are more problematic, since selection of the change-point value t introduces bias into the estimation of ξ which may be exacerbated by a biased estimate of ρ. In Section 3 we consider a different method for estimating the parameters ξ. The maximum Z value is about 4.60, with ρ estimated to be 0.16, and occurs in 1969.
The appropriate (two-sided) p-value is about 6 × 10 −5 . See Figure 1 , where it appears that there might also be earlier change-points. Although the mercury thermometer was relatively young in 1659, those early years are interesting since they appear to involve the "little ice age." We return to these data below for our discussion of segmentation.
Meanwhile it may be interesting to note that an approximate 95% confidence region based on the detected change in 1969 is (1895, 1990) , which reflects the local shape of the plot in Figure 1 and suggests the possibility that an increase in temperature may have actually begun during the late 19th century. 2.3 Estimation of σ 2 and ρ: A Simulated Example.
Although we have assumed the variance known in our theoretical calculations, in our numerical studies we have estimated it by the residual mean square under the null hypothesis.
Other possibilities when the data are independent are sums of squares of first or second
, which may be preferable. First order differences remove most of the effect of changing mean values, and second order differences also mitigate the effect of slope changes.
Regarding correlation of the observations, a useful model is one that helps to control false positive errors with minimal loss of power; and to that end our working model is a first order autoregression. Using a value of ρ that is too small can lead to an increase in false positives. But when there are change-points, the maximum likelihood estimator of ρ under the null hypothesis that ξ = 0 can be upwardly biased, leading to a loss of power, as we can easily see from the simulated data described below.
The simulated example in Table 1 illustrates the problem of estimating ρ and the effect on the power to detect a change.
For the data in the third and fourth rows of the Table, if we assume it known that ρ = 0 and use second order differences to estimate σ 2 , the estimated value is 0.89 and the maximum Z value is 6.77. In all cases, the maximizing value of t is 94. In most other cases as well the maximizing value of t provides a reasonable estimator of the true value.
It is clear from these examples and others not shown that the existence of a change-point causes the estimator of ρ to be upwardly biased and results in a loss of power compared to using the true value of ρ. This is usually worse if there are multiple change-points. If a plot of the data indicates long stretches without change-points, one might use an estimator of ρ from that segment. For the data in the next to last row of Table 1 , an estimate of ρ based on 13 A second consideration in estimation of ρ is robustness of our procedure against more complex forms of dependence. Without exploring this question in detail, consider again the next to last row of Table 1 A different approach (e.g., Robbins, Gallagher, and Lund (2016) ) for dealing with dependence in the standard change-point model (i.e., f (s) = I{s > 0}) is to assume that the estimated residuals from a least squares fit of the null model form a locally dependent stationary process and estimate the autocovariances accordingly. Weak convergence arguments typically indicate that a rescaled process converges weakly to a Gaussian process, often related to a Brownian Bridge, which can be studied to look for change-points.
Although this method has some advantages, we prefer our approach for facilitating quantitative statements in the form of hypothesis tests and confidence intervals after a minimal amount of "data snooping" to arrive at a useful model. If we are successful in detecting the signal locations, our model for the other parameters becomes a standard linear model.
In Section 3.1 we discuss analysis of this linear model as a test of our approach.
We now turn to the problem of segmentation when there may be several change-points. It is helpful conceptually to think of two different cases. In the first the signals exhibit no particular pattern. In the second two or more changes are expected to produce a definite signature. For the problem of changes in levels of the process, paired changes frequently take the form of an abrupt departure from followed by a return to a baseline level.
Here we initially concentrate on the case of signals having no apparent pattern and consider versions of the two methods suggested in Fang, Li and Siegmund (2018) for detection of jump changes. First we consider a pseudo-sequential procedure, called Seq below, where it is convenient to write Z(t, T ) to denote the statistic Z(t) when the interval of observation
We also use a minimum sample size m 0 and a minimum lag n 0 , both of which we often set equal to 5. Let t 0 = m 0 and let T be the smallest integer exceeding t 0 + n 0 such that for some t 0 < t < T − n 0 , the value of Z(t, T ) exceeds a threshold b to be determined by a constraint on the probability of false positive error. The first change-point is taken to be t 1 : equal to arg max |Z t |, or the smallest or largest value of t for that smallest T . (In numerical experiments we have been unable to find a consistent preference.) The process is then iterated starting from t 1 . The nuisance parameter ρ is assumed to be known. In practice it is estimated from the data or more often from a subset of the data to mitigate the effect of the bias discussed above. The other nuisance parameters, are estimated from the current, evolving interval of observations, so that the nuisance parameters associated with one change-point do not confound detection of another. It would also be possible to estimate ρ from the currently studied subset of the data, but this estimator appears to be unstable.
To control the global false positive error rate, we want an approximation for
where m is the number of observations and m 0 , n 0 represent minimal sample sizes, both of which we frequently take to be 5. To emphasize that T is a variable quantity, we also write
A derivation of the approximation (15) begins from the initial approximation of the probability of interest by a sum over T of the sum over t ≤ T − n 0 of the integral over x from 0 to ∞ of the product of three factors:
. If we choose δ, so that δb is very small, the integral is dominated by small values of x. The third factor is approximately P{max j S j ≤ −x}, for a suitable random walk. For small x, it is approximately P{max j S j < 0}. Putting these two observations together leads to the approximation given by (15).
A slightly more complex calculation leads to a different and apparently preferable approximation, which is based on a similar analysis but incorporates in addition the conditional probability, given
T ) 2 and hence is an independent χ 2 1 random variable. This leads to an expression similar to (15), but with the integrand in (15) replaced by
In spite of its somewhat different appearance, this second approximation is slightly less conservative than the first, but gives very similar numerical values. For m 0 = n 0 = 3 and m = 250, the first approximation gives a 0.05 threshold of 4.08, while the second gives 3.97. Since these probabilities vary substantially for small values of m 0 , n 0 , and neighboring changes in slope seem difficult to detect and rarely to occur in the data we have studied, we typically use the values m 0 = n 0 = 5, which for this example would reduce the 0.05 threshold given by (16) from 3.97 to 3.89.
The approximation (16) When this process detects a change in slope, for the smallest value of T there typically is a cluster of values of t where |Z(t, T )| exceeds the chosen threshold. In practice we might reasonably choose the smallest such t, or the value that maximizes |Z(t, T )|. We then iterate the process starting from the selected t. Upon iteration the background linear regression is estimated using the data beginning with that t. It may not be clear which 18 choice is better in a particular application.
Remarks. (i) Although we use this "pseudo-sequential" method with a given amount of data, it can also be used for actual sequential detection, as in the problem of monitoring kidney transplants discussed in Smith and Cook (1980) .
(ii) Although sample paths of the process Z(t, T ) as a function of t are smooth, as a function of T they behave locally like a random walk. A consequence is that to obtain a mathematically precise derivation of the approximation (15), local perturbations of the paramter T must be scaled by a small factor, say ∆, which converges to 0 at the rate of Fang, Li and Siegmund (2018) for the case of abrupt changes in the mean. A similar remark applies to (17) given below.
A second method of segmentation, called MS below (for Maximum Score Statistic), uses
An asymptotic approximation to the null probability that this expression exceeds b can be computed by a modification of the argument behind (16), and we obtain the approximating expression
where β 0 (t, T 0 ) is defined similarly to β(t, T ). The evaluation of (17) can be simplified by a summation by parts and the observation that the various functions of three variables, (T 0 , t, T 1 ), occurring in (17) are in fact functions of two variables:
These methods determine a list of putative local signals together with generally overlapping backgrounds. We consider different algorithms for selecting a set of local signals and backgrounds from this list, with a preference for the additional constraint that no background overlap two local signals. In Fang, Li and Siegmund (2018) examples of particular interest were the shortest background, (cf. also Baranowski, Chen, and Fryzlewicz (2019)), and the largest Z value. For the second of these methods we have paid in advance for false positive errors, so we can also consider other methods for searching the list, e.g., subjective
consistency with a plot of the data. We also find it convenient computationally to follow the suggestion of Fryzlewicz (2014) by searching a random subset of intervals, which seems to work very well. When it appears that the number of changes is small, as in the examples considered in this paper, one might also use a combination of methods, e.g., searching a small random subset of intervals at a low threshold to generate candidates, then sorting those by hand with tests to detect at least one change at the appropriate higher threshold.
The methods discussed so far are "bottom up" methods in the sense that they attempt to identify one local signal at a time against a background appropriate for that signal. A popular "top down" method in the change-point literature scans all the data looking for a pair of change-points (e.g., Olshen et al. (2004) ). In some cases, especially if there is only one change to be detected, the method will "detect" two changes, as required, but it will put one of them near an end of the data sequence, where it can be recognized and ignored.
After an initial discovery of one or two change-points, the sequence is broken into two or three parts by those change-points, and the process is iterated as long as new change-points are identified. Here we also consider the statistic max s<t−h U s,t , where
and h is a parameter that represents a minimum distance between changes that we find interesting (usually taken to be 5 or 10 in the examples below). An appropriate threshold may be determined from the approximation (closely related to (9))
This approximation is approximately linear in the length of the sequence searched, so iterating the process leads in a final step to roughly the same probability of a false positive as is present in the first step. If we assume that in early iterations, only true positives are detected, then the iterative process does not create a multiple comparisons issue. As we will see in the examples below, this method has very good power properties, although its use requires somewhat more thoughtful analysis than the other methods. Since the process Z t can have relatively large local correlations, when there is only one change to be detected in an interval, and hence we expect to find a second "fake" change near an end-point of the interval searched, that "fake" change may be somewhat distant from the end-point, so it is not clear without additional analysis whether it should be ignored. Use of a linear analysis as in the following section can be very helpful.
A search algorithm based on the maximum of (18) If we use (18) The plot of Z t (Figure 2 ) contains a broad peak and a local maximum in the late 20th century, which suggests the possibility of more than one change. A 95% confidence interval based on the theory suggested above is [1866, 1972] , which leaves open the possibility of a second change. Seq and MS confirm a 19th century change, but they fail to detect one in the 20th century even if the autocorrelation is set to 0. With the estimated autocorrelation of 0.24, the statistic (18) detects changes in both 1897 and 1980. A linear analysis with these two change points yields and R 2 of 0.33 and an autocorrelation of 0.03, which does not test to be different from 0. Temperature anomalies of several other countries in Europe exhibit similar, but not identical behavior. In these cases the test to detect a single change will sometimes pick out a local maximum in the 19th century and sometimes in the 20th, while the statistic (18) 
We consider again the model (2), now with µ u an unknown constant and f (u) = 1{u > 0}, so the t k are change-points in the level of the observed process. For some applications it may seem more appropriate to use f (u) = 1{0 < u/τ ≤ 1}, so the length of a local interval where the mean changes is an unknown constant τ , which may vary from one interval to another when there are multiple changepoints. Since the case ρ = 0 has been widely studied, here we concentrate on time series, where the first order autoregressive dependence discussed above may be pertinent. Again it turns out that change-points introduce bias into the estimation of ρ, so we consider estimators based on subsets of the data that appear to be free of change-points.
We consider primarily methods taken from Fang, Li and Siegmund (2018), which, as noted above, motivate the structure of the bottom up segmentation methods suggested in this paper.
We also consider the top down method of circular binary segmentation (CBS) (Olshen et al. (2004) ), which may have advantages when a change that increases (decreases) the mean value tends to be followed by a change that decreases (increases) the mean, since it directly models such paired change-points. And finally we consider the analogue of (18), which was not studied in Fang, Li and Siegmund (2018) . To that end we take ρ = 0 for simplicity and let S t = t 1 Y u . Then V t = (tS m /m − S t ) , so U s,t is defined by the equation (18) for the approriately modified covariance matrix with entries s(1 − t/m) for s ≤ t. In this section we continue to refer to this as the statistic (18). Two interesting climate related time series are AMO (Atlantic Multidecadal Oscillation), for which there are data beginning in 1856, and PDO (Pacific Decadal Oscillation) with data from 1900. In both cases the data have been recorded monthly, but we average 27 the monthly data to get annual time series with 163 and 118 values, respectively.
A visual inspection of the AMO data suggests that there is no change-point before about 1900, so we estimate ρ from the first 45 observations to get the value 0.38. All methods agree that there is a decrease in the mean at 46 (1901) , an increase at 70 (1925), a decrease about at 108 (1963) , and the latest increase at 139 (1995) . A linear analysis with these change-points returns an R 2 of about 0.72 and an auto-correlation of 0. It may be interesting to observe that the changes in slope detected in the Northern Hemisphere ocean temperature anomalies (Example 6 in Section 3.1) fall very close to the mid-points of these three intervals.
The PDO data are similar, but changes appear to be more frequent, as the name suggests, and a cumulative sum plot appears very noisy. For these data level changes seem to alternate positive and negative directions, so CBS and related methods that look for paired changes appear to have an advantage. Although the pseudo-sequential method detected the same changes in the PDO data as CBS, its detections at 76 and 114 were borderline.
It is relatively straightforward to extend the analysis given above to generalized linear models, at least for independent observations. For example, suppose that Y u has the log likelihood
where
The efficient score to test the hypothesis ξ = 0 or segment the sequence is much as before, although computation of Ψ(t) is more complicated. Fortunately at least for the following examples, it seems plausible to assume independent (quasi)Poisson observations.
A perhaps more flexible model than (20) is an asymptotic quasi-likelihood model, e.g., Heyde (1997) , where there is the possibility of different models for the mean and variance. A very similar model to the one we developed in Section 2 is obtained by setting ψ(θ u ) = θ 2 u /2 in (20), but we interpret θ u as the conditional expectation of Y u given the observations up to u − 1, so the efficient score to detect a local signal is a sum of martingale differences. It is then straightforward to extend the methods developed above. We have studied several score statistics to detect local signals in the form of changes of level, slope, or (in the online Supplement) the autoregressive coefficient. To segment the observations, we consider two "bottom up" methods patterned after the methods of multiple change-point segmentation in Fang, Li and Siegmund (2018) and one "top down" method, defined in (18).
Estimation of the nuisance parameters, ρ and σ 2 , pose special problems. For the bottom up methods, in examining an interval of observations for possible change-points, intercept, level, slope and variance of the process are estimated locally, i.e., using only the data from the interval under consideration. For the autoregressive coefficient, our theory suggests using the (global) maximum likelihood estimator estimator under the hypothesis of no change. Since this estimator can be badly biased and result in a loss of power when there are change-points, we also consider ad hoc methods based on examination of different parts of the data. For use of (18), which estimates of σ 2 based on the (usually long) interval being searched, a possible method for dealing with variance heterogeneity arising from the early part of the data is to discard initial segments of varying length and consider the stability of the results obtained.
The multiple regression analsis suggested in Section 3.1, based on the assumption that detected break-points are correct, allows us to see if our segmentation is reasonable, estimate the magnitude of the detected signals, and reconsider our chosen value of ρ to see if we have used one that is too small. In most examples, it appears that we have used a larger value than necessary, which may result in a loss of power in the case of a large discrepancy.
A different approach would be to use initially the value ρ = 0, or some other small value to detect change points, apply the multiple regression analysis to discard changes that seem to be superfluous, and then iterate the process with the value of ρ suggested by the regression analysis. We have not chosen this path because our goal has been to use a value of ρ obtained with at most a small amount of "data snooping," in order to claim that due to minimal re-use of the data the false positive error rate has been adequately controlled.
Because our emphasis has been on detection of local signals, we have admittedly been For problems involving sparse, bump like changes as discussed briefly the Supplement, some superficial analysis suggests that our methods work well since most of the data are consistent with the null model, and hence global estimation of nuisance parameters does not pose a serious problem. In view of the ubiquity of genomic applications where "bumps" have different signatures suggested by a combination of science and experimental technique, and where the background may involve use of a control group it seems worthwhile to pursue a more systematic study.
Here we have considered signals in one-dimensional processes, where the number of possible "shapes" of the signals is relatively small. In view of the much larger variety of possible multi-dimensional signal shapes and the variety of approaches already existing in the literature a systematic comparative study of that problem may be valuable.
In this appendix we record basic formulas we have used in the paper. Assume the model is given by (2) 
or alternatively the solution of (3) with Y 0 = 0 is
Hence, under the null hypothesis that ξ = 0 we have
For the special case of broken line regression in continuous time the first and second coor-
The matrix A −1 = E( θ θ ) is straightforward to compute and somewhat tedious to invert. A very useful result is
It follows that the numerator of Z t , has covariance function
which does not depend on α, β, nor on ρ. In particular
To derive (15) and (17) (see below), we must consider T as variable and study the quantities σ 2 , Ψ, etc., as functions of both t and T . It is obvious from (25) that ∂σ 2 (t, T )/∂T does not depend on nuisance parameters. Since Ψ(t, T ) depends on the nuisance parameters (α, β, ρ) only in its third coordinate, if we first differentiate (23) with respect to T , then multiply on the right by Ψ(t, T ) it follows from (23) that neither (∂Ψ/∂T ) A T Ψ(t, T ) nor Ψ(t, T ) (∂A T /∂T )Ψ(t, T ) depends on nuisance paramters. The approximation (17) requires that we introduce T 0 < t as a third parameter. For broken line regression we also require
To study the behavior of Z t when ξ is different from 0, and to help with the interpretation of a plot of Z t , we put t = (t 1 , . . . , t k ), ξ = (ξ 1 , . . . , ξ k ) and use the notation E t,ξ (·). By combining the results given above we find that Figure 6 ). Since estimates of ρ based on segments of the data that appear not to involve changes were no more than 0.9, we tried using Seq with that value of suggests there may be a slope decrease in 1943, which together with the two other changes detected increases R 2 to 0.997. This change may be interesting insofar as it appears to coincide roughly with the hiatus in temperature increases in Greenland shown in Example 3 and the dip in northern hemisphere ocean anomalies mentioned in Example 10.
Appendix C: Bump Hunting: fixed shape, variable amplitude and scale.
In a variety of problems, the local signal to be detected occurs in a departure from, followed by a return to, a baseline value. In many cases the shape of the local signal arises naturally from the scientific context. An important example is inherited copy number variation, where there is an abrupt increase or decrease from the baseline value of two, which is followed quickly by a return to the baseline (e.g., Zhang et al. (2010) ). Other examples are ChIP-Seq, where a shape to expect is roughly triangular or double exponential (e.g., Shin et al. (2013) ) , or differential methylation, where a normal probability density function or an integral thereof might be reasonable (cf. ).
We modify (2) to include a scale paramter τ , to obtain the log likelihood function
Here f is a positive, symmetric integrable function, e.g., the square root of (i) a standard normal density function, (ii) a triangular probability density on [−1, 1], (iii) a double exponential probability density, or (iv) a uniform probability density function on [−1/2, 1/2].
Consider the case of at most one bump and assume that the search interval is long enough relative to the scale parameter that it is reasonable to ignore end effects. We define our standardized statistic as above, but we now write it as Z t,τ to reflect the unknown location t and width τ of the bump. An approximation to its false positive error probability is given in display (9). (For (iv) above a different approximation is required because of the discontinuities in f . See Fang, Li and Siegmund (2018) .)
Often bumps are relatively sparse, which makes estimation of ρ and σ 2 relatively easy, so one can see intervals of the data that can be used for estimating autocorrelation and other nuisance parameters without serious bias.
Although a serious study of these methods is warranted by the variety of potential applications, we do not discuss these in detail here, which would require a lengthy investigation of the specific scientific features of those applications. In both the South Korean and the China data, the net slope change is close to zero. It may be interesting to consider a (quasi)Poisson log linear analysis, which usually agrees with the linear least squares analysis. The dispersion parameter often indicates large overdispersion. We have also considered a detection strategy using a square root transformation of the dependent variable, which would be expected to stabilize Poisson like variability.
Simulations suggest that this technique can be useful, but also that heteroscedasticity is usually not a serious problem.
As a final illustration of our methods we consider a simple case of Threshold Autoregression, which has been studied in a number of papers by Tong and colleagues.
For a model assume
To test ξ = 0, in the notation of Section 2
and asymptotically under the null hypothesis we have
and A −1 is the 2 × 2 matrix with entries a 11 = T, a 12 = T E(Y ), a 22 = T E(Y 2 ), where Y denotes the stationary distribution of Y u under the null hypothesis (and the stationarity assumption that |ρ| < 1). Hence
tions show that the covariance of the standardized score statistics Z s and Z t are given by [G(min(s, t) ) − Ψ (s)AΨ(t)]/σ(s)σ(t). Hence locally for small δ
where the functions G, Ψ andĠ are evaluated at t, and slightly more generally Cov[(Z t+δ 1 , Z t+δ 2 )|Z t ] ≈ min(δ 1 , δ 2 )[G − Ψ AΨ] −1Ġ . This allows us to calculate an approximation to P{max t |Z t | ≥ b}. The methods of, e.g., Woodroofe (1976) or Yakir (2013) lead after some calculation to
We have evaluated (31) by numerical integration, and the p-values given below come from that evaluation with the observed mean value and standard deviation. We have also implemented an empirical version of this approximation, where the entire computation is based on the appropriately estimated quantities. These two approximations are in rough agreement for our examples, where the p-value is very small. Simulations suggest that the Type I error control of (31) is adequate under the model The power also seems reasonable.
In principle, one should learn something from the maximizing point of the statistic Z s relative to the estimated values of µ and of σ, but in simulations the maximizing value of shows more variability than we can easily interpret. A calculation similar to that giving (30) indicates that in the case that ξ = 0 the expected value of Z t is approximately ξT 1/2 σ(t), as it would be for likelihood theory with standard regularity conditions. This approximation seems relatively stable and suggests a rough approximation for ξ by equating the approximation to the observed value of max s Z s (or min s Z s ).
With a slight modification in principle, but substantially more detailed calculation in application, one can consider higher order autoregressions. Suppose, for example, the model is
Now X t is a two-dimensional vector, while G(t), Ψ(t), and σ 2 (t) = Σ t , say, are 2 × 2
matrices. An appropriate test statistic has the form max t V t Σ −1 t V t . Putting Z t = Σ −1/2 V t , we can control the false positive probability based on the approximation
where e = (cos(ω), sin(ω)) , and where π −1 times the integral over ω equals the trace of Σ −1 tĠ (t).
Simulations suggest that the Type I error control of this procedure is adequate under the model. The power also seems reasonable. In principle, one should learn something from the maximizing point of the statistic Z t relative to the estimated values of µ and of σ, but in simulations the maximizing value of t shows more variability than we can easily interpret.
Examples. As illustrative applications, we follow in considering the Canadian lynx data and Nicholson's blowfly data. We do not, however, try to give a complete discussion of these well studied data.
The Canadian lynx data consist of annual counts over a 114 year period, which are considered a surrogate for the size of the lynx population. The scientific reasoning behind a two-phase model is the hypothesis that when the lynx population is small, it finds a more 48 than adequate food supply and increases until it reaches a size where the food supply is inadequate, when it then decreases until the cycle begins again. Empirical observations suggest that the decrease in times of food shortage is more rapid than the increase in times of abundance. If we use a first order auto regressive model, these fluctuations suggest a positive autoregressive parameter in times of food abundance, which decreases when there is a food shortage. Although there are clear outliers in the data, which has lead others to consider transformations to shorten the tails, we analyze the original data. Employing the first order autoregressive model with a possible change in the autoregressive parameter, as suggested above yields null estimators of E(Y ) and Var(Y ) of 1538 and 1560, respectively.
The maximum value of Z t is 3.89, which occurs at about t = 3500, and the resulting (two-sided) p-value based on (31) is approximately 0.004.
Nicholson's blowfly data (Brillinger (et al. 1980 ) and the web site www.stat.berkeley.edu/ brill/blowfly97I.html contains 361 observations in four columns, labeled respectively "births," "nonemerging," "emerging," and "deaths." The first, third, and fourth columns vary from numbers close to 0 to numbers in the several thousands. The second column typically involves substantially smaller numbers. For the first column. the estimated mean value, autocorrelation, and standard deviation under the null model are 1438, 0.6 and 1670, respectively. The maximum Z value is about 4.73, which occurs for t ≈ 5300. The p-value for the hypothesis of no change is approximately 2 × 10 −4 . For "emerging" the maximum Z-value is about 5.22, which gives a p-value of about 6 × 10 −5 . For these data the change in the autocorrelation is positive.