<jupyter_start><jupyter_text>Introduction to STL ForecastingIn this notebook we present an example [STL forecasting](https://www.statsmodels.org/stable/generated/statsmodels.tsa.forecasting.stl.STLForecast.html) (Seasonal and Trend decomposition using Loess) with [statsmodels](https://www.statsmodels.org/). From the documentation:> *This is a model-based forecasting using STL to remove seasonality. Forecasts are produced by first subtracting the seasonality estimated using STL, then forecasting the deseasonalized data using a time-series model, for example, ARIMA.*For more details see [Forecasting: Principles and Practice, Section 3.6 STL Decomposition](https://otexts.com/fpp3/stl.html).  Prepare Notebook<jupyter_code># Main data packages. 
import numpy as np
import pandas as pd
# Data Viz. 
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_style(
    style='darkgrid', 
    rc={'axes.facecolor': '.9', 'grid.color': '.8'}
)
sns.set_palette(palette='deep')
sns_c = sns.color_palette(palette='deep')
%matplotlib inline
from pandas.plotting import register_matplotlib_converters
register_matplotlib_converters()

plt.rcParams['figure.figsize'] = [12, 6]
plt.rcParams['figure.dpi'] = 100<jupyter_output><empty_output><jupyter_text>Read DataWe will use the Basel temperature data set.<jupyter_code>from pathlib import Path

file_path = Path('../../data/basel_weather.csv')

raw_df = pd.read_csv(file_path)

data_df = raw_df.copy()<jupyter_output><empty_output><jupyter_text>EDA<jupyter_code>data_df = data_df.rename(columns={
    'Basel Temperature [2 m elevation corrected]': 'temperature', 
    'Basel Precipitation Total': 'precipitation', 
    'Basel Wind Speed [10 m]': 'wind_speed', 
    'Basel Wind Direction [10 m]': 'wind_direction'
    }
)

data_df = data_df.assign(
    timestamp = lambda x: pd.to_datetime(x['timestamp']), 
    date = lambda x: x['timestamp'].dt.date,
    year = lambda x: x['timestamp'].dt.year,
    month = lambda x: x['timestamp'].dt.month,
    day = lambda x: x['timestamp'].dt.day,
    dayofyear = lambda x: x['timestamp'].dt.dayofyear,
    hour = lambda x: x['timestamp'].dt.hour,
)

# Next, we aggregate by day. 
daily_data_df = data_df \
    .groupby(['date', 'year', 'month', 'day', 'dayofyear'], as_index=False)\
    .agg({'temperature': np.mean}) \
    .assign(date = lambda x: pd.to_datetime(x['date']))

daily_data_df.head()
# Plot daily data. 
fig, ax = plt.subplots()
sns.lineplot(x='date', y='temperature', data=daily_data_df, ax=ax)
ax.set(title='Basel Temperature (Daily)', ylabel=r'$^\circ$C');<jupyter_output><empty_output><jupyter_text>Train-Test Split<jupyter_code>train_test_cut_date = '2019-01-01'
df_train = daily_data_df.query(f'date < "{train_test_cut_date}"')
df_test = daily_data_df.query(f'date >= "{train_test_cut_date}"')

fig, ax = plt.subplots()
sns.lineplot(x='date', y='temperature', data=df_train, color=sns_c[0], label='y_train', ax=ax)
sns.lineplot(x='date', y='temperature', data=df_test, color=sns_c[1], label='y_test', ax=ax)
ax.axvline(x=train_test_cut_date, linestyle='--', color=sns_c[3], label='train-test-split')
ax.set(title='Basel Temperature (Daily)', ylabel=r'$^\circ$C');
y_train = df_train['temperature'].to_numpy()
t_train = df_train['date'].to_numpy()
y_test =df_test['temperature'].to_numpy()
t_test = df_test['date'].to_numpy()<jupyter_output><empty_output><jupyter_text>Model FitWe fit an exponential smoothing and an ARIMA model.<jupyter_code>from statsmodels.tsa.forecasting.stl import STLForecast
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.statespace.exponential_smoothing import ExponentialSmoothing

arima_model = STLForecast(
    endog=y_train,
    period=365, 
    model=ARIMA,
    model_kwargs={'order': (2, 1, 1)},
    trend_deg=1
)

es_model = STLForecast(
    endog=y_train,
    period=365, 
    model=ExponentialSmoothing,
    model_kwargs={'trend': True}, 
    trend_deg=1
)

arima_result = arima_model.fit()
es_result = es_model.fit()<jupyter_output><empty_output><jupyter_text>Generate Forecast<jupyter_code>y_pred_arima = arima_result.forecast(y_test.size)
y_pred_es = es_result.forecast(y_test.size)
from sklearn.metrics import mean_squared_error

mse_arima = mean_squared_error(y_true=y_test, y_pred=y_pred_arima)
mse_es = mean_squared_error(y_true=y_test, y_pred=y_pred_es)

fig, ax = plt.subplots()
sns.lineplot(x=t_test, y=y_test, color=sns_c[1], label='y_test', ax=ax)
sns.lineplot(x=t_test, y=y_pred_arima, color=sns_c[2], label=f'ARIMA pred (mse = {mse_arima: 0.3f})', ax=ax)
sns.lineplot(x=t_test, y=y_pred_es, color=sns_c[4], label=f'ES pred (mse = {mse_es: 0.3f})', ax=ax)
ax.set(title='Basel Temperature (Daily)', ylabel=r'$^\circ$C');<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14733276883996613, 'fraction_numerical': 0.009525825571549535, 'mean_word_length': 6.224770642201835, 'pattern_counts': {'":': 0, '<': 26, '<?xml version=': 0, '>': 27, 'https://': 3, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '475016', 'license': 'MIT', 'n_tokens_mistral': 1848, 'n_tokens_neox': 1774, 'n_words': 341, 'path': 'python/intro_forecasting/basel_stl_forecasting.ipynb', 'repo_name': 'datenzauberai/btsa', 'size': '4724'}