# Group by ‘unique_id’ and calculate the length of each group
lengths = y_cl4.groupby('unique_id').agg(pl.count().alias('length'))

# Count the occurrences of each length
counts = lengths.groupby('length').agg(pl.count().alias('count'))
counts = counts.sort('length')
 
pl.Config.set_tbl_rows(200)  


print(counts)  length ┆ count │
│ ---    ┆ ---   │
│ u32    ┆ u32   │
╞════════╪═══════╡
│ 1      ┆ 1942  │
│ 2      ┆ 357   │
│ 3      ┆ 157   │
│ 4      ┆ 107   │
│ 5      ┆ 74    │
│ 6      ┆ 40    │
│ 7      ┆ 48    │
│ 8      ┆ 37    │
│ 9      ┆ 39    │
│ 10     ┆ 47    │
│ 11     ┆ 54    │
│ 12     ┆ 36    │
│ 13     ┆ 35    │
│ 14     ┆ 43    │
│ 15     ┆ 47    │
│ 16     ┆ 45    │
│ 17     ┆ 36    │
│ 18     ┆ 37    │
│ 19     ┆ 51    │
│ 20     ┆ 35    │
│ 21     ┆ 41    │
│ 22     ┆ 29    │
│ 23     ┆ 26    │
│ 24     ┆ 33    │
│ 25     ┆ 35    │
│ 26     ┆ 41    │
│ 27     ┆ 39    │
│ 28     ┆ 34    │
│ 29     ┆ 37    │
│ 30     ┆ 31    │
│ 31     ┆ 32    │
│ 32     ┆ 26    │
│ 33     ┆ 30    │
│ 34     ┆ 22    │
│ 35     ┆ 39    │
│ 36     ┆ 32    │
│ 37     ┆ 32    │
│ 38     ┆ 33    │
│ 39     ┆ 37    │
│ 40     ┆ 34    │
│ 41     ┆ 24    │
│ 42     ┆ 22    │
│ 43     ┆ 17    │
│ 44     ┆ 18    │
│ 45     ┆ 13    │
│ 46     ┆ 10    │
│ 47     ┆ 18    │
│ 48     ┆ 15    │
│ 49     ┆ 17    │
│ 50     ┆ 12    │
│ 51     ┆ 15    │
│ 52     ┆ 10    │
│ 53     ┆ 11    │
│ 54     ┆ 6     │
│ 55     ┆ 9     │
│ 56     ┆ 7     │
│ 57     ┆ 11    │
│ 58     ┆ 11    │
│ 59     ┆ 9     │
│ 60     ┆ 13    │
│ 61     ┆ 14    │
│ 62     ┆ 7     │
│ 63     ┆ 5     │
│ 64     ┆ 3     │
│ 65     ┆ 6     │
│ 66     ┆ 6     │
│ 67     ┆ 5     │
│ 68     ┆ 11    │
│ 69     ┆ 7     │
│ 70     ┆ 4     │
│ 71     ┆ 2     │
│ 72     ┆ 4     │
│ 73     ┆ 3     │
│ 74     ┆ 1     │ 4276 series, now my model 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=12),
    AutoETS(season_length=12),
    AutoCES(season_length=12),
    DynamicOptimizedTheta(season_length=12)
]

# 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[['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_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('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'),
])

# Calculate the ensemble forecast
ensemble_forecast = forecasts_df[['AutoARIMA', 'AutoETS', 'CES', '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), I had to filtered out all series having less than 15 weeks for the cross validation to work, and to fit and do 24 month forecast, it is this dataset # For dataset with length > 8
lengths_1 = y_cl4.groupby('unique_id').agg(pl.count().alias('length'))
lengths_fit_1 = lengths_1.filter(pl.col('length') > 8)
y_cl4_fit_1 = y_cl4.join(
    lengths_fit_1.select(pl.col('unique_id')), 
    on='unique_id', 
    how='semi'  
) # Sort, group, and count the lengths of data sequences per unique_id
y_cl4_fit_1 = y_cl4_fit_1.sort("ds")
lengths_1 = y_cl4_fit_1.groupby('unique_id').agg(pl.count().alias('length'))
counts_1 = lengths_1.groupby('length').agg(pl.count().alias('count'))
counts_1 = counts_1.sort('length')
pl.Config.set_tbl_rows(200)  
print(counts_1) has to be at least 9 data points, now how do we fit the series with only 8 data points or less? especially series with only 1 data point? assume the ensemble model is already good satisfying the accuracy, bias requirements