CA3 - Seasonal Passenger Forecast for the Economy Sector of Dublin Airport¶

Roberto Tavares de Oliveira Neto | Student ID: 2025444¶

Module: Strategic Thinking - Higher Diploma in Data Analytics
Lecturer: Taufique Ahmed
CCT College Dublin | 2025

Table of Contents¶

  1. Exploration Data Analysis - EDA
    1.1 Libraries
    1.2 Data Loading via API
    1.3 Descriptive and Univariate Analysis
    1.4 Data Preparation
    1.5 Feature Understanding and Historical Charts
    1.6 Seasonal Decomposition
    1.7 Outliers Analysis
    1.8 Relational Analysis
  2. Stationarity Testing - ADF Test
  3. Baseline Models
  4. Model 1 - ARIMA
    4.1 Setup A - With COVID
    4.2 Setup B - COVID Excluded
  5. Model 2 - SARIMA
    5.1 Setup A - With COVID
    5.2 Setup B - COVID Excluded + Cross-Validation
  6. Model 3 - Holt-Winters Exponential Smoothing
    6.1 Setup A - With COVID
    6.2 Setup B - Post-COVID Only (2022+)
  7. Model Comparison - Setup A vs Setup B
  8. 12-Month Forecast
  9. Model Interpretability

1. Exploration Data Analysis - EDA¶

1.1 Libraries¶

In [2]:
import warnings
warnings.filterwarnings('ignore')

import logging
logging.captureWarnings(True)
logging.getLogger('py.warnings').setLevel(logging.ERROR)

from statsmodels.tools.sm_exceptions import ConvergenceWarning, ValueWarning
warnings.simplefilter('ignore', ConvergenceWarning)
warnings.simplefilter('ignore', ValueWarning)

import pandas as pd
import numpy as np
from collections import OrderedDict

import matplotlib.pyplot as plt
import seaborn as sns

from statsmodels.tsa.stattools import adfuller
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.tsa.holtwinters import ExponentialSmoothing as HWES

from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score

1.2 Data Loading via API¶

The dataset ASM01 - Air and Sea Travel is published by the Central Statistics Office (CSO) of Ireland. It is accessed via a public REST API and returns data in JSON-stat format, covering monthly arrivals and departures across all Irish airports and seaports since January 2010.

In [3]:
import requests
from pyjstat import pyjstat

url = "https://ws.cso.ie/public/api.restful/PxStat.Data.Cube_API.ReadDataset/ASM01/JSON-stat/2.0/en"

r = requests.get(url)
r.raise_for_status()

json_stat_obj = r.json(object_pairs_hook=OrderedDict)
df_list = pyjstat.from_json_stat(json_stat_obj)
df = pd.DataFrame(df_list[0])

df
Out[3]:
STATISTIC Month Direction Mode value
0 Air and Sea Travel 2010 January Arrivals All modes of transport 773.2
1 Air and Sea Travel 2010 January Arrivals Air 728.4
2 Air and Sea Travel 2010 January Arrivals Sea 44.8
3 Air and Sea Travel 2010 January Departures All modes of transport 791.7
4 Air and Sea Travel 2010 January Departures Air 732.4
... ... ... ... ... ...
1165 Air and Sea Travel 2026 March Arrivals Air 1690.2
1166 Air and Sea Travel 2026 March Arrivals Sea 58.2
1167 Air and Sea Travel 2026 March Departures All modes of transport 1758.9
1168 Air and Sea Travel 2026 March Departures Air 1706.5
1169 Air and Sea Travel 2026 March Departures Sea 52.4

1170 rows × 5 columns

1.3 Descriptive and Univariate Analysis¶

In [4]:
print(f"Total Variables   : {df.shape[1]}")
print(f"Total Observations: {df.shape[0]}")

Total Variables   : 5
Total Observations: 1170
In [5]:
df.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 1170 entries, 0 to 1169
Data columns (total 5 columns):
 #   Column     Non-Null Count  Dtype  
---  ------     --------------  -----  
 0   STATISTIC  1170 non-null   object 
 1   Month      1170 non-null   object 
 2   Direction  1170 non-null   object 
 3   Mode       1170 non-null   object 
 4   value      1170 non-null   float64
dtypes: float64(1), object(4)
memory usage: 45.8+ KB
In [6]:
df.describe()
Out[6]:
value
count 1170.000000
mean 872.161026
std 695.686830
min 4.000000
25% 109.600000
50% 942.100000
75% 1425.000000
max 2513.900000
In [7]:
df.describe(include=object)
Out[7]:
STATISTIC Month Direction Mode
count 1170 1170 1170 1170
unique 1 195 2 3
top Air and Sea Travel 2010 January Arrivals All modes of transport
freq 1170 6 585 390
In [8]:
df.duplicated().sum()
Out[8]:
0
In [9]:
df.isnull().sum()
Out[9]:
STATISTIC    0
Month        0
Direction    0
Mode         0
value        0
dtype: int64
In [10]:
df.columns
Out[10]:
Index(['STATISTIC', 'Month', 'Direction', 'Mode', 'value'], dtype='object')
In [11]:
# Initial chart to inspect the overall passenger volume across all modes and directions
plt.figure(figsize=(20, 6))
sns.lineplot(x=df['Month'], y=df['value'], color='steelblue')
plt.xlabel('Date')
plt.ylabel('Passengers (thousands)')
plt.title('Total Passenger Movements - Ireland (Air + Sea, 2010-2024)')
plt.xticks(rotation=45)
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
No description has been provided for this image
In [12]:
# Identify the temporal range of the dataset
first = pd.to_datetime(df['Month']).dt.date.min()
last  = pd.to_datetime(df['Month']).dt.date.max()
print(f"First observation : {first}")
print(f"Last observation  : {last}")

First observation : 2010-01-01
Last observation  : 2026-03-01

1.4 Data Preparation¶

In [13]:
# The STATISTIC column contains a single constant value and adds no analytical value
df = df.drop(columns=['STATISTIC'])
df.head()
Out[13]:
Month Direction Mode value
0 2010 January Arrivals All modes of transport 773.2
1 2010 January Arrivals Air 728.4
2 2010 January Arrivals Sea 44.8
3 2010 January Departures All modes of transport 791.7
4 2010 January Departures Air 732.4
In [14]:
# Filter for air transport only - the focus of this project is Dublin Airport
df_air = df[df["Mode"] == "Air"].copy()

# Convert the Month string to a proper datetime and extract year and month components
df_air['Date']      = pd.to_datetime(df_air['Month'])
df_air['Year']      = pd.DatetimeIndex(df_air['Date']).year
df_air['Month_num'] = pd.DatetimeIndex(df_air['Date']).month

df_air.head()
Out[14]:
Month Direction Mode value Date Year Month_num
1 2010 January Arrivals Air 728.4 2010-01-01 2010 1
4 2010 January Departures Air 732.4 2010-01-01 2010 1
7 2010 February Arrivals Air 757.2 2010-02-01 2010 2
10 2010 February Departures Air 757.2 2010-02-01 2010 2
13 2010 March Arrivals Air 906.6 2010-03-01 2010 3
In [15]:
# Reorganise columns and derive the Dublin Airport volume
# Dublin Airport handles approximately 83% of all Irish air passengers (CSO, 2025)
df_air = df_air[['Direction', 'Date', 'Year', 'Month_num', 'value']]
df_air['Dublin_airport'] = df_air['value'] * 0.83

df_air = df_air.sort_values('Date').reset_index(drop=True)
df_air.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 390 entries, 0 to 389
Data columns (total 6 columns):
 #   Column          Non-Null Count  Dtype         
---  ------          --------------  -----         
 0   Direction       390 non-null    object        
 1   Date            390 non-null    datetime64[ns]
 2   Year            390 non-null    int32         
 3   Month_num       390 non-null    int32         
 4   value           390 non-null    float64       
 5   Dublin_airport  390 non-null    float64       
dtypes: datetime64[ns](1), float64(2), int32(2), object(1)
memory usage: 15.4+ KB
In [16]:
# Visualise annual passenger volume to observe trend and the COVID disruption
plt.figure(figsize=(18, 5))
plt.bar(df_air['Year'], df_air['Dublin_airport'], color='steelblue', alpha=0.8)
plt.grid(True, axis='y', linestyle='--', alpha=0.4)
plt.xlabel('Year')
plt.ylabel('Passengers (thousands)')
plt.title('Total Passengers per Year - Dublin Airport')
plt.tight_layout()
plt.show()
No description has been provided for this image
In [17]:
# Calculate month-over-month growth rate per direction
df_air['Growth_rate'] = df_air.groupby('Direction')['Dublin_airport'].pct_change()

# Separate arrivals and departures for independent analysis
df_arr = df_air[df_air["Direction"] == "Arrivals"].drop(columns=["Direction"]).reset_index(drop=True)
df_dep = df_air[df_air["Direction"] == "Departures"].drop(columns=["Direction"]).reset_index(drop=True)

print(f"Arrivals   shape: {df_arr.shape}")
print(f"Departures shape: {df_dep.shape}")

Arrivals   shape: (195, 6)
Departures shape: (195, 6)
In [18]:
display(df_arr.style.format({'Growth_rate': '{:.1%}'}))
  Date Year Month_num value Dublin_airport Growth_rate
0 2010-01-01 00:00:00 2010 1 728.400000 604.572000 nan%
1 2010-02-01 00:00:00 2010 2 757.200000 628.476000 4.0%
2 2010-03-01 00:00:00 2010 3 906.600000 752.478000 19.7%
3 2010-04-01 00:00:00 2010 4 735.800000 610.714000 -18.8%
4 2010-05-01 00:00:00 2010 5 968.500000 803.855000 31.6%
5 2010-06-01 00:00:00 2010 6 1087.800000 902.874000 12.3%
6 2010-07-01 00:00:00 2010 7 1286.800000 1068.044000 18.3%
7 2010-08-01 00:00:00 2010 8 1300.400000 1079.332000 1.1%
8 2010-09-01 00:00:00 2010 9 1091.200000 905.696000 -16.1%
9 2010-10-01 00:00:00 2010 10 1059.000000 878.970000 -3.0%
10 2010-11-01 00:00:00 2010 11 792.200000 657.526000 -25.2%
11 2010-12-01 00:00:00 2010 12 686.700000 569.961000 -13.3%
12 2011-01-01 00:00:00 2011 1 715.800000 594.114000 4.2%
13 2011-02-01 00:00:00 2011 2 722.100000 599.343000 0.9%
14 2011-03-01 00:00:00 2011 3 880.000000 730.400000 21.9%
15 2011-04-01 00:00:00 2011 4 963.200000 799.456000 9.5%
16 2011-05-01 00:00:00 2011 5 1080.900000 897.147000 12.2%
17 2011-06-01 00:00:00 2011 6 1114.200000 924.786000 3.1%
18 2011-07-01 00:00:00 2011 7 1302.300000 1080.909000 16.9%
19 2011-08-01 00:00:00 2011 8 1285.400000 1066.882000 -1.3%
20 2011-09-01 00:00:00 2011 9 1073.900000 891.337000 -16.5%
21 2011-10-01 00:00:00 2011 10 971.200000 806.096000 -9.6%
22 2011-11-01 00:00:00 2011 11 785.800000 652.214000 -19.1%
23 2011-12-01 00:00:00 2011 12 752.600000 624.658000 -4.2%
24 2012-01-01 00:00:00 2012 1 683.300000 567.139000 -9.2%
25 2012-02-01 00:00:00 2012 2 714.100000 592.703000 4.5%
26 2012-03-01 00:00:00 2012 3 867.600000 720.108000 21.5%
27 2012-04-01 00:00:00 2012 4 996.000000 826.680000 14.8%
28 2012-05-01 00:00:00 2012 5 1039.900000 863.117000 4.4%
29 2012-06-01 00:00:00 2012 6 1149.100000 953.753000 10.5%
30 2012-07-01 00:00:00 2012 7 1283.800000 1065.554000 11.7%
31 2012-08-01 00:00:00 2012 8 1309.300000 1086.719000 2.0%
32 2012-09-01 00:00:00 2012 9 1120.500000 930.015000 -14.4%
33 2012-10-01 00:00:00 2012 10 1035.100000 859.133000 -7.6%
34 2012-11-01 00:00:00 2012 11 812.400000 674.292000 -21.5%
35 2012-12-01 00:00:00 2012 12 795.500000 660.265000 -2.1%
36 2013-01-01 00:00:00 2013 1 684.400000 568.052000 -14.0%
37 2013-02-01 00:00:00 2013 2 704.600000 584.818000 3.0%
38 2013-03-01 00:00:00 2013 3 923.200000 766.256000 31.0%
39 2013-04-01 00:00:00 2013 4 996.500000 827.095000 7.9%
40 2013-05-01 00:00:00 2013 5 1112.900000 923.707000 11.7%
41 2013-06-01 00:00:00 2013 6 1194.400000 991.352000 7.3%
42 2013-07-01 00:00:00 2013 7 1374.900000 1141.167000 15.1%
43 2013-08-01 00:00:00 2013 8 1366.900000 1134.527000 -0.6%
44 2013-09-01 00:00:00 2013 9 1156.500000 959.895000 -15.4%
45 2013-10-01 00:00:00 2013 10 1082.600000 898.558000 -6.4%
46 2013-11-01 00:00:00 2013 11 844.900000 701.267000 -22.0%
47 2013-12-01 00:00:00 2013 12 846.000000 702.180000 0.1%
48 2014-01-01 00:00:00 2014 1 738.500000 612.955000 -12.7%
49 2014-02-01 00:00:00 2014 2 749.500000 622.085000 1.5%
50 2014-03-01 00:00:00 2014 3 905.400000 751.482000 20.8%
51 2014-04-01 00:00:00 2014 4 1105.800000 917.814000 22.1%
52 2014-05-01 00:00:00 2014 5 1188.200000 986.206000 7.5%
53 2014-06-01 00:00:00 2014 6 1268.600000 1052.938000 6.8%
54 2014-07-01 00:00:00 2014 7 1451.000000 1204.330000 14.4%
55 2014-08-01 00:00:00 2014 8 1457.600000 1209.808000 0.5%
56 2014-09-01 00:00:00 2014 9 1229.300000 1020.319000 -15.7%
57 2014-10-01 00:00:00 2014 10 1185.700000 984.131000 -3.5%
58 2014-11-01 00:00:00 2014 11 934.400000 775.552000 -21.2%
59 2014-12-01 00:00:00 2014 12 924.800000 767.584000 -1.0%
60 2015-01-01 00:00:00 2015 1 833.600000 691.888000 -9.9%
61 2015-02-01 00:00:00 2015 2 863.500000 716.705000 3.6%
62 2015-03-01 00:00:00 2015 3 1048.500000 870.255000 21.4%
63 2015-04-01 00:00:00 2015 4 1208.900000 1003.387000 15.3%
64 2015-05-01 00:00:00 2015 5 1308.500000 1086.055000 8.2%
65 2015-06-01 00:00:00 2015 6 1450.300000 1203.749000 10.8%
66 2015-07-01 00:00:00 2015 7 1637.800000 1359.374000 12.9%
67 2015-08-01 00:00:00 2015 8 1597.900000 1326.257000 -2.4%
68 2015-09-01 00:00:00 2015 9 1368.500000 1135.855000 -14.4%
69 2015-10-01 00:00:00 2015 10 1346.300000 1117.429000 -1.6%
70 2015-11-01 00:00:00 2015 11 1060.700000 880.381000 -21.2%
71 2015-12-01 00:00:00 2015 12 1020.200000 846.766000 -3.8%
72 2016-01-01 00:00:00 2016 1 963.800000 799.954000 -5.5%
73 2016-02-01 00:00:00 2016 2 992.000000 823.360000 2.9%
74 2016-03-01 00:00:00 2016 3 1222.800000 1014.924000 23.3%
75 2016-04-01 00:00:00 2016 4 1300.100000 1079.083000 6.3%
76 2016-05-01 00:00:00 2016 5 1453.300000 1206.239000 11.8%
77 2016-06-01 00:00:00 2016 6 1565.100000 1299.033000 7.7%
78 2016-07-01 00:00:00 2016 7 1764.000000 1464.120000 12.7%
79 2016-08-01 00:00:00 2016 8 1741.500000 1445.445000 -1.3%
80 2016-09-01 00:00:00 2016 9 1502.300000 1246.909000 -13.7%
81 2016-10-01 00:00:00 2016 10 1434.000000 1190.220000 -4.5%
82 2016-11-01 00:00:00 2016 11 1163.100000 965.373000 -18.9%
83 2016-12-01 00:00:00 2016 12 1130.400000 938.232000 -2.8%
84 2017-01-01 00:00:00 2017 1 1050.900000 872.247000 -7.0%
85 2017-02-01 00:00:00 2017 2 1023.100000 849.173000 -2.6%
86 2017-03-01 00:00:00 2017 3 1253.300000 1040.239000 22.5%
87 2017-04-01 00:00:00 2017 4 1430.500000 1187.315000 14.1%
88 2017-05-01 00:00:00 2017 5 1500.600000 1245.498000 4.9%
89 2017-06-01 00:00:00 2017 6 1646.300000 1366.429000 9.7%
90 2017-07-01 00:00:00 2017 7 1845.200000 1531.516000 12.1%
91 2017-08-01 00:00:00 2017 8 1831.900000 1520.477000 -0.7%
92 2017-09-01 00:00:00 2017 9 1587.000000 1317.210000 -13.4%
93 2017-10-01 00:00:00 2017 10 1482.800000 1230.724000 -6.6%
94 2017-11-01 00:00:00 2017 11 1234.200000 1024.386000 -16.8%
95 2017-12-01 00:00:00 2017 12 1190.300000 987.949000 -3.6%
96 2018-01-01 00:00:00 2018 1 1110.500000 921.715000 -6.7%
97 2018-02-01 00:00:00 2018 2 1057.900000 878.057000 -4.7%
98 2018-03-01 00:00:00 2018 3 1262.400000 1047.792000 19.3%
99 2018-04-01 00:00:00 2018 4 1531.300000 1270.979000 21.3%
100 2018-05-01 00:00:00 2018 5 1632.200000 1354.726000 6.6%
101 2018-06-01 00:00:00 2018 6 1762.500000 1462.875000 8.0%
102 2018-07-01 00:00:00 2018 7 1953.000000 1620.990000 10.8%
103 2018-08-01 00:00:00 2018 8 1927.300000 1599.659000 -1.3%
104 2018-09-01 00:00:00 2018 9 1696.500000 1408.095000 -12.0%
105 2018-10-01 00:00:00 2018 10 1598.500000 1326.755000 -5.8%
106 2018-11-01 00:00:00 2018 11 1305.300000 1083.399000 -18.3%
107 2018-12-01 00:00:00 2018 12 1264.800000 1049.784000 -3.1%
108 2019-01-01 00:00:00 2019 1 1174.600000 974.918000 -7.1%
109 2019-02-01 00:00:00 2019 2 1142.400000 948.192000 -2.7%
110 2019-03-01 00:00:00 2019 3 1400.500000 1162.415000 22.6%
111 2019-04-01 00:00:00 2019 4 1596.900000 1325.427000 14.0%
112 2019-05-01 00:00:00 2019 5 1712.000000 1420.960000 7.2%
113 2019-06-01 00:00:00 2019 6 1823.700000 1513.671000 6.5%
114 2019-07-01 00:00:00 2019 7 2017.700000 1674.691000 10.6%
115 2019-08-01 00:00:00 2019 8 2033.600000 1687.888000 0.8%
116 2019-09-01 00:00:00 2019 9 1764.400000 1464.452000 -13.2%
117 2019-10-01 00:00:00 2019 10 1646.300000 1366.429000 -6.7%
118 2019-11-01 00:00:00 2019 11 1305.100000 1083.233000 -20.7%
119 2019-12-01 00:00:00 2019 12 1290.400000 1071.032000 -1.1%
120 2020-01-01 00:00:00 2020 1 1195.200000 992.016000 -7.4%
121 2020-02-01 00:00:00 2020 2 1171.700000 972.511000 -2.0%
122 2020-03-01 00:00:00 2020 3 603.600000 500.988000 -48.5%
123 2020-04-01 00:00:00 2020 4 12.100000 10.043000 -98.0%
124 2020-05-01 00:00:00 2020 5 19.900000 16.517000 64.5%
125 2020-06-01 00:00:00 2020 6 40.200000 33.366000 102.0%
126 2020-07-01 00:00:00 2020 7 188.100000 156.123000 367.9%
127 2020-08-01 00:00:00 2020 8 311.600000 258.628000 65.7%
128 2020-09-01 00:00:00 2020 9 224.000000 185.920000 -28.1%
129 2020-10-01 00:00:00 2020 10 142.500000 118.275000 -36.4%
130 2020-11-01 00:00:00 2020 11 90.900000 75.447000 -36.2%
131 2020-12-01 00:00:00 2020 12 135.900000 112.797000 49.5%
132 2021-01-01 00:00:00 2021 1 101.100000 83.913000 -25.6%
133 2021-02-01 00:00:00 2021 2 50.500000 41.915000 -50.0%
134 2021-03-01 00:00:00 2021 3 55.100000 45.733000 9.1%
135 2021-04-01 00:00:00 2021 4 61.400000 50.962000 11.4%
136 2021-05-01 00:00:00 2021 5 73.200000 60.756000 19.2%
137 2021-06-01 00:00:00 2021 6 136.100000 112.963000 85.9%
138 2021-07-01 00:00:00 2021 7 326.600000 271.078000 140.0%
139 2021-08-01 00:00:00 2021 8 720.400000 597.932000 120.6%
140 2021-09-01 00:00:00 2021 9 743.000000 616.690000 3.1%
141 2021-10-01 00:00:00 2021 10 853.200000 708.156000 14.8%
142 2021-11-01 00:00:00 2021 11 754.500000 626.235000 -11.6%
143 2021-12-01 00:00:00 2021 12 641.500000 532.445000 -15.0%
144 2022-01-01 00:00:00 2022 1 563.400000 467.622000 -12.2%
145 2022-02-01 00:00:00 2022 2 754.600000 626.318000 33.9%
146 2022-03-01 00:00:00 2022 3 1056.800000 877.144000 40.0%
147 2022-04-01 00:00:00 2022 4 1377.100000 1142.993000 30.3%
148 2022-05-01 00:00:00 2022 5 1501.900000 1246.577000 9.1%
149 2022-06-01 00:00:00 2022 6 1605.700000 1332.731000 6.9%
150 2022-07-01 00:00:00 2022 7 1783.600000 1480.388000 11.1%
151 2022-08-01 00:00:00 2022 8 1836.400000 1524.212000 3.0%
152 2022-09-01 00:00:00 2022 9 1626.800000 1350.244000 -11.4%
153 2022-10-01 00:00:00 2022 10 1557.200000 1292.476000 -4.3%
154 2022-11-01 00:00:00 2022 11 1292.800000 1073.024000 -17.0%
155 2022-12-01 00:00:00 2022 12 1237.000000 1026.710000 -4.3%
156 2023-01-01 00:00:00 2023 1 1216.600000 1009.778000 -1.6%
157 2023-02-01 00:00:00 2023 2 1184.200000 982.886000 -2.7%
158 2023-03-01 00:00:00 2023 3 1414.000000 1173.620000 19.4%
159 2023-04-01 00:00:00 2023 4 1676.600000 1391.578000 18.6%
160 2023-05-01 00:00:00 2023 5 1770.500000 1469.515000 5.6%
161 2023-06-01 00:00:00 2023 6 1846.200000 1532.346000 4.3%
162 2023-07-01 00:00:00 2023 7 2030.800000 1685.564000 10.0%
163 2023-08-01 00:00:00 2023 8 2070.200000 1718.266000 1.9%
164 2023-09-01 00:00:00 2023 9 1814.800000 1506.284000 -12.3%
165 2023-10-01 00:00:00 2023 10 1742.100000 1445.943000 -4.0%
166 2023-11-01 00:00:00 2023 11 1352.800000 1122.824000 -22.3%
167 2023-12-01 00:00:00 2023 12 1362.300000 1130.709000 0.7%
168 2024-01-01 00:00:00 2024 1 1271.300000 1055.179000 -6.7%
169 2024-02-01 00:00:00 2024 2 1268.400000 1052.772000 -0.2%
170 2024-03-01 00:00:00 2024 3 1550.200000 1286.666000 22.2%
171 2024-04-01 00:00:00 2024 4 1712.200000 1421.126000 10.5%
172 2024-05-01 00:00:00 2024 5 1889.000000 1567.870000 10.3%
173 2024-06-01 00:00:00 2024 6 1938.600000 1609.038000 2.6%
174 2024-07-01 00:00:00 2024 7 2078.500000 1725.155000 7.2%
175 2024-08-01 00:00:00 2024 8 2205.500000 1830.565000 6.1%
176 2024-09-01 00:00:00 2024 9 1905.900000 1581.897000 -13.6%
177 2024-10-01 00:00:00 2024 10 1764.200000 1464.286000 -7.4%
178 2024-11-01 00:00:00 2024 11 1387.700000 1151.791000 -21.3%
179 2024-12-01 00:00:00 2024 12 1371.600000 1138.428000 -1.2%
180 2025-01-01 00:00:00 2025 1 1283.100000 1064.973000 -6.5%
181 2025-02-01 00:00:00 2025 2 1276.600000 1059.578000 -0.5%
182 2025-03-01 00:00:00 2025 3 1550.600000 1286.998000 21.5%
183 2025-04-01 00:00:00 2025 4 1830.200000 1519.066000 18.0%
184 2025-05-01 00:00:00 2025 5 1996.500000 1657.095000 9.1%
185 2025-06-01 00:00:00 2025 6 2072.900000 1720.507000 3.8%
186 2025-07-01 00:00:00 2025 7 2261.800000 1877.294000 9.1%
187 2025-08-01 00:00:00 2025 8 2318.200000 1924.106000 2.5%
188 2025-09-01 00:00:00 2025 9 2001.400000 1661.162000 -13.7%
189 2025-10-01 00:00:00 2025 10 1920.400000 1593.932000 -4.0%
190 2025-11-01 00:00:00 2025 11 1541.300000 1279.279000 -19.7%
191 2025-12-01 00:00:00 2025 12 1538.000000 1276.540000 -0.2%
192 2026-01-01 00:00:00 2026 1 1442.300000 1197.109000 -6.2%
193 2026-02-01 00:00:00 2026 2 1388.600000 1152.538000 -3.7%
194 2026-03-01 00:00:00 2026 3 1690.200000 1402.866000 21.7%

1.5 Feature Understanding and Historical Charts¶

In [19]:
# Arrivals time series - the primary series used for modelling
plt.figure(figsize=(20, 6))
sns.lineplot(x=df_arr['Date'], y=df_arr['Dublin_airport'], color='crimson', linestyle='--')
plt.xlabel('Date')
plt.ylabel('Passengers (thousands)')
plt.title('Monthly Arrivals - Dublin Airport (2010-2024)')
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
No description has been provided for this image
In [20]:
# Departures time series for comparison
plt.figure(figsize=(20, 6))
sns.lineplot(x=df_dep['Date'], y=df_dep['Dublin_airport'], color='steelblue')
plt.xlabel('Date')
plt.ylabel('Passengers (thousands)')
plt.title('Monthly Departures - Dublin Airport (2010-2024)')
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
No description has been provided for this image
In [21]:
# Overlay arrivals and departures to observe relative magnitude and the COVID impact
plt.figure(figsize=(20, 6))
sns.lineplot(x=df_arr['Date'], y=df_arr['Dublin_airport'],
             color='crimson', linestyle='--', label='Arrivals')
sns.lineplot(x=df_dep['Date'], y=df_dep['Dublin_airport'],
             color='steelblue', label='Departures')
plt.axvspan(pd.Timestamp('2020-01-01'), pd.Timestamp('2022-01-31'),
            alpha=0.12, color='grey', label='COVID-19 Period (excluded from models)')
plt.xlabel('Date')
plt.ylabel('Passengers (thousands)')
plt.title('Arrivals vs Departures - Dublin Airport (2010-2024)')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
No description has been provided for this image
In [22]:
# Average monthly passengers post-COVID (2022-2024)
# This confirms the strong and consistent seasonal pattern driven by summer travel
df_post = df_arr[df_arr['Year'] >= 2022].copy()
monthly_avg = df_post.groupby('Month_num')['Dublin_airport'].mean()
months_lbl  = ['Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec']

plt.figure(figsize=(12, 5))
plt.bar(monthly_avg.index, monthly_avg.values, color='steelblue', alpha=0.85)
plt.xticks(range(1, 13), months_lbl)
plt.ylabel('Average Passengers (thousands)')
plt.title('Average Monthly Passengers - Post-COVID Recovery (2022-2024)')
plt.grid(True, alpha=0.3, axis='y')
plt.tight_layout()
plt.show()
No description has been provided for this image

1.6 Seasonal Decomposition¶

The time series is decomposed into three components: Trend, Seasonality, and Residuals. A logarithmic transformation is applied before decomposition to stabilise the variance, which grows alongside the level of the series.

In [23]:
# Build the full time-indexed series for decomposition
ts_full = df_arr.set_index('Date')['Dublin_airport'].copy()
ts_full = ts_full.asfreq('MS')

# Log transformation stabilises variance; additive decomposition on log scale
log_ts = np.log(ts_full)
log_ts.interpolate(inplace=True)

decomposition = seasonal_decompose(log_ts, model='additive')
fig = decomposition.plot()
fig.set_size_inches(16, 10)
plt.suptitle('Seasonal Decomposition of Log(Arrivals) - Dublin Airport', y=1.01)
plt.tight_layout()
plt.show()
No description has been provided for this image

1.7 Outliers Analysis¶

The COVID-19 pandemic generated extreme low values between January 2020 and January 2022. These observations do not reflect the structural behaviour of the series and would significantly bias any forecasting model. For this reason, this period is excluded from all modelling datasets. The full series is preserved only for contextual visualisation.

In [24]:
fig, ax = plt.subplots(figsize=(7, 5))
sns.boxplot(y=df_air['Dublin_airport'], ax=ax, color='steelblue')
ax.set_title('Boxplot - Dublin Airport Monthly Passengers (All Years)')
ax.set_ylabel('Passengers (thousands)')
plt.tight_layout()
plt.show()
No description has been provided for this image
In [25]:
# Scatter plot to visually isolate the COVID anomaly period
plt.figure(figsize=(18, 5))
sns.scatterplot(x='Date', y='Dublin_airport', data=df_arr, hue='Year', legend=False)
plt.axvspan(pd.Timestamp('2020-01-01'), pd.Timestamp('2022-01-31'),
            alpha=0.12, color='red', label='COVID Outlier Period')
plt.legend()
plt.title('Monthly Arrivals Scatter - COVID Period Highlighted')
plt.xlabel('Date')
plt.ylabel('Passengers (thousands)')
plt.tight_layout()
plt.show()
No description has been provided for this image
In [26]:
# Monthly distribution post-recovery to confirm consistent seasonal structure
df_post2 = df_air[df_air['Year'] >= 2022].copy()
df_post2['Month_name'] = pd.Categorical(
    df_post2['Date'].dt.strftime('%b'),
    categories=['Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec'],
    ordered=True)

plt.figure(figsize=(14, 5))
sns.boxplot(data=df_post2, x='Month_name', y='Dublin_airport', palette='Blues')
plt.title('Passenger Distribution by Month - Post-COVID (2022-2024)')
plt.xlabel('Month')
plt.ylabel('Passengers (thousands)')
plt.grid(True, alpha=0.3, axis='y')
plt.tight_layout()
plt.show()
No description has been provided for this image

1.8 Relational Analysis¶

In [27]:
# Correlation between year, month number and passenger volume
df_corr = df_arr[['Year', 'Month_num', 'Dublin_airport']].corr()

plt.figure(figsize=(7, 5))
sns.heatmap(df_corr, annot=True, fmt='.2f', cmap='coolwarm', linewidths=0.5)
plt.title('Correlation Matrix - Arrivals DataFrame')
plt.tight_layout()
plt.show()
No description has been provided for this image

2. Stationarity Testing - Augmented Dickey-Fuller (ADF)¶

ARIMA and SARIMA require a stationary input series. The Augmented Dickey-Fuller (ADF) test determines whether differencing is needed. The null hypothesis states the series contains a unit root (non-stationary); a p-value below 0.05 allows us to reject it.

Note on data scope: The COVID period (January 2020 to January 2022) is excluded from all modelling. This decision is motivated by the extreme distortion caused by travel restrictions - a structural break that does not represent the long-term passenger dynamics of Dublin Airport.

In [28]:
# Build the clean modelling series: Arrivals with COVID period removed
ts_full = df_arr.set_index('Date')['Dublin_airport'].copy()
ts_full = ts_full.asfreq('MS')

# Exclude January 2020 to January 2022 - anomalous COVID disruption
ts = ts_full[(ts_full.index < '2020-01-01') | (ts_full.index > '2022-02-01')].copy()

print(f"Modelling series: {len(ts)} monthly observations")
print(f"  Pre-COVID : {ts.index[0].date()} to 2019-12-01")
print(f"  Post-COVID: 2022-02-01 to {ts.index[-1].date()}")
print(f"  Removed   : {len(ts_full) - len(ts)} months (COVID disruption)")

Modelling series: 169 monthly observations
  Pre-COVID : 2010-01-01 to 2019-12-01
  Post-COVID: 2022-02-01 to 2026-03-01
  Removed   : 26 months (COVID disruption)
In [29]:
# ADF test on the raw clean series
ADF_result = adfuller(ts)
print(f"ADF Statistic : {ADF_result[0]:.4f}")
print(f"p-value       : {ADF_result[1]:.4f}")
print(f"Critical Values: {ADF_result[4]}")
if ADF_result[1] < 0.05:
    print("\n Stationary series (p < 0.05)")
else:
    print("\n Non-stationary series - first-order differencing required")

ADF Statistic : 0.1719
p-value       : 0.9706
Critical Values: {'1%': -3.4729792729247793, '5%': -2.8802520918255534, '10%': -2.5767469641683105}

 Non-stationary series - first-order differencing required
In [30]:
# Apply first-order differencing and re-test
ts_diff = ts.diff(1).dropna()

ADF_diff = adfuller(ts_diff)
print(f"ADF Statistic (1st diff): {ADF_diff[0]:.4f}")
print(f"p-value (1st diff)       : {ADF_diff[1]:.4f}")
if ADF_diff[1] < 0.05:
    print("\n Differenced series is STATIONARY - d=1 is appropriate")

ADF Statistic (1st diff): -4.9172
p-value (1st diff)       : 0.0000

 Differenced series is STATIONARY - d=1 is appropriate
In [31]:
# Plot raw vs differenced series
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(16, 8), sharex=False)

ax1.plot(ts, color='steelblue')
ax1.set_title('Original Arrivals Series - Dublin Airport (COVID period excluded)')
ax1.set_ylabel('Passengers (thousands)')
ax1.grid(True, alpha=0.3)

ax2.plot(ts_diff, color='crimson')
ax2.set_title('First-Order Differenced Series')
ax2.set_ylabel('Delta Passengers (thousands)')
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()
No description has been provided for this image
In [32]:
# ACF guides the MA(q) order selection; PACF guides the AR(p) order selection
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(14, 8))
plot_acf(ts_diff,  lags=36, ax=ax1)
ax1.set_title('Autocorrelation Function (ACF) - Differenced Series')
plot_pacf(ts_diff, lags=36, ax=ax2)
ax2.set_title('Partial Autocorrelation Function (PACF) - Differenced Series')
plt.tight_layout()
plt.show()
No description has been provided for this image

3. Baseline Models¶

Four simple forecasting baselines are established before applying statistical models. These provide a performance floor - any model that cannot beat them offers no real predictive value. MAPE (Mean Absolute Percentage Error) is the primary evaluation metric here.

In [33]:
# Train / Test split - last 12 months reserved for testing
train = ts[:-12]
test  = ts[-12:]

print(f"Train: {len(train)} months  ({train.index[0].date()} to {train.index[-1].date()})")
print(f"Test : {len(test)} months  ({test.index[0].date()} to {test.index[-1].date()})")

Train: 157 months  (2010-01-01 to 2025-03-01)
Test : 12 months  (2025-04-01 to 2026-03-01)
In [34]:
def mape(y_true, y_pred):
    """Mean Absolute Percentage Error"""
    return np.mean(np.abs((np.array(y_true) - np.array(y_pred)) / np.array(y_true))) * 100
In [35]:
# Baseline 1 - Historical mean of the entire training period
hist_mean = np.mean(train)
test_df = test.to_frame(name='Actual')
test_df['pred_mean'] = hist_mean
mape_hist = mape(test_df['Actual'], test_df['pred_mean'])
print(f"Historical Mean = {hist_mean:.1f}k | MAPE = {mape_hist:.2f}%")

Historical Mean = 1079.3k | MAPE = 27.10%
In [36]:
# Baseline 2 - Mean of the last 12 months in the training set
last_yr_mean = np.mean(train[-12:])
test_df['pred_last_yr_mean'] = last_yr_mean
mape_last_yr = mape(test_df['Actual'], test_df['pred_last_yr_mean'])
print(f"Last Year Mean = {last_yr_mean:.1f}k | MAPE = {mape_last_yr:.2f}%")

Last Year Mean = 1408.5k | MAPE = 14.98%
In [37]:
# Baseline 3 - Last known observation repeated for all test periods
last_val = train.iloc[-1]
test_df['pred_last'] = last_val
mape_last = mape(test_df['Actual'], test_df['pred_last'])
print(f"Last Known Value = {last_val:.1f}k | MAPE = {mape_last:.2f}%")

Last Known Value = 1287.0k | MAPE = 16.50%
In [38]:
# Baseline 4 - Naive seasonal: repeat the last full 12-month seasonal cycle
test_df['pred_naive_seasonal'] = train[-12:].values
mape_naive = mape(test_df['Actual'], test_df['pred_naive_seasonal'])
print(f"Naive Seasonal MAPE = {mape_naive:.2f}%")

Naive Seasonal MAPE = 7.69%
In [39]:
fig, ax = plt.subplots(figsize=(14, 6))
ax.plot(train.index[-24:], train[-24:], 'g-.', label='Train (last 24 months)')
ax.plot(test_df.index, test_df['Actual'],             'b-',  label='Actual')
ax.plot(test_df.index, test_df['pred_mean'],           'r--', label=f'Hist. Mean  (MAPE={mape_hist:.1f}%)')
ax.plot(test_df.index, test_df['pred_last_yr_mean'],   'm-.', label=f'Last Yr Mean (MAPE={mape_last_yr:.1f}%)')
ax.plot(test_df.index, test_df['pred_last'],           'y:',  label=f'Last Value  (MAPE={mape_last:.1f}%)')
ax.plot(test_df.index, test_df['pred_naive_seasonal'], 'k--', label=f'Naive Seasonal (MAPE={mape_naive:.1f}%)')
ax.set_xlabel('Date')
ax.set_ylabel('Passengers (thousands)')
ax.set_title('Baseline Forecasts vs Actual - Dublin Airport Arrivals')
ax.legend(loc='upper left')
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
No description has been provided for this image
In [40]:
fig, ax = plt.subplots(figsize=(9, 4))
labels = ['Hist. Mean', 'Last Yr Mean', 'Last Value', 'Naive Seasonal']
mapes  = [mape_hist, mape_last_yr, mape_last, mape_naive]
bars = ax.bar(labels, mapes, width=0.5, color=['#e74c3c','#e67e22','#f1c40f','#2ecc71'])
for b, v in zip(bars, mapes):
    ax.text(b.get_x()+b.get_width()/2, v+0.5, f'{v:.1f}%', ha='center', fontsize=10)
ax.set_ylabel('MAPE (%)')
ax.set_title('Baseline Model Comparison - MAPE')
ax.set_ylim(0, max(mapes)*1.25)
ax.grid(True, alpha=0.3, axis='y')
plt.tight_layout()
plt.show()
print(f"\nBest baseline: Naive Seasonal MAPE = {mape_naive:.2f}%")
No description has been provided for this image 
Best baseline: Naive Seasonal MAPE = 7.69%

Setup A vs Setup B – Methodology Note¶

To measure the impact of the COVID-19 pandemic on model accuracy, each model is evaluated under two configurations:

  • Setup A – Full Series: January 2010 to February 2026, including the pandemic period. Included for reference to demonstrate the impact of the structural break.
  • Setup B – COVID Excluded (Primary): Pandemic period removed. Primary setup for business forecasting.
  • Holt-Winters exception: Setup B uses only post-COVID data (2022+) because the model's exponential smoothing mechanism is heavily influenced by the recovery gradient. Training on the full clean series (2010-2019 + 2022-2026) causes the model to project an unrealistic upward trend.

All metrics are computed against the same 12-month test set: March 2025 – February 2026.

4. Model 1 - ARIMA¶

ARIMA (Autoregressive Integrated Moving Average) combines three components: p (autoregressive order), d (degree of differencing), and q (moving average order). A grid search over (p,d,q) combinations is performed, selecting the model with the lowest AIC (Akaike Information Criterion).

4.1 Setup A – ARIMA With COVID¶

Reference only. ARIMA trained on the full series including the pandemic structural break.

In [41]:
# Grid search over (p, d, q) to find the combination with the lowest AIC
aicVal = []

for d in range(1, 3):
    for p in range(0, 4):
        for q in range(0, 4):
            try:
                model  = ARIMA(train, order=(p, d, q))
                fitted = model.fit()
                aicVal.append([p, d, q, fitted.aic])
            except Exception:
                pass

aic_df = pd.DataFrame(aicVal, columns=['p','d','q','AIC']).sort_values('AIC').reset_index(drop=True)
print('Top 10 ARIMA models by AIC:');
print(aic_df.head(10).to_string(index=False));

Top 10 ARIMA models by AIC:
 p  d  q         AIC
 2  1  3 1804.901514
 3  1  3 1807.080379
 3  1  2 1816.398922
 2  1  2 1837.634554
 3  1  1 1871.016888
 1  1  3 1914.353393
 1  2  3 1927.351011
 0  2  3 1930.480630
 3  1  0 1934.736349
 3  2  2 1934.784497
In [42]:
best   = aic_df.iloc[0]
p_best = int(best['p'])
d_best = int(best['d'])
q_best = int(best['q'])
print(f"Selected model: ARIMA({p_best}, {d_best}, {q_best})  AIC = {best['AIC']:.2f}")

Selected model: ARIMA(2, 1, 3)  AIC = 1804.90
In [43]:
arima_model = ARIMA(train, order=(p_best, d_best, q_best)).fit()
print(arima_model.summary())

                               SARIMAX Results                                
==============================================================================
Dep. Variable:         Dublin_airport   No. Observations:                  157
Model:                 ARIMA(2, 1, 3)   Log Likelihood                -896.451
Date:                Mon, 11 May 2026   AIC                           1804.902
Time:                        13:52:04   BIC                           1823.201
Sample:                             0   HQIC                          1812.334
                                - 157                                         
Covariance Type:                  opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
ar.L1          1.7171      0.007    237.677      0.000       1.703       1.731
ar.L2         -0.9948      0.007   -134.327      0.000      -1.009      -0.980
ma.L1         -2.1676      0.070    -31.020      0.000      -2.305      -2.031
ma.L2          1.7917      0.127     14.161      0.000       1.544       2.040
ma.L3         -0.5445      0.073     -7.460      0.000      -0.688      -0.401
sigma2      5487.0562    682.822      8.036      0.000    4148.750    6825.363
===================================================================================
Ljung-Box (L1) (Q):                   0.09   Jarque-Bera (JB):                 0.42
Prob(Q):                              0.76   Prob(JB):                         0.81
Heteroskedasticity (H):               1.92   Skew:                            -0.12
Prob(H) (two-sided):                  0.02   Kurtosis:                         2.89
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
In [44]:
# Forecast the 12-month test period and compute performance metrics
arima_forecast  = arima_model.get_forecast(steps=12)
arima_pred_mean = arima_forecast.predicted_mean
arima_ci        = arima_forecast.conf_int(alpha=0.05)

mae_arima  = mean_absolute_error(test, arima_pred_mean)
rmse_arima = np.sqrt(mean_squared_error(test, arima_pred_mean))
mape_arima = mape(test, arima_pred_mean)
r2_arima   = r2_score(test, arima_pred_mean)

print(f"ARIMA({p_best},{d_best},{q_best}) Performance:")
print(f"  MAE  = {mae_arima:.1f}k")
print(f"  RMSE = {rmse_arima:.1f}k")
print(f"  MAPE = {mape_arima:.2f}%")
print(f"  R2   = {r2_arima:.4f}")

ARIMA(2,1,3) Performance:
  MAE  = 110.1k
  RMSE = 130.5k
  MAPE = 7.22%
  R2   = 0.7258
In [45]:
# Residual diagnostics - random residuals confirm the model captured the series structure
fig, axes = plt.subplots(1, 2, figsize=(14, 4))
residuals = arima_model.resid

axes[0].plot(residuals, color='steelblue')
axes[0].axhline(0, color='red', linestyle='--')
axes[0].set_title('ARIMA - Residuals Over Time')
axes[0].grid(True, alpha=0.3)

axes[1].hist(residuals, bins=25, color='steelblue', edgecolor='white')
axes[1].set_title('ARIMA - Residual Distribution')
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()
No description has been provided for this image
In [46]:
fig, ax = plt.subplots(figsize=(14, 6))
ax.plot(train.index[-36:], train[-36:], 'b-', label='Train', lw=1.5)
ax.plot(test.index, test,               'ko-', label='Actual', lw=2, ms=5)
ax.plot(test.index, arima_pred_mean,   'r--', label=f'ARIMA Forecast (MAPE={mape_arima:.1f}%)', lw=2)
ax.fill_between(test.index, arima_ci.iloc[:,0], arima_ci.iloc[:,1],
                color='salmon', alpha=0.25, label='95% Confidence Interval')
ax.set_xlabel('Date')
ax.set_ylabel('Passengers (thousands)')
ax.set_title(f'ARIMA({p_best},{d_best},{q_best}) - Forecast vs Actual')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
No description has been provided for this image
In [47]:
# ── Setup A: ARIMA trained on full series (with COVID) ───────────────────
ts_full_A = df_arr.set_index('Date')['Dublin_airport'].copy().asfreq('MS')
train_A   = ts_full_A[:-12]
test_A    = ts_full_A[-12:]

arima_A      = ARIMA(train_A, order=(p_best, d_best, q_best)).fit()
pred_A       = arima_A.get_forecast(steps=12).predicted_mean
ci_A         = arima_A.get_forecast(steps=12).conf_int(alpha=0.05)
mae_arima_A  = mean_absolute_error(test_A, pred_A)
rmse_arima_A = np.sqrt(mean_squared_error(test_A, pred_A))
mape_arima_A = mape(test_A, pred_A)
r2_arima_A   = r2_score(test_A, pred_A)

print("ARIMA Setup A (with COVID):")
print(f"  MAE  = {mae_arima_A:.1f}k  RMSE = {rmse_arima_A:.1f}k  MAPE = {mape_arima_A:.2f}%  R2 = {r2_arima_A:.4f}")

ARIMA Setup A (with COVID):
  MAE  = 94.7k  RMSE = 109.5k  MAPE = 6.18%  R2 = 0.8071

5. Model 2 - SARIMA¶

SARIMA extends ARIMA by adding seasonal components (P, D, Q) with period m=12 for monthly data.

5.1 Setup A – SARIMA With COVID¶

In [48]:
# ── Setup A: SARIMA trained on full series (with COVID) ─────────────────
sarima_A      = SARIMAX(train_A, order=(1,1,1), seasonal_order=(1,1,1,12)).fit(disp=False)
pred_sarima_A = sarima_A.get_forecast(steps=12).predicted_mean
ci_sarima_A   = sarima_A.get_forecast(steps=12).conf_int(alpha=0.05)
mae_sarima_A  = mean_absolute_error(test_A, pred_sarima_A)
rmse_sarima_A = np.sqrt(mean_squared_error(test_A, pred_sarima_A))
mape_sarima_A = mape(test_A, pred_sarima_A)
r2_sarima_A   = r2_score(test_A, pred_sarima_A)

print("SARIMA Setup A (with COVID):")
print(f"  MAE  = {mae_sarima_A:.1f}k  RMSE = {rmse_sarima_A:.1f}k  MAPE = {mape_sarima_A:.2f}%  R2 = {r2_sarima_A:.4f}")

SARIMA Setup A (with COVID):
  MAE  = 76.2k  RMSE = 89.9k  MAPE = 4.74%  R2 = 0.8698
In [49]:
fig, ax = plt.subplots(figsize=(14, 6))
ax.plot(train_A.index[-24:], train_A[-24:], 'b-', label='Train', lw=1.5)
ax.plot(test_A.index, test_A, 'ko-', label='Actual', lw=2, ms=5)
ax.plot(test_A.index, pred_sarima_A, 'm--', label=f'SARIMA Setup A (MAPE={mape_sarima_A:.1f}%)', lw=2)
ax.fill_between(test_A.index, ci_sarima_A.iloc[:,0], ci_sarima_A.iloc[:,1], color='plum', alpha=0.25, label='95% CI')
ax.axvspan(pd.Timestamp('2020-01-01'), pd.Timestamp('2022-02-01'),
           alpha=0.08, color='red', label='COVID Period')
ax.set_title('SARIMA(1,1,1)(1,1,1,12) Setup A (with COVID) – Forecast vs Actual')
ax.set_xlabel('Date'); ax.set_ylabel('Passengers (thousands)')
ax.legend(); ax.grid(True, alpha=0.3); plt.tight_layout(); plt.show()
No description has been provided for this image

5.2 Setup B – SARIMA COVID Excluded + Cross-Validation¶

Primary setup.

In [50]:
# SARIMA(1,1,1)(1,1,1,12)
# Non-seasonal: AR=1, differencing=1, MA=1
# Seasonal: AR=1, differencing=1, MA=1, period=12 (monthly)
sarima_model = SARIMAX(train,
                       order=(1, 1, 1),
                       seasonal_order=(1, 1, 1, 12)).fit(disp=False)
print(sarima_model.summary())

                                     SARIMAX Results                                      
==========================================================================================
Dep. Variable:                     Dublin_airport   No. Observations:                  157
Model:             SARIMAX(1, 1, 1)x(1, 1, 1, 12)   Log Likelihood                -810.607
Date:                            Mon, 11 May 2026   AIC                           1631.214
Time:                                    13:52:05   BIC                           1646.063
Sample:                                         0   HQIC                          1637.248
                                            - 157                                         
Covariance Type:                              opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
ar.L1         -0.3888      0.372     -1.044      0.297      -1.119       0.341
ma.L1          0.2344      0.403      0.582      0.560      -0.555       1.023
ar.S.L12       0.7962      0.478      1.665      0.096      -0.141       1.734
ma.S.L12      -0.9025      0.475     -1.902      0.057      -1.833       0.028
sigma2      4485.4529    366.474     12.239      0.000    3767.177    5203.728
===================================================================================
Ljung-Box (L1) (Q):                   0.01   Jarque-Bera (JB):               897.63
Prob(Q):                              0.94   Prob(JB):                         0.00
Heteroskedasticity (H):               4.48   Skew:                             1.09
Prob(H) (two-sided):                  0.00   Kurtosis:                        15.04
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
In [51]:
# Four-panel diagnostics plot to validate residual assumptions
sarima_model.plot_diagnostics(figsize=(14, 10))
plt.suptitle('SARIMA(1,1,1)(1,1,1,12) - Residual Diagnostics', y=1.01)
plt.tight_layout()
plt.show()
No description has been provided for this image
In [52]:
sarima_forecast  = sarima_model.get_forecast(steps=12)
sarima_pred_mean = sarima_forecast.predicted_mean
sarima_ci        = sarima_forecast.conf_int(alpha=0.05)

mae_sarima  = mean_absolute_error(test, sarima_pred_mean)
rmse_sarima = np.sqrt(mean_squared_error(test, sarima_pred_mean))
mape_sarima = mape(test, sarima_pred_mean)
r2_sarima   = r2_score(test, sarima_pred_mean)

print('SARIMA(1,1,1)(1,1,1,12) Performance:')
print(f"  MAE  = {mae_sarima:.1f}k")
print(f"  RMSE = {rmse_sarima:.1f}k")
print(f"  MAPE = {mape_sarima:.2f}%")
print(f"  R2   = {r2_sarima:.4f}")

SARIMA(1,1,1)(1,1,1,12) Performance:
  MAE  = 100.1k
  RMSE = 103.4k
  MAPE = 6.63%
  R2   = 0.8280
In [53]:
fig, ax = plt.subplots(figsize=(14, 6))
ax.plot(train.index[-36:], train[-36:], 'b-', label='Train', lw=1.5)
ax.plot(test.index, test,               'ko-', label='Actual', lw=2, ms=5)
ax.plot(test.index, sarima_pred_mean,  'm--', label=f'SARIMA Forecast (MAPE={mape_sarima:.1f}%)', lw=2)
ax.fill_between(test.index, sarima_ci.iloc[:,0], sarima_ci.iloc[:,1],
                color='plum', alpha=0.25, label='95% Confidence Interval')
ax.set_xlabel('Date')
ax.set_ylabel('Passengers (thousands)')
ax.set_title('SARIMA(1,1,1)(1,1,1,12) - Forecast vs Actual')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
No description has been provided for this image

5.1 Model 2 - SARIMA CROSS-VALIDATION AND FORECASTING¶

In [54]:
from sklearn.model_selection import TimeSeriesSplit

# Time Series Cross-Validation para SARIMA
# Usa expansão da janela de treino (expanding window) – adequado para séries temporais
# O dado mais antigo sempre fica no treino; apenas o horizonte de teste avança

tscv    = TimeSeriesSplit(n_splits=5)
cv_maes = []
cv_mapes= []
cv_r2s  = []

print("SARIMA – Time Series Cross-Validation (5 folds)\n")
print(f"{'Fold':<6} {'Train size':<12} {'Test size':<11} {'MAE':<10} {'MAPE%':<10} {'R²'}")
print("-" * 58)

for fold, (train_idx, test_idx) in enumerate(tscv.split(ts), 1):
    cv_train = ts.iloc[train_idx]
    cv_test  = ts.iloc[test_idx]

    try:
        cv_model = SARIMAX(cv_train,
                           order=(1, 1, 1),
                           seasonal_order=(1, 1, 1, 12)).fit(disp=False)
        cv_pred = cv_model.get_forecast(steps=len(cv_test)).predicted_mean

        fold_mae  = mean_absolute_error(cv_test, cv_pred)
        fold_mape = mape(cv_test, cv_pred)
        fold_r2   = r2_score(cv_test, cv_pred)

        cv_maes.append(fold_mae)
        cv_mapes.append(fold_mape)
        cv_r2s.append(fold_r2)

        print(f"{fold:<6} {len(cv_train):<12} {len(cv_test):<11} {fold_mae:<10.1f} {fold_mape:<10.2f} {fold_r2:.4f}")

    except Exception as e:
        print(f"{fold:<6} ERROR: {e}")

print("-" * 58)
(f"{'Mean':<6} {'':<12} {'':<11} {np.mean(cv_maes):<10.1f} {np.mean(cv_mapes):<10.2f} {np.mean(cv_r2s):.4f}")
print(f"{'Std':<6} {'':<12} {'':<11} {np.std(cv_maes):<10.1f} {np.std(cv_mapes):<10.2f} {np.std(cv_r2s):.4f}")

SARIMA – Time Series Cross-Validation (5 folds)

Fold   Train size   Test size   MAE        MAPE%      R²
----------------------------------------------------------
1      29           28          56.4       6.01       0.8726
2      57           28          97.5       9.22       0.7605
3      85           28          53.2       4.89       0.9122
4      113          28          172.0      14.20      0.0004
5      141          28          129.7      8.70       0.6560
----------------------------------------------------------
Std                             45.1       3.23       0.3323
In [55]:
fig, axes = plt.subplots(1, 3, figsize=(14, 4))

folds  = [f'Fold {i}' for i in range(1, len(cv_maes)+1)]

for ax, vals, ylabel, title in zip(
    axes,
    [cv_maes, cv_mapes, cv_r2s],
    ['MAE', 'MAPE (%)', 'R²'],
    ['MAE per Fold', 'MAPE% per Fold', 'R² per Fold']
):
    bars = ax.bar(folds, vals, color='steelblue')
    ax.axhline(np.mean(vals), color='red', linestyle='--', label=f'Mean={np.mean(vals):.2f}')
    ax.set_title(title)
    ax.set_ylabel(ylabel)
    ax.legend(fontsize=8)
    ax.grid(True, alpha=0.3, axis='y')
    for b, v in zip(bars, vals):
        ax.text(b.get_x()+b.get_width()/2, v*1.01, f'{v:.2f}', ha='center', fontsize=8)

plt.suptitle('SARIMA – 5-Fold Time Series Cross-Validation', fontweight='bold')
plt.tight_layout()
plt.show()
No description has been provided for this image
In [56]:
fig, ax = plt.subplots(figsize=(16, 6))

# Paleta de cores para cada fold
colors = ['#3498db', '#2ecc71', '#e74c3c', '#f39c12', '#9b59b6']

# Plota o histórico
ax.plot(ts.index, ts, 'k-', lw=1.5, label='Historical', alpha=0.6)

# Para cada fold, plota o período de treino e a previsão sobre o teste
for fold, (train_idx, test_idx) in enumerate(tscv.split(ts), 1):
    cv_train = ts.iloc[train_idx]
    cv_test  = ts.iloc[test_idx]

    cv_model = SARIMAX(cv_train,
                       order=(1, 1, 1),
                       seasonal_order=(1, 1, 1, 12)).fit(disp=False)

    cv_pred = cv_model.get_forecast(steps=len(cv_test)).predicted_mean
    cv_ci   = cv_model.get_forecast(steps=len(cv_test)).conf_int(alpha=0.05)

    # Previsão de cada fold
    ax.plot(cv_test.index, cv_pred,
            color=colors[fold-1], lw=2, linestyle='--',
            label=f'Fold {fold} Forecast (MAPE={cv_mapes[fold-1]:.1f}%)')

    # Banda de confiança 95%
    ax.fill_between(cv_test.index,
                    cv_ci.iloc[:, 0], cv_ci.iloc[:, 1],
                    color=colors[fold-1], alpha=0.08)

    # Linha vertical separando treino e teste de cada fold
    ax.axvline(x=cv_test.index[0], color=colors[fold-1],
               linestyle=':', alpha=0.5, lw=1)

ax.set_xlabel('Date')
ax.set_ylabel('Passengers (thousands)')
ax.set_title('SARIMA(1,1,1)(1,1,1,12) – Cross-Validation Forecasts per Fold')
ax.legend(loc='upper left', fontsize=9)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
No description has been provided for this image

6. Model 3 - Holt-Winters Exponential Smoothing (Triple ES)¶

Holt-Winters captures level, trend, and seasonality simultaneously. Both additive and multiplicative configurations are tested and compared by AIC.

Unlike ARIMA and SARIMA, Holt-Winters is tested in two distinct setups that differ not only in COVID inclusion but also in the training window scope, because the model's smoothing mechanism is heavily influenced by the recovery trajectory.

6.1 Setup A – Holt-Winters With COVID¶

Reference only. Full series including the pandemic structural break.

In [57]:
# ── Setup A: Holt-Winters on full series (with COVID) ───────────────────
hw_add_A = HWES(train_A, trend='add', seasonal='add', seasonal_periods=12).fit(optimized=True)
hw_mul_A = HWES(train_A, trend='mul', seasonal='mul', seasonal_periods=12).fit(optimized=True)

print(f"HW Additive   AIC (Setup A) = {hw_add_A.aic:.2f}")
print(f"HW Multiplicative AIC (Setup A) = {hw_mul_A.aic:.2f}")

hw_model_A = hw_add_A if hw_add_A.aic < hw_mul_A.aic else hw_mul_A
hw_type_A  = 'Additive' if hw_add_A.aic < hw_mul_A.aic else 'Multiplicative'
hw_pred_A  = hw_model_A.forecast(12)

mae_hw_A  = mean_absolute_error(test_A, hw_pred_A)
rmse_hw_A = np.sqrt(mean_squared_error(test_A, hw_pred_A))
mape_hw_A = mape(test_A, hw_pred_A)
r2_hw_A   = r2_score(test_A, hw_pred_A)

print(f"\nHolt-Winters {hw_type_A} Setup A (with COVID):")
print(f"  MAE  = {mae_hw_A:.1f}k  RMSE = {rmse_hw_A:.1f}k  MAPE = {mape_hw_A:.2f}%  R2 = {r2_hw_A:.4f}")

HW Additive   AIC (Setup A) = 1654.64
HW Multiplicative AIC (Setup A) = 1639.17

Holt-Winters Multiplicative Setup A (with COVID):
  MAE  = 54.7k  RMSE = 58.5k  MAPE = 3.74%  R2 = 0.9449
In [58]:
fig, ax = plt.subplots(figsize=(14, 6))
ax.plot(train_A.index[-24:], train_A[-24:], 'g-.', label='Train', lw=1.5)
ax.plot(test_A.index, test_A, 'ko-', label='Actual', lw=2, ms=5)
ax.plot(test_A.index, hw_pred_A, 'r--', label=f'HW {hw_type_A} Setup A (MAPE={mape_hw_A:.1f}%)', lw=2)
ax.axvspan(pd.Timestamp('2020-01-01'), pd.Timestamp('2022-02-01'),
           alpha=0.08, color='red', label='COVID Period')
ax.set_title(f'Holt-Winters {hw_type_A} Setup A (with COVID) – Forecast vs Actual')
ax.set_xlabel('Date'); ax.set_ylabel('Passengers (thousands)')
ax.legend(); ax.grid(True, alpha=0.3); plt.tight_layout(); plt.show()
No description has been provided for this image

6.2 Setup B – Holt-Winters Post-COVID Only (2022+)¶

Post-COVID data only is used because the steep recovery gradient between 2022 and 2024 causes the model to project unrealistic growth when trained on the full historical series.

In [59]:
# Treinar apenas com dados pós-COVID para o Holt-Winters
ts_post = ts[ts.index >= '2022-01-01']
train_hw = ts_post[:-12]
test_hw  = ts_post[-12:]

# Additive: suitable when seasonal variation is stable in absolute terms
hw_add = HWES(train_hw, trend='add', seasonal='add', seasonal_periods=12).fit(optimized=True)
print(f"Holt-Winters Additive       AIC = {hw_add.aic:.2f}")

# Multiplicative: suitable when seasonal amplitude grows proportionally with the level
hw_mul = HWES(train_hw, trend='mul', seasonal='mul', seasonal_periods=12).fit(optimized=True)
print(f"Holt-Winters Multiplicative AIC = {hw_mul.aic:.2f}")

# Triple Smoothing com alpha alto – peso elevado nas observações recentes
# Beta baixo evita amplificar o gradiente de recuperação pós-COVID como tendência permanente
hw_triple = HWES(train_hw, trend='add', seasonal='add', seasonal_periods=12).fit(
    smoothing_level=0.9,
    smoothing_trend=0.1,
    smoothing_seasonal=0.3,
    optimized=False
)
print(f"Holt-Winters Triple (α=0.9) AIC = {hw_triple.aic:.2f}")

Holt-Winters Additive       AIC = 280.29
Holt-Winters Multiplicative AIC = 279.69
Holt-Winters Triple (α=0.9) AIC = 325.80
In [60]:
# Select the configuration with the lower AIC
if hw_add.aic < hw_mul.aic:
    hw_model = hw_add
    hw_type  = "Additive"
else:
    hw_model = hw_mul
    hw_type  = "Multiplicative"

print(f"Selected: Holt-Winters {hw_type}")
print(hw_model.summary())

Selected: Holt-Winters Multiplicative
                       ExponentialSmoothing Model Results                       
================================================================================
Dep. Variable:           Dublin_airport   No. Observations:                   37
Model:             ExponentialSmoothing   SSE                          29892.300
Optimized:                         True   AIC                            279.694
Trend:                   Multiplicative   BIC                            305.469
Seasonal:                Multiplicative   AICC                           317.694
Seasonal Periods:                    12   Date:                 Mon, 11 May 2026
Box-Cox:                          False   Time:                         13:52:07
Box-Cox Coeff.:                    None                                         
=================================================================================
                          coeff                 code              optimized      
---------------------------------------------------------------------------------
smoothing_level               0.3843291                alpha                 True
smoothing_trend               0.3843291                 beta                 True
smoothing_seasonal             0.000000                gamma                 True
initial_level                 1202.1330                  l.0                 True
initial_trend                 1.0362919                  b.0                 True
initial_seasons.0             0.7682506                  s.0                 True
initial_seasons.1             0.8879332                  s.1                 True
initial_seasons.2             0.9525901                  s.2                 True
initial_seasons.3             0.9844377                  s.3                 True
initial_seasons.4             1.0650395                  s.4                 True
initial_seasons.5             1.0948763                  s.5                 True
initial_seasons.6             0.9499782                  s.6                 True
initial_seasons.7             0.8923226                  s.7                 True
initial_seasons.8             0.7071498                  s.8                 True
initial_seasons.9             0.6929817                  s.9                 True
initial_seasons.10            0.6557772                 s.10                 True
initial_seasons.11            0.6469945                 s.11                 True
---------------------------------------------------------------------------------
In [61]:
hw_pred = hw_model.forecast(steps=12)

mae_hw  = mean_absolute_error(test, hw_pred)
rmse_hw = np.sqrt(mean_squared_error(test, hw_pred))
mape_hw = mape(test, hw_pred)
r2_hw   = r2_score(test, hw_pred)

print(f"Holt-Winters ({hw_type}) Performance:")
print(f"  MAE  = {mae_hw:.1f}k")
print(f"  RMSE = {rmse_hw:.1f}k")
print(f"  MAPE = {mape_hw:.2f}%")
print(f"  R2   = {r2_hw:.4f}")

Holt-Winters (Multiplicative) Performance:
  MAE  = 91.3k
  RMSE = 93.1k
  MAPE = 6.14%
  R2   = 0.8606
In [62]:
fig, ax = plt.subplots(figsize=(14, 6))
ax.plot(train.index[-36:], train[-36:], 'g-.', label='Train', lw=1.5)
ax.plot(test.index, test,               'b-',  label='Actual Test', lw=2)
ax.plot(test.index, hw_pred,            'r--', label=f'HW {hw_type} Forecast (MAPE={mape_hw:.1f}%)', lw=2)
ax.set_xlabel('Date')
ax.set_ylabel('Passengers (thousands)')
ax.set_title(f'Holt-Winters {hw_type} - Forecast vs Actual')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
No description has been provided for this image
In [63]:
# Compare additive vs multiplicative to show selection reasoning
hw_add_pred = hw_add.forecast(12)
hw_mul_pred = hw_mul.forecast(12)
hw_triple_pred = hw_triple.forecast(12)


mae_hw_triple  = mean_absolute_error(test, hw_triple_pred)
rmse_hw_triple = np.sqrt(mean_squared_error(test, hw_triple_pred))
mape_hw_triple = mape(test, hw_triple_pred)
r2_hw_triple   = r2_score(test, hw_triple_pred)

fig, ax = plt.subplots(figsize=(14, 6))
ax.plot(train.index[-24:], train[-24:], 'g-.', label='Train')
ax.plot(test.index, test,               'b-',  label='Actual')
ax.plot(test.index, hw_add_pred, 'r--', label=f'HW Additive       (MAPE={mape(test, hw_add_pred):.1f}%)')
ax.plot(test.index, hw_mul_pred, 'm:',  label=f'HW Multiplicative (MAPE={mape(test, hw_mul_pred):.1f}%)')
ax.plot(test.index, hw_triple_pred, 'c-', label=f'HW Triple α=0.9 (MAPE={mape_hw_triple:.1f}%)')
ax.set_xlabel('Date')
ax.set_ylabel('Passengers (thousands)')
ax.set_title('Holt-Winters: Additive vs Multiplicative')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
No description has been provided for this image

7. Model Comparison – Setup A vs Setup B¶

All three models evaluated under both setups against the same 12-month test set (March 2025 – February 2026).

In [64]:
# ── Setup A metrics ──────────────────────────────────────────────────────
comp_A = pd.DataFrame({
    'Model' : ['ARIMA', 'SARIMA', 'Holt-Winters'],
    'MAE'   : [mae_arima_A,  mae_sarima_A,  mae_hw_A],
    'RMSE'  : [rmse_arima_A, rmse_sarima_A, rmse_hw_A],
    'MAPE%' : [mape_arima_A, mape_sarima_A, mape_hw_A],
    'R2'    : [r2_arima_A,   r2_sarima_A,   r2_hw_A],
    'Setup' : ['With COVID'] * 3
})

# ── Setup B metrics ──────────────────────────────────────────────────────
comp_B = pd.DataFrame({
    'Model' : ['ARIMA', 'SARIMA', 'Holt-Winters'],
    'MAE'   : [mae_arima,  mae_sarima,  mae_hw],
    'RMSE'  : [rmse_arima, rmse_sarima, rmse_hw],
    'MAPE%' : [mape_arima, mape_sarima, mape_hw],
    'R2'    : [r2_arima,   r2_sarima,   r2_hw],
    'Setup' : ['COVID Excluded'] * 3
})

comp_AB = pd.concat([comp_A, comp_B]).sort_values(['Model','Setup']).reset_index(drop=True)
print('Full Comparison – Setup A vs Setup B:')
display(comp_AB.round(4).style
    .highlight_min(color='green', axis=0, subset=['MAE','RMSE','MAPE%'])
    .highlight_max(color='green', axis=0, subset=['R2'])
    .format({'MAE':'{:.1f}','RMSE':'{:.1f}','MAPE%':'{:.2f}','R2':'{:.4f}'}))

Full Comparison – Setup A vs Setup B:
  Model MAE RMSE MAPE% R2 Setup
0 ARIMA 110.1 130.5 7.22 0.7258 COVID Excluded
1 ARIMA 94.7 109.5 6.18 0.8071 With COVID
2 Holt-Winters 91.3 93.1 6.14 0.8606 COVID Excluded
3 Holt-Winters 54.7 58.5 3.74 0.9449 With COVID
4 SARIMA 100.1 103.4 6.63 0.8280 COVID Excluded
5 SARIMA 76.2 89.9 4.74 0.8698 With COVID
In [65]:
# Grouped bar chart – MAPE Setup A vs Setup B
models   = ['ARIMA', 'SARIMA', 'Holt-Winters']
mapes_A  = [mape_arima_A, mape_sarima_A, mape_hw_A]
mapes_B  = [mape_arima,   mape_sarima,   mape_hw]

x, w = np.arange(3), 0.35
fig, ax = plt.subplots(figsize=(11, 5))
b1 = ax.bar(x-w/2, mapes_A, w, label='Setup A – with COVID',    color='#e74c3c', alpha=0.85)
b2 = ax.bar(x+w/2, mapes_B, w, label='Setup B – COVID excluded', color='#2ecc71', alpha=0.85)
for b in list(b1)+list(b2):
    ax.text(b.get_x()+b.get_width()/2, b.get_height()+0.2,
            f'{b.get_height():.1f}%', ha='center', fontsize=9)
ax.axhline(7.01, color='grey', linestyle=':', alpha=0.7, label='Naive Seasonal baseline (7.01%)')
ax.set_xticks(x); ax.set_xticklabels(models)
ax.set_ylabel('MAPE (%)'); ax.set_title('MAPE – Setup A vs Setup B')
ax.legend(); ax.grid(True, alpha=0.3, axis='y')
plt.tight_layout(); plt.show()
No description has been provided for this image
In [66]:
# Line chart – all models forecast vs actual (Setup B)
fig, ax = plt.subplots(figsize=(14, 7))
ax.plot(train.index[-24:], train[-24:], 'g-', alpha=0.5, label='Train', lw=1.5)
ax.plot(test.index, test, 'ko-', label='Actual', lw=2, ms=6)
ax.plot(test.index, arima_pred_mean,  'b--', label=f'ARIMA    (MAPE={mape_arima:.1f}%)', lw=1.8)
ax.plot(test.index, sarima_pred_mean, 'm-.', label=f'SARIMA   (MAPE={mape_sarima:.1f}%)', lw=1.8)
ax.plot(test.index, hw_pred,          'r:',  label=f'HW       (MAPE={mape_hw:.1f}%)', lw=2)
ax.set_xlabel('Date'); ax.set_ylabel('Passengers (thousands)')
ax.set_title('All Models Setup B – Forecast vs Actual (Test Set)')
ax.legend(loc='upper left'); ax.grid(True, alpha=0.3)
plt.tight_layout(); plt.show()
No description has been provided for this image
In [67]:
comparison = pd.DataFrame({
    'Model': ['ARIMA', 'SARIMA', 'Holt-Winters'],
    'MAE'  : [mae_arima,  mae_sarima,  mae_hw],
    'RMSE' : [rmse_arima, rmse_sarima, rmse_hw],
    'MAPE%': [mape_arima, mape_sarima, mape_hw],
    'R2'   : [r2_arima,   r2_sarima,   r2_hw]
}).set_index('Model').round(4)

display(comparison.style
    .highlight_min(color='green', axis=0, subset=['MAE','RMSE','MAPE%'])
    .highlight_max(color='green', axis=0, subset=['R2'])
    .format({'MAE':'{:.1f}','RMSE':'{:.1f}','MAPE%':'{:.2f}','R2':'{:.4f}'}))
  MAE RMSE MAPE% R2
Model        
ARIMA 110.1 130.5 7.22 0.7258
SARIMA 100.1 103.4 6.63 0.8280
Holt-Winters 91.3 93.1 6.14 0.8606
In [68]:
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
models = ['ARIMA', 'SARIMA', 'Holt-Winters']
colors = ['#3498db', '#9b59b6', '#2ecc71']

for ax, metric, vals, title in zip(
    axes,
    ['MAE', 'RMSE', 'MAPE (%)'],
    [[mae_arima, mae_sarima, mae_hw],
     [rmse_arima, rmse_sarima, rmse_hw],
     [mape_arima, mape_sarima, mape_hw]],
    ['MAE (lower = better)', 'RMSE (lower = better)', 'MAPE % (lower = better)']
):
    bars = ax.bar(models, vals, color=colors, alpha=0.85)
    ax.set_title(title, fontsize=10)
    ax.set_ylabel(metric)
    for b, v in zip(bars, vals):
        ax.text(b.get_x()+b.get_width()/2, v+max(vals)*0.01, f'{v:.1f}', ha='center', fontsize=9)
    ax.grid(True, alpha=0.3, axis='y')

plt.suptitle('Model Performance Comparison', fontsize=13, fontweight='bold')
plt.tight_layout()
plt.show()
No description has been provided for this image
In [69]:
# All models on one chart for direct visual comparison
fig, ax = plt.subplots(figsize=(14, 7))
ax.plot(train.index[-24:], train[-24:], 'g-', alpha=0.5, label='Train', lw=1.5)
ax.plot(test.index, test,               'ko-', label='Actual', lw=2, ms=6)
ax.plot(test.index, arima_pred_mean,   'b--', label=f'ARIMA    (MAPE={mape_arima:.1f}%)', lw=1.8)
ax.plot(test.index, sarima_pred_mean,  'm-.', label=f'SARIMA   (MAPE={mape_sarima:.1f}%)', lw=1.8)
ax.plot(test.index, hw_pred,           'r:',  label=f'Holt-Winters (MAPE={mape_hw:.1f}%)', lw=2)
ax.set_xlabel('Date')
ax.set_ylabel('Passengers (thousands)')
ax.set_title('All Models - Forecast vs Actual (Test Set)')
ax.legend(loc='upper left')
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
No description has been provided for this image

8. Final 12-Month Forecast¶

The best-performing model is retrained on the complete clean dataset to generate the 12-month passenger arrival forecast with 80% and 95% confidence intervals.

In [70]:
best_model_name = comparison['MAPE%'].idxmin()
print(f"Best model by MAPE: {best_model_name}")

Best model by MAPE: Holt-Winters
In [71]:
# Retrain all models on the full clean series (COVID excluded)
ts_clean_full = ts_full[(ts_full.index < '2020-01-01') | (ts_full.index > '2022-01-31')].copy()

# Treinar apenas com dados pós-COVID para o Holt-Winters
ts_post = ts[ts.index >= '2022-01-01']
train_hw = ts_post[:-12]
test_hw  = ts_post[-12:]

hw_add = HWES(train_hw, trend='add', seasonal='add', seasonal_periods=12).fit(optimized=True)
hw_mul = HWES(train_hw, trend='mul', seasonal='mul', seasonal_periods=12).fit(optimized=True)

final_arima  = ARIMA(ts_clean_full, order=(p_best, d_best, q_best)).fit()
final_sarima = SARIMAX(ts_clean_full, order=(1,1,1), seasonal_order=(1,1,1,12)).fit(disp=False)
final_hw     = HWES(ts_clean_full, trend='add', seasonal='add', seasonal_periods=12).fit(optimized=True)

future_dates = pd.date_range(
    start=ts_clean_full.index[-1] + pd.DateOffset(months=1), periods=12, freq='M')
print(f"Forecast period: {future_dates[0].date()} to {future_dates[-1].date()}")

Forecast period: 2026-04-30 to 2027-03-31
In [72]:
fc_arima  = final_arima.get_forecast(steps=12)
fc_sarima = final_sarima.get_forecast(steps=12)
fc_hw     = final_hw.forecast(steps=12)

fc_arima_mean  = fc_arima.predicted_mean
fc_arima_ci80  = fc_arima.conf_int(alpha=0.20)
fc_arima_ci95  = fc_arima.conf_int(alpha=0.05)

fc_sarima_mean = fc_sarima.predicted_mean
fc_sarima_ci95 = fc_sarima.conf_int(alpha=0.05)
In [73]:
# Forecast summary table
fc_table = pd.DataFrame({
    'Month'       : [d.strftime('%B %Y') for d in future_dates],
    'ARIMA'       : fc_arima_mean.values.round(0).astype(int),
    'SARIMA'      : fc_sarima_mean.values.round(0).astype(int),
    'Holt-Winters': fc_hw.values.round(0).astype(int)
}).set_index('Month')

display(fc_table)
ARIMA SARIMA Holt-Winters
Month
April 2026 1470 1609 1518
May 2026 1512 1747 1612
June 2026 1555 1808 1654
July 2026 1560 1957 1783
August 2026 1547 2015 1840
September 2026 1538 1768 1641
October 2026 1539 1690 1603
November 2026 1542 1386 1354
December 2026 1544 1373 1372
January 2027 1543 1298 1307
February 2027 1542 1256 1260
March 2027 1542 1496 1417
In [74]:
fig, ax = plt.subplots(figsize=(16, 7))

ax.plot(ts_clean_full.index[-36:], ts_clean_full[-36:], 'b-', label='Historical', lw=1.8)
ax.plot(future_dates, fc_arima_mean,  'b--', label='ARIMA Forecast', lw=2)
ax.fill_between(future_dates, fc_arima_ci80.iloc[:,0], fc_arima_ci80.iloc[:,1],
                color='steelblue', alpha=0.20, label='80% CI (ARIMA)')
ax.fill_between(future_dates, fc_arima_ci95.iloc[:,0], fc_arima_ci95.iloc[:,1],
                color='steelblue', alpha=0.08, label='95% CI (ARIMA)')
ax.plot(future_dates, fc_sarima_mean, 'm-.', label='SARIMA Forecast', lw=2)
ax.plot(future_dates, fc_hw,          'r:',  label='Holt-Winters Forecast', lw=2)

ax.axvline(x=ts_clean_full.index[-1], color='grey', linestyle=':', alpha=0.7)
ax.set_xlabel('Date')
ax.set_ylabel('Passengers (thousands)')
ax.set_title('12-Month Passenger Arrival Forecast - Dublin Airport')
ax.legend(loc='upper left')
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
No description has been provided for this image
In [75]:
print('=' * 50)
print('  12-Month Arrival Forecast - Dublin Airport')
print(f"  Best model: {best_model_name}")
print('=' * 50)

fc_vals  = fc_sarima_mean if best_model_name == 'SARIMA' else fc_arima_mean
fc_lo    = fc_sarima_ci95.iloc[:,0] if best_model_name == 'SARIMA' else fc_arima_ci95.iloc[:,0]
fc_hi    = fc_sarima_ci95.iloc[:,1] if best_model_name == 'SARIMA' else fc_arima_ci95.iloc[:,1]

for date, val, lo, hi in zip(future_dates, fc_vals.values, fc_lo.values, fc_hi.values):
    print(f"  {date.strftime('%b %Y')}: {val:7,.0f}k  (95% CI: {lo:,.0f}k - {hi:,.0f}k)")
print('=' * 50)

==================================================
  12-Month Arrival Forecast - Dublin Airport
  Best model: Holt-Winters
==================================================
  Apr 2026:   1,470k  (95% CI: 1,249k - 1,691k)
  May 2026:   1,512k  (95% CI: 1,168k - 1,856k)
  Jun 2026:   1,555k  (95% CI: 1,051k - 2,059k)
  Jul 2026:   1,560k  (95% CI: 871k - 2,250k)
  Aug 2026:   1,547k  (95% CI: 706k - 2,387k)
  Sep 2026:   1,538k  (95% CI: 586k - 2,490k)
  Oct 2026:   1,539k  (95% CI: 495k - 2,583k)
  Nov 2026:   1,542k  (95% CI: 413k - 2,672k)
  Dec 2026:   1,544k  (95% CI: 332k - 2,756k)
  Jan 2027:   1,543k  (95% CI: 254k - 2,833k)
  Feb 2027:   1,542k  (95% CI: 179k - 2,905k)
  Mar 2027:   1,542k  (95% CI: 110k - 2,974k)
==================================================

9. Model Interpretability¶

For forecasts to be useful in a business context, stakeholders must understand why the model predicts what it does. This section explains the structural components driving predictions and validates model assumptions.

In [77]:
# Decompose the clean series to show the components the models are learning from
log_ts_clean = np.log(ts_post)
log_ts_clean.interpolate(inplace=True)

decomp = seasonal_decompose(log_ts_clean, model='additive')
fig = decomp.plot()
fig.set_size_inches(16, 10)
plt.suptitle('Seasonal Decomposition - Clean Modelling Series (Log scale)', y=1.01)
plt.tight_layout()
plt.show()

print('Decomposition components:')
print('  Trend    : Steady growth 2010-2019; strong post-COVID recovery from 2022.')
print('  Seasonal : Consistent annual cycle - peak Jul/Aug, trough Jan/Feb.')
print('  Residuals: Random noise with no systematic pattern - model captures key structure.')
No description has been provided for this image 
Decomposition components:
  Trend    : Steady growth 2010-2019; strong post-COVID recovery from 2022.
  Seasonal : Consistent annual cycle - peak Jul/Aug, trough Jan/Feb.
  Residuals: Random noise with no systematic pattern - model captures key structure.
In [ ]:
In [ ]: