import pandas as pd
import numpy as np
from sktime.forecasting.model_selection import (
ForecastingOptunaSearchCV,
ExpandingWindowSplitter,
temporal_train_test_split
)
from sktime.forecasting.base import ForecastingHorizon
from sktime.performance_metrics.forecasting import MeanAbsolutePercentageError
from sktime.forecasting.compose import TransformedTargetForecaster
from sktime.forecasting.statsforecast import (
StatsForecastMSTL,
StatsForecastAutoETS,
StatsForecastAutoARIMA,
StatsForecastAutoTheta
)
from sktime.transformations.series.detrend import Detrender
from sktime.transformations.series.deseasonalize import Deseasonalizer
import optuna
import warnings
warnings.filterwarnings('ignore')
# Load your time series data
# Ensure 'pivot_table' is defined and contains the 'PAN4_PIBPMG4' series
y = pivot_table['PAN4_PIBPMG4']
# Split the data into train and test sets
y_train, y_test = temporal_train_test_split(y, test_size=8)
# Define the forecasting horizon
fh = ForecastingHorizon(np.arange(1, 9), is_relative=True)
# Set up cross-validation with an expanding window splitter
cv = ExpandingWindowSplitter(fh=fh, initial_window=len(y_train) - 8)
# Define the parameter space for tuning
param_distributions = {
'forecaster**season_length': optuna.distributions.CategoricalDistribution([(4,), (8,)]),
'forecaster**trend_forecaster': optuna.distributions.CategoricalDistribution([
StatsForecastAutoETS(model="ZZZ"),
StatsForecastAutoARIMA(seasonal=True),
StatsForecastAutoTheta()
]),
'forecaster\_\_stl_kwargs': {
'robust': optuna.distributions.CategoricalDistribution([True, False]),
'period': optuna.distributions.IntUniformDistribution(4, 8)
}
}
# Initialize the MSTL forecaster
mstl_forecaster = StatsForecastMSTL()
# Create a pipeline with optional transformations
forecaster = TransformedTargetForecaster(steps=[
("detrender", Detrender()),
("deseasonalizer", Deseasonalizer()),
("mstl_forecaster", mstl_forecaster)
])
# Set up the OptunaSearchCV
optuna_search = ForecastingOptunaSearchCV(
forecaster=forecaster,
cv=cv,
param_distributions=param_distributions,
scoring=MeanAbsolutePercentageError(symmetric=True),
n_trials=100,
random_state=42
)
# Fit the model
optuna_search.fit(y_train)
# Predict
y_pred = optuna_search.predict(fh)
# Evaluate
mape = MeanAbsolutePercentageError(symmetric=True)
final_mape = mape(y_test, y_pred)
print(f"Final sMAPE: {final_mape:.2f}")
# Plot results
import matplotlib.pyplot as plt
plt.figure(figsize=(15, 7))
plt.plot(y_train.index, y_train.values, label='Training Data', color='blue')
plt.plot(y_test.index, y_test.values, label='Test Data', color='green')
plt.plot(y_pred.index, y_pred.values, label='Predictions', color='red', linestyle='--')
plt.title('MSTL Forecast Results with Optuna Optimization')
plt.legend()
plt.grid(True)
plt.show()
# Save the best model
from joblib import dump
dump(optuna*search.best_forecaster*, 'best_mstl_model_optuna.joblib')
print("\nBest model saved as 'best_mstl_model_optuna.joblib'")
# Print additional optimization results
print("\nOptimization Results:")
print("="\*50)
print(f"Number of completed trials: {len(optuna*search.cv_results*)}")
print(f"Best trial number: {optuna*search.best_index*}")
print(f"Best sMAPE achieved during optimization: {optuna*search.best_score*:.2f}")
# Print best parameters
print("\nBest Parameters Found:")
print("="\*50)
for param, value in optuna*search.best_params*.items():
print(f"{param}: {value}")
