尝试在R中使用stl和decompose函数时出现错误

18

我做了一个简单的时间序列,给正弦函数添加了一点噪声,并尝试使用R中的“stl”和“decompose”函数对其进行分解。尽管我的序列明显具有多个周期且是周期性的,但这两个函数都给出了以下错误:

x
  [1]  1.4537365796  2.7185844368  2.8394728999  3.8926989923  4.3405508086  5.1959080871
  [7]  5.6602505790  5.4829985648  5.6357660330  4.6084976233  4.6617322922  4.0286486832
 [13]  3.3641752333  1.7408063182  0.8815147612  0.2895139342 -0.5402768515 -1.5612641107
 [19] -2.1584502547 -2.9878043526 -3.5545638149 -4.0530074199 -4.0748538612 -4.7581704662
 [25] -4.6555349052 -4.0726206240 -3.1646413472 -2.6934453823 -2.2364605277 -1.2643569882
 [31] -0.1202011946  1.1136371449  2.2504199271  3.0313528996  3.5384449109  4.5176211013
 [37]  5.4013172839  5.4252837451  5.4768196692  5.8979709077  5.6698285659  4.5133489450
 [43]  4.2702602998  3.5180837069  2.2652913344  1.1975595698  0.5412697849 -0.5966162032
 [49] -1.0827728340 -1.8488242277 -3.4118061838 -3.9009752140 -3.9102671954 -4.3486102172
 [55] -4.7481017993 -4.0097598695 -3.9078554267 -3.8070416888 -2.5968567322 -2.2567568949
 [61] -1.1423907008  0.0002492447  0.4338279080  1.2431986797  2.3216397323  3.3235925116
 [67]  4.1591487169  4.9028481873  5.4535861470  5.0579349546  5.1548777627  4.9707124992
 [73]  5.4496833187  4.4563072487  4.1301372986  2.4594352788  1.7253019929  0.6961453965
 [79]  0.4281167695 -1.3152944759 -1.8645880957 -2.5764132038 -3.7681528112 -4.3731672862
 [85] -3.9940201611 -4.5497596299 -4.9496796983 -4.1233093447 -3.7759837204 -3.3359027749
 [91] -2.3518009102 -1.7488933797 -0.7225148838  0.5395759836  1.0496249652  2.0383715782
 [97]  3.2357449979  3.8028316517  5.0346866280  5.2154265148

fit<- stl(x, t.window=15, s.window="per", robust=TRUE)
Error in stl(x, t.window = 15, s.window = "per", robust = TRUE) :series is not periodic or has less      than two periods

fit<- decompose(x,type="multiplicative")
Error in decompose(x, type = "multiplicative") :time series has no or less than 2 periods

有人能帮我解决这个问题吗?


2
一如既往,第一条规则是:绘制数据图以确保其符合您的预期。 - Carl Witthoft
3
尝试运行 stl(ts(x,freq=10), t.window=15, s.window="per", robust=TRUE)。其中,ts(x,freq=10) 表示输入数据是一个时间序列,采样频率为 10;t.window=15 表示趋势分量的窗口长度为 15;s.window="per" 表示季节分量的窗口长度根据数据周期自适应;robust=TRUE 表示使用鲁棒性较强的方法进行拟合。 - Ben Bolker
1个回答

18

从错误信息中不是显而易见吗:

time series has no or less than 2 periods

?R并不是告诉你你的数据没有周期性,而是你传递给它的数据没有表明它们具有周期性。

从R的角度来看,你的时间序列x没有周期性。看起来你忘记告诉ts()周期是多少,或者在告知时出现了错误。由于你没有展示x是如何创建的,所以我们除了建议你回去创建一个至少包含两个周期的x之外,没有更多可做的了。

这里的关键是,单独使用R不能推断出每个单位时间内观察的频率。你必须告诉ts()这个信息。你可以用以下几种方式实现:

frequency:每个时间单位的观测次数。

deltat: 每个连续观测之间的采样时间分数;例如月度数据为1/12。只需提供frequencydeltat中的一个。

如果你没有提供其中之一,ts()会使用默认值frequency = 1deltat = 1,这将指示一个每个单位时间的观测(例如每年一个观测)的时间序列。

stl()需要一个“ts”类的对象——如果你没有提供这样的对象,它将通过as.ts()强制将输入数据转换为“ts”对象。此函数将使用我上面描述的默认值。

我认为在这里发生的情况是,你没有意识到stl()需要一个“ts”类对象,也没有意识到当你只提供观测向量时,它会为你的数据创建一个不合适的对象。

解决方法是通过ts()显式创建“ts”类的对象,并指定frequencydeltat之一。

例如:

dat <- cumsum(rnorm(12*4))
x <- ts(dat)
stl(x, "periodic")
xx <- ts(dat, frequency = 12)
stl(xx, "periodic")

我在这里使用 frequency = 12 表示每个单位时间内有12个观测值,例如月度数据。

以上内容是给我看的。

R> stl(x, "periodic")
Error in stl(x, "periodic") : 
  series is not periodic or has less than two periods

R> stl(xx, "periodic")
 Call:
 stl(x = xx, s.window = "periodic")

Components
       seasonal    trend remainder
Jan 1 -0.103529  0.55245  -0.44301
Feb 1  0.001333  0.56981   0.86135
Mar 1 -0.382075  0.58717   1.11162
Apr 1  0.010552  0.59891  -1.04966
....

针对您的数据,我怀疑您希望frequency = 10,鉴于时间序列的长度; 这意味着您每年有十个观测值。如果该系列每年有更多的观测值,比如说月度数据的12个观测值,但是你没有最后两个或者前两个月的数据(即没有缺失数据,NAs),那么你只是在这一年的后期(早期)开始(结束)了,因此在系列的一个或两个端点上没有完整的一年数据。


2
抱歉,如果这很明显的话,我是R语言的新手,但是我的假设是周期性序列是常规序列,在其傅里叶变换中包含某些频率,我的序列是使用Sin函数创建的,因此它具有这种特性。 - Etan A Ehsanfar
2
是的,但就计算机而言,它只是一组数字向量,没有指示这些数字观察的时间点。 - Gavin Simpson
1
回答很好且详细,但是当我尝试分析小时数据(xxx <- ts(x, freq=24))时出现了错误信息。也许有时间序列长度或每个周期观测次数的限制吗?编辑:我有缺失值,并且我已经尝试使用na.action=na.passna.action=na.exclude两种方法。 - rbatt
1
@rbatt 那个检查是在NA被移除之后进行的,因此很可能你有超过2 * 24个观测值,直到NA被移除后,一旦它们被移除,你的数据就不再满足长度阈值。看起来na.action参数实际上没有起到任何作用。 - Gavin Simpson
@GavinSimpson 我有大约119个时间段。我刚刚尝试将其子集化为前6个时间段---这些值中没有一个是NA,而且它也正常工作。然后我对整个时间序列进行了线性插值,也成功了。看起来stl不能处理NA值!感谢您的想法。 - rbatt
显示剩余3条评论

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接