I'm performing a time series analysis with multiple breakpoints in R.
I managed to identify three breakpoints using the procedure suggested in strucchange package but I'm struggling to get the significance (p-value) for these break points.
Here there is a dummy dataset and the code I was working with. the dataset:
x=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,
29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,
54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,
79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94)
z=c(128,103,29,117,53,49,84,67,76,111,81,38,36,-35,-12,21,121,38,84,173,153,99,
91,110,69,50,15,-50,15,-97,-2,13,107,47,137,25,-19,54,4,87,72,58,32,-4,75,50,
80,65,124,56,58,-3,30,42,55,212,245,18,106,128,88,216,205,234,120,171,195,230,
237,143,225,253,202,218,283,227,291,192,179,197,337,259,261,215,290,293,255,
316,355,312,337,341,388,338)
df=data.frame(z,x)
plot(x,z)
the code:
# https://www.marinedatascience.co/blog/2019/09/28/comparison-of-change-point-detection-methods/
library(strucchange)
library(sandwich)
library(fxregime)
# get best model
opt_bpts <- function(x) {
#x = bpts_sum$RSS["BIC",]
n <- length(x)
lowest <- vector("logical", length = n-1)
lowest[1] <- FALSE
for (i in 2:n) {
lowest[i] <- x[i] < x[i-1] & x[i] < x[i+1]
}
out <- as.integer(names(x)[lowest])
return(out)
}
#################################################################
#marinedatascience.co/blog/2019/09/28/comparison-of-change-point-detection-methods/
#Zeileis, A., Leisch, F., Hornik, K. & Kleiber, C. (2002), strucchange: An R Package for Testing for Structural Change in Linear Regression Models. J Stat Softw 7(2), 38p., doi: 10.18637/jss.v007.i02↩
z_ts <- as.ts(df$z) #crate time series
bpts <- breakpoints(z ~ x, data = df)
plot(bpts)
bpts_sum <- summary(bpts)
opt_brks <- opt_bpts(bpts_sum$RSS["BIC",])
opt_brks
# Nested syntax with 3 breaks:
ci=confint(bpts,breaks = 3)#, level = 0.99)
bpts <- breakpoints(breakpoints(z ~ x, data = df), breaks = 3)
Fst=Fstats(z_ts~1)
#here I get a p-value for the analysis with three breakpoints
plot(Fst)
sctest(Fst)
I get as Fst output:
supF test
data: Fst
sup.F = 203.23, p-value < 2.2e-16
I would like to obtain (if it's possible) the p-value of each breakpoint. Something like this:
F test
data: Fst
breakpoint1:p-value < brk1.pvalue
breakpoint2:p-value < brk2.pvalue
breakpoint3:p-value < brk3.pvalue