# Finding the change point before a significant increase

I have the following time series:

df <- structure(c(288L, 259L, 265L, 293L, 271L, 278L, 300L, 286L, 278L,
275L, 282L, 285L, 290L, 296L, 296L, 279L, 270L, 292L, 283L, 289L,
280L, 269L, 289L, 290L, 287L, 271L, 280L, 299L, 278L, 287L, 293L,
286L, 297L, 281L, 285L, 305L, 288L, 295L, 277L, 292L, 286L, 281L,
287L, 302L, 292L, 297L, 292L, 279L, 281L, 291L), .Tsp = c(1961,
2010, 1), class = "ts")


This is the plot:

Questions/Problems

1. I would like to know, which point in the time series is the start of a significant trend. Applying a Mann-Kendall trend test for the whole time series gives a significant trend. However,by visual inspection, the time series does not have an increasing trend at least over the fisrt half.

2. This is like a changepoint problem, but I am not sure, which detection method is appropriate for this. I just want to determine the point in the time series before the most significant increasing trend.

3. I am thinking of a moving Mann-Kendall trend test and plotting the corresponding tau statistic. But I dont know how to set a threshold value and test the significance of the tau statistics.

Any suggestions on how to do this correctly in R?

I'll appreciate any help.

• Do you know why and when the change point occurs? Jul 18, 2018 at 2:02
• I actually dont know when the change point occurrs. I just want to know at which point in the series i can see the most significant change before an increasing trend. I'll modify my question above. Jul 18, 2018 at 2:09

You could use the EnvCpt package in R which fits mean, trend, AR and changepoint models. It then gives you the best fit of all models and you can use AIC (or another metric) to choose the best model.

library('EnvCpt')
out=envcpt(df)
out$summary # gives the fit and number of parameters for each model plot(out,type='aic') # plots the aic values, trendcpt+AR2 model is clearly the best out$trendar2cpt # gives the fit for trendcpt+AR2


The above gives the best model as a Trend+AR2+cpt model (assuming Normal errors) with a changepoint after 6 observations. I'm not sure you would want to fit a Trend+AR2 model to 6 observations though - but that is your call.

The next lowest AIC value is the Trend+AR2 model. Thus indicating that there is not a clear changepoint in this data. This could be that there is not enough data to be confident there is a changepoint at observation 6, or the change is too small (relative to the noise), or a combination.

• Thank you for this. I was thinking of another solution like a moving Mann-Kendall test. Like in a sliding correlation, the significance is tested using a t-statistic (fixed threshold, if the correlation exceeds this value then it is considered significant). Do you have any idea if its also possible to do the same in Mann-Kendall test? Jul 19, 2018 at 1:53
• You can do a moving window test but then you have to choose the window size that you use. Potentially you could get different changes with different window sizes and then how would you make inference? Jul 20, 2018 at 11:47
• The approach above considers the entire data with no choice of window size and the ability to incorporate a minimum distance between the changes through the minseglen argument. The idea is to make it as data driven as possible. Jul 20, 2018 at 11:48

Use mcp if you (1) want to quantify uncertainty about the location of the change point, and (2) want to specify a more informed model structure, e.g., that the first segment is a plateau. I arranged the data so that it is a regular data frame. Then fit an AR(1) model with a plateau + joined slope:

model = list(
y ~ 1 + ar(1),
~ 0 + x
)

library(mcp)
fit = mcp(model, data = df)


You can use summary(fit), plot(fit), and plot_pars(fit) to see the change point. Here's plot:

The distribution of the change point (bottom) is very broad because there is very little information about a change point here (no clear change and quite few data). This is not a weakness to the method - this should be an accurate inference, and other change point methods often completely ignore uncertainty.

You can help it by specifying more informed priors for the change point location and slope strengths in segment 2.

My bias is that there is no silver bullet here and the results apparently depend on the particular models or tools chosen. Both EnvCpt and mcp are great tools. Here are some additional perspectives from another changepoint and trend analysis package Rbeas in R (https://github.com/zhaokg/Rbeast):

df <- structure(c(288L, 259L, 265L, 293L, 271L, 278L, 300L, 286L, 278L,
275L, 282L, 285L, 290L, 296L, 296L, 279L, 270L, 292L, 283L, 289L,
280L, 269L, 289L, 290L, 287L, 271L, 280L, 299L, 278L, 287L, 293L,
286L, 297L, 281L, 285L, 305L, 288L, 295L, 277L, 292L, 286L, 281L,
287L, 302L, 292L, 297L, 292L, 279L, 281L, 291L), .Tsp = c(1961,
2010, 1), class = "ts")

library(Rbeast)
out = beast(df)
plot(out)
print(out)

#####################################################################
#                      Trend  Changepoints                          #
#####################################################################
.-------------------------------------------------------------------.
| Ascii plot of probability distribution for number of chgpts (ncp) |
.-------------------------------------------------------------------.
|Pr(ncp = 0 )=0.345|*****************************************       |
|Pr(ncp = 1 )=0.397|*********************************************** |
|Pr(ncp = 2 )=0.162|********************                            |
|Pr(ncp = 3 )=0.064|********                                        |
|Pr(ncp = 4 )=0.023|***                                             |
|Pr(ncp = 5 )=0.006|*                                               |
|Pr(ncp = 6 )=0.002|*                                               |
|Pr(ncp = 7 )=0.000|*                                               |
|Pr(ncp = 8 )=0.000|*                                               |
|Pr(ncp = 9 )=0.000|*                                               |
|Pr(ncp = 10)=0.000|*                                               |
.-------------------------------------------------------------------.
|    Summary for number of Trend ChangePoints (tcp)                 |
.-------------------------------------------------------------------.
|ncp_max    = 10   | MaxTrendKnotNum: A parameter you set           |
|ncp_mode   = 1    | Pr(ncp= 1)=0.40: There is a 39.7% probability  |
|                  | that the trend component has  1 changepoint(s).|
|ncp_mean   = 1.05 | Sum{ncp*Pr(ncp)} for ncp = 0,...,10            |
|ncp_pct10  = 0.00 | 10% percentile for number of changepoints      |
|ncp_median = 1.00 | 50% percentile: Median number of changepoints  |
|ncp_pct90  = 2.00 | 90% percentile for number of changepoints      |
.-------------------------------------------------------------------.
| List of probable trend changepoints ranked by probability of      |
| occurrence: Please combine the ncp reported above to determine    |
| which changepoints below are  practically meaningful              |
'-------------------------------------------------------------------'
|tcp#              |time (cp)                  |prob(cpPr)          |
|------------------|---------------------------|--------------------|
|1                 |1967.000000                |0.13171             |
|2                 |1986.000000                |0.09733             |
|3                 |1996.000000                |0.09254             |
|4                 |1977.000000                |0.08875             |
|5                 |2004.000000                |0.08242             |
.-------------------------------------------------------------------.


On average, there is 1 changepoint detected; its most likely location is 1967 with a posterior probability of 0.1317 (not too strong evidence). The prob subplot below shows the probability of changepoint occurrence over time. The sgnSlp subplot shows the tim-varying probabilities of slope being positive (red part), zero(green part), and negative (blue part). The three parts sum up to 1.0. The green part is the largest, that is, the slope is most likely to be zero. In addition, the read part is larger than the blue part, that is, if there is a non-zero slope/trend, it is more likely to be increasing than decreasing.