I have time series data that I need to transform into a stationary series. It is a monthly series so according to this graph, I think the break point is in September, 2007. However, I am not sure if I should go with my hunch or find a module in R to get the accurate break point.

enter image description here


There are a lot of R packages. Here is a non-exhaustive overview with worked examples. Below, I'll show a solution using the mcp package.

Let's start by simulating some data that look a bit like yours:

df = data.frame(x = 1:150)
df$y = c(arima.sim(n = 75, list(ar = 0.8), sd = 2) + 0, 
         arima.sim(n = 75, list(ar = 0.5), sd = 2) + 7)

You can fit a model consisting of a two plateaus with 1st-order autoregressive residuals.

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

Maybe you just want one autoregressive coefficient that is shared for both segments. In that case, just delete ar(1) in the formula for the second segment - that will make the first ar(1) carry over.

Now we fit the model:

fit = mcp(model, data = df, par_x = "x")

Sometimes you want do more and faster sampling via mcp(..., iter = 10000, cores = 3).

Let's see the result visually:


enter image description here

It seems to capture the data OK. The black dots are the raw data. The grey lines are samples from the posterior distribution. The blue distribution is the inferred change point. In this case, the change in intercept is probably the clearest "signal", but this article on time-series analysis with mcp shows examples without such abrupt changes in the mean.

You can see posteriors for each parameter using plot_pars(fit). Let's also see if it recovered the parameters well:

> summary(fit)
Population-level parameters:
    name  mean lower upper Rhat n.eff
   ar1_1  0.74  0.56  0.91  1.0  7131
   ar1_2  0.61  0.43  0.80  1.0  6133
    cp_1 75.48 72.05 76.00  1.2   147
   int_1  0.64 -1.11  2.50  1.0  6288
   int_2  6.56  5.22  7.97  1.0  4948
 sigma_1  2.09  1.86  2.35  1.0  6743

We simualted the data from ar1_1 = 0.8, ar1_2 = 0.5, cp_1 = 75, int_1 = 0, int_2 = 7, and sigma_1 = 2. So it looks OK.

This was a very simple model with intercepts-only and a single change point. Several other packages are suitable here as well, including segmented and EnvCpt. Strengths of mcp include (1) it quantifies uncertainty of the change points, (2) it can do per-segment regression modeling, and (3) it can handle more complex autoregressive structures, among other things. The main limitations are (1) it's slower - but not unreasonably slow, (2) it requires JAGS (which is free).


The heretofore referenced "changepoint package in R" assumes that there is no arima memory i.e. the data is independent and normally distributed see https://www.jstatsoft.org/article/view/v058i03 section 4.1 . If the data is not independent then the underlying arima model needs to be identified but no provision is made for the identification of this critical component just requiring the user to pre-specify it.

Furthermore if the user has some notion about causal predictors and possibly their contemporary and lag effects that might need to be conditioned for … notsomuch … which limits it's general usefulness.

Breakpoints can occur in the arima component over time and certainly error variance change points might also need to be conditioned for see http://docplayer.net/12080848-Outliers-level-shifts-and-variance-changes-in-time-series.html . One needs to KNOW the assumptions under which each and every proposed solution is valid. All that glitters is not gold !


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.