R拟合和预测每日时间序列

Question

我正在处理每日时间序列，我需要根据我的历史记录构建 90 天(或更多)的预测 - 当前时间序列大约有 298 个数据点。

我遇到的问题是最终预测中著名的平线 - 是的，我可能没有季节性，但我正在努力解决这个问题。另一个问题是如何找到最佳模型并从现在开始针对这种行为进行调整。

我创建了一个测试用例来进一步调查此问题，我们将不胜感激。

谢谢，

开始

x <- day_data  # My time serie
z <- 90        # Days to forecast

low_bound_date <- as.POSIXlt(min(x$time), format = "%m/%d/%Y") # oldest date in the DF.

> low_bound_date
[1] "2015-12-21 PST"

low_bound_date$yday 
> low_bound_date$yday  # Day in Julian
[1] 354

lbyear <- as.numeric(substr(low_bound_date, 1, 4))
> lbyear
[1] 2015

这是我的时间系列内容

> ts
Time Series:
Start = c(2065, 4) 
End = c(2107, 7) 
Frequency = 7 
  [2] 20.73 26.19 27.51 26.11 26.28 27.58 26.84 27.00 26.30 28.75 28.43 39.03 41.36 45.42 44.80 45.33 47.79 44.70 45.17
 [20] 34.90 32.54 32.75 33.35 34.76 34.11 33.59 33.60 38.08 30.45 29.66 31.09 31.36 31.96 29.30 30.04 30.85 31.13 25.09
 [39] 17.88 23.73 25.31 31.30 35.18 34.13 34.96 35.12 27.36 38.33 38.59 38.14 38.54 41.72 37.15 35.92 37.37 32.39 30.64
 [58] 30.57 30.66 31.16 31.50 30.68 32.21 32.27 32.55 33.61 34.80 33.53 33.09 20.90  6.91  7.82 15.78  7.25  6.19  6.38
 [77] 38.06 39.82 35.53 38.63 41.91 39.76 37.26 38.79 37.74 35.61 39.70 35.79 35.36 29.63 22.07 35.39 35.99 37.35 38.82
 [96] 25.80 21.31 18.85  9.52 20.75 36.83 44.12 37.79 34.45 36.05 16.39 21.84 31.39 34.26 31.50 30.87 28.88 42.83 41.52
[115] 42.34 47.35 44.47 44.10 44.49 26.89 18.17 40.44 43.93 41.56 39.98 40.31 40.59 40.17 40.22 40.50 32.68 35.89 36.06
[134] 34.30 22.67 12.56 13.29 12.34 28.00 35.27 36.57 33.78 32.15 33.58 34.62 30.96 32.06 33.05 30.66 32.47 30.42 32.83
[153] 31.74 29.39 22.39 12.58 16.46  5.36  4.01 15.32 32.79 31.66 32.02 27.60 31.47 31.61 34.96 27.77 31.91 33.94 33.43
[172] 26.94 28.38 21.42 24.51 23.82 31.71 26.64 27.96 29.29 29.25 28.70 27.02 27.62 30.90 27.46 27.37 26.46 27.77 13.61
[191]  5.87 12.18  5.68  4.15  4.35  4.42 16.42 25.18 26.06 27.39 27.57 28.86 15.18  5.19  5.61  8.28  7.78  5.13  4.90
[210]  5.02  5.27 16.31 25.01 26.19 25.96 24.93 25.53 25.56 26.39 26.80 26.73 26.00 25.61 25.90 25.89 13.80  6.66  6.41
[229]  5.28  5.64  5.71  5.38  5.76  7.20  7.27  5.55  5.31  5.94  5.75  5.93  5.77  6.57  5.52  5.51  5.47  5.69 19.75
[248] 29.22 30.75 29.63 30.49 29.48 31.83 30.42 29.27 30.40 29.91 32.00 30.09 28.93 14.54  7.75  5.63 17.17 22.27 24.93
[267] 35.94 37.42 33.13 25.88 24.27 37.64 37.42 38.33 35.20 21.32  7.32  4.81  5.17 17.49 23.77 23.36 27.60 26.53 24.99
[286] 24.22 23.76 24.10 24.22 27.06 25.53 23.40 37.07 26.52 25.19 28.02 28.53 26.67

第一步，我在 ts 中获取我的数据

day_data_ts <- ts(x$avg_day, start = c(lbyear,low_bound_date$yday), frequency=7)

plot(day_data_ts)

plot_ts

acf(day_data_ts)

acf_ts

第二步，我在 msts 中获取我的数据

day_data_msts <- msts(x$avg_day, seasonal.periods=c(7,365.25), start = c(lbyear,low_bound_date$yday))

plot(day_data_msts)

acf(day_data_msts)

我进行了几次拟合迭代，试图找出最佳拟合和预测模型。

第一个拟合测试仅使用 ts。

fit1 <- HoltWinters(day_data_ts)
> fit1
    Holt-Winters exponential smoothing with trend and additive seasonal component.
    Call: HoltWinters(x = day_data_ts)
    Smoothing parameters: alpha: 1   beta : 0.006757112  gamma: 0

    Coefficients:
             [,1]
    a  28.0922449
    b   0.1652477
    s1  0.6241837
    s2  1.9084694
    s3  0.9913265
    s4  0.8198980
    s5 -1.7015306
    s6 -1.2201020
    s7 -1.4222449


fit2 <- tbats(day_data_ts)
> fit2
    BATS(1, {0,0}, 0.8, -)
    Parameters:   Alpha: 1.309966     Beta: -0.3011143    Damping Parameter: 0.800001
    Seed States:
              [,1]
    [1,] 15.282259
    [2,]  2.177787
    Sigma: 5.501356     AIC: 2723.911


fit3 <- ets(day_data_ts)
> fit3
    ETS(A,N,N) 
      Smoothing parameters: alpha = 0.9999 
      Initial states:       l = 25.2275 
      sigma:  5.8506
         AIC     AICc      BIC 
    2756.597 2756.678 2767.688 


fit4 <- auto.arima(day_data_ts)
> fit4
    ARIMA(1,1,2)                    
    Coefficients:
             ar1      ma1      ma2
          0.7396  -0.6897  -0.2769
    s.e.  0.0545   0.0690   0.0621
    sigma^2 estimated as 30.47:  log likelihood=-927.9
    AIC=1863.81   AICc=1863.94   BIC=1878.58

第二个测试使用 msts。我还将 ets 模型更改为 MAM。

fit5 <- tbats(day_data_msts)
> fit5
    BATS(1, {0,0}, 0.8, -)
    Parameters:   Alpha: 1.309966     Beta: -0.3011143    Damping Parameter: 0.800001
    Seed States:
              [,1]
    [1,] 15.282259
    [2,]  2.177787
    Sigma: 5.501356     AIC: 2723.911


fit6 <- ets(day_data_msts, model="MAN")
> fit6
    ETS(M,A,N) 
      Smoothing parameters:     alpha = 0.9999      beta  = 9e-04 
      Initial states:           l = 52.8658         b = 3.9184 
      sigma:  0.3459
         AIC     AICc      BIC 
    3042.744 3042.949 3061.229 


fit7 <- auto.arima(day_data_msts)
> fit7
    ARIMA(1,1,2)                    
    Coefficients:
             ar1      ma1      ma2
          0.7396  -0.6897  -0.2769
    s.e.  0.0545   0.0690   0.0621
    sigma^2 estimated as 30.47:  log likelihood=-927.9
    AIC=1863.81   AICc=1863.94   BIC=1878.58

Answer 1

您可以按如下方式对先前估计的模型进行预测(使用内置时间序列 LakeHuron):

library(forecast)
y <- LakeHuron
tsdisplay(y)
# estimate ARMA(1,1)
mod_2 <- Arima(y, order = c(1, 0, 1))
#make forecast for 5 periods (years in this case)
fHuron <- forecast(mod_2, h = 5)
#show results in table
fHuron
#plot results
plot(fHuron)

这会给你: 请注意，ARIMA 模型的预测基于先前的值，因此如果我们对多个周期进行预测，该模型将使用已经预测的值来预测下一个。这会降低准确性。

要拟合最佳 ARIMA 模型，请使用此函数:

library(R.utils) #for the function 'withTimeout'
fitARIMA<-function(timeseriesObject, timout)
{
    final.aic <- Inf
    final.order <- c(0,0,0)
    for (p in 0:5) for (q in 0:5) {
        if ( p == 0 && q == 0) {
        next
        }

        arimaFit = tryCatch( 
        withTimeout(arima(timeseriesObject
                            ,order=c(p, 0, q))
                    ,timeout = timeout)
        ,error=function( err ) FALSE
        ,warning=function( err ) FALSE )

        if( !is.logical( arimaFit ) ) {
        current.aic <- AIC(arimaFit)
        if (current.aic < final.aic) {
            final.aic <- current.aic
            final.order <- c(p, 0, q)
            final.arima <- arima(timeseriesObject, order=final.order)
        }
        } else {
        next
        }
    }
    final.order<-c(final.order,final.aic)
    final.order
}