首先,我会利用您的时间序列的差异来消除它,它具有相当大的随机趋势。除了时间序列稳步上升之外,最明显的迹象是 ACF 成分需要很长时间才能消亡。虽然您可以使用漂移项将模型拟合到数据中,但当趋势不存在时,更容易阅读 ACF 和 PACF 图。
tf <- structure(c(58082, 48500, 45723, 53662, 46723, 45070, 49782, 55437,
57672, 61121, 43857, 49819, 50750, 53589, 53812, 53575, 52339, 51115,
56529, 61498, 58757, 72876, 55999, 58374, 63885, 63287, 60027, 65795,
62850, 61908, 68108, 72639, 77105, 84996, 65488, 62178, 74750, 66085,
59711, 69304, 68357, 67133, 74545, 73623, 82071, 89533, 72117, 69004,
72681, 80214, 80744, 81643, 87599, 86213, 97495, 97841, 104953, 110353,
90415, 83875, 93160, 89539, 85021, 91314, 87036, 83731, 91047, 94552,
105628, 94743, 84954, 72535, 77898, 68418, 60609, 73703, 67298, 64375,
73550, 76887, 77538, 92233, 73267, 77779, 80634, 72736, 81475, 87595,
87386, 88874, 95145, 96991, 95186, 106122, 81173, 77941, 88576, 86372,
77850, 91188, 90547, 87803, 95264, 90054, 100544, 96302, 82402, 78297,
91847, 86235, 87557, 91139, 93116, 93855, 94172, 100003, 97051, 86785,
84849, 81682, 88273, 85645, 80121, 92187, 96409, 97609, 94971, 111356,
102049, 110838, 97596, 88747, 100882, 97801, 99312, 100163, 112241, 101667,
122227, 127548, 123216, 131987, 112248, 118140, 128127, 114529, 151671,
135476, 148513, 141155, 142314, 142144, 139774, 142715, 124773, 111401,
129554, 140624, 128378, 130208, 141051, 132299, 145779, 152341, 146552,
150930, 139732, 133423, 154363, 148374, 137392, 149258, 160086, 154738,
159570, 164496, 166885, 188369, 144181, 148121, 169758, 158890, 159699,
161691, 165828, 175617, 181875, 182883), .Tsp = c(1, 188, 1), class = "ts")
par(mfrow=c(3, 1), mar=c(3, 3, 2, 1), mgp=c(2, 0.6, 0), oma=c(0, 0, 0, 0))
plot(tf); acf(tf, main=""); pacf(tf, main="")
图1
tf.d <- diff(tf)
plot(tf.d); acf(tf.d, main=""); pacf(tf.d, main="")
图2
现在我们可以开始解释绘图了,从 ACF 开始。我们可以看到滞后 1 处的相关性是负的,这经常发生在我们取一个序列的差异时,被称为“过度差异”。从具有与过去数据点过度相似的数据点的系列,现在它们与过去数据点过度不同,增加了周期 1 的波动。从图中我们可以看到水平大约在 -0.3,而可能不足以引起问题,值得牢记。传统智慧还告诉我们,在 1 滞后后像这样死亡的 ACF 应该配备 MA(1) 系数,我的经验是最好从 AR 系数开始,因为它们可能会使 MA 系数变得多余。但情况也可能相反,因此 MA(1) 应该在候选模型中。
查看 PACF 图,我们在滞后 1 和 2 处看到两个重要分量,这表明是 AR(2) 模型。
进一步观察 ACF 和 PACF 图,我们可以看到类似波浪的特征和滞后 12 处的正峰值,这告诉我这是每月或每两小时的数据,并且我们有季节性成分。
找出季节性成分与找出非季节性成分并没有太大区别,我们将仅根据季节性/周期性滞后而不是样本滞后来描述相关性。
tf.d12 <- ts(tf.d, f=12)
plot(tf.d12); acf(tf.d12, main="", lag.max=12*4); pacf(tf.d12, main="", lag.max=12*4)
图3
查看 ACF 图,我们看到滞后 1、2、3、4... 周期的强分量,有点让人想起第一个 ACF 图,这意味着我们有一个一阶季节性差分分量,但不是再次转换数据,我们将在现在已成为 SARIMA 模型的情况下将 D 设置为 1。如开头所述,非平稳性会掩盖 ACF 图中 MA 过程的迹象,因此我们必须等待,看看是否需要。
在 PACF 图中,我们可以在 1、2 甚至 3 个周期滞后处确定重要分量,但本着简约的精神,我们假设 SAR(2) 就足够了。
下一步是拟合所有模型并评估残差的 ACF 和 PACF。如果我们在选择候选模型方面很聪明,那么这里不应该太多。对数据的性质一无所知,这三个人或多或少是随意选择的。
ari1 <- arima(tf.d12, order=c(2, 0, 0), seasonal=c(2, 1, 0))
ari2 <- arima(tf.d12, order=c(1, 0, 1), seasonal=c(2, 1, 0))
ari3 <- arima(tf.d12, order=c(1, 0, 1), seasonal=c(2, 1, 1))
par(mfcol=c(3, 2), mar=c(3, 3, 2, 1), mgp=c(2, 0.6, 0), oma=c(0, 0, 1.5, 0))
acf(residuals(ari1), main="", ylim=c(-0.2, 1), lag.max=12*4)
mtext(ari1$call, 3, cex=0.8)
acf(residuals(ari2), main="", ylim=c(-0.2, 1), lag.max=12*4)
mtext(ari2$call, 3, cex=0.8)
acf(residuals(ari3), main="", ylim=c(-0.2, 1), lag.max=12*4)
mtext(ari3$call, 3, cex=0.8)
pacf(residuals(ari1), main="", ylim=c(-0.2, 0.44), lag.max=12*4)
mtext(ari1$call, 3, cex=0.8)
pacf(residuals(ari2), main="", ylim=c(-0.2, 0.44), lag.max=12*4)
mtext(ari2$call, 3, cex=0.8)
pacf(residuals(ari3), main="", ylim=c(-0.2, 0.44), lag.max=12*4)
mtext(ari3$call, 3, cex=0.8)
mtext("residual values", 3, outer=TRUE, cex=1.3)
图4
此类图表、关于数据的经验和理论知识、信息标准(如 AICc)的应用和验证(如 CV),将引导您找到合适的模型。盲目相信auto.arima() 是不行的。
一些进一步的说明:
par(mfrow=c(1, 1), mar=c(3, 3, 2, 1), mgp=c(2, 0.6, 0), oma=c(0, 0, 0, 0))
plot(stl(tf.d12, "periodic"))
图5
如果我们分解tf.d12,我们会看到数据中存在轻微趋势。向模型添加非季节性差异可能是合适的:
arima(tf.d12, order=c(2, 1, 1), seasonal=c(2, 1, 1))
分解还揭示了其余部分看起来像是一个时间分量,数量级似乎随着时间的推移而增加。我们的模型不处理这个问题。
tsoutliers::locate.outliers 表示在索引 145 处有一个加性异常值,在 71 和 149 处有几个时间变化。
抱歉,这有点啰嗦,我开始研究数据,无法停止。最后这整件事可能更适合Cross Validated,那里也有很多更有知识的人可以提供第二个意见。