forecast x13 special #151
Comments
|
Hi, the forecasts do work internally since library(feasts)
#> Loading required package: fabletools
fit_36 <- AirPassengers %>%
tsibble::as_tsibble() %>%
model(
seats = X_13ARIMA_SEATS(value ~ transform(`function` = "log") + forecast(maxlead = 36,save = "fct"))
)Forecasts exist in model object: fit_36$seats[[1]]$fit$fit$series$fct
#> forecast lowerci upperci
#> Jan 1961 444.2964 418.1670 472.0585
#> Feb 1961 413.5093 381.4051 448.3158
#> Mar 1961 465.5498 422.4688 513.0240
#> Apr 1961 497.5508 445.3449 555.8766
#> May 1961 498.3558 440.6152 563.6630
#> Jun 1961 575.2956 503.0633 657.8994
#> Jul 1961 670.7166 580.5151 774.9337
#> Aug 1961 654.9312 561.3962 764.0501
#> Sep 1961 556.7305 472.9611 655.3369
#> Oct 1961 490.9876 413.5452 582.9322
#> Nov 1961 422.5169 352.9779 505.7555
#> Dec 1961 478.1585 396.3657 576.8297
#> Jan 1962 481.0627 391.6394 590.9040
#> Feb 1962 452.3742 362.7822 564.0917
#> Mar 1962 514.5911 406.8467 650.8694
#> Apr 1962 538.7240 420.2625 690.5767
#> May 1962 545.1952 419.9158 707.8509
#> Jun 1962 635.8976 483.9242 835.5973
#> Jul 1962 726.2196 546.2771 965.4347
#> Aug 1962 716.4868 532.9609 963.2101
#> Sep 1962 615.3769 452.8749 836.1884
#> Oct 1962 531.6177 387.1698 729.9571
#> Nov 1962 462.2284 333.2868 641.0547
#> Dec 1962 523.0996 373.5303 732.5597
#> Jan 1963 526.2768 369.1401 750.3038
#> Feb 1963 494.8919 341.6832 716.7984
#> Mar 1963 568.7985 386.7674 836.5021
#> Apr 1963 583.3043 390.8428 870.5389
#> May 1963 596.4370 394.0242 902.8305
#> Jun 1963 702.8835 458.0716 1078.5327
#> Jul 1963 786.3155 505.6679 1222.7235
#> Aug 1963 791.9620 502.8166 1247.3808
#> Sep 1963 666.3004 417.7790 1062.6581
#> Oct 1963 581.5833 360.2171 938.9870
#> Nov 1963 510.9198 312.7342 834.6995
#> Dec 1963 566.3870 342.6866 936.1157
As far as I’m aware, the forecasts from X-13ARIMA-SEATS are from the regARIMA model it uses, so be mindful that the forecasts won’t use all options specified for the decomposition. Note that even though x11() is specified for seasonal adjustment, the forecasts do not change. fit_12 <- AirPassengers %>%
tsibble::as_tsibble() %>%
model(
seats = X_13ARIMA_SEATS(value ~ transform(`function` = "log") + x11() + forecast(maxlead = 12,save = "fct"))
)
waldo::compare(
fit_12$seats[[1]]$fit$fit$series$fct[1:12,],
fit_36$seats[[1]]$fit$fit$series$fct[1:12,]
)
#> ✔ No differencesCreated on 2022-08-25 by the reprex package (v2.0.1) I would like to add a |
I've been struggling to figure out how to incorporate forecasts from X_13ARIMA_SEATS into a tsibble/feasts workflow, e.g. to use the feasts graphics methods on a forecast from X_13ARIMA_SEATS. In leiu of a more integrated approach (i'd love the I'm thinking some pivoting and rbinding of the original data and the forecast data back into a new tsibble might do the trick, but is there something more elegant than that? Thanks!! |
|
I'm working on writing a As for using X_13ARIMA_SEATS forecasts in tsibble/feasts, I guess it depends what you're trying to do. If you're trying to plot the forecasts, you'll want to convert the mts object shown above (forecast, lowerci, upperci) into a fable with a distribution column for the forecasted values. Then you can plot these forecasts with |
|
I've pushed my work in progress code to the Here's how it works, and how the intervals currently slightly differ for library(feasts)
#> Loading required package: fabletools
library(dplyr)
#>
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#>
#> filter, lag
#> The following objects are masked from 'package:base':
#>
#> intersect, setdiff, setequal, union
fit_36 <- AirPassengers %>%
tsibble::as_tsibble() %>%
model(
seats = X_13ARIMA_SEATS(value ~ forecast(maxlead = 36,lognormal = "no", save = c("ftr", "fct", "fvr")))
)
fit_36 %>%
forecast() %>%
mutate(hilo(value, 95), median(value))
#> # A fable: 24 x 6 [1M]
#> # Key: .model [1]
#> .model index value .mean `hilo(value, 95)` `median(value)`
#> <chr> <mth> <dist> <dbl> <hilo> <dbl>
#> 1 seats 1961 Jan lN(6.1, 0.00096) 445. [418.1670, 472.0585]95 444.
#> 2 seats 1961 Feb lN(6, 0.0017) 418. [384.8105, 452.3187]95 417.
#> 3 seats 1961 Mar lN(6.1, 0.0025) 466. [422.4688, 513.0240]95 466.
#> 4 seats 1961 Apr lN(6.2, 0.0032) 498. [445.3449, 555.8766]95 498.
#> 5 seats 1961 May lN(6.2, 0.0039) 499. [440.6152, 563.6630]95 498.
#> 6 seats 1961 Jun lN(6.4, 0.0047) 577. [503.0633, 657.8994]95 575.
#> 7 seats 1961 Jul lN(6.5, 0.0054) 673. [580.5151, 774.9337]95 671.
#> 8 seats 1961 Aug lN(6.5, 0.0062) 657. [561.3962, 764.0501]95 655.
#> 9 seats 1961 Sep lN(6.3, 0.0069) 559. [472.9611, 655.3369]95 557.
#> 10 seats 1961 Oct lN(6.2, 0.0077) 493. [413.5452, 582.9322]95 491.
#> # … with 14 more rows
fit_36$seats[[1]]$fit$fit$series$fct[1:12,]
#> forecast lowerci upperci
#> [1,] 444.2964 418.1670 472.0585
#> [2,] 413.5093 381.4051 448.3158
#> [3,] 465.5498 422.4688 513.0240
#> [4,] 497.5508 445.3449 555.8766
#> [5,] 498.3558 440.6152 563.6630
#> [6,] 575.2956 503.0633 657.8994
#> [7,] 670.7166 580.5151 774.9337
#> [8,] 654.9312 561.3962 764.0501
#> [9,] 556.7305 472.9611 655.3369
#> [10,] 490.9876 413.5452 582.9322
#> [11,] 422.5169 352.9779 505.7555
#> [12,] 478.1585 396.3657 576.8297Created on 2022-08-27 by the reprex package (v2.0.1) |
|
1st & foremost 🤓: this is on my radar and appreciate the response. plan to go through the examples in more depth and provide back more code. I think the end goal - in my mind is how do we swap between nsa/sa models & forecasts e.g. components() %>% forecast() %>% reCompose() I think the seasonality forecasts get us a portion of the way there thoughts? |
|
I'm not quite sure what you mean here, are you trying to produce forecasts of seasonal and seasonally adjusted data? If so, you would typically do: Also do you have more familiarity with X13-ARIMA-SEATS, and can explain why the forecasts differ for |
|
I can take a look on the 2nd comment! on the first, I think we want to treat it like a hierarchical forecast, where we forecast the seasonal component and the SA series separately and then recombine into the NSA forecast and historical contribution. does that help? |
|
For a hierarchical reconciliation of decomposed data, you can set this up by manipulating the components data: library(fable)
#> Loading required package: fabletools
library(feasts)
library(tidyr)
library(dplyr)
#>
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#>
#> filter, lag
#> The following objects are masked from 'package:base':
#>
#> intersect, setdiff, setequal, union
# Prepare/tidy dataset
as_tsibble(USAccDeaths) %>%
# Specify and estimate model for decomposition
model(STL(value)) %>%
# Extract decomposed data
components() %>%
# Convert components of interest into long form
select(index, season_year, season_adjust) %>%
pivot_longer(c(season_year, season_adjust), names_to = "component") %>%
# Compute aggregate (the original data) from decomposition
aggregate_key(component, value = sum(value)) %>%
# Model each series
model(ets = ETS(value)) %>%
# Impose reconciliation constraint
reconcile(ets_dcmp_rec = min_trace(ets)) %>%
# Produce forecasts
forecast(h = "2 years") %>%
# Optionally, focus only on original series
filter(is_aggregated(component)) %>%
# Plot and compare forecasts
# (seems reconciled forecasts have smaller intervals but similar means)
autoplot()Created on 2022-08-31 by the reprex package (v2.0.1) I'd be very interested to hear about what you think of the |
|
@mitchelloharawild I connected with one of my colleagues and we believe this is due to leap-year adjustments in the x13 model. This would line up with the forecasts differing in February / |
|
If I can intervene in the conversation, this certainly come from the pre-leap-year adjustment. Indeed, by default, when a log-transformation is used the series is pre-adjusted from leap-year instead of adding an external regressor in the regarima model. In that case the month of February is renormalise to corresponds to a month of length 28.5: the values are multiplied by 28.25/28 in leap years and 28.25/29 otherwise. This correction is not performed with additive models (and as you see no difference in that case for library(feasts)
library(dplyr)
fit_36 <- AirPassengers %>%
tsibble::as_tsibble() %>%
model(
seats = X_13ARIMA_SEATS(value ~ forecast(maxlead = 36,lognormal = "no", save = c("ftr", "fct", "fvr")))
)
fit_36$seats[[1]]$fit$fit # no leap year regressor
#>
#> Call:
#> NULL
#>
#> Coefficients:
#> Weekday Easter[1] AO1951.May MA-Nonseasonal-01
#> -0.00295 0.01777 0.10016 0.11562
#> MA-Seasonal-12
#> 0.49736
fit_36 %>%
forecast() %>%
mutate(hilo(value, 95), median(value))
#> # A fable: 24 x 6 [1M]
#> # Key: .model [1]
#> .model index value .mean `hilo(value, 95)` median(valu…¹
#> <chr> <mth> <dist> <dbl> <hilo> <dbl>
#> 1 seats 1961 janv. lN(6.1, 0.00096) 445. [418.1670, 472.0585]95 444.
#> 2 seats 1961 févr. lN(6, 0.0017) 418. [384.8105, 452.3187]95 417.
#> 3 seats 1961 mars lN(6.1, 0.0025) 466. [422.4688, 513.0240]95 466.
#> 4 seats 1961 avr. lN(6.2, 0.0032) 498. [445.3449, 555.8766]95 498.
#> 5 seats 1961 mai lN(6.2, 0.0039) 499. [440.6152, 563.6630]95 498.
#> 6 seats 1961 juin lN(6.4, 0.0047) 577. [503.0633, 657.8994]95 575.
#> 7 seats 1961 juil. lN(6.5, 0.0054) 673. [580.5151, 774.9337]95 671.
#> 8 seats 1961 août lN(6.5, 0.0062) 657. [561.3962, 764.0501]95 655.
#> 9 seats 1961 sept. lN(6.3, 0.0069) 559. [472.9611, 655.3369]95 557.
#> 10 seats 1961 oct. lN(6.2, 0.0077) 493. [413.5452, 582.9322]95 491.
#> # … with 14 more rows, and abbreviated variable name ¹`median(value)`
# 1961 is not a leap year
diag(c(1, 28.25/28, rep(1,10))) %*%
fit_36$seats[[1]]$fit$fit$series$fct[1:12,] # same value
#> forecast lowerci upperci
#> [1,] 444.2964 418.1670 472.0585
#> [2,] 417.2013 384.8105 452.3187
#> [3,] 465.5498 422.4688 513.0240
#> [4,] 497.5508 445.3449 555.8766
#> [5,] 498.3558 440.6152 563.6630
#> [6,] 575.2956 503.0633 657.8994
#> [7,] 670.7166 580.5151 774.9337
#> [8,] 654.9312 561.3962 764.0501
#> [9,] 556.7305 472.9611 655.3369
#> [10,] 490.9876 413.5452 582.9322
#> [11,] 422.5169 352.9779 505.7555
#> [12,] 478.1585 396.3657 576.8297
fit_36_add <- AirPassengers %>%
tsibble::as_tsibble() %>%
model(
seats = X_13ARIMA_SEATS(value ~ forecast(maxlead = 36,lognormal = "no", save = c("ftr", "fct", "fvr")) + transform(`function` = "none"))
)
fit_36_add$seats[[1]]$fit$fit # leap year regressor
#>
#> Call:
#> NULL
#>
#> Coefficients:
#> Constant Leap Year Weekday Easter[1]
#> 30.6208 11.3210 -0.9036 6.8937
#> AR-Nonseasonal-01
#> 0.8193
fit_36_add %>%
forecast() %>%
mutate(hilo(value, 95), median(value))
#> # A fable: 24 x 6 [1M]
#> # Key: .model [1]
#> .model index value .mean `hilo(value, 95)` `median(value)`
#> <chr> <mth> <dist> <dbl> <hilo> <dbl>
#> 1 seats 1961 janv. N(439, 100) 439. [419.3149, 458.4862]95 439.
#> 2 seats 1961 févr. N(407, 178) 407. [380.5100, 432.7901]95 407.
#> 3 seats 1961 mars N(446, 212) 446. [417.3459, 474.4348]95 446.
#> 4 seats 1961 avr. N(492, 244) 492. [461.1230, 522.3311]95 492.
#> 5 seats 1961 mai N(497, 266) 497. [464.9678, 528.9404]95 497.
#> 6 seats 1961 juin N(564, 281) 564. [530.7221, 596.4164]95 564.
#> 7 seats 1961 juil. N(651, 291) 651. [617.4773, 684.4027]95 651.
#> 8 seats 1961 août N(635, 299) 635. [601.3547, 669.1328]95 635.
#> 9 seats 1961 sept. N(541, 305) 541. [506.4526, 574.8578]95 541.
#> 10 seats 1961 oct. N(488, 309) 488. [453.1060, 521.9616]95 488.
#> # … with 14 more rows
fit_36_add$seats[[1]]$fit$fit$series$fct[1:12,]
#> forecast lowerci upperci
#> [1,] 438.9006 419.3149 458.4862
#> [2,] 406.6501 380.5100 432.7901
#> [3,] 445.8903 417.3459 474.4348
#> [4,] 491.7271 461.1230 522.3311
#> [5,] 496.9541 464.9678 528.9404
#> [6,] 563.5693 530.7221 596.4164
#> [7,] 650.9400 617.4773 684.4027
#> [8,] 635.2437 601.3547 669.1328
#> [9,] 540.6552 506.4526 574.8578
#> [10,] 487.5338 453.1060 521.9616
#> [11,] 419.8635 385.3110 454.4159
#> [12,] 465.1630 430.4826 499.8433
Note that the forecasts have an impact in the decomposition, especially with the X-11 decomposition: the decomposes series is a series extended with forecasts.
This is normal because the forecasts only depends on the pre-adjustment process and SEATS/X-11 is the decomposition algorithm. |
Is the intended behavior of the
forecastspecial within the x13 model to match that of the seasonal package?By way of example, with the
forecast.saveandforecast.maxleadarguments the return behavior is different than the seasonal example (could be user error).Is it intended that this would be performed within a downstream model and forecast or is there a method for extracting these from the feasts special?
I appreciate the help and guidance in advance.
Created on 2022-08-06 by the reprex package (v2.0.1)
The text was updated successfully, but these errors were encountered: