In this blog post, I will show you how you can quickly and easily forecast a univariate time series.
I am going to use data from the EU Open Data Portal on air passenger transport. You can find the
data here. I downloaded
the data in the TSV format for Luxembourg Airport, but you could repeat the analysis for any airport.
Once you have the data, load some of the package we are going to need:
library(tidyverse)
library(lubridate)
library(forecast)
library(tsibble)
library(brotools)
and define the following function:
ihs %
gather(date, passengers, -`unit,tra_meas,airp_pr\\time`)
The first line makes it possible to only select the columns that contain the string “20”, so
selecting columns from 2000 onward. Then, using gather, I go from long to wide. The data looks
like this now:
## # A tibble: 117,070 x 3
## `unit,tra_meas,airp_pr\\time` date passengers
## ## 1 FLIGHT,CAF_PAS,LU_ELLX_AT_LOWW 2018Q1 511 ## 2 FLIGHT,CAF_PAS,LU_ELLX_BE_EBBR 2018Q1 : ## 3 FLIGHT,CAF_PAS,LU_ELLX_CH_LSGG 2018Q1 : ## 4 FLIGHT,CAF_PAS,LU_ELLX_CH_LSZH 2018Q1 485 ## 5 FLIGHT,CAF_PAS,LU_ELLX_DE_EDDF 2018Q1 834 ## 6 FLIGHT,CAF_PAS,LU_ELLX_DE_EDDI 2018Q1 : ## 7 FLIGHT,CAF_PAS,LU_ELLX_DE_EDDM 2018Q1 1095 ## 8 FLIGHT,CAF_PAS,LU_ELLX_DE_EDDR 2018Q1 : ## 9 FLIGHT,CAF_PAS,LU_ELLX_DE_EDDT 2018Q1 : ## 10 FLIGHT,CAF_PAS,LU_ELLX_DK_EKCH 2018Q1 : ## # ... with 117,060 more rows
Now, let’s separate the first column into 3 columns:
avia %__%
select(&quot;unit,tra_meas,airp_pr\\time&quot;, contains(&quot;20&quot;)) %__%
gather(date, passengers, -`unit,tra_meas,airp_pr\\time`) %__%
separate(col = `unit,tra_meas,airp_pr\\time`, into = c(&quot;unit&quot;, &quot;tra_meas&quot;, &quot;air_pr\\time&quot;), sep = &quot;,&quot;)
This separates the first column into 3 new columns, “unit”, “tra_meas” and “air_pr\time”. This step
is not necessary for the rest of the analysis, but might as well do it. The data looks like this now:
## # A tibble: 117,070 x 5
## unit tra_meas `air_pr\\time` date passengers
## ## 1 FLIGHT CAF_PAS LU_ELLX_AT_LOWW 2018Q1 511 ## 2 FLIGHT CAF_PAS LU_ELLX_BE_EBBR 2018Q1 : ## 3 FLIGHT CAF_PAS LU_ELLX_CH_LSGG 2018Q1 : ## 4 FLIGHT CAF_PAS LU_ELLX_CH_LSZH 2018Q1 485 ## 5 FLIGHT CAF_PAS LU_ELLX_DE_EDDF 2018Q1 834 ## 6 FLIGHT CAF_PAS LU_ELLX_DE_EDDI 2018Q1 : ## 7 FLIGHT CAF_PAS LU_ELLX_DE_EDDM 2018Q1 1095 ## 8 FLIGHT CAF_PAS LU_ELLX_DE_EDDR 2018Q1 : ## 9 FLIGHT CAF_PAS LU_ELLX_DE_EDDT 2018Q1 : ## 10 FLIGHT CAF_PAS LU_ELLX_DK_EKCH 2018Q1 : ## # ... with 117,060 more rows
The next steps are simple renamings. I have copy-pasted the information from the Eurostat page
where you can view the data online.
If you click here:
you will be able to select the variables you want displayed in the table, as well as the dictionary
of the variables. I simply copy pasted it and recoded the variables. You can take a look at the
whole cleaning workflow by clicking “Click to expand” below:
Click here to take a look at the whole cleaning workflow
avia_clean %
select(&quot;unit,tra_meas,airp_pr\\time&quot;, contains(&quot;20&quot;)) %__%
gather(date, passengers, -`unit,tra_meas,airp_pr\\time`) %__%
separate(col = `unit,tra_meas,airp_pr\\time`, into = c(&quot;unit&quot;, &quot;tra_meas&quot;, &quot;air_pr\\time&quot;), sep = &quot;,&quot;) %__%
mutate(tra_meas = fct_recode(tra_meas,
`Passengers on board` = &quot;PAS_BRD&quot;,
`Passengers on board (arrivals)` = &quot;PAS_BRD_ARR&quot;,
`Passengers on board (departures)` = &quot;PAS_BRD_DEP&quot;,
`Passengers carried` = &quot;PAS_CRD&quot;,
`Passengers carried (arrival)` = &quot;PAS_CRD_ARR&quot;,
`Passengers carried (departures)` = &quot;PAS_CRD_DEP&quot;,
`Passengers seats available` = &quot;ST_PAS&quot;,
`Passengers seats available (arrivals)` = &quot;ST_PAS_ARR&quot;,
`Passengers seats available (departures)` = &quot;ST_PAS_DEP&quot;,
`Commercial passenger air flights` = &quot;CAF_PAS&quot;,
`Commercial passenger air flights (arrivals)` = &quot;CAF_PAS_ARR&quot;,
`Commercial passenger air flights (departures)` = &quot;CAF_PAS_DEP&quot;)) %__%
mutate(unit = fct_recode(unit,
Passenger = &quot;PAS&quot;,
Flight = &quot;FLIGHT&quot;,
`Seats and berths` = &quot;SEAT&quot;)) %__%
mutate(destination = fct_recode(`air_pr\\time`,
`WIEN-SCHWECHAT` = &quot;LU_ELLX_AT_LOWW&quot;,
`BRUSSELS` = &quot;LU_ELLX_BE_EBBR&quot;,
`GENEVA` = &quot;LU_ELLX_CH_LSGG&quot;,
`ZURICH` = &quot;LU_ELLX_CH_LSZH&quot;,
`FRANKFURT/MAIN` = &quot;LU_ELLX_DE_EDDF&quot;,
`HAMBURG` = &quot;LU_ELLX_DE_EDDH&quot;,
`BERLIN-TEMPELHOF` = &quot;LU_ELLX_DE_EDDI&quot;,
`MUENCHEN` = &quot;LU_ELLX_DE_EDDM&quot;,
`SAARBRUECKEN` = &quot;LU_ELLX_DE_EDDR&quot;,
`BERLIN-TEGEL` = &quot;LU_ELLX_DE_EDDT&quot;,
`KOBENHAVN/KASTRUP` = &quot;LU_ELLX_DK_EKCH&quot;,
`HURGHADA / INTL` = &quot;LU_ELLX_EG_HEGN&quot;,
`IRAKLION/NIKOS KAZANTZAKIS` = &quot;LU_ELLX_EL_LGIR&quot;,
`FUERTEVENTURA` = &quot;LU_ELLX_ES_GCFV&quot;,
`GRAN CANARIA` = &quot;LU_ELLX_ES_GCLP&quot;,
`LANZAROTE` = &quot;LU_ELLX_ES_GCRR&quot;,
`TENERIFE SUR/REINA SOFIA` = &quot;LU_ELLX_ES_GCTS&quot;,
`BARCELONA/EL PRAT` = &quot;LU_ELLX_ES_LEBL&quot;,
`ADOLFO SUAREZ MADRID-BARAJAS` = &quot;LU_ELLX_ES_LEMD&quot;,
`MALAGA/COSTA DEL SOL` = &quot;LU_ELLX_ES_LEMG&quot;,
`PALMA DE MALLORCA` = &quot;LU_ELLX_ES_LEPA&quot;,
`SYSTEM - PARIS` = &quot;LU_ELLX_FR_LF90&quot;,
`NICE-COTE D'AZUR` = &quot;LU_ELLX_FR_LFMN&quot;,
`PARIS-CHARLES DE GAULLE` = &quot;LU_ELLX_FR_LFPG&quot;,
`STRASBOURG-ENTZHEIM` = &quot;LU_ELLX_FR_LFST&quot;,
`KEFLAVIK` = &quot;LU_ELLX_IS_BIKF&quot;,
`MILANO/MALPENSA` = &quot;LU_ELLX_IT_LIMC&quot;,
`BERGAMO/ORIO AL SERIO` = &quot;LU_ELLX_IT_LIME&quot;,
`ROMA/FIUMICINO` = &quot;LU_ELLX_IT_LIRF&quot;,
`AGADIR/AL MASSIRA` = &quot;LU_ELLX_MA_GMAD&quot;,
`AMSTERDAM/SCHIPHOL` = &quot;LU_ELLX_NL_EHAM&quot;,
`WARSZAWA/CHOPINA` = &quot;LU_ELLX_PL_EPWA&quot;,
`PORTO` = &quot;LU_ELLX_PT_LPPR&quot;,
`LISBOA` = &quot;LU_ELLX_PT_LPPT&quot;,
`STOCKHOLM/ARLANDA` = &quot;LU_ELLX_SE_ESSA&quot;,
`MONASTIR/HABIB BOURGUIBA` = &quot;LU_ELLX_TN_DTMB&quot;,
`ENFIDHA-HAMMAMET INTERNATIONAL` = &quot;LU_ELLX_TN_DTNH&quot;,
`ENFIDHA ZINE EL ABIDINE BEN ALI` = &quot;LU_ELLX_TN_DTNZ&quot;,
`DJERBA/ZARZIS` = &quot;LU_ELLX_TN_DTTJ&quot;,
`ANTALYA (MIL-CIV)` = &quot;LU_ELLX_TR_LTAI&quot;,
`ISTANBUL/ATATURK` = &quot;LU_ELLX_TR_LTBA&quot;,
`SYSTEM - LONDON` = &quot;LU_ELLX_UK_EG90&quot;,
`MANCHESTER` = &quot;LU_ELLX_UK_EGCC&quot;,
`LONDON GATWICK` = &quot;LU_ELLX_UK_EGKK&quot;,
`LONDON/CITY` = &quot;LU_ELLX_UK_EGLC&quot;,
`LONDON HEATHROW` = &quot;LU_ELLX_UK_EGLL&quot;,
`LONDON STANSTED` = &quot;LU_ELLX_UK_EGSS&quot;,
`NEWARK LIBERTY INTERNATIONAL, NJ.` = &quot;LU_ELLX_US_KEWR&quot;,
`O.R TAMBO INTERNATIONAL` = &quot;LU_ELLX_ZA_FAJS&quot;)) %__%
mutate(passengers = as.numeric(passengers)) %__%
select(unit, tra_meas, destination, date, passengers)
## Warning in evalq(as.numeric(passengers), ): NAs introduced by
## coercion
There is quarterly data and monthly data. Let’s separate the two:
avia_clean_quarterly %
filter(tra_meas == &quot;Passengers on board (arrivals)&quot;,
!is.na(passengers)) %__%
filter(str_detect(date, &quot;Q&quot;)) %__%
mutate(date = yq(date))
In the “date” column, I detect the observations with “Q” in their name, indicating that it is quarterly data.
I do the same for monthly data, but I have to add the string “01” to the dates. This transforms
a date that looks like this “2018M1” to this “2018M101”. “2018M101” can then be converted into a
date by using the ymd() function from lubridate. yq() was used for the quarterly data.
avia_clean_monthly %
filter(tra_meas == &quot;Passengers on board (arrivals)&quot;,
!is.na(passengers)) %__%
filter(str_detect(date, &quot;M&quot;)) %__%
mutate(date = paste0(date, &quot;01&quot;)) %__%
mutate(date = ymd(date)) %__%
select(destination, date, passengers)
Time for some plots. Let’s start with the raw data:
avia_clean_monthly %__%
group_by(date) %__%
summarise(total = sum(passengers)) %__%
ggplot() +
ggtitle(&quot;Raw data&quot;) +
geom_line(aes(y = total, x = date), colour = &quot;#82518c&quot;) +
scale_x_date(date_breaks = &quot;1 year&quot;, date_labels = &quot;%m-%Y&quot;) + theme_blog()
And now with the logged data (or rather, the data transformed using the inverted hyperbolic sine
transformation):
avia_clean_monthly %__%
group_by(date) %__%
summarise(total = sum(passengers)) %__%
mutate(total_ihs = ihs(total)) %__%
ggplot() +
ggtitle(&quot;Logged data&quot;) +
geom_line(aes(y = total_ihs, x = date), colour = &quot;#82518c&quot;) +
scale_x_date(date_breaks = &quot;1 year&quot;, date_labels = &quot;%m-%Y&quot;) + theme_blog()
We clearly see a seasonal pattern in the data. There is also an upward trend. We will have to deal
with these two problems if we want to do some forecasting. For this, let’s limit ourselves to data
from before 2015, and convert the “passengers” column from the data to a time series object, using
the ts() function:
avia_clean_train %
select(date, passengers) %__%
filter(year(date) __ 2015) %__%
group_by(date) %__%
summarise(total_passengers = sum(passengers)) %__%
pull(total_passengers) %__%
ts(., frequency = 12, start = c(2005, 1))
We will try to pseudo-forecast the data from 2015 to the last point available, March 2018.
First, let’s tranform the data:
logged_data %
rename(date = index,
upper80 = `80%`,
upper95 = `95%`)
lower %
as_tsibble() %__%
spread(key, value) %__%
rename(date = index,
lower80 = `80%`,
lower95 = `95%`)
estimated_data %
group_by(date) %__%
summarise(total = sum(passengers)) %__%
mutate(total_ihs = ihs(total)) %__%
ggplot() +
ggtitle(&quot;Logged data&quot;) +
geom_line(aes(y = total_ihs, x = date), colour = &quot;#82518c&quot;) +
scale_x_date(date_breaks = &quot;1 year&quot;, date_labels = &quot;%m-%Y&quot;) +
geom_ribbon(data = estimated_data, aes(x = date, ymin = lower95, ymax = upper95), fill = &quot;#666018&quot;, alpha = 0.2) +
geom_line(data = estimated_data, aes(x = date, y = point_estimate), linetype = 2, colour = &quot;#8e9d98&quot;) +
theme_blog()
The pseudo-forecast (the dashed line) is not very far from the truth, only overestimating the
seasonal peaks, but the true line is within the 95% confidence interval, which is good!
If you found this blog post useful, you might want to follow me on twitter for blog post updates.