|
90 | 90 | gasoline %>% stlf() %>% autoplot()
|
91 | 91 |
|
92 | 92 | # An alternative approach is to use a dynamic harmonic regression model. Example using gasoline data again.
|
93 |
| -gasoline_dreg.best <- list(aicc=Inf) |
94 |
| -
|
95 |
| -for(K in seq(26)){ |
96 |
| - gasoline_dreg <- auto.arima( |
97 |
| - #substitute seasonal component with Fourier terms. |
| 93 | +# Use parallel computing to choose the number of Fourier pairs that yields smallest AIC value. |
| 94 | +gasoline_aiccs <- foreach( |
| 95 | + i = 1:26, |
| 96 | + .packages = 'fpp2' |
| 97 | + ) %dopar% |
| 98 | + auto.arima( |
98 | 99 | gasoline,
|
99 |
| - xreg=fourier(gasoline, K=K), |
| 100 | + xreg=fourier(gasoline, K=i), |
100 | 101 | seasonal=FALSE
|
101 |
| - ) |
102 |
| - |
103 |
| - if(gasoline_dreg$aicc < gasoline_dreg.best$aicc) |
104 |
| - { |
105 |
| - gasoline_dreg.best <- gasoline_dreg |
106 |
| - gasoline_K.best <- K |
107 |
| - } |
108 |
| -} |
| 102 | + )$aic %>% |
| 103 | + unlist() |
| 104 | +
|
| 105 | +gasoline_K.best <- which( |
| 106 | + gasoline_aiccs == min(gasoline_aiccs) |
| 107 | +) |
| 108 | +
|
| 109 | +# Using for-loop to get the number of Fourier pairs. |
| 110 | +#gasoline_dreg.best <- list(aicc=Inf) |
| 111 | +#for(K in seq(26)){ |
| 112 | +# gasoline_dreg <- auto.arima( |
| 113 | +# #substitute seasonal component with Fourier terms. |
| 114 | +# gasoline, |
| 115 | +# xreg=fourier(gasoline, K=K), |
| 116 | +# seasonal=FALSE |
| 117 | +# ) |
| 118 | +# |
| 119 | +# if(gasoline_dreg$aicc < gasoline_dreg.best$aicc) |
| 120 | +# { |
| 121 | +# gasoline_dreg.best <- gasoline_dreg |
| 122 | +# gasoline_K.best <- K |
| 123 | +# } |
| 124 | +#} |
109 | 125 |
|
110 | 126 | # forecast the next 2 years of data. If I assume that 1 year is about 52 weeks, forecast horizon(h) is 104.
|
111 | 127 | fc_gasoline_dreg.best <- forecast(
|
@@ -211,12 +227,23 @@ fc_gas_ets <- forecast(gas_ets, h=6)
|
211 | 227 | # Simulate 10000 future sample paths
|
212 | 228 | nsim <- 10000
|
213 | 229 | h <- 6
|
214 |
| -sim <- numeric(nsim) |
215 | 230 |
|
216 |
| -for(i in seq_len(nsim)){ |
217 |
| - # for each sample path, add 6 months' forecasts. |
218 |
| - sim[i] <- sum(simulate(gas_ets, future=TRUE, nsim=h)) |
219 |
| -} |
| 231 | +# Use parallel computing to get sample paths and the sums of their forecasts. |
| 232 | +sim <- foreach( |
| 233 | + i = 1:nsim, |
| 234 | + .packages = 'forecast' |
| 235 | + ) %dopar% |
| 236 | + sum( |
| 237 | + simulate(gas_ets, future = TRUE, nsim = h) |
| 238 | + ) %>% |
| 239 | + unlist() |
| 240 | + |
| 241 | +# Use for-loop to get sample paths and the sums of their forecasts. |
| 242 | +#sim <- numeric(nsim) |
| 243 | +#for(i in seq_len(nsim)){ |
| 244 | +# # for each sample path, add 6 months' forecasts. |
| 245 | +# sim[i] <- sum(simulate(gas_ets, future=TRUE, nsim=h)) |
| 246 | +#} |
220 | 247 |
|
221 | 248 | # get final aggregated forecast.
|
222 | 249 | gas_ets.sim.meanagg <- mean(sim)
|
|
0 commit comments