# Identify two major change points in time series data and summarise for several replicates in R

I'm working with temporal simulations (forward in time) and I would like to find the best way to detect the two major change points in this time series.

I have several replicates of the same simulation, then my idea was to collect the two change points identified for each replicate and average 1) the first change point among replicates; and 2) the second change point among replicates.

I found the package tsoutliers but does not find what I would expect.

This is an example of one replicate. The data are not continuous in the sense that I'm collecting data with an interval time of 50.

    data<-c(0.01516377,0.01426401,0.01293092,0.01197672,0.01777328,0.01513655,0.01464312,0.01134877,0.01394351,
0.01841415,0.01465141,0.01662577,0.01431855,0.01709126,0.01427008,0.01633598,0.01592693,0.01428456,
0.01478259,0.01409589,0.01477847,0.02530393,0.03750182,0.03212164,0.02343575,0.03001415,0.04005386,
0.03490057,0.03762527,0.03798278,0.02948298,0.03047519,0.03823259,0.03836810,0.04254081,0.04268210,
0.03375728,0.03724411,0.03615817,0.04067189,0.05058244,0.04777934,0.04200034,0.04249145,0.05403044,
0.04748773,0.04695485,0.04226534,0.03882807,0.04164517,0.04562888,0.03819429,0.04346850,0.04188714,
0.04042580,0.03863569,0.02794996,0.03403659,0.04118528,0.03850582,0.04560576,0.02849807,0.04056327,
0.04398593,0.03483359,0.03847835,0.03761654,0.03682544,0.03955815,0.03892210,0.04148052,0.04137503,
0.04147201,0.03277047,0.03959029,0.04737270,0.03709050,0.04513550,0.05063302,0.05353608,0.04153208,
0.04001323,0.04014528,0.05278923,0.04249564,0.04494758,0.04896572,0.04406259,0.04368667,0.04159133,
0.04218148,0.04632765,0.04816140,0.03725149,0.05110241,0.04445239,0.04348772,0.03297161,0.03249867,
0.04658435,0.04358190)

dat.ts<- ts(data,frequency=1)
library(tsoutliers)
data.ts.outliers <- tso(dat.ts, types = c("AO", "LS", "TC"))
data.ts.outliers
plot(data.ts.outliers)


Do you have any suggestion on which method would be the best?

The other option I was considering is to obtain the mean and variance for each time point and make pairwise statistical test for consecutive time points.

My first suggestion would be to use a package for changepoint analysis rather than outlier detection. An example would be the changepoint package in R. Here is the code for the analysis of the above data:

library(changepoint)
out=cpt.meanvar(dat.ts)
plot(out)
cpts(out)


The mean and variance of your data is changing so it makes sense to use the cpt.meanvar function.

In terms of a multivariate comparison, I would advise using a multivariate changepoint detection approach. If you want the changes in the same place in all series then you can use the InspectChangepoint R package which appears to work well.

• I followed your suggestion and I used cpt.meanvar. It find TWO change points, corresponding to the points where I see big differences. I tried to use also the InspectChangepoint R package, but it did not work as good. It finds a lot change points. To overcome this issue my idea was to run cpt.meanvar() on each replicate independently and then average the time value of each change point. At the end I have an average of the first change point (value 1) and an average of the second change point (value 2). What do you think? – CafféSospeso Aug 31 '17 at 17:23
• It is unfortunate that the InspectChangepoint package didn't work out of the box, have you tried changing the lambda value to be larger? As lambda increases you get less changes. You also might want to set M=100 or something similar as M=0 is the default which uses the Binary Segmentation algorithm whereas M>0 uses the Wild Binary Segmentation algorithm which is preferred. – adunaic Sep 2 '17 at 18:32
• Oh and I suppose you can average the first and second changes if you want to but i'm not sure what that tells you as the average probably isn't optimal in any sense. – adunaic Sep 2 '17 at 18:33
• I realised that with cpt.meanvar I can find the changing point but it is not so sensitive on detecting when the "changing phase" stops. In other terms, this function give me the point in which the mean variance change, but it doesn t give the point corresponding to the end of the "changing phase". Now i'm looking to some method that allows me to extract change rates, may be with some fitting and extraction of the parameters obtained from the fitting. What do you think? – CafféSospeso Sep 3 '17 at 13:41