# Easy time-series prediction with R: a tutorial with air traffic data from Lux Airport

November 13, 2018
By

(This article was first published on Econometrics and Free Software, and kindly contributed to R-bloggers)

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
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 <- function(x){
log(x + sqrt(x**2 + 1))
}``````

This function, the inverse hyperbolic sine, is useful to transform data in a manner that is very
close to logging it, but that allows for 0’s. The data from Eurostat is not complete for some reason,
so there are some 0 sometimes. To avoid having to log 0, which in R yields `-Inf`, I use this
transformation.

``avia <- read_tsv("avia_par_lu.tsv")``
``````## Parsed with column specification:
## cols(
##   .default = col_character()
## )``````
``## See spec(...) for full column specifications.``

Let’s take a look at the data:

``head(avia)``
``````## # A tibble: 6 x 238
##   `unit,tra_meas,… `2018Q1` `2018M03` `2018M02` `2018M01` `2017Q4` `2017Q3`
##
## 1 FLIGHT,CAF_PAS,… 511      172       161       178       502      475
## 2 FLIGHT,CAF_PAS,… :        :         :         :         :        :
## 3 FLIGHT,CAF_PAS,… :        :         :         :         399      306
## 4 FLIGHT,CAF_PAS,… 485      167       151       167       493      497
## 5 FLIGHT,CAF_PAS,… 834      293       267       274       790      728
## 6 FLIGHT,CAF_PAS,… :        :         :         :         :        :
## # ... with 231 more variables: `2017Q2` , `2017Q1` ,
## #   `2017M12` , `2017M11` , `2017M10` , `2017M09` ,
## #   `2017M08` , `2017M07` , `2017M06` , `2017M05` ,
## #   `2017M04` , `2017M03` , `2017M02` , `2017M01` ,
## #   `2017` , `2016Q4` , `2016Q3` , `2016Q2` ,
## #   `2016Q1` , `2016M12` , `2016M11` , `2016M10` ,
## #   `2016M09` , `2016M08` , `2016M07` , `2016M06` ,
## #   `2016M05` , `2016M04` , `2016M03` , `2016M02` ,
## #   `2016M01` , `2016` , `2015Q4` , `2015Q3` ,
## #   `2015Q2` , `2015Q1` , `2015M12` , `2015M11` ,
## #   `2015M10` , `2015M09` , `2015M08` , `2015M07` ,
## #   `2015M06` , `2015M05` , `2015M04` , `2015M03` ,
## #   `2015M02` , `2015M01` , `2015` , `2014Q4` ,
## #   `2014Q3` , `2014Q2` , `2014Q1` , `2014M12` ,
## #   `2014M11` , `2014M10` , `2014M09` , `2014M08` ,
## #   `2014M07` , `2014M06` , `2014M05` , `2014M04` ,
## #   `2014M03` , `2014M02` , `2014M01` , `2014` ,
## #   `2013Q4` , `2013Q3` , `2013Q2` , `2013Q1` ,
## #   `2013M12` , `2013M11` , `2013M10` , `2013M09` ,
## #   `2013M08` , `2013M07` , `2013M06` , `2013M05` ,
## #   `2013M04` , `2013M03` , `2013M02` , `2013M01` ,
## #   `2013` , `2012Q4` , `2012Q3` , `2012Q2` ,
## #   `2012Q1` , `2012M12` , `2012M11` , `2012M10` ,
## #   `2012M09` , `2012M08` , `2012M07` , `2012M06` ,
## #   `2012M05` , `2012M04` , `2012M03` , `2012M02` ,
## #   `2012M01` , `2012` , …``````

So yeah, useless in that state. The first column actually is composed of 3 columns, merged together,
and instead of having one column with the date, and another with the value, we have one column
per date. Some cleaning is necessary before using this data.

Let’s start with going from a wide to a long data set:

``````avia %>%
select("unit,tra_meas,airp_pr\\time", contains("20")) %>%
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("unit,tra_meas,airp_pr\\time", contains("20")) %>%
gather(date, passengers, -`unit,tra_meas,airp_pr\\time`) %>%
separate(col = `unit,tra_meas,airp_pr\\time`, into = c("unit", "tra_meas", "air_pr\\time"), sep = ",")``````

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.

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:

``````avia_clean <- avia %>%
select("unit,tra_meas,airp_pr\\time", contains("20")) %>%
gather(date, passengers, -`unit,tra_meas,airp_pr\\time`) %>%
separate(col = `unit,tra_meas,airp_pr\\time`, into = c("unit", "tra_meas", "air_pr\\time"), sep = ",") %>%
mutate(tra_meas = fct_recode(tra_meas,
`Passengers on board` = "PAS_BRD",
`Passengers on board (arrivals)` = "PAS_BRD_ARR",
`Passengers on board (departures)` = "PAS_BRD_DEP",
`Passengers carried` = "PAS_CRD",
`Passengers carried (arrival)` = "PAS_CRD_ARR",
`Passengers carried (departures)` = "PAS_CRD_DEP",
`Passengers seats available` = "ST_PAS",
`Passengers seats available (arrivals)` = "ST_PAS_ARR",
`Passengers seats available (departures)` = "ST_PAS_DEP",
`Commercial passenger air flights` = "CAF_PAS",
`Commercial passenger air flights (arrivals)` = "CAF_PAS_ARR",
`Commercial passenger air flights (departures)` = "CAF_PAS_DEP")) %>%
mutate(unit = fct_recode(unit,
Passenger = "PAS",
Flight = "FLIGHT",
`Seats and berths` = "SEAT")) %>%
mutate(destination = fct_recode(`air_pr\\time`,
`WIEN-SCHWECHAT` = "LU_ELLX_AT_LOWW",
`BRUSSELS` = "LU_ELLX_BE_EBBR",
`GENEVA` = "LU_ELLX_CH_LSGG",
`ZURICH` = "LU_ELLX_CH_LSZH",
`FRANKFURT/MAIN` = "LU_ELLX_DE_EDDF",
`HAMBURG` = "LU_ELLX_DE_EDDH",
`BERLIN-TEMPELHOF` = "LU_ELLX_DE_EDDI",
`MUENCHEN` = "LU_ELLX_DE_EDDM",
`SAARBRUECKEN` = "LU_ELLX_DE_EDDR",
`BERLIN-TEGEL` = "LU_ELLX_DE_EDDT",
`KOBENHAVN/KASTRUP` = "LU_ELLX_DK_EKCH",
`IRAKLION/NIKOS KAZANTZAKIS` = "LU_ELLX_EL_LGIR",
`FUERTEVENTURA` = "LU_ELLX_ES_GCFV",
`GRAN CANARIA` = "LU_ELLX_ES_GCLP",
`LANZAROTE` = "LU_ELLX_ES_GCRR",
`TENERIFE SUR/REINA SOFIA` = "LU_ELLX_ES_GCTS",
`BARCELONA/EL PRAT` = "LU_ELLX_ES_LEBL",
`MALAGA/COSTA DEL SOL` = "LU_ELLX_ES_LEMG",
`PALMA DE MALLORCA` = "LU_ELLX_ES_LEPA",
`SYSTEM - PARIS` = "LU_ELLX_FR_LF90",
`NICE-COTE D'AZUR` = "LU_ELLX_FR_LFMN",
`PARIS-CHARLES DE GAULLE` = "LU_ELLX_FR_LFPG",
`STRASBOURG-ENTZHEIM` = "LU_ELLX_FR_LFST",
`KEFLAVIK` = "LU_ELLX_IS_BIKF",
`MILANO/MALPENSA` = "LU_ELLX_IT_LIMC",
`BERGAMO/ORIO AL SERIO` = "LU_ELLX_IT_LIME",
`ROMA/FIUMICINO` = "LU_ELLX_IT_LIRF",
`AMSTERDAM/SCHIPHOL` = "LU_ELLX_NL_EHAM",
`WARSZAWA/CHOPINA` = "LU_ELLX_PL_EPWA",
`PORTO` = "LU_ELLX_PT_LPPR",
`LISBOA` = "LU_ELLX_PT_LPPT",
`STOCKHOLM/ARLANDA` = "LU_ELLX_SE_ESSA",
`MONASTIR/HABIB BOURGUIBA` = "LU_ELLX_TN_DTMB",
`ENFIDHA-HAMMAMET INTERNATIONAL` = "LU_ELLX_TN_DTNH",
`ENFIDHA ZINE EL ABIDINE BEN ALI` = "LU_ELLX_TN_DTNZ",
`DJERBA/ZARZIS` = "LU_ELLX_TN_DTTJ",
`ANTALYA (MIL-CIV)` = "LU_ELLX_TR_LTAI",
`ISTANBUL/ATATURK` = "LU_ELLX_TR_LTBA",
`SYSTEM - LONDON` = "LU_ELLX_UK_EG90",
`MANCHESTER` = "LU_ELLX_UK_EGCC",
`LONDON GATWICK` = "LU_ELLX_UK_EGKK",
`LONDON/CITY` = "LU_ELLX_UK_EGLC",
`LONDON HEATHROW` = "LU_ELLX_UK_EGLL",
`LONDON STANSTED` = "LU_ELLX_UK_EGSS",
`NEWARK LIBERTY INTERNATIONAL, NJ.` = "LU_ELLX_US_KEWR",
`O.R TAMBO INTERNATIONAL` = "LU_ELLX_ZA_FAJS")) %>%
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 <- avia_clean %>%
filter(tra_meas == "Passengers on board (arrivals)",
!is.na(passengers)) %>%
filter(str_detect(date, "Q")) %>%
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 <- avia_clean %>%
filter(tra_meas == "Passengers on board (arrivals)",
!is.na(passengers)) %>%
filter(str_detect(date, "M")) %>%
mutate(date = paste0(date, "01")) %>%
mutate(date = ymd(date)) %>%
select(destination, date, passengers)``````

``````avia_clean_monthly %>%
group_by(date) %>%
summarise(total = sum(passengers)) %>%
ggplot() +
ggtitle("Raw data") +
geom_line(aes(y = total, x = date), colour = "#82518c") +
scale_x_date(date_breaks = "1 year", date_labels = "%m-%Y") +
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("Logged data") +
geom_line(aes(y = total_ihs, x = date), colour = "#82518c") +
scale_x_date(date_breaks = "1 year", date_labels = "%m-%Y") +
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 <- avia_clean_monthly %>%
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 <- ihs(avia_clean_train)``

Taking the log, or ihs of the data deals with stabilizing the variance of the time series.

There might also be a need to difference the data. Computing the differences between consecutive
observations makes the time-series stationary. This will be taken care of by the `auto.arima()`
function, if needed. The `auto.arima()` function returns the best ARIMA model according to different
statistical criterions, sach us the AIC, AICc or BIC.

``(model_fit <- auto.arima(logged_data))``
``````## Series: logged_data
## ARIMA(2,1,1)(2,1,0)[12]
##
## Coefficients:
##           ar1      ar2      ma1     sar1     sar2
##       -0.4061  -0.2431  -0.3562  -0.5590  -0.3282
## s.e.   0.2003   0.1432   0.1994   0.0911   0.0871
##
## sigma^2 estimated as 0.004503:  log likelihood=137.11
## AIC=-262.21   AICc=-261.37   BIC=-246.17``````

`auto.arima()` found that the best model would be an \(ARIMA(2, 1, 1)(2, 1, 0)_{12}\). This is an
seasonal autoregressive model, with p = 2, d = 1, q = 1, P = 2 and D = 1.

``model_forecast <- forecast(model_fit, h = 39)``

I can now forecast the model for the next 39 months (which correspond to the data available).

To plot the forecast, one could do a simple call to the plot function. But the resulting plot
is not very aesthetic. To plot my own, I have to grab the data that was forecast, and do some
munging again:

``````point_estimate <- model_forecast\$mean %>%
as_tsibble() %>%
rename(point_estimate = value,
date = index)

upper <- model_forecast\$upper %>%
as_tsibble() %>%
rename(date = index,
upper80 = `80%`,
upper95 = `95%`)

lower <- model_forecast\$lower %>%
as_tsibble() %>%
rename(date = index,
lower80 = `80%`,
lower95 = `95%`)

estimated_data <- reduce(list(point_estimate, upper, lower), full_join, by = "date")``````

`as_tibble()` is a function from the `{tsibble}` package that converts objects that are time-series aware
to time-aware tibbles. If you are not familiar with `ts_tibble()`, I urge you to run the above lines
one by one, and especially to compare `as_tsibble()` with the standard `as_tibble()` from the `{tibble}`
package.

This is how `estimated_data` looks:

``head(estimated_data)``
``````## # A tsibble: 6 x 6 [1M]
##       date point_estimate upper80 upper95 lower80 lower95
##
## 1 2015 Jan           11.9    12.0    12.1    11.8    11.8
## 2 2015 Feb           11.9    12.0    12.0    11.8    11.7
## 3 2015 Mar           12.1    12.2    12.3    12.0    12.0
## 4 2015 Apr           12.2    12.3    12.4    12.1    12.1
## 5 2015 May           12.3    12.4    12.4    12.2    12.1
## 6 2015 Jun           12.3    12.4    12.5    12.2    12.1``````

We can now plot the data, with the forecast, and with the 95% confidence interval:

``````avia_clean_monthly %>%
group_by(date) %>%
summarise(total = sum(passengers)) %>%
mutate(total_ihs = ihs(total)) %>%
ggplot() +
ggtitle("Logged data") +
geom_line(aes(y = total_ihs, x = date), colour = "#82518c") +
scale_x_date(date_breaks = "1 year", date_labels = "%m-%Y") +
geom_ribbon(data = estimated_data, aes(x = date, ymin = lower95, ymax = upper95), fill = "#666018", alpha = 0.2) +
geom_line(data = estimated_data, aes(x = date, y = point_estimate), linetype = 2, colour = "#8e9d98") +
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.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...