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)

# Calculate the ensemble mean
ensemble = crossvalidation_df[['AutoARIMA', 'AutoETS', 'CES', 'DynamicOptimizedTheta']].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)

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

# Calculate the ensemble forecast
ensemble_forecast = forecasts_df.select(
    [
      pl.when(pl.col('AutoARIMA') < 0).then(0).otherwise(pl.col('AutoARIMA')).alias('AutoARIMA'),
      pl.when(pl.col('AutoETS') < 0).then(0).otherwise(pl.col('AutoETS')).alias('AutoETS'),
      pl.when(pl.col('CES') < 0).then(0).otherwise(pl.col('CES')).alias('CES'),
      pl.when(pl.col('DynamicOptimizedTheta') < 0).then(0).otherwise(pl.col('DynamicOptimizedTheta')).alias('DynamicOptimizedTheta'),
    ]
).mean(axis=1)

# Calculate the lower and upper prediction intervals for the ensemble forecast
ensemble_lo_95 = forecasts_df.select(
    [
        pl.when(pl.col('AutoARIMA-lo-95') < 0).then(0).otherwise(pl.col('AutoARIMA-lo-95')).alias('AutoARIMA-lo-95'),
        pl.when(pl.col('AutoETS-lo-95') < 0).then(0).otherwise(pl.col('AutoETS-lo-95')).alias('AutoETS-lo-95'),
        pl.when(pl.col('CES-lo-95') < 0).then(0).otherwise(pl.col('CES-lo-95')).alias('CES-lo-95'),
        pl.when(pl.col('DynamicOptimizedTheta-lo-95') < 0).then(0).otherwise(pl.col('DynamicOptimizedTheta-lo-95')).alias('DynamicOptimizedTheta-lo-95'),
   ]
).mean(axis=1)
ensemble_hi_95 = forecasts_df[['AutoARIMA-hi-95', 'AutoETS-hi-95', 'CES-hi-95', 'DynamicOptimizedTheta-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",
    "AutoARIMA", 
    "AutoARIMA-lo-95", 
    "AutoARIMA-hi-95", 
    "AutoETS", 
    "AutoETS-lo-95", 
    "AutoETS-hi-95", 
    "CES", 
    "CES-lo-95", 
    "CES-hi-95",
    "DynamicOptimizedTheta", 
    "DynamicOptimizedTheta-lo-95", 
    "DynamicOptimizedTheta-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)
 output {'unique_id': '6922332', 'ds': datetime.datetime(2025, 6, 16, 0, 0), 'EnsembleForecast': 19444, 'Ensemble-lo-95': 16762, 'Ensemble-hi-95': 22125, 'AutoARIMA': 21242.763671875, 'AutoARIMA-lo-95': 19512.095703125, 'AutoARIMA-hi-95': 22973.431640625, 'AutoETS': 21279.15625, 'AutoETS-lo-95': 19626.21484375, 'AutoETS-hi-95': 22932.09765625, 'CES': 14011.3857421875, 'CES-lo-95': 8402.5244140625, 'CES-hi-95': 19620.248046875, 'DynamicOptimizedTheta': 21241.310546875, 'DynamicOptimizedTheta-lo-95': 19507.21875, 'DynamicOptimizedTheta-hi-95': 22975.40234375}
{'unique_id': '6922332', 'ds': datetime.datetime(2025, 6, 23, 0, 0), 'EnsembleForecast': 23127, 'Ensemble-lo-95': 20446, 'Ensemble-hi-95': 25809, 'AutoARIMA': 21242.763671875, 'AutoARIMA-lo-95': 19512.095703125, 'AutoARIMA-hi-95': 22973.431640625, 'AutoETS': 21279.15625, 'AutoETS-lo-95': 19626.21484375, 'AutoETS-hi-95': 22932.09765625, 'CES': 28746.3984375, 'CES-lo-95': 23137.537109375, 'CES-hi-95': 34355.26171875, 'DynamicOptimizedTheta': 21241.310546875, 'DynamicOptimizedTheta-lo-95': 19507.21875, 'DynamicOptimizedTheta-hi-95': 22975.40234375}
{'unique_id': '6922332', 'ds': datetime.datetime(2025, 6, 30, 0, 0), 'EnsembleForecast': 22671, 'Ensemble-lo-95': 19990, 'Ensemble-hi-95': 25353, 'AutoARIMA': 21242.763671875, 'AutoARIMA-lo-95': 19512.095703125, 'AutoARIMA-hi-95': 22973.431640625, 'AutoETS': 21279.15625, 'AutoETS-lo-95': 19626.21484375, 'AutoETS-hi-95': 22932.09765625, 'CES': 26922.3515625, 'CES-lo-95': 21313.490234375, 'CES-hi-95': 32531.21484375, 'DynamicOptimizedTheta': 21241.310546875, 'DynamicOptimizedTheta-lo-95': 19507.21875, 'DynamicOptimizedTheta-hi-95': 22975.40234375}
{'unique_id': '6922332', 'ds': datetime.datetime(2025, 7, 7, 0, 0), 'EnsembleForecast': 19052, 'Ensemble-lo-95': 16370, 'Ensemble-hi-95': 21733, 'AutoARIMA': 21242.763671875, 'AutoARIMA-lo-95': 19512.095703125, 'AutoARIMA-hi-95': 22973.431640625, 'AutoETS': 21279.15625, 'AutoETS-lo-95': 19626.21484375, 'AutoETS-hi-95': 22932.09765625, 'CES': 12443.15234375, 'CES-lo-95': 6834.291015625, 'CES-hi-95': 18052.013671875, 'DynamicOptimizedTheta': 21241.310546875, 'DynamicOptimizedTheta-lo-95': 19507.21875, 'DynamicOptimizedTheta-hi-95': 22975.40234375}
{'unique_id': '6922332', 'ds': datetime.datetime(2025, 7, 14, 0, 0), 'EnsembleForecast': 21149, 'Ensemble-lo-95': 18467, 'Ensemble-hi-95': 23831, 'AutoARIMA': 21242.763671875, 'AutoARIMA-lo-95': 19512.095703125, 'AutoARIMA-hi-95': 22973.431640625, 'AutoETS': 21279.15625, 'AutoETS-lo-95': 19626.21484375, 'AutoETS-hi-95': 22932.09765625, 'CES': 20832.93359375, 'CES-lo-95': 15224.0712890625, 'CES-hi-95': 26441.794921875, 'DynamicOptimizedTheta': 21241.310546875, 'DynamicOptimizedTheta-lo-95': 19507.21875, 'DynamicOptimizedTheta-hi-95': 22975.40234375}
{'unique_id': '6922332', 'ds': datetime.datetime(2025, 7, 21, 0, 0), 'EnsembleForecast': 20599, 'Ensemble-lo-95': 17918, 'Ensemble-hi-95': 23281, 'AutoARIMA': 21242.763671875, 'AutoARIMA-lo-95': 19512.095703125, 'AutoARIMA-hi-95': 22973.431640625, 'AutoETS': 21279.15625, 'AutoETS-lo-95': 19626.21484375, 'AutoETS-hi-95': 22932.09765625, 'CES': 18633.75390625, 'CES-lo-95': 13024.892578125, 'CES-hi-95': 24242.615234375, 'DynamicOptimizedTheta': 21241.310546875, 'DynamicOptimizedTheta-lo-95': 19507.21875, 'DynamicOptimizedTheta-hi-95': 22975.40234375}
{'unique_id': '6922332', 'ds': datetime.datetime(2025, 7, 28, 0, 0), 'EnsembleForecast': 20781, 'Ensemble-lo-95': 18099, 'Ensemble-hi-95': 23462, 'AutoARIMA': 21242.763671875, 'AutoARIMA-lo-95': 19512.095703125, 'AutoARIMA-hi-95': 22973.431640625, 'AutoETS': 21279.15625, 'AutoETS-lo-95': 19626.21484375, 'AutoETS-hi-95': 22932.09765625, 'CES': 19359.54296875, 'CES-lo-95': 13750.6826171875, 'CES-hi-95': 24968.40625, 'DynamicOptimizedTheta': 21241.310546875, 'DynamicOptimizedTheta-lo-95': 19507.21875, 'DynamicOptimizedTheta-hi-95': 22975.40234375}
{'unique_id': '6922332', 'ds': datetime.datetime(2025, 8, 4, 0, 0), 'EnsembleForecast': 18442, 'Ensemble-lo-95': 15761, 'Ensemble-hi-95': 21124, 'AutoARIMA': 21242.763671875, 'AutoARIMA-lo-95': 19512.095703125, 'AutoARIMA-hi-95': 22973.431640625, 'AutoETS': 21279.15625, 'AutoETS-lo-95': 19626.21484375, 'AutoETS-hi-95': 22932.09765625, 'CES': 10005.392578125, 'CES-lo-95': 4396.53125, 'CES-hi-95': 15614.2548828125, 'DynamicOptimizedTheta': 21241.310546875, 'DynamicOptimizedTheta-lo-95': 19507.21875, 'DynamicOptimizedTheta-hi-95': 22975.40234375}
{'unique_id': '6922332', 'ds': datetime.datetime(2025, 8, 11, 0, 0), 'EnsembleForecast': 20543, 'Ensemble-lo-95': 17862, 'Ensemble-hi-95': 23225, 'AutoARIMA': 21242.763671875, 'AutoARIMA-lo-95': 19512.095703125, 'AutoARIMA-hi-95': 22973.431640625, 'AutoETS': 21279.15625, 'AutoETS-lo-95': 19626.21484375, 'AutoETS-hi-95': 22932.09765625, 'CES': 18410.65625, 'CES-lo-95': 12801.794921875, 'CES-hi-95': 24019.517578125, 'DynamicOptimizedTheta': 21241.310546875, 'DynamicOptimizedTheta-lo-95': 19507.21875, 'DynamicOptimizedTheta-hi-95': 22975.40234375}
{'unique_id': '6922332', 'ds': datetime.datetime(2025, 8, 18, 0, 0), 'EnsembleForecast': 19830, 'Ensemble-lo-95': 17148, 'Ensemble-hi-95': 22511, 'AutoARIMA': 21242.763671875, 'AutoARIMA-lo-95': 19512.095703125, 'AutoARIMA-hi-95': 22973.431640625, 'AutoETS': 21279.15625, 'AutoETS-lo-95': 19626.21484375, 'AutoETS-hi-95': 22932.09765625, 'CES': 15556.056640625, 'CES-lo-95': 9947.1953125, 'CES-hi-95': 21164.91796875, 'DynamicOptimizedTheta': 21241.310546875, 'DynamicOptimizedTheta-lo-95': 19507.21875, 'DynamicOptimizedTheta-hi-95': 22975.40234375}
{'unique_id': '6922332', 'ds': datetime.datetime(2025, 8, 25, 0, 0), 'EnsembleForecast': 19568, 'Ensemble-lo-95': 16886, 'Ensemble-hi-95': 22249, 'AutoARIMA': 21242.763671875, 'AutoARIMA-lo-95': 19512.095703125, 'AutoARIMA-hi-95': 22973.431640625, 'AutoETS': 21279.15625, 'AutoETS-lo-95': 19626.21484375, 'AutoETS-hi-95': 22932.09765625, 'CES': 14507.791015625, 'CES-lo-95': 8898.9296875, 'CES-hi-95': 20116.65234375, 'DynamicOptimizedTheta': 21241.310546875, 'DynamicOptimizedTheta-lo-95': 19507.21875, 'DynamicOptimizedTheta-hi-95': 22975.40234375}
{'unique_id': '6922332', 'ds': datetime.datetime(2025, 9, 1, 0, 0), 'EnsembleForecast': 18793, 'Ensemble-lo-95': 16111, 'Ensemble-hi-95': 21474, 'AutoARIMA': 21242.763671875, 'AutoARIMA-lo-95': 19512.095703125, 'AutoARIMA-hi-95': 22973.431640625, 'AutoETS': 21279.15625, 'AutoETS-lo-95': 19626.21484375, 'AutoETS-hi-95': 22932.09765625, 'CES': 11407.0390625, 'CES-lo-95': 5798.177734375, 'CES-hi-95': 17015.900390625, 'DynamicOptimizedTheta': 21241.310546875, 'DynamicOptimizedTheta-lo-95': 19507.21875, 'DynamicOptimizedTheta-hi-95': 22975.40234375}
{'unique_id': '6922332', 'ds': datetime.datetime(2025, 9, 8, 0, 0), 'EnsembleForecast': 17686, 'Ensemble-lo-95': 15004, 'Ensemble-hi-95': 20367, 'AutoARIMA': 21242.763671875, 'AutoARIMA-lo-95': 19512.095703125, 'AutoARIMA-hi-95': 22973.431640625, 'AutoETS': 21279.15625, 'AutoETS-lo-95': 19626.21484375, 'AutoETS-hi-95': 22932.09765625, 'CES': 6979.3349609375, 'CES-lo-95': 1370.473388671875, 'CES-hi-95': 12588.1962890625, 'DynamicOptimizedTheta': 21241.310546875, 'DynamicOptimizedTheta-lo-95': 19507.21875, 'DynamicOptimizedTheta-hi-95': 22975.40234375}
{'unique_id': '6922332', 'ds': datetime.datetime(2025, 9, 15, 0, 0), 'EnsembleForecast': 21029, 'Ensemble-lo-95': 18347, 'Ensemble-hi-95': 23711, 'AutoARIMA': 21242.763671875, 'AutoARIMA-lo-95': 19512.095703125, 'AutoARIMA-hi-95': 22973.431640625, 'AutoETS': 21279.15625, 'AutoETS-lo-95': 19626.21484375, 'AutoETS-hi-95': 22932.09765625, 'CES': 20352.7578125, 'CES-lo-95': 14743.896484375, 'CES-hi-95': 25961.62109375, 'DynamicOptimizedTheta': 21241.310546875, 'DynamicOptimizedTheta-lo-95': 19507.21875, 'DynamicOptimizedTheta-hi-95': 22975.40234375}
{'unique_id': '6922332', 'ds': datetime.datetime(2025, 9, 22, 0, 0), 'EnsembleForecast': 20262, 'Ensemble-lo-95': 17580, 'Ensemble-hi-95': 22944, 'AutoARIMA': 21242.763671875, 'AutoARIMA-lo-95': 19512.095703125, 'AutoARIMA-hi-95': 22973.431640625, 'AutoETS': 21279.15625, 'AutoETS-lo-95': 19626.21484375, 'AutoETS-hi-95': 22932.09765625, 'CES': 17284.330078125, 'CES-lo-95': 11675.4677734375, 'CES-hi-95': 22893.19140625, 'DynamicOptimizedTheta': 21241.310546875, 'DynamicOptimizedTheta-lo-95': 19507.21875, 'DynamicOptimizedTheta-hi-95': 22975.40234375}
{'unique_id': '6922332', 'ds': datetime.datetime(2025, 9, 29, 0, 0), 'EnsembleForecast': 20816, 'Ensemble-lo-95': 18134, 'Ensemble-hi-95': 23498, 'AutoARIMA': 21242.763671875, 'AutoARIMA-lo-95': 19512.095703125, 'AutoARIMA-hi-95': 22973.431640625, 'AutoETS': 21279.15625, 'AutoETS-lo-95': 19626.21484375, 'AutoETS-hi-95': 22932.09765625, 'CES': 19500.451171875, 'CES-lo-95': 13891.58984375, 'CES-hi-95': 25109.3125, 'DynamicOptimizedTheta': 21241.310546875, 'DynamicOptimizedTheta-lo-95': 19507.21875, 'DynamicOptimizedTheta-hi-95': 22975.40234375} unique_count = y_cl4['unique_id'].unique().shape[0]
print('Unique count of unique_id:', unique_count)Unique count of unique_id: 113 y_cl4.head() unique_id	ds	y
str	datetime[μs]	f64
"6665053"	2020-11-30 00:00:00	5484.0
"5338737"	2020-11-30 00:00:00	52310.0
"5339184"	2020-11-30 00:00:00	97457.0
"1131041"	2020-11-30 00:00:00	15202.588
"1131030"	2020-11-30 00:00:00	17838.0 based on all the information I have given you, is there a way to see the y in y_cl4 dataframe with the ensemble forecasted values? Normally if I just do 1 series, I could create a plot and see if the forecasted values would make sense given historical data, but there are 113 series here