read these wbsites and learn from them carefully, first one https://github.com/Nixtla/statsforecast/blob/main/experiments/m3/src/experiment.py, second one https://nixtlaverse.nixtla.io/statsforecast/docs/models/autoces.html, third one https://nixtlaverse.nixtla.io/statsforecast/docs/models/autoets.html, fourth one https://nixtlaverse.nixtla.io/statsforecast/docs/models/dynamicoptimizedtheta.html, is my ensemble model code correct? from statsforecast import StatsForecast
from statsforecast.models import AutoARIMA, AutoETS, AutoCES, DynamicOptimizedTheta
from statsforecast.utils import ConformalIntervals
import numpy as np
import polars as pl

# Polars option to display all rows
pl.Config.set_tbl_rows(None)

# Initialize the models
models = [
    AutoARIMA(season_length=52),
    AutoETS(season_length=52),
    AutoCES(season_length=52),
    DynamicOptimizedTheta(season_length=52)
]

# Initialize the StatsForecast model
sf = StatsForecast(models=models, freq='1w', n_jobs=-1)

# Perform cross-validation with a step size of 1 to mimic an expanding window
crossvalidation_df = sf.cross_validation(df=y_cl4, h=16, step_size=1, n_windows=18, sort_df=True)

# Create an ensemble out-of-sample predictions by taking the median of the predictions from all models
crossvalidation_df = crossvalidation_df.with_columns(
    pl.col(['AutoARIMA', 'AutoETS', 'CES', 'DynamicOp']).mean().alias("Ensemble")
)

def wmape(y_true, y_pred):
    return np.abs(y_true - y_pred).sum() / np.abs(y_true).sum()

# Calculate the WMAPE for the ensemble model
wmape_value = wmape(crossvalidation_df['y'], crossvalidation_df['Ensemble'])
print('Average WMAPE for Ensemble: ', round(wmape_value, 4))

# Calculate the errors for the ensemble model
errors = crossvalidation_df['y'] - crossvalidation_df['Ensemble']

# For an individual forecast
individual_accuracy = 1 - (abs(crossvalidation_df['y'] - crossvalidation_df['Ensemble']) / crossvalidation_df['y'])
individual_bias = (crossvalidation_df['Ensemble'] / crossvalidation_df['y']) - 1

# Add these calculations as new columns to your DataFrame
crossvalidation_df = crossvalidation_df.with_columns([
    individual_accuracy.alias("individual_accuracy"),
    individual_bias.alias("individual_bias")
])

# Print the individual accuracy and bias for each week
for row in crossvalidation_df.to_dicts():
    id = row['unique_id']
    date = row['ds'] 
    accuracy = row['individual_accuracy']
    bias = row['individual_bias']
    print(f"{id}, {date}, Individual Accuracy: {accuracy:.4f}, Individual Bias: {bias:.4f}")

# For groups of forecasts
group_accuracy = 1 - (errors.abs().sum() / crossvalidation_df['y'].sum())
group_bias = (crossvalidation_df['Ensemble'].sum() / crossvalidation_df['y'].sum()) - 1

# Print the average group accuracy and group bias over all folds for the ensemble model
print('Average Group Accuracy: ', round(group_accuracy, 4))
print('Average Group Bias: ', round(group_bias, 4))


# Fit the models on the entire dataset
sf.fit(y_cl4)

# Instantiate the ConformalIntervals class
prediction_intervals = ConformalIntervals()

# Generate 24 months forecasts
forecasts_df = sf.forecast(h=4*24, prediction_intervals=prediction_intervals, level=[95], id_col='unique_id', sort_df=True)

# Create an ensemble forecast by taking the mean of the forecasts from all models
forecasts_df = forecasts_df.with_columns(
    pl.col(['AutoARIMA', 'AutoETS', 'CES', 'DynamicOp']).mean().alias("Ensemble")
)

forecasts_df = forecasts_df.with_columns([
    pl.col("Ensemble").round().alias("Ensemble"),
    pl.col("Ensemble-lo-95").round().alias("Ensemble-lo-95"),
    pl.col("Ensemble-hi-95").round().alias("Ensemble-hi-95")
])

# Create an empty list
forecasts_list = []

# Append each row to the list
for row in forecasts_df.to_dicts():
    forecasts_list.append(row)

# Print the list
for forecast in forecasts_list:
    print(forecast)