检测时间序列中的变化点需要构建一个强大的全局 ARIMA 模型（在您的情况下，模型变化和参数随时间变化肯定存在缺陷），然后识别该模型参数中最重要的变化点。使用您的 509 值，最显着的变化点是在 353 周期左右。我使用了 AUTOBOX（我帮助开发）中可用的一些专有算法，这些算法可能会被许可用于您的定制应用程序。基本思想是将数据分成两部分，并在找到最重要的变化点后，分别重新分析两个时间范围中的每一个（1-352；353-509），以确定两组中的每一个内的进一步变化点。重复此操作，直到您有 k 个子集。我已经附上了使用这种方法的第一步。

我对这个问题的解释是，OP 正在寻找适合所提供示例形状的方法，而不是 HAC 残差。此外，需要不需要大量人工或分析师干预的自动化例程。Box-Jenkins 可能不合适，尽管他们在这个线程中强调，因为他们确实需要大量的分析师参与。

R 模块用于这种类型的非基于矩的模式匹配。排列分布聚类是由马克斯普朗克研究所的科学家开发的一种模式匹配技术，符合您概述的标准。它的应用是时间序列数据，但不限于此。这是已开发的 R 模块的引用：

pdc：Andreas Brandmaier 的基于复杂性的时间序列聚类的 R 包

除了 PDC，还有机器学习，由加州大学欧文分校的 Eamon Keogh 开发的 iSax 例程也值得比较。

最后，还有这篇关于数据粉碎的论文：Uncovering Lurking Order in Data查托帕迪耶和利普森。除了聪明的标题之外，还有一个严肃的目的在起作用。以下是摘要：“从自动语音识别到发现不寻常的星星，几乎所有自动发现任务的基础是能够相互比较和对比数据流，识别连接和发现异常值。然而，尽管数据很普遍，但自动化方法没有跟上步伐。一个关键瓶颈是当今大多数数据比较算法依赖于人类专家来指定数据的哪些“特征”与比较相关。在这里，我们提出了一种估计任意来源之间相似性的新原则数据流，既不使用领域知识也不使用学习。我们展示了将这一原则应用于分析来自许多现实世界具有挑战性问题的数据，包括消除与癫痫发作有关的脑电图模式的歧义，从心音记录中检测异常心脏活动以及从原始光度法对天文物体进行分类。在所有这些情况下，在没有获得任何领域知识的情况下，我们展示了与领域专家设计的专业算法和启发式算法所达到的准确性相当的性能。我们建议，数据粉碎原则可能会为理解日益复杂的观察打开大门，尤其是当专家不知道该寻找什么时。” 在所有这些情况下，在没有获得任何领域知识的情况下，我们展示了与领域专家设计的专业算法和启发式算法所达到的准确性相当的性能。我们建议，数据粉碎原则可能会为理解日益复杂的观察打开大门，尤其是当专家不知道该寻找什么时。” 在所有这些情况下，在没有获得任何领域知识的情况下，我们展示了与领域专家设计的专业算法和启发式算法所达到的准确性相当的性能。我们建议，数据粉碎原则可能会为理解日益复杂的观察打开大门，尤其是当专家不知道该寻找什么时。”

这种方法远远超出了曲线拟合。值得一试。

我认为该线程的标题具有误导性：您不是要比较密度函数，而是实际上是在寻找时间序列中的结构中断。但是，您并没有指定这些结构性中断是应该在滚动时间窗口中发现还是通过查看时间序列的总历史而在事后发现。从这个意义上说，您的问题实际上与此重复：什么方法可以检测时间序列的结构中断？

正如 Rob Hyndman 在此链接中所提到的，R 为此目的提供了 strucchange 包。我玩弄了你的数据，但我必须说结果令人失望[第一个数据点真的是 4 还是应该是 54？]：

raw = c(54,53,53,53,53,58,56,52,49,52,56,51,44,39,39,39,37,33,27,21,18,12,19,30,45,66,92,118,135,148,153,160,168,174,181,187,191,190,191,192,194,194,194,193,193,201,200,199,199,199,197,193,190,187,176,162,157,154,144,126,110,87,74,57,46,44,51,60,65,66,90,106,99,87,84,85,83,91,95,99,101,102,102,103,105,110,107,108,135,171,171,141,120,78,42,44,52,54,103,128,82,103,46,27,73,123,125,77,24,30,27,36,42,49,32,55,20,16,21,31,78,140,116,99,58,139,70,22,44,7,48,32,18,16,25,16,17,35,29,11,13,8,8,18,14,0,10,18,2,1,4,0,61,87,91,2,0,2,9,40,21,2,14,5,9,49,116,100,114,115,62,41,119,191,190,164,156,109,37,15,0,5,1,0,0,2,4,2,0,48,129,168,112,98,95,119,125,191,241,209,229,230,231,246,249,240,99,32,0,0,2,13,28,39,15,15,19,31,47,61,92,91,99,108,114,118,121,125,129,129,125,125,131,135,138,142,147,141,149,153,152,153,159,161,158,158,162,167,171,173,174,176,178,184,190,190,185,190,200,199,189,196,197,197,196,199,200,195,187,191,192,190,186,184,184,179,173,171,170,164,156,155,156,151,141,141,139,143,143,140,146,145,130,126,127,127,125,122,122,127,131,134,140,150,160,166,175,192,208,243,251,255,255,255,249,221,190,181,181,181,181,179,173,165,159,153,162,169,165,154,144,142,145,136,134,131,130,128,124,119,115,103,78,54,40,25,8,2,7,12,25,13,22,15,33,34,57,71,48,16,1,2,0,2,21,112,174,191,190,152,153,161,159,153,71,16,28,3,4,0,14,26,30,26,15,12,19,21,18,53,89,125,139,140,142,141,135,136,140,159,170,173,176,184,180,170,167,168,170,167,161,163,170,164,161,160,163,163,160,160,163,169,166,161,156,155,156,158,160,150,149,149,151,154,156,156,156,151,149,150,153,154,151,146,144,149,150,151,152,151,150,148,147,144,141,137,133,130,128,128,128,136,143,159,180,196,205,212,218,222,225,227,227,225,223,222,222,221,220,220,220,220,221,222,223,221,223,225,226,227,228,232,235,234,236,238,240,241,240,239,237,238,240,240,237,236,239,238,235)
raw = log(raw+1)
d = as.ts(raw,frequency = 12)
dd = ts.intersect(d = d, d1 = lag(d, -1),d2 = lag(d, -2),d3 = lag(d, -3),d4 = lag(d, -4),d5 = lag(d, -5),d6 = lag(d, -6),d7 = lag(d, -7),d8 = lag(d, -8),d9 = lag(d, -9),d10 = lag(d, -10),d11 = lag(d, -11),d12 = lag(d, -12))

(breakpoints(d ~d1 + d2+ d3+ d4+ d5+ d6+ d7+ d8+ d9+ d10+ d11+ d12, data = dd))
>Breakpoints at observation number:
>151 
>Corresponding to breakdates:
>163 

(breakpoints(d ~d1 + d2, data = dd))
>Breakpoints at observation number:
>95 178 
>Corresponding to breakdates:
>107 190 

我不是该软件包的常规用户。如您所见，这取决于您适合数据的模型。你可以尝试

library(forecast)
auto.arima(raw)

为您提供最合适的 ARIMA 模型。