R 的 decompose/stl 函数无法从每日时间序列中完全提取年度季节性
R's decompose/stl functions unable to completely extract annual seasonality from daily time series
上下文优先，问题在底部。
我有 10 年的每日降水数据，这些数据显示出年度季节性，我正在尝试使用 ARMA 方法对其进行建模，然后进行预测。数据here，时间序列对象创建如下。
我知道常见的 R 包和函数很难处理每日时间序列。例如，Forecast 的 arima() 函数不接受高于 350 的频率，而 ts() 不接受频率的非整数值（这两者都会有用，因为一年的平均天数是 365.25）。
显然，Forecast 的 msts() 函数可以接受 seasonal.perdiods 参数的非整数值，因此我创建了这样的时间序列：
the_ts <- msts(data$PRCP, start=c(2007, 10), end=c(2017, 9), seasonal.periods=c(365.25))
Time Series:
Start = c(2007, 10)
End = c(2017, 9)
Frequency = 365
[1] 0.09 0.75 1.63 0.06 0.36 0.63 0.76 0.43 0.13 0.00 0.00 0.02 0.31 1.80 0.03 0.19 0.25 0.01 0.00 0.52 0.01 0.00 0.00 0.00 0.00 ... etc
--
plot(the_ts)
这个级数是平稳的，所以不需要差分。
然后我希望分解序列，提取季节性和趋势，留下残差，如果成功，应该近似于白噪声。
下面是分解图。我也尝试过使用 decompose() 并进行了大量的参数调整。
the_ts_decomp = stl(the_ts, s.window = "periodic")
plot(the_ts_decomp)
显然，残差中遗留了某种类型的季节性，因为数据中季节性趋势的形状在残差中是平行的。让我们删除已识别的季节性成分并检查假定的去季节性数据：
the_ts_deseasonal <- seasadj(the_ts_decomp)
plot(the_ts_deseasonal)
对我来说仍然很季节性。 ACF 和 PACF（未显示）确认正在发生自相关。
Acf(the_ts_deseasonal, lag.max=1000)
Pacf(the_ts_deseasonal, lag.max=1000)
显然，可以通过将傅立叶变换作为外生协变量传递到 ARMA 模型中来在预测时考虑序列的季节性成分，如预测包 here and here 的作者所述，但我不清楚这与在拟合模型之前分解和去季节化序列的需求究竟有何关系。
我能够根据上述看似分解不充分的数据使用该方法生成以下预测。预测看起来不像是家 运行，残差证实有问题：
the_ts.fit <- auto.arima(the_ts, seasonal=FALSE, xreg=fourier(the_ts, K=5))
plot(forecast(the_ts.fit, h=365, xreg=fourier(the_ts, K=5, h=365)))
tsdisplay(residuals(the_ts.fit), lag.max=1000)
我不太了解 Hyndman 的 exreg=fourier 解决方案，可能是我不知道如何将其正确应用于我的情况。
问题 1：我不需要能够将数据分解为无趋势、无季节的白噪声，以便为预测重建它吗？使用 exreg=fourier 解时怎么样？
问题 2：为什么上述代码无法提取系列的季节性成分，我该如何解决？
问题 3：可以使用什么包、函数或技术来指定 365.25 的年度季节性？
我将尝试回答以下问题：
问题 1：我不需要能够将数据分解为无趋势、无季节的白噪声，以便为预测重建它吗？使用 exreg=fourier 解时怎么样？
傅里叶解考虑了确定性的季节性，即它假设季节性模式不随时间变化。一般来说，我不会说趋势和季节性调整的系列确实需要是白噪声，因为系列中可能会保留一些短期模式。
问题 2：为什么上述代码无法提取系列的季节性成分，我该如何解决？
当您在 stl 函数中指定 s.window = "periodic" 时，您基本上假设季节性模式不会随时间改变。这可能是问题的根源之一。
问题 3：可以使用什么包、函数或技术来指定 365.25 的年度季节性？
dsa包中的dsa函数，需要将时间序列转换为xts格式。例如像这样：
data = rnorm(365.25*10, 100, 1)
data_xts <- xts::xts(data, seq.Date(as.Date("2007-01-01"), by="days", length.out = length(data)))
sa = dsa::dsa(data_xts, fourier_number = 24) # the fourier_number is used to model monthly recurring seasonal patterns in the regARIMA part
data_adjusted <- sa$output[,1]
上下文优先，问题在底部。
我有 10 年的每日降水数据，这些数据显示出年度季节性，我正在尝试使用 ARMA 方法对其进行建模，然后进行预测。数据here，时间序列对象创建如下。
我知道常见的 R 包和函数很难处理每日时间序列。例如，Forecast 的 arima() 函数不接受高于 350 的频率，而 ts() 不接受频率的非整数值（这两者都会有用，因为一年的平均天数是 365.25）。
显然，Forecast 的 msts() 函数可以接受 seasonal.perdiods 参数的非整数值，因此我创建了这样的时间序列：
the_ts <- msts(data$PRCP, start=c(2007, 10), end=c(2017, 9), seasonal.periods=c(365.25))
Time Series:
Start = c(2007, 10)
End = c(2017, 9)
Frequency = 365
[1] 0.09 0.75 1.63 0.06 0.36 0.63 0.76 0.43 0.13 0.00 0.00 0.02 0.31 1.80 0.03 0.19 0.25 0.01 0.00 0.52 0.01 0.00 0.00 0.00 0.00 ... etc
--
plot(the_ts)
这个级数是平稳的，所以不需要差分。
然后我希望分解序列，提取季节性和趋势，留下残差，如果成功，应该近似于白噪声。
下面是分解图。我也尝试过使用 decompose() 并进行了大量的参数调整。
the_ts_decomp = stl(the_ts, s.window = "periodic")
plot(the_ts_decomp)
显然，残差中遗留了某种类型的季节性，因为数据中季节性趋势的形状在残差中是平行的。让我们删除已识别的季节性成分并检查假定的去季节性数据：
the_ts_deseasonal <- seasadj(the_ts_decomp)
plot(the_ts_deseasonal)
对我来说仍然很季节性。 ACF 和 PACF（未显示）确认正在发生自相关。
Acf(the_ts_deseasonal, lag.max=1000)
Pacf(the_ts_deseasonal, lag.max=1000)
显然，可以通过将傅立叶变换作为外生协变量传递到 ARMA 模型中来在预测时考虑序列的季节性成分，如预测包 here and here 的作者所述，但我不清楚这与在拟合模型之前分解和去季节化序列的需求究竟有何关系。
我能够根据上述看似分解不充分的数据使用该方法生成以下预测。预测看起来不像是家 运行，残差证实有问题：
the_ts.fit <- auto.arima(the_ts, seasonal=FALSE, xreg=fourier(the_ts, K=5))
plot(forecast(the_ts.fit, h=365, xreg=fourier(the_ts, K=5, h=365)))
tsdisplay(residuals(the_ts.fit), lag.max=1000)
我不太了解 Hyndman 的 exreg=fourier 解决方案，可能是我不知道如何将其正确应用于我的情况。
问题 1：我不需要能够将数据分解为无趋势、无季节的白噪声，以便为预测重建它吗？使用 exreg=fourier 解时怎么样？
问题 2：为什么上述代码无法提取系列的季节性成分，我该如何解决？
问题 3：可以使用什么包、函数或技术来指定 365.25 的年度季节性？
我将尝试回答以下问题：
问题 1：我不需要能够将数据分解为无趋势、无季节的白噪声，以便为预测重建它吗？使用 exreg=fourier 解时怎么样？
傅里叶解考虑了确定性的季节性，即它假设季节性模式不随时间变化。一般来说，我不会说趋势和季节性调整的系列确实需要是白噪声，因为系列中可能会保留一些短期模式。
问题 2：为什么上述代码无法提取系列的季节性成分，我该如何解决？
当您在 stl 函数中指定 s.window = "periodic" 时，您基本上假设季节性模式不会随时间改变。这可能是问题的根源之一。
问题 3：可以使用什么包、函数或技术来指定 365.25 的年度季节性？
dsa包中的dsa函数，需要将时间序列转换为xts格式。例如像这样：
data = rnorm(365.25*10, 100, 1)
data_xts <- xts::xts(data, seq.Date(as.Date("2007-01-01"), by="days", length.out = length(data)))
sa = dsa::dsa(data_xts, fourier_number = 24) # the fourier_number is used to model monthly recurring seasonal patterns in the regARIMA part
data_adjusted <- sa$output[,1]