library(tidyverse)
library(fpp3)
library(fable.prophet)Clase 7
lax_passengers <- read.csv("https://raw.githubusercontent.com/mitchelloharawild/fable.prophet/master/data-raw/lax_passengers.csv")
lax_passengers <- lax_passengers |>
mutate(datetime = mdy_hms(ReportPeriod)) |>
group_by(month = yearmonth(datetime), type = Domestic_International) |>
summarise(passengers = sum(Passenger_Count), .groups = "drop") |>
as_tsibble(index = month, key = type)lax_passengers# A tsibble: 318 x 3 [1M]
# Key: type [2]
month type passengers
<mth> <chr> <int>
1 2006 Jan Domestic 3315937
2 2006 Feb Domestic 3045963
3 2006 Mar Domestic 3711057
4 2006 Apr Domestic 3683040
5 2006 May Domestic 3743141
6 2006 Jun Domestic 3994798
7 2006 Jul Domestic 4236581
8 2006 Aug Domestic 4128492
9 2006 Sep Domestic 3386421
10 2006 Oct Domestic 3598030
# ℹ 308 more rowslax_passengers |>
autoplot(passengers) +
theme(legend.position = "top")
train/test splits
Vamos a hacer un pronóstico a dos años.
lax_pass_train <- lax_passengers |>
filter_index(. ~ "2017 Mar")
lax_pass_train# A tsibble: 270 x 3 [1M]
# Key: type [2]
month type passengers
<mth> <chr> <int>
1 2006 Jan Domestic 3315937
2 2006 Feb Domestic 3045963
3 2006 Mar Domestic 3711057
4 2006 Apr Domestic 3683040
5 2006 May Domestic 3743141
6 2006 Jun Domestic 3994798
7 2006 Jul Domestic 4236581
8 2006 Aug Domestic 4128492
9 2006 Sep Domestic 3386421
10 2006 Oct Domestic 3598030
# ℹ 260 more rowsModel fitting
lax_pass_fit <- lax_pass_train |>
model(
snaive = SNAIVE(passengers),
snaive_drift = decomposition_model(
STL(passengers, robust = TRUE),
RW(season_adjust ~ drift())
),
ets = ETS(passengers),
sarima = ARIMA(passengers),
reg_lin = TSLM(passengers ~ season() + trend()),
reg_x_partes = TSLM(passengers ~ season() + trend(knots = as.Date("2010-01-01"))),
harmonic = ARIMA(passengers ~ fourier(K = 2) + PDQ(0,0,0)),
prophet = prophet(passengers ~ growth("linear") + season("year", type = "multiplicative")),
prophet_auto = prophet(passengers)
)
lax_pass_fit# A mable: 2 x 10
# Key: type [2]
type snaive snaive_drift ets
<chr> <model> <model> <model>
1 Domestic <SNAIVE> <STL decomposition model> <ETS(M,Ad,M)>
2 International <SNAIVE> <STL decomposition model> <ETS(M,A,M)>
# ℹ 6 more variables: sarima <model>, reg_lin <model>, reg_x_partes <model>,
# harmonic <model>, prophet <model>, prophet_auto <model>accuracy(lax_pass_fit) |>
arrange(type, MAPE) |>
select(type, .model, .type, MAPE, MASE, MAE, RMSE)# A tibble: 18 × 7
type .model .type MAPE MASE MAE RMSE
<chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
1 Domestic prophet Training 1.59 0.326 61442. 74259.
2 Domestic ets Training 1.60 0.330 62046. 80544.
3 Domestic sarima Training 1.60 0.329 61892. 86413.
4 Domestic prophet_auto Training 1.65 0.341 64238. 81226.
5 Domestic snaive_drift Training 1.75 0.361 67990. 88424.
6 Domestic harmonic Training 4.08 0.850 160007. 196711.
7 Domestic snaive Training 4.72 1 188278. 221068.
8 Domestic reg_lin Training 5.28 1.07 201989. 229489.
9 Domestic reg_x_partes Training 5.28 1.07 201989. 229489.
10 International ets Training 1.71 0.269 25276. 33315.
11 International snaive_drift Training 1.80 0.292 27412. 39392.
12 International sarima Training 1.89 0.297 27904. 38573.
13 International prophet Training 2.04 0.316 29699. 38842.
14 International prophet_auto Training 2.53 0.400 37553. 50182.
15 International harmonic Training 4.27 0.667 62669. 80863.
16 International snaive Training 6.19 1 93910. 115513.
17 International reg_lin Training 7.58 1.19 111985. 126628.
18 International reg_x_partes Training 7.58 1.19 111985. 126628.Forecast on the test set
lax_pass_fc <- lax_pass_fit |>
forecast(h = "2 years")Warning: There were 2 warnings in `mutate()`.
The first warning was:
ℹ In argument: `reg_x_partes = (function (object, ...) ...`.
Caused by warning:
! prediction from a rank-deficient fit may be misleading
ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.lax_pass_fc |>
autoplot(lax_passengers |> filter_index("2014 Jan." ~ .), level = NULL)
p <- lax_pass_fc |>
ggplot(aes(x = month, y = .mean)) +
geom_line(data = lax_passengers |> filter_index("2014 Jan." ~ .), aes(y = passengers)) +
geom_line(aes(color = .model)) +
facet_wrap(~ type, scales = "free_y", ncol = 1)Warning: The output of `fortify(<fable>)` has changed to better suit usage with the ggdist package.
If you're using it to extract intervals, consider using `hilo()` to compute intervals, and `unpack_hilo()` to obtain values.plotly::ggplotly(p, dynamicTicks = TRUE)Vamos a combinar los modelos “ets”, “prophet_auto” y “sarima”
lax_pass_fit <- lax_pass_fit |>
mutate(combinado = (ets + prophet_auto + sarima)/3)
lax_pass_fit# A mable: 2 x 11
# Key: type [2]
type snaive snaive_drift ets
<chr> <model> <model> <model>
1 Domestic <SNAIVE> <STL decomposition model> <ETS(M,Ad,M)>
2 International <SNAIVE> <STL decomposition model> <ETS(M,A,M)>
# ℹ 7 more variables: sarima <model>, reg_lin <model>, reg_x_partes <model>,
# harmonic <model>, prophet <model>, prophet_auto <model>, combinado <model>Corremos los pronósticos nuevamente
lax_pass_fc <- lax_pass_fit |>
forecast(h = "2 years")Warning: There were 2 warnings in `mutate()`.
The first warning was:
ℹ In argument: `reg_x_partes = (function (object, ...) ...`.
Caused by warning:
! prediction from a rank-deficient fit may be misleading
ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.p <- lax_pass_fc |>
ggplot(aes(x = month, y = .mean)) +
geom_line(data = lax_passengers |> filter_index("2014 Jan." ~ .), aes(y = passengers)) +
geom_line(aes(color = .model)) +
facet_wrap(~ type, scales = "free_y", ncol = 1)Warning: The output of `fortify(<fable>)` has changed to better suit usage with the ggdist package.
If you're using it to extract intervals, consider using `hilo()` to compute intervals, and `unpack_hilo()` to obtain values.plotly::ggplotly(p, dynamicTicks = TRUE)lax_pass_fc |>
accuracy(lax_passengers) |>
arrange(type, MAPE) |>
select(type, .model, .type, MAPE, MASE, MAE, RMSE)# A tibble: 20 × 7
type .model .type MAPE MASE MAE RMSE
<chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
1 Domestic combinado Test 1.07 0.286 53817. 70465.
2 Domestic sarima Test 1.46 0.395 74318. 88136.
3 Domestic snaive_drift Test 1.77 0.492 92694. 111599.
4 Domestic ets Test 2.40 0.636 119747. 139337.
5 Domestic prophet_auto Test 2.73 0.707 133059. 172953.
6 Domestic prophet Test 3.52 0.965 181761. 208428.
7 Domestic snaive Test 4.45 1.20 225134. 251961.
8 Domestic harmonic Test 6.49 1.80 337965. 396988.
9 Domestic reg_lin Test 7.17 1.96 368297. 378755.
10 Domestic reg_x_partes Test 7.17 1.96 368297. 378755.
11 International prophet Test 3.07 0.681 63942. 87985.
12 International sarima Test 3.28 0.748 70206. 82727.
13 International combinado Test 3.49 0.801 75265. 87059.
14 International prophet_auto Test 4.09 0.914 85807. 112381.
15 International ets Test 4.85 1.11 104617. 115124.
16 International snaive_drift Test 7.02 1.65 154942. 172200.
17 International snaive Test 8.93 2.04 191874. 208335.
18 International harmonic Test 9.51 2.22 208069. 227662.
19 International reg_lin Test 15.4 3.54 332155. 342893.
20 International reg_x_partes Test 15.4 3.54 332155. 342893.lax_pass_fit <- lax_pass_fit |>
mutate(combinado2 = (combinado + prophet)/2,
combinado3 = (prophet + prophet_auto + ets + sarima)/4)
lax_pass_fc <- lax_pass_fit |>
select(combinado, combinado2, combinado3, prophet, prophet_auto, sarima, ets) |>
forecast(h = "2 years")
p <- lax_pass_fc |>
ggplot(aes(x = month, y = .mean)) +
geom_line(data = lax_passengers, aes(y = passengers)) +
geom_line(aes(color = .model)) +
facet_wrap(~ type, scales = "free_y", ncol = 1)Warning: The output of `fortify(<fable>)` has changed to better suit usage with the ggdist package.
If you're using it to extract intervals, consider using `hilo()` to compute intervals, and `unpack_hilo()` to obtain values.plotly::ggplotly(p, dynamicTicks = TRUE) |>
plotly::rangeslider()lax_pass_fc |>
accuracy(lax_passengers) |>
arrange(type, MAPE) |>
select(type, .model, .type, MAPE, MASE, MAE, RMSE)# A tibble: 14 × 7
type .model .type MAPE MASE MAE RMSE
<chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
1 Domestic combinado Test 1.07 0.286 53841. 70415.
2 Domestic combinado3 Test 1.22 0.326 61380. 81799.
3 Domestic sarima Test 1.46 0.395 74318. 88136.
4 Domestic combinado2 Test 1.89 0.513 96607. 117079.
5 Domestic ets Test 2.40 0.636 119747. 139337.
6 Domestic prophet_auto Test 2.75 0.713 134284. 174305.
7 Domestic prophet Test 3.53 0.969 182422. 209294.
8 International combinado2 Test 2.89 0.645 60586. 75587.
9 International prophet Test 3.05 0.676 63494. 87306.
10 International combinado3 Test 3.07 0.696 65357. 78713.
11 International sarima Test 3.28 0.748 70206. 82727.
12 International combinado Test 3.48 0.801 75240. 87041.
13 International prophet_auto Test 4.10 0.914 85864. 112493.
14 International ets Test 4.85 1.11 104617. 115124.Quarto
Quarto enables you to weave together content and executable code into a finished document. To learn more about Quarto see https://quarto.org.
Running Code
When you click the Render button a document will be generated that includes both content and the output of embedded code. You can embed code like this:
1 + 1[1] 2You can add options to executable code like this
[1] 4The echo: false option disables the printing of code (only output is displayed).