read thru this website carefully and learn from it https://nixtlaverse.nixtla.io/statsforecast/docs/tutorials/intermittentdata.html, see if I have implemented the code correctly from the website for my ensemble model from statsforecast import StatsForecast
from statsforecast.models import ADIDA, CrostonClassic, IMAPA, TSB
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 = [
    ADIDA(), 
    CrostonClassic(), 
    IMAPA(), 
    TSB(alpha_d = 0.2, alpha_p = 0.2)
]

# 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_filtered, h=2, step_size=1, n_windows=8, sort_df=True)

# Calculate the ensemble mean
ensemble = crossvalidation_df[['ADIDA', 'CrostonClassic', 'IMAPA', 'TSB']].mean(axis=1)

# Create a Series for the ensemble mean
ensemble_series = pl.Series('Ensemble', ensemble)

# Add the ensemble mean as a new column to the DataFrame
crossvalidation_df = crossvalidation_df.with_columns(ensemble_series)

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 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_fit_1)

# Instantiate the ConformalIntervals class
prediction_intervals = ConformalIntervals()

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

# Apply the non-negative constraint to the forecasts of individual models
forecasts_df = forecasts_df.with_columns([
    pl.when(pl.col('ADIDA') < 0).then(0).otherwise(pl.col('ADIDA')).alias('ADIDA'),
    pl.when(pl.col('CrostonClassic') < 0).then(0).otherwise(pl.col('CrostonClassic')).alias('CrostonClassic'),
    pl.when(pl.col('IMAPA') < 0).then(0).otherwise(pl.col('IMAPA')).alias('IMAPA'),
    pl.when(pl.col('TSB') < 0).then(0).otherwise(pl.col('TSB')).alias('TSB'),
])

# Calculate the ensemble forecast
ensemble_forecast = forecasts_df[['ADIDA', 'CrostonClassic', 'IMAPA', 'TSB']].mean(axis=1)

# Calculate the lower and upper prediction intervals for the ensemble forecast
ensemble_lo_95 = forecasts_df.select(
    [
        pl.when(pl.col('ADIDA-lo-95') < 0).then(0).otherwise(pl.col('ADIDA-lo-95')).alias('ADIDA-lo-95'),
        pl.when(pl.col('CrostonClassic-lo-95') < 0).then(0).otherwise(pl.col('CrostonClassic-lo-95')).alias('CrostonClassic-lo-95'),
        pl.when(pl.col('IMAPA-lo-95') < 0).then(0).otherwise(pl.col('IMAPA-lo-95')).alias('IMAPA-lo-95'),
        pl.when(pl.col('TSB-lo-95') < 0).then(0).otherwise(pl.col('TSB-lo-95')).alias('TSB-lo-95'),
   ]
).mean(axis=1)
ensemble_hi_95 = forecasts_df[['ADIDA-hi-95', 'CrostonClassic-hi-95', 'IMAPA-hi-95', 'TSB-hi-95']].mean(axis=1)

# Create Series for the ensemble forecast and its prediction intervals
ensemble_forecast_series = pl.Series('EnsembleForecast', ensemble_forecast)
ensemble_lo_95_series = pl.Series('Ensemble-lo-95', ensemble_lo_95)
ensemble_hi_95_series = pl.Series('Ensemble-hi-95', ensemble_hi_95)

# Add the ensemble forecast and its prediction intervals as new columns to the DataFrame
forecasts_df = forecasts_df.with_columns([ensemble_forecast_series, ensemble_lo_95_series, ensemble_hi_95_series])

# Round the ensemble forecast and prediction intervals and convert to integer
forecasts_df = forecasts_df.with_columns([
    pl.col("EnsembleForecast").round().cast(pl.Int32),
    pl.col("Ensemble-lo-95").round().cast(pl.Int32),
    pl.col("Ensemble-hi-95").round().cast(pl.Int32)
])

# Reorder the columns
forecasts_df = forecasts_df.select([
    "unique_id", 
    "ds", 
    "EnsembleForecast", 
    "Ensemble-lo-95", 
    "Ensemble-hi-95",
    "ADIDA", 
    "ADIDA-lo-95", 
    "ADIDA-hi-95", 
    "CrostonClassic", 
    "CrostonClassic-lo-95", 
    "CrostonClassic-hi-95", 
    "IMAPA", 
    "IMAPA-lo-95", 
    "IMAPA-hi-95",
    "TSB", 
    "TSB-lo-95", 
    "TSB-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)