Skip to contents

All mosum code is commented out due to [#4].

All mcp code is commented out due to link 1 and link 2.

All bcp code is commented out due to [#5]

for (
  package in c(
    "ggplot2", "mvtnorm",
    # "bcp",
    "breakfast", "changepoint", "cpm", "CptNonPar", "ecp", "fpop", "gfpop",
    "InspectChangepoint", "jointseg", "mcp", "mosum", "not", "Rbeast",
    "segmented", "stepR", "strucchange", "VARDetect", "wbs"
  )
) {
  if (!requireNamespace(package, quietly = TRUE)) utils::install.packages(
    package, repos = "https://cloud.r-project.org", quiet = TRUE
  )
}

Data setup

Univariate mean change

# Univariate mean change
set.seed(1)
p <- 1
mean_data_1 <- rbind(
  mvtnorm::rmvnorm(300, mean = rep(0, p), sigma = diag(100, p)),
  mvtnorm::rmvnorm(400, mean = rep(50, p), sigma = diag(100, p)),
  mvtnorm::rmvnorm(300, mean = rep(2, p), sigma = diag(100, p))
)

plot.ts(mean_data_1)

Univariate mean and/or variance change

# Univariate mean and/or variance change
set.seed(1)
p <- 1
mv_data_1 <- rbind(
  mvtnorm::rmvnorm(300, mean = rep(0, p), sigma = diag(1, p)),
  mvtnorm::rmvnorm(400, mean = rep(10, p), sigma = diag(1, p)),
  mvtnorm::rmvnorm(300, mean = rep(0, p), sigma = diag(50, p)),
  mvtnorm::rmvnorm(300, mean = rep(0, p), sigma = diag(1, p)),
  mvtnorm::rmvnorm(400, mean = rep(10, p), sigma = diag(1, p)),
  mvtnorm::rmvnorm(300, mean = rep(10, p), sigma = diag(50, p))
)

plot.ts(mv_data_1)

Multivariate mean change

# Multivariate mean change
set.seed(1)
p <- 3
mean_data_3 <- rbind(
  mvtnorm::rmvnorm(300, mean = rep(0, p), sigma = diag(100, p)),
  mvtnorm::rmvnorm(400, mean = rep(50, p), sigma = diag(100, p)),
  mvtnorm::rmvnorm(300, mean = rep(2, p), sigma = diag(100, p))
)

plot.ts(mean_data_3)

Multivariate mean and/or variance change

# Multivariate mean and/or variance change
set.seed(1)
p <- 3
mv_data_3 <- rbind(
  mvtnorm::rmvnorm(300, mean = rep(0, p), sigma = diag(1, p)),
  mvtnorm::rmvnorm(400, mean = rep(10, p), sigma = diag(1, p)),
  mvtnorm::rmvnorm(300, mean = rep(0, p), sigma = diag(50, p)),
  mvtnorm::rmvnorm(300, mean = rep(0, p), sigma = diag(1, p)),
  mvtnorm::rmvnorm(400, mean = rep(10, p), sigma = diag(1, p)),
  mvtnorm::rmvnorm(300, mean = rep(10, p), sigma = diag(50, p))
)

plot.ts(mv_data_3)

Linear regression

# Linear regression
set.seed(1)
n <- 300
p <- 4
x <- mvtnorm::rmvnorm(n, rep(0, p), diag(p))
theta_0 <- rbind(c(1, 3.2, -1, 0), c(-1, -0.5, 2.5, -2), c(0.8, 0, 1, 2))
y <- c(
  x[1:100, ] %*% theta_0[1, ] + rnorm(100, 0, 3),
  x[101:200, ] %*% theta_0[2, ] + rnorm(100, 0, 3),
  x[201:n, ] %*% theta_0[3, ] + rnorm(100, 0, 3)
)
lm_data <- data.frame(y = y, x = x)

plot.ts(lm_data)

Logistic regression

# Logistic regression
set.seed(1)
n <- 500
p <- 4
x <- mvtnorm::rmvnorm(n, rep(0, p), diag(p))
theta <- rbind(rnorm(p, 0, 1), rnorm(p, 2, 1))
y <- c(
  rbinom(300, 1, 1 / (1 + exp(-x[1:300, ] %*% theta[1, ]))),
  rbinom(200, 1, 1 / (1 + exp(-x[301:n, ] %*% theta[2, ])))
)
binomial_data <- data.frame(y = y, x = x)

plot.ts(binomial_data)

Poisson regression

# Poisson regression
set.seed(1)
n <- 1100
p <- 3
x <- mvtnorm::rmvnorm(n, rep(0, p), diag(p))
delta <- rnorm(p)
theta_0 <- c(1, 0.3, -1)
y <- c(
  rpois(500, exp(x[1:500, ] %*% theta_0)),
  rpois(300, exp(x[501:800, ] %*% (theta_0 + delta))),
  rpois(200, exp(x[801:1000, ] %*% theta_0)),
  rpois(100, exp(x[1001:1100, ] %*% (theta_0 - delta)))
)
poisson_data <- data.frame(y = y, x = x)

plot.ts(log(poisson_data$y))

plot.ts(poisson_data[, -1])

Lasso

# Lasso
set.seed(1)
n <- 480
p_true <- 6
p <- 50
x <- mvtnorm::rmvnorm(n, rep(0, p), diag(p))
theta_0 <- rbind(
  runif(p_true, -5, -2),
  runif(p_true, -3, 3),
  runif(p_true, 2, 5),
  runif(p_true, -5, 5)
)
theta_0 <- cbind(theta_0, matrix(0, ncol = p - p_true, nrow = 4))
y <- c(
  x[1:80, ] %*% theta_0[1, ] + rnorm(80, 0, 1),
  x[81:200, ] %*% theta_0[2, ] + rnorm(120, 0, 1),
  x[201:320, ] %*% theta_0[3, ] + rnorm(120, 0, 1),
  x[321:n, ] %*% theta_0[4, ] + rnorm(160, 0, 1)
)
lasso_data <- data.frame(y = y, x = x)

plot.ts(lasso_data[, seq_len(p_true + 1)])

AR(3)

# AR(3)
set.seed(1)
n <- 1000
x <- rep(0, n + 3)
for (i in 1:600) {
  x[i + 3] <- 0.6 * x[i + 2] - 0.2 * x[i + 1] + 0.1 * x[i] + rnorm(1, 0, 3)
}
for (i in 601:1000) {
  x[i + 3] <- 0.3 * x[i + 2] + 0.4 * x[i + 1] + 0.2 * x[i] + rnorm(1, 0, 3)
}
ar_data <- x[-seq_len(3)]

plot.ts(ar_data)

GARCH(1, 1)

# GARCH(1, 1)
set.seed(1)
n <- 400
sigma_2 <- rep(1, n + 1)
x <- rep(0, n + 1)
for (i in seq_len(200)) {
  sigma_2[i + 1] <- 20 + 0.5 * x[i]^2 + 0.1 * sigma_2[i]
  x[i + 1] <- rnorm(1, 0, sqrt(sigma_2[i + 1]))
}
for (i in 201:400) {
  sigma_2[i + 1] <- 1 + 0.1 * x[i]^2 + 0.5 * sigma_2[i]
  x[i + 1] <- rnorm(1, 0, sqrt(sigma_2[i + 1]))
}
garch_data <- x[-1]

plot.ts(garch_data)

VAR(2)

# VAR(2)
set.seed(1)
n <- 800
p <- 2
theta_1 <- matrix(c(-0.3, 0.6, -0.5, 0.4, 0.2, 0.2, 0.2, -0.2), nrow = p)
theta_2 <- matrix(c(0.3, -0.4, 0.1, -0.5, -0.5, -0.2, -0.5, 0.2), nrow = p)
x <- matrix(0, n + 2, p)
for (i in 1:500) {
  x[i + 2, ] <- theta_1 %*% c(x[i + 1, ], x[i, ]) + rnorm(p, 0, 1)
}
for (i in 501:n) {
  x[i + 2, ] <- theta_2 %*% c(x[i + 1, ], x[i, ]) + rnorm(p, 0, 1)
}
var_data <- x[-seq_len(2), ]

plot.ts(var_data)

Univariate mean change

The true change points are 300 and 700. Some methods are plotted due to the un-retrievable change points.

fastcpd::fastcpd.mean(mean_data_1, r.progress = FALSE)@cp_set
#> [1] 300 700
CptNonPar::np.mojo(mean_data_1, G = floor(length(mean_data_1) / 6))$cpts
#> [1] 300 700
# Slow
strucchange::breakpoints(y ~ 1, data = data.frame(y = mean_data_1))$breakpoints
#> [1] 300 700
# Slower
ecp::e.divisive(mean_data_1)$estimates
#> [1]    1  301  701 1001
# Data need to be processed
changepoint::cpt.mean(c(mean_data_1))@cpts
#> [1]  300 1000
breakfast::breakfast(mean_data_1)$cptmodel.list[[6]]$cpts
#> [1] 300 700
wbs::wbs(mean_data_1)$cpt$cpt.ic$mbic.penalty
#> [1] 300 700
# Data need to be processed. `G` is selected based on the example
if (interactive()) {
  mosum::mosum(c(mean_data_1), G = 40)$cpts.info$cpts
}
fpop::Fpop(mean_data_1, nrow(mean_data_1))$t.est
#> [1]  300  700 1000
gfpop::gfpop(
  data = mean_data_1,
  mygraph = gfpop::graph(
    penalty = 2 * log(nrow(mean_data_1)) * gfpop::sdDiff(mean_data_1) ^ 2,
    type = "updown"
  ),
  type = "mean"
)$changepoints
#> [1]  300  700 1000
invisible(
  suppressMessages(
    capture.output(
      result_InspectChangepoint <- InspectChangepoint::inspect(
        t(mean_data_1),
        threshold = InspectChangepoint::compute.threshold(
          nrow(mean_data_1), ncol(mean_data_1)
        )
      )
    )
  )
)
result_InspectChangepoint$changepoints[, "location"]
#> [1] 300 700
jointseg::jointSeg(mean_data_1, K = 2)$bestBkp
#> [1] 300 700
Rbeast::beast(
  mean_data_1, season = "none", print.progress = FALSE, quiet = TRUE
)$trend$cp
#>  [1] 301 701 NaN NaN NaN NaN NaN NaN NaN NaN
stepR::stepFit(mean_data_1, alpha = 0.5)$rightEnd
#> [1]  300  700 1000
cpm::processStream(mean_data_1, cpmType = "Student")$changePoints
#> [1] 299 699
segmented::segmented(
  lm(y ~ 1 + x, data.frame(y = mean_data_1, x = seq_len(nrow(mean_data_1)))),
  seg.Z = ~ x
)$psi[, "Est."]
#> [1] 495
# Slowest
if (interactive()) {
  plot(
    mcp::mcp(
      list(y ~ 1, ~ 1, ~ 1),
      data = data.frame(y = mean_data_1, x = seq_len(nrow(mean_data_1))),
      par_x = "x"
    )
  )
}
plot(not::not(mean_data_1, contrast = "pcwsConstMean"))

# plot(bcp::bcp(mean_data_1))

Univariate mean and/or variance change

The true change points are 300, 700, 1000, 1300 and 1700. Some methods are plotted due to the un-retrievable change points.

fastcpd::fastcpd.mv(mv_data_1, r.progress = FALSE)@cp_set
#> [1]  300  700 1000 1300 1700
# Slow
ecp::e.divisive(mv_data_1)$estimates
#> [1]    1  301  701 1001 1301 1701 2001
# Data need to be processed
changepoint::cpt.meanvar(c(mv_data_1))@cpts
#> [1]  300 2000
CptNonPar::np.mojo(mv_data_1, G = floor(length(mv_data_1) / 6))$cpts
#> [1]  333  700 1300
cpm::processStream(mv_data_1, cpmType = "GLR")$changePoints
#>  [1]  293  300  403  408  618  621  696 1000 1021 1024 1293 1300 1417 1693 1981
invisible(
  suppressMessages(
    capture.output(
      result_InspectChangepoint <- InspectChangepoint::inspect(
        t(mv_data_1),
        threshold = InspectChangepoint::compute.threshold(
          nrow(mv_data_1), ncol(mv_data_1)
        )
      )
    )
  )
)
result_InspectChangepoint$changepoints[, "location"]
#>   [1]  300  700  702  704  707  708  712  715  716  717  718  721  722  723  726
#>  [16]  727  729  731  732  734  736  740  742  744  746  748  750  753  755  756
#>  [31]  757  759  760  762  764  765  768  771  772  774  776  777  784  785  786
#>  [46]  789  791  792  794  798  799  801  802  803  807  809  810  813  815  817
#>  [61]  819  826  827  828  829  831  833  835  836  837  838  840  841  842  843
#>  [76]  845  848  849  852  854  860  862  864  866  868  870  872  875  879  881
#>  [91]  884  886  887  888  889  896  897  898  899  901  903  904  905  906  909
#> [106]  912  913  915  917  919  921  922  923  927  928  932  934  936  937  940
#> [121]  944  945  947  948  949  958  959  961  962  963  964  966  967  968  972
#> [136]  974  976  978  979 1300 1700 1702 1703 1704 1705 1708 1710 1712 1714 1716
#> [151] 1717 1718 1721 1723 1725 1726 1731 1733 1735 1736 1737 1739 1742 1745 1747
#> [166] 1748 1752 1754 1756 1758 1759 1768 1770 1771 1778 1785 1790 1793 1795 1796
#> [181] 1797 1799 1800 1803 1804 1805 1806 1807 1808 1821 1828 1829 1833 1835 1837
#> [196] 1840 1841 1842 1848 1849 1851 1852 1854 1855 1857 1862 1863 1865 1867 1868
#> [211] 1876 1878 1879 1880 1882 1883 1884 1886 1887 1889 1894 1898 1899 1905 1906
#> [226] 1907 1908 1909 1912 1919 1926 1927 1928 1930 1933 1934 1935 1936 1938 1940
#> [241] 1952 1955 1956 1960 1962 1963 1966 1967 1969 1970 1974 1976 1977 1978 1980
#> [256] 1985 1987 1990 1997 1998
Rbeast::beast(
  mv_data_1, season = "none", print.progress = FALSE, quiet = TRUE
)$trend$cp
#>  [1] 1301 1986  301 1794 1855  703 1981 1769  709  713
# Slower
if (interactive()) {
  plot(
    mcp::mcp(
      list(y ~ 1, ~ 1, ~ 1, ~ 1, ~ 1, ~ 1),
      data = data.frame(y = mv_data_1, x = seq_len(nrow(mv_data_1))),
      par_x = "x"
    )
  )
}
plot(not::not(mv_data_1, contrast = "pcwsConstMeanVar"))

#> Press [enter] to continue

Multivariate mean change

The true change points are 300 and 700. Some methods are plotted due to the un-retrievable change points.

fastcpd::fastcpd.mean(mean_data_3, r.progress = FALSE)@cp_set
#> [1] 300 700
CptNonPar::np.mojo(mean_data_3, G = floor(nrow(mean_data_3) / 6))$cpts
#> [1] 300 700
invisible(
  suppressMessages(
    capture.output(
      result_InspectChangepoint <- InspectChangepoint::inspect(
        t(mean_data_3),
        threshold = InspectChangepoint::compute.threshold(
          nrow(mean_data_3), ncol(mean_data_3)
        )
      )
    )
  )
)
result_InspectChangepoint$changepoints[, "location"]
#> [1] 300 700
jointseg::jointSeg(mean_data_3, K = 2)$bestBkp
#> [1] 300 700
invisible(
  capture.output(
    result_Rbeast <- Rbeast::beast123(
      mean_data_3,
      metadata = list(whichDimIsTime = 1),
      season = "none"
    )
  )
)
result_Rbeast$trend$cp
#>       [,1] [,2] [,3]
#>  [1,]  301  701  701
#>  [2,]  701  301  301
#>  [3,]  142  685  182
#>  [4,]  770  NaN  760
#>  [5,]  291  NaN  707
#>  [6,]  NaN  NaN  NaN
#>  [7,]  NaN  NaN  NaN
#>  [8,]  NaN  NaN  NaN
#>  [9,]  NaN  NaN  NaN
#> [10,]  NaN  NaN  NaN
# Slow
strucchange::breakpoints(
  cbind(y.1, y.2, y.3) ~ 1, data = data.frame(y = mean_data_3)
)$breakpoints
#> [1] 300 700
# Slower
ecp::e.divisive(mean_data_3)$estimates
#> [1]    1  301  701 1001
# plot(bcp::bcp(mean_data_3))

Multivariate mean and/or variance change

The true change points are 300, 700, 1000, 1300 and 1700. Some methods are plotted due to the un-retrievable change points.

fastcpd::fastcpd.mv(mv_data_3, r.progress = FALSE)@cp_set
#> [1]  300  700 1000 1300 1700
# Slow
ecp::e.divisive(mv_data_3)$estimates
#> [1]    1  301  701 1001 1301 1701 2001
invisible(
  suppressMessages(
    capture.output(
      result_InspectChangepoint <- InspectChangepoint::inspect(
        t(mv_data_3),
        threshold = InspectChangepoint::compute.threshold(
          nrow(mv_data_3), ncol(mv_data_3)
        )
      )
    )
  )
)
result_InspectChangepoint$changepoints[, "location"]
#>   [1]  300  700  701  703  704  705  707  708  710  711  712  714  715  716  717
#>  [16]  719  721  722  723  725  726  730  731  732  734  735  736  738  739  740
#>  [31]  741  742  744  745  747  748  749  752  753  754  756  757  759  761  762
#>  [46]  763  764  766  768  769  770  772  773  774  776  777  778  780  781  783
#>  [61]  784  786  788  789  791  793  795  797  799  800  801  803  805  806  808
#>  [76]  809  811  813  814  816  817  819  821  822  824  825  826  828  829  830
#>  [91]  831  832  834  835  837  838  839  841  842  843  845  846  847  852  853
#> [106]  855  857  860  861  862  864  865  866  868  869  870  871  873  876  877
#> [121]  878  880  882  883  884  886  888  891  893  895  896  898  900  901  903
#> [136]  904  905  906  908  909  910  911  914  915  916  919  920  921  923  924
#> [151]  926  927  928  930  931  933  934  936  938  939  940  942  944  945  946
#> [166]  948  949  951  952  954  956  957  959  960  961  962  964  965  966  968
#> [181]  969  970  971  973  974  976  977  979  980  982  983  984  985  986  988
#> [196]  989  990  992  993  995  996  998 1000 1300 1700 1701 1702 1705 1708 1710
#> [211] 1711 1712 1713 1715 1716 1717 1718 1720 1721 1724 1725 1726 1728 1729 1730
#> [226] 1731 1733 1734 1735 1737 1738 1740 1742 1743 1745 1747 1749 1750 1752 1753
#> [241] 1754 1755 1756 1758 1759 1760 1762 1763 1764 1766 1767 1769 1770 1771 1773
#> [256] 1774 1776 1777 1779 1780 1782 1783 1785 1787 1788 1790 1791 1792 1793 1794
#> [271] 1796 1798 1800 1801 1803 1805 1809 1811 1812 1813 1814 1816 1817 1819 1820
#> [286] 1821 1823 1824 1826 1828 1830 1832 1833 1834 1836 1838 1839 1841 1842 1843
#> [301] 1846 1847 1850 1851 1853 1854 1856 1857 1859 1861 1862 1864 1865 1866 1869
#> [316] 1870 1872 1874 1876 1878 1879 1881 1882 1883 1885 1887 1888 1889 1890 1892
#> [331] 1894 1896 1898 1899 1900 1902 1904 1905 1907 1908 1909 1911 1913 1916 1918
#> [346] 1919 1920 1922 1923 1925 1926 1928 1930 1931 1933 1935 1936 1937 1939 1940
#> [361] 1942 1943 1945 1947 1948 1950 1951 1952 1953 1955 1956 1957 1959 1960 1961
#> [376] 1963 1964 1966 1967 1968 1970 1972 1974 1976 1977 1978 1982 1983 1985 1987
#> [391] 1990 1992 1993 1994 1996 1997
invisible(
  capture.output(
    result_Rbeast <- Rbeast::beast123(
      mv_data_3,
      metadata = list(whichDimIsTime = 1),
      season = "none"
    )
  )
)
result_Rbeast$trend$cp
#>       [,1] [,2] [,3]
#>  [1,]  301  301  702
#>  [2,] 1301  701 1301
#>  [3,]  701 1301  301
#>  [4,] 1978  830 1924
#>  [5,] 1994  807  833
#>  [6,]  877  989  845
#>  [7,]  820  881  769
#>  [8,]  861 1960 1971
#>  [9,]  814  874 1880
#> [10,] 1905 1970  775

Linear regression

The true change points are 100 and 200.

fastcpd::fastcpd.lm(lm_data, r.progress = FALSE)@cp_set
#> [1]  97 201
strucchange::breakpoints(y ~ . - 1, data = lm_data)$breakpoints
#> [1] 100 201
segmented::segmented(
  lm(
    y ~ . - 1,
    data.frame(y = lm_data$y, x = lm_data[, -1], index = seq_len(nrow(lm_data)))
  ),
  seg.Z = ~ index
)$psi[, "Est."]
#> [1] 233

Logistic regression

The true change point is 300.

fastcpd::fastcpd.binomial(binomial_data, r.progress = FALSE)@cp_set
#> [1] 301
strucchange::breakpoints(y ~ . - 1, data = binomial_data)$breakpoints
#> [1] 297

Poisson regression

The true change points are 500, 800 and 1000.

fastcpd::fastcpd.poisson(poisson_data, r.progress = FALSE)@cp_set
#> [1]  505  800 1003
# Slow
strucchange::breakpoints(y ~ . - 1, data = poisson_data)$breakpoints
#> [1] 935

Lasso

The true change points are 80, 200 and 320.

fastcpd::fastcpd.lasso(lasso_data, r.progress = FALSE)@cp_set
#> [1]  79 200 325
# Slow
strucchange::breakpoints(y ~ . - 1, data = lasso_data)$breakpoints
#> [1]  80 200 321

AR(3)

The true change point is 600. Some methods are plotted due to the un-retrievable change points.

fastcpd::fastcpd.ar(ar_data, 3, r.progress = FALSE)@cp_set
#> [1] 614
CptNonPar::np.mojo(ar_data, G = floor(length(ar_data) / 6))$cpts
#> integer(0)
segmented::segmented(
  lm(
    y ~ x + 1, data.frame(y = ar_data, x = seq_along(ar_data))
  ),
  seg.Z = ~ x
)$psi[, "Est."]
#> [1] 690
# Slow
if (interactive()) {
  plot(
    mcp::mcp(
      list(y ~ 1 + ar(3), ~ 0 + ar(3)),
      data = data.frame(y = ar_data, x = seq_along(ar_data)),
      par_x = "x"
    )
  )
}

GARCH(1, 1)

The true change point is 200.

fastcpd::fastcpd.garch(garch_data, c(1, 1), r.progress = FALSE)@cp_set
#> Registered S3 method overwritten by 'quantmod':
#>   method            from
#>   as.zoo.data.frame zoo
#> Residual calculation failed.
#> [1] 205
CptNonPar::np.mojo(garch_data, G = floor(length(garch_data) / 6))$cpts
#> [1] 206
strucchange::breakpoints(x ~ 1, data = data.frame(x = garch_data))$breakpoints
#> [1] NA

VAR(2)

The true change points is 500.

fastcpd::fastcpd.var(
  var_data, 2, cost_adjustment = NULL, r.progress = FALSE
)@cp_set
#> [1] 516
VARDetect::tbss(var_data)
#> [1] "first.brk.points:"
#>  [1] 140 196 252 308 364 420 476 532 588 644 700
#> [1] "selected lambda1:"
#> [1] 0.2107081
#> [1] "selected lambda2:"
#> [1] 0.02943525
#> [1] "second.brk.points:"
#> [1] 308 364 588
#> [1] "second.brk.points:"
#> [1] 308 476 588
#> [1] "second.brk.points:"
#> [1] 476 532
#> [1] "second.brk.points:"
#> [1] 476 532
#> [1] "second.brk.points:"
#> [1] 476 532
#> no refit for the parameter estimation
#> Estimated change points are: 501

Appendix: all code snippets

knitr::opts_chunk$set(
  collapse = TRUE, comment = "#>", eval = TRUE, cache = FALSE,
  warning = FALSE, fig.width = 8, fig.height = 5
)
for (
  package in c(
    "ggplot2", "mvtnorm",
    # "bcp",
    "breakfast", "changepoint", "cpm", "CptNonPar", "ecp", "fpop", "gfpop",
    "InspectChangepoint", "jointseg", "mcp", "mosum", "not", "Rbeast",
    "segmented", "stepR", "strucchange", "VARDetect", "wbs"
  )
) {
  if (!requireNamespace(package, quietly = TRUE)) utils::install.packages(
    package, repos = "https://cloud.r-project.org", quiet = TRUE
  )
}
# Univariate mean change
set.seed(1)
p <- 1
mean_data_1 <- rbind(
  mvtnorm::rmvnorm(300, mean = rep(0, p), sigma = diag(100, p)),
  mvtnorm::rmvnorm(400, mean = rep(50, p), sigma = diag(100, p)),
  mvtnorm::rmvnorm(300, mean = rep(2, p), sigma = diag(100, p))
)

plot.ts(mean_data_1)
# Univariate mean and/or variance change
set.seed(1)
p <- 1
mv_data_1 <- rbind(
  mvtnorm::rmvnorm(300, mean = rep(0, p), sigma = diag(1, p)),
  mvtnorm::rmvnorm(400, mean = rep(10, p), sigma = diag(1, p)),
  mvtnorm::rmvnorm(300, mean = rep(0, p), sigma = diag(50, p)),
  mvtnorm::rmvnorm(300, mean = rep(0, p), sigma = diag(1, p)),
  mvtnorm::rmvnorm(400, mean = rep(10, p), sigma = diag(1, p)),
  mvtnorm::rmvnorm(300, mean = rep(10, p), sigma = diag(50, p))
)

plot.ts(mv_data_1)
# Multivariate mean change
set.seed(1)
p <- 3
mean_data_3 <- rbind(
  mvtnorm::rmvnorm(300, mean = rep(0, p), sigma = diag(100, p)),
  mvtnorm::rmvnorm(400, mean = rep(50, p), sigma = diag(100, p)),
  mvtnorm::rmvnorm(300, mean = rep(2, p), sigma = diag(100, p))
)

plot.ts(mean_data_3)
# Multivariate mean and/or variance change
set.seed(1)
p <- 3
mv_data_3 <- rbind(
  mvtnorm::rmvnorm(300, mean = rep(0, p), sigma = diag(1, p)),
  mvtnorm::rmvnorm(400, mean = rep(10, p), sigma = diag(1, p)),
  mvtnorm::rmvnorm(300, mean = rep(0, p), sigma = diag(50, p)),
  mvtnorm::rmvnorm(300, mean = rep(0, p), sigma = diag(1, p)),
  mvtnorm::rmvnorm(400, mean = rep(10, p), sigma = diag(1, p)),
  mvtnorm::rmvnorm(300, mean = rep(10, p), sigma = diag(50, p))
)

plot.ts(mv_data_3)
# Linear regression
set.seed(1)
n <- 300
p <- 4
x <- mvtnorm::rmvnorm(n, rep(0, p), diag(p))
theta_0 <- rbind(c(1, 3.2, -1, 0), c(-1, -0.5, 2.5, -2), c(0.8, 0, 1, 2))
y <- c(
  x[1:100, ] %*% theta_0[1, ] + rnorm(100, 0, 3),
  x[101:200, ] %*% theta_0[2, ] + rnorm(100, 0, 3),
  x[201:n, ] %*% theta_0[3, ] + rnorm(100, 0, 3)
)
lm_data <- data.frame(y = y, x = x)

plot.ts(lm_data)
# Logistic regression
set.seed(1)
n <- 500
p <- 4
x <- mvtnorm::rmvnorm(n, rep(0, p), diag(p))
theta <- rbind(rnorm(p, 0, 1), rnorm(p, 2, 1))
y <- c(
  rbinom(300, 1, 1 / (1 + exp(-x[1:300, ] %*% theta[1, ]))),
  rbinom(200, 1, 1 / (1 + exp(-x[301:n, ] %*% theta[2, ])))
)
binomial_data <- data.frame(y = y, x = x)

plot.ts(binomial_data)
# Poisson regression
set.seed(1)
n <- 1100
p <- 3
x <- mvtnorm::rmvnorm(n, rep(0, p), diag(p))
delta <- rnorm(p)
theta_0 <- c(1, 0.3, -1)
y <- c(
  rpois(500, exp(x[1:500, ] %*% theta_0)),
  rpois(300, exp(x[501:800, ] %*% (theta_0 + delta))),
  rpois(200, exp(x[801:1000, ] %*% theta_0)),
  rpois(100, exp(x[1001:1100, ] %*% (theta_0 - delta)))
)
poisson_data <- data.frame(y = y, x = x)

plot.ts(log(poisson_data$y))
plot.ts(poisson_data[, -1])
# Lasso
set.seed(1)
n <- 480
p_true <- 6
p <- 50
x <- mvtnorm::rmvnorm(n, rep(0, p), diag(p))
theta_0 <- rbind(
  runif(p_true, -5, -2),
  runif(p_true, -3, 3),
  runif(p_true, 2, 5),
  runif(p_true, -5, 5)
)
theta_0 <- cbind(theta_0, matrix(0, ncol = p - p_true, nrow = 4))
y <- c(
  x[1:80, ] %*% theta_0[1, ] + rnorm(80, 0, 1),
  x[81:200, ] %*% theta_0[2, ] + rnorm(120, 0, 1),
  x[201:320, ] %*% theta_0[3, ] + rnorm(120, 0, 1),
  x[321:n, ] %*% theta_0[4, ] + rnorm(160, 0, 1)
)
lasso_data <- data.frame(y = y, x = x)

plot.ts(lasso_data[, seq_len(p_true + 1)])
# AR(3)
set.seed(1)
n <- 1000
x <- rep(0, n + 3)
for (i in 1:600) {
  x[i + 3] <- 0.6 * x[i + 2] - 0.2 * x[i + 1] + 0.1 * x[i] + rnorm(1, 0, 3)
}
for (i in 601:1000) {
  x[i + 3] <- 0.3 * x[i + 2] + 0.4 * x[i + 1] + 0.2 * x[i] + rnorm(1, 0, 3)
}
ar_data <- x[-seq_len(3)]

plot.ts(ar_data)
# GARCH(1, 1)
set.seed(1)
n <- 400
sigma_2 <- rep(1, n + 1)
x <- rep(0, n + 1)
for (i in seq_len(200)) {
  sigma_2[i + 1] <- 20 + 0.5 * x[i]^2 + 0.1 * sigma_2[i]
  x[i + 1] <- rnorm(1, 0, sqrt(sigma_2[i + 1]))
}
for (i in 201:400) {
  sigma_2[i + 1] <- 1 + 0.1 * x[i]^2 + 0.5 * sigma_2[i]
  x[i + 1] <- rnorm(1, 0, sqrt(sigma_2[i + 1]))
}
garch_data <- x[-1]

plot.ts(garch_data)
# VAR(2)
set.seed(1)
n <- 800
p <- 2
theta_1 <- matrix(c(-0.3, 0.6, -0.5, 0.4, 0.2, 0.2, 0.2, -0.2), nrow = p)
theta_2 <- matrix(c(0.3, -0.4, 0.1, -0.5, -0.5, -0.2, -0.5, 0.2), nrow = p)
x <- matrix(0, n + 2, p)
for (i in 1:500) {
  x[i + 2, ] <- theta_1 %*% c(x[i + 1, ], x[i, ]) + rnorm(p, 0, 1)
}
for (i in 501:n) {
  x[i + 2, ] <- theta_2 %*% c(x[i + 1, ], x[i, ]) + rnorm(p, 0, 1)
}
var_data <- x[-seq_len(2), ]

plot.ts(var_data)
fastcpd::fastcpd.mean(mean_data_1, r.progress = FALSE)@cp_set
CptNonPar::np.mojo(mean_data_1, G = floor(length(mean_data_1) / 6))$cpts
# Slow
strucchange::breakpoints(y ~ 1, data = data.frame(y = mean_data_1))$breakpoints
# Slower
ecp::e.divisive(mean_data_1)$estimates
# Data need to be processed
changepoint::cpt.mean(c(mean_data_1))@cpts
breakfast::breakfast(mean_data_1)$cptmodel.list[[6]]$cpts
wbs::wbs(mean_data_1)$cpt$cpt.ic$mbic.penalty
# Data need to be processed. `G` is selected based on the example
if (interactive()) {
  mosum::mosum(c(mean_data_1), G = 40)$cpts.info$cpts
}
fpop::Fpop(mean_data_1, nrow(mean_data_1))$t.est
gfpop::gfpop(
  data = mean_data_1,
  mygraph = gfpop::graph(
    penalty = 2 * log(nrow(mean_data_1)) * gfpop::sdDiff(mean_data_1) ^ 2,
    type = "updown"
  ),
  type = "mean"
)$changepoints
invisible(
  suppressMessages(
    capture.output(
      result_InspectChangepoint <- InspectChangepoint::inspect(
        t(mean_data_1),
        threshold = InspectChangepoint::compute.threshold(
          nrow(mean_data_1), ncol(mean_data_1)
        )
      )
    )
  )
)
result_InspectChangepoint$changepoints[, "location"]
jointseg::jointSeg(mean_data_1, K = 2)$bestBkp
Rbeast::beast(
  mean_data_1, season = "none", print.progress = FALSE, quiet = TRUE
)$trend$cp
stepR::stepFit(mean_data_1, alpha = 0.5)$rightEnd
cpm::processStream(mean_data_1, cpmType = "Student")$changePoints
segmented::segmented(
  lm(y ~ 1 + x, data.frame(y = mean_data_1, x = seq_len(nrow(mean_data_1)))),
  seg.Z = ~ x
)$psi[, "Est."]
# Slowest
if (interactive()) {
  plot(
    mcp::mcp(
      list(y ~ 1, ~ 1, ~ 1),
      data = data.frame(y = mean_data_1, x = seq_len(nrow(mean_data_1))),
      par_x = "x"
    )
  )
}
plot(not::not(mean_data_1, contrast = "pcwsConstMean"))
# plot(bcp::bcp(mean_data_1))
fastcpd::fastcpd.mv(mv_data_1, r.progress = FALSE)@cp_set
# Slow
ecp::e.divisive(mv_data_1)$estimates
# Data need to be processed
changepoint::cpt.meanvar(c(mv_data_1))@cpts
CptNonPar::np.mojo(mv_data_1, G = floor(length(mv_data_1) / 6))$cpts
cpm::processStream(mv_data_1, cpmType = "GLR")$changePoints
invisible(
  suppressMessages(
    capture.output(
      result_InspectChangepoint <- InspectChangepoint::inspect(
        t(mv_data_1),
        threshold = InspectChangepoint::compute.threshold(
          nrow(mv_data_1), ncol(mv_data_1)
        )
      )
    )
  )
)
result_InspectChangepoint$changepoints[, "location"]
Rbeast::beast(
  mv_data_1, season = "none", print.progress = FALSE, quiet = TRUE
)$trend$cp
# Slower
if (interactive()) {
  plot(
    mcp::mcp(
      list(y ~ 1, ~ 1, ~ 1, ~ 1, ~ 1, ~ 1),
      data = data.frame(y = mv_data_1, x = seq_len(nrow(mv_data_1))),
      par_x = "x"
    )
  )
}
plot(not::not(mv_data_1, contrast = "pcwsConstMeanVar"))
fastcpd::fastcpd.mean(mean_data_3, r.progress = FALSE)@cp_set
CptNonPar::np.mojo(mean_data_3, G = floor(nrow(mean_data_3) / 6))$cpts
invisible(
  suppressMessages(
    capture.output(
      result_InspectChangepoint <- InspectChangepoint::inspect(
        t(mean_data_3),
        threshold = InspectChangepoint::compute.threshold(
          nrow(mean_data_3), ncol(mean_data_3)
        )
      )
    )
  )
)
result_InspectChangepoint$changepoints[, "location"]
jointseg::jointSeg(mean_data_3, K = 2)$bestBkp
invisible(
  capture.output(
    result_Rbeast <- Rbeast::beast123(
      mean_data_3,
      metadata = list(whichDimIsTime = 1),
      season = "none"
    )
  )
)
result_Rbeast$trend$cp
# Slow
strucchange::breakpoints(
  cbind(y.1, y.2, y.3) ~ 1, data = data.frame(y = mean_data_3)
)$breakpoints
# Slower
ecp::e.divisive(mean_data_3)$estimates
# plot(bcp::bcp(mean_data_3))
fastcpd::fastcpd.mv(mv_data_3, r.progress = FALSE)@cp_set
# Slow
ecp::e.divisive(mv_data_3)$estimates
invisible(
  suppressMessages(
    capture.output(
      result_InspectChangepoint <- InspectChangepoint::inspect(
        t(mv_data_3),
        threshold = InspectChangepoint::compute.threshold(
          nrow(mv_data_3), ncol(mv_data_3)
        )
      )
    )
  )
)
result_InspectChangepoint$changepoints[, "location"]
invisible(
  capture.output(
    result_Rbeast <- Rbeast::beast123(
      mv_data_3,
      metadata = list(whichDimIsTime = 1),
      season = "none"
    )
  )
)
result_Rbeast$trend$cp
fastcpd::fastcpd.lm(lm_data, r.progress = FALSE)@cp_set
strucchange::breakpoints(y ~ . - 1, data = lm_data)$breakpoints
segmented::segmented(
  lm(
    y ~ . - 1,
    data.frame(y = lm_data$y, x = lm_data[, -1], index = seq_len(nrow(lm_data)))
  ),
  seg.Z = ~ index
)$psi[, "Est."]
fastcpd::fastcpd.binomial(binomial_data, r.progress = FALSE)@cp_set
strucchange::breakpoints(y ~ . - 1, data = binomial_data)$breakpoints
fastcpd::fastcpd.poisson(poisson_data, r.progress = FALSE)@cp_set
# Slow
strucchange::breakpoints(y ~ . - 1, data = poisson_data)$breakpoints
fastcpd::fastcpd.lasso(lasso_data, r.progress = FALSE)@cp_set
# Slow
strucchange::breakpoints(y ~ . - 1, data = lasso_data)$breakpoints
fastcpd::fastcpd.ar(ar_data, 3, r.progress = FALSE)@cp_set
CptNonPar::np.mojo(ar_data, G = floor(length(ar_data) / 6))$cpts
segmented::segmented(
  lm(
    y ~ x + 1, data.frame(y = ar_data, x = seq_along(ar_data))
  ),
  seg.Z = ~ x
)$psi[, "Est."]
# Slow
if (interactive()) {
  plot(
    mcp::mcp(
      list(y ~ 1 + ar(3), ~ 0 + ar(3)),
      data = data.frame(y = ar_data, x = seq_along(ar_data)),
      par_x = "x"
    )
  )
}
fastcpd::fastcpd.garch(garch_data, c(1, 1), r.progress = FALSE)@cp_set
CptNonPar::np.mojo(garch_data, G = floor(length(garch_data) / 6))$cpts
strucchange::breakpoints(x ~ 1, data = data.frame(x = garch_data))$breakpoints
fastcpd::fastcpd.var(
  var_data, 2, cost_adjustment = NULL, r.progress = FALSE
)@cp_set
VARDetect::tbss(var_data)