Time series Analysis - Formula 1
A simple time series forecast using search trends for Formula 1
JPitt
1) Introduction
Time series analysis is the study of observations that are measured sequentially over time, with the aim of identifying patterns and forecasting future observations. Time series data may present trends, seasonality, and autocorrelation, which can be assessed using decomposition and statistical modelling techniques. Models such as Autoregressive Integrated Moving Average (ARIMA) and Exponential Smoothing (ETS) are examples of techniques that can be used to create forecasts and predictive evaluations.
Formula 1 is a worldwide racing competition and is considered to be the pinnacle of Motorsport. A Formula 1 season consists of a series of races that take place in multiple countries and continents. Therefore, it is hypothesised that Google search interest in Formula 1 will demonstrate an increasing longitudinal trend from a global perspective. Consequently, the aim of the report is to perform time series analysis and forecasting procedures on the Google search trends for Formula 1.
Research Question: How has Google search interest for Formula 1 changed over time, and how accurate can extrapolation procedures predict future events using a range of time-series models?
2) Load Libraries
First of all, we load all required packages for time series
modelling. Use the install.packages() function if you do
not have them installed.
library(tidymodels)
library(tidyverse)
library(timetk)
library(fpp3)
library(forecast)
library(feasts)
library(fabletools)
library(modeltime)
library(rsample)
library(f1dataR)3) Import Data
Data was sourced from Google trends using the search term “Formula
1”. Observations range from 2004 to 2026 and the search location was set
to “Global”. This was exported as a .csv file so it needs to be loaded
using the read_csv() function.
4) Exploratory Data Analysis
Let’s plot the data to assess the search trends over time. Use ggplot() to create a line graph to observe fluctuations in the data.
f1 |>
ggplot(aes(x = time, y = formula_1)) +
geom_line(colour = "steelblue") +
geom_point(size = 0.8) +
labs(
title = "Formula 1 Search Trends",
x = "Time",
y = "Search Frequency"
) +
theme_classic()
Seasonal Diagnostics
Plot the seasonal diagnostics using the timetk package
to provide a more granular understanding of the data by assessing the
trends across each month, quarter, and year.
Anomaly Detection Diagnostics
The Anomaly diagnostics plot from timetk is used to
identify time points that deviate from the expected pattern. From this,
several anomalies were observed. These anomalies may reflect external
events; thus, further research could assess why the search frequencies
deviate from the expected patterns. For example, the 2020 observation
could be due to the COVID-19 pandemic, significantly reducing the
observations during this time.
ACF & PACF Plots
ACF: The Autocorrelation Function (ACF) plots the correlations of the time series with different lags to identify patterns such as trend and seasonality. For example, a Lag 1 indicates correlation between observations at time t and t-1. Lag 2 is t and t-2 etc. Patterns of significant autocorrelation across multiple lags can suggest evidence of trend, seasonality, and possible non-stationarity.
PACF: The Partial Autocorrelation Function (PACF) plot measures the correlation between two observations at a specific lag after accounting for the shorter lags. For example, the partial correlation for Lag 3 represents the relationship at time t and t-3 after accounting for the relationships at Lag 1 and Lag 2.
Below is a visual representation of these plots using the
forecast package.
Interpretation: Based on the ACF plot, the Google search trends in Formula 1 are related to previous observations due to the positive correlations across multiple lags. Evidently, a large spike at lag 12 suggests a 12 month seasonality component. In comparison, the PACF plot shows a large positive partial correlation at lag 1, indicating that the previous month search interest has an important relationship with current observations. The increase in partial correlation at lag 12 also suggests an element of seasonality.
STL Decomposition
Below is a Seasonal Trend and Loess (STL) decomposition plot that is
used to assess the raw data, overall trend, seasonal patterns, and
residuals in a faceted graphic. First, convert the original data into a
tsibble() and plot the STL component using
autoplot().
f1_ts <- f1 |>
mutate(month = yearmonth(time)) |>
select(month, formula_1) |>
as_tsibble(index = month)
f1_stl <- f1_ts |>
model(STL(formula_1))
f1_stl |>
components() |>
autoplot()
Interpretation: STL Decomposition identified a seasonal component and non-linear trend in search interest. Evidently, after the 2020 time point, the trend significantly increases and the magnitude of the troughs within the seasonal pattern become larger.
5) Split Data
Split the data using an 80/20 split and assign this to a train and
test data frame using the rsample package. You can
visualise the splits using the timetk package through
tk_time_series_cv_plan() to prepare the splits and
plot_time_series_cv_plan() to visualise the
data.
splits <- rsample::initial_time_split(f1, prop = 0.8)
train_f1 <- rsample::training(splits)
test_f1 <- rsample::testing(splits)
splits |>
timetk::tk_time_series_cv_plan() |>
timetk::plot_time_series_cv_plan(.date_var = time, .value = formula_1)6) Model Fit
Based on the ACF & PACF interpretation, a range of models can be used that fit the aforementioned assumptions. A list of models can be found using the following link: https://www.tidymodels.org/find/parsnip/
arima_fit <- arima_reg() |>
parsnip::set_engine("auto_arima") |>
parsnip::fit(formula_1 ~ time, data = train_f1)
ets_fit <- exp_smoothing() |>
parsnip::set_engine("ets") |>
parsnip::fit(formula_1 ~ time, data = train_f1)
prophet_fit <- prophet_reg() |>
parsnip::set_engine("prophet") |>
parsnip::fit(formula_1 ~ time, data = train_f1)
snaive_fit <- naive_reg() |>
parsnip::set_engine("snaive") |>
parsnip::fit(formula_1 ~ time, data = train_f1)
stlm_ets_fit <- seasonal_reg(seasonal_period_1 = 12) |>
parsnip::set_engine("stlm_ets") |>
parsnip::fit(formula_1 ~ time, data = train_f1)
stlm_arima_fit <- seasonal_reg(seasonal_period_1 = 12) |>
parsnip::set_engine("stlm_arima") |>
parsnip::fit(formula_1 ~ time, data = train_f1)7) Model Table
Create a model table using the modeltime package so we
can view the performance later.
mod_table <- modeltime::modeltime_table(snaive_fit, arima_fit, prophet_fit, ets_fit, stlm_arima_fit, stlm_ets_fit)8) Calibration
use the modeltime package to calibrate the models using
the test dataset.
9) Forecast
Generate forecasts using the modeltime package to
visually compare model performance to the test observations. Assess the
model performance visually by examining the trend, seasonal patterns,
peaks and toughs, and magnitude of fluctuations.
calibrate |>
modeltime::modeltime_forecast(actual_data = f1, new_data = test_f1) |>
modeltime::plot_modeltime_forecast()Interpretation: Each model demonstrated an ability to capture the seasonal pattern. The SNAIVE, Prophet, ETS, and STLM ARIMA models produced stable forecasting levels. Conversely, the auto ARIMA and STLM ETS showed an annual increasing level of search interest over time. Further, the models underestimated the magnitude of some peaks and troughs. Therefore, to determine the best model, forecasting accuracy will be used to compare each model to the test data.
10) Forecast Results
Check the results of each model using R-squared and error values.
## # A tibble: 6 × 9
## .model_id .model_desc .type mae mape mase smape rmse rsq
## <int> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 SNAIVE [12] Test 14.7 34.6 0.738 25.6 21.5 0.0982
## 2 2 ARIMA(0,1,2)(0,1,1)[12] Test 21.6 55.0 1.09 36.6 26.1 0.151
## 3 3 PROPHET Test 13.5 27.0 0.681 25.8 16.7 0.504
## 4 4 ETS(M,N,A) Test 12.6 29.4 0.635 24.8 15.3 0.586
## 5 5 SEASONAL DECOMP: ARIMA(0… Test 12.3 27.7 0.618 23.5 15.2 0.543
## 6 6 SEASONAL DECOMP: ETS(M,A… Test 25.5 62.7 1.28 41.2 29.7 0.130Interpretation: The STLM ARIMA produced the lowest MAE (12.26) and RMSE (15.23) while producing the second highest R-squared value (0.54). In comparison, the ETS model produced a MAE of 12.61, RMSE of 15.33 and the highest R-squared of 0.58. A visual depiction of each result is shown below.
accuracy <- calibrate |>
modeltime::modeltime_accuracy()
accuracy <- accuracy |>
dplyr::select(.model_id, .model_desc, rmse, rsq, mae) |>
tidyr::pivot_longer(cols = c(rmse, rsq, mae),
names_to = "metric",
values_to = "value")
accuracy |>
ggplot(aes(x = .model_id, y = value, fill = .model_desc)) +
geom_col() +
facet_wrap(~metric, scales = "free_y") +
geom_text(aes(label = round(value, 2)), vjust = -0.5, size = 2, fontface = "bold") +
scale_x_continuous(breaks = seq(1,6, 1)) +
labs(
x = "Model",
y = "Value",
title = "Forecast Performance by Model") +
theme_bw() +
theme(legend.title = element_blank(), legend.position = "bottom")
11) 3 year Forecast
Similar to the previous forecast, we will use extrapolation, where observations are forecast 3 years into the future.
refit_table <- calibrate |>
modeltime::modeltime_refit(data = f1)
refit_table |>
modeltime::modeltime_forecast(h = "3 years", actual_data = f1) |>
modeltime::plot_modeltime_forecast()12) Summary
The STL decomposition showed an increase in trend from approximately 2020 onward and a clear seasonal pattern which was supported by the strong autocorrelation found in the ACF and PACF analysis.
Several anomalies were identified, suggesting external events may have caused these deviations.
The forecasting models reproduced the annual seasonal patterns but differed in their ability to capture the magnitude of the peaks and troughs.
Forecast accuracy differed between models with the STLM-ARIMA demonstrating the greatest predictive accuracy, producing the lowest MAE and RMSE. This was followed closely by the ETS model that produced a low MAE and RMSE, with the highest R-squared.
Overall, the time series forecasting methodologies can be implemented across a range of longitudinal data sets with different modelling techniques selected according to the characteristics of the data.