my model some of the individual accuracy and group accuracy output 90084053_US01_8_6067465, 2023-08-28 00:00:00, Individual Accuracy: 0.9278, Individual Bias: -0.0722 90084053_US01_8_6067465, 2023-09-18 00:00:00, Individual Accuracy: 0.7643, Individual Bias: 0.2357 90084053_US01_8_6067465, 2023-09-25 00:00:00, Individual Accuracy: 0.8620, Individual Bias: -0.1380 90084053_US01_8_6067465, 2023-10-16 00:00:00, Individual Accuracy: 0.8762, Individual Bias: -0.1238 90084053_US01_8_6067465, 2023-09-18 00:00:00, Individual Accuracy: 0.8725, Individual Bias: 0.1275 90084053_US01_8_6067465, 2023-09-25 00:00:00, Individual Accuracy: 0.8697, Individual Bias: -0.1303 90084053_US01_8_6067465, 2023-10-16 00:00:00, Individual Accuracy: 0.9412, Individual Bias: -0.0588 90084053_US01_8_6067465, 2023-11-06 00:00:00, Individual Accuracy: 0.8918, Individual Bias: 0.1082 Average Group Accuracy: 0.6598 Average Group Bias: -0.0194. screenshot of the excel file shows an quick example of how to do lag forecast in supply chain. Now I need to do 12 weeks out for lag forecast starting first week of 2023. would I do this given my model and output already generated. I also have the forecast output from the model as well. basically we have to forecast 12 weeks out from first week 2023 cutoff as lag 0 then 12 weeks out forecast form 2nd week 2023 cutoff as lag 1. am I getting this right from the excel file?, this is my mdoel code from statsforecast import StatsForecast
from statsforecast.models import AutoARIMA
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)]

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

# Perform cross-validation with a step size of 1 to do an expanding window
crossvalidation_df = sf.cross_validation(df=y_cl4_over_10, h=4, step_size=1, n_windows=5, sort_df=True)

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

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

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

# For an individual forecast
individual_accuracy = 1 - (abs(crossvalidation_df['y'] - crossvalidation_df['AutoARIMA']) / crossvalidation_df['y'])
individual_bias = (crossvalidation_df['AutoARIMA'] / 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}")

# Filter individual accuracy and individual bias to be within absolute value of 15
filtered_df = crossvalidation_df.filter(
    (crossvalidation_df["individual_accuracy"].abs() <= 5) 
    & (crossvalidation_df["individual_bias"].abs() <= 5)
)

# Calculate errors for the filtered DataFrame
errors_filtered = filtered_df['y'] - filtered_df['AutoARIMA']

# For groups of forecasts
group_accuracy = 1 - (errors_filtered.abs().sum() / filtered_df['y'].sum())
group_bias = (filtered_df['AutoARIMA'].sum() / filtered_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))


# Instantiate the ConformalIntervals class
prediction_intervals = ConformalIntervals()

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

# Apply the non-negative constraint to the forecasts of the model
forecasts_df = forecasts_df.with_columns([
    pl.when(pl.col('AutoARIMA') < 0).then(0).otherwise(pl.col('AutoARIMA')).alias('AutoARIMA'),
    pl.when(pl.col('AutoARIMA-lo-95') < 0).then(0).otherwise(pl.col('AutoARIMA-lo-95')).alias('AutoARIMA-lo-95')
])

forecasts_df = forecasts_df.with_columns([
    pl.col("AutoARIMA").round().cast(pl.Int32),
    pl.col("AutoARIMA-lo-95").round().cast(pl.Int32),
    pl.col("AutoARIMA-hi-95").round().cast(pl.Int32)
])

# Split the unique_id concat into the original columns
def split_unique_id(unique_id):
    parts = unique_id.split('_')
    return parts if len(parts) >= 4 else (parts + [None] * (4 - len(parts)))

forecasts_df = (
    forecasts_df
    .with_columns([
        pl.col('unique_id').apply(lambda uid: split_unique_id(uid)[0]).alias('MaterialID'),
        pl.col('unique_id').apply(lambda uid: split_unique_id(uid)[1]).alias('SalesOrg'),
        pl.col('unique_id').apply(lambda uid: split_unique_id(uid)[2]).alias('DistrChan'),
        pl.col('unique_id').apply(lambda uid: split_unique_id(uid)[3]).alias('CL4'),
    ])
    .drop('unique_id')
)

# Rename ‘ds’ to ‘WeekDate’
forecasts_df = forecasts_df.rename({'ds': 'WeekDate'})

# Reorder the columns
forecasts_df = forecasts_df.select([
    "MaterialID",
    "SalesOrg",
    "DistrChan", 
    "CL4",
    "WeekDate",
    "AutoARIMA", 
    "AutoARIMA-lo-95", 
    "AutoARIMA-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), forecast output {'MaterialID': '90075618', 'SalesOrg': 'US01', 'DistrChan': '8', 'CL4': '6067435', 'WeekDate': datetime.datetime(2024, 3, 18, 0, 0), 'AutoARIMA': 4, 'AutoARIMA-lo-95': 2, 'AutoARIMA-hi-95': 6}
{'MaterialID': '90075618', 'SalesOrg': 'US01', 'DistrChan': '8', 'CL4': '6067435', 'WeekDate': datetime.datetime(2024, 3, 25, 0, 0), 'AutoARIMA': 4, 'AutoARIMA-lo-95': 2, 'AutoARIMA-hi-95': 6}
{'MaterialID': '90075618', 'SalesOrg': 'US01', 'DistrChan': '8', 'CL4': '6067435', 'WeekDate': datetime.datetime(2024, 4, 1, 0, 0), 'AutoARIMA': 4, 'AutoARIMA-lo-95': 2, 'AutoARIMA-hi-95': 6}
{'MaterialID': '90075618', 'SalesOrg': 'US01', 'DistrChan': '8', 'CL4': '6067435', 'WeekDate': datetime.datetime(2024, 4, 8, 0, 0), 'AutoARIMA': 3, 'AutoARIMA-lo-95': 1, 'AutoARIMA-hi-95': 5}
{'MaterialID': '90075618', 'SalesOrg': 'US01', 'DistrChan': '8', 'CL4': '6067435', 'WeekDate': datetime.datetime(2024, 4, 15, 0, 0), 'AutoARIMA': 4, 'AutoARIMA-lo-95': 2, 'AutoARIMA-hi-95': 6}
{'MaterialID': '90075618', 'SalesOrg': 'US01', 'DistrChan': '8', 'CL4': '6067435', 'WeekDate': datetime.datetime(2024, 4, 22, 0, 0), 'AutoARIMA': 4, 'AutoARIMA-lo-95': 2, 'AutoARIMA-hi-95': 6} , my code attempt to solve this problem from datetime import datetime, timedelta
from collections import defaultdict

# Function to extract forecasts and actuals for a given lag and cutoff date
def extract_data_for_lag(lag, cutoff_date):
    data = defaultdict(lambda: {"forecasts": [], "actuals": []})
    start_date = cutoff_date
    end_date = cutoff_date + timedelta(weeks=12)
    for row in forecasts_list:
        date = row['WeekDate']
        if start_date <= date < end_date:
            key = (row['MaterialID'], row['SalesOrg'], row['DistrChan'], row['CL4'])
            data[key]["forecasts"].append(row['AutoARIMA'])
            actual_row = crossvalidation_df[(crossvalidation_df['ds'] == date) & (crossvalidation_df['unique_id'] == f"{key[0]}_{key[1]}_{key[2]}_{key[3]}")]
            if not actual_row.empty:
                data[key]["actuals"].append(actual_row['y'].values[0])
    return data

# Function to calculate group accuracy and bias for a given lag
def calculate_lag_accuracy_bias(lag, cutoff_date):
    data = extract_data_for_lag(lag, cutoff_date)
    lag_accuracy = defaultdict(float)
    lag_bias = defaultdict(float)
    for key, values in data.items():
        forecasts = values["forecasts"]
        actuals = values["actuals"]
        abs_errors = [abs(forecast - actual) for forecast, actual in zip(forecasts, actuals)]
        sum_abs_errors = sum(abs_errors)
        sum_actuals = sum(actuals)
        sum_forecasts = sum(forecasts)
        if sum_actuals > 0:
            lag_accuracy[key] = 1 - (sum_abs_errors / sum_actuals)
            lag_bias[key] = (sum_forecasts / sum_actuals) - 1
    return lag_accuracy, lag_bias

# Set the cutoff date for lag 0 (e.g., first week of 2023)
cutoff_date = datetime(2023, 1, 2)

# Calculate accuracy and bias for each lag
lag_results = {}
for lag in range(12):
    lag_accuracy, lag_bias = calculate_lag_accuracy_bias(lag, cutoff_date)
    lag_results[f'lag{lag}'] = {'accuracy': lag_accuracy, 'bias': lag_bias}
    cutoff_date += timedelta(weeks=1)

# Print the results along with the 12-week out forecasts
for lag, result in lag_results.items():
    print(f"{lag}:")
    for key, accuracy in result['accuracy'].items():
        forecasts, _ = extract_data_for_lag(lag, cutoff_date - timedelta(weeks=12))
        print(f"  {key}: Accuracy={accuracy:.4f}, Bias={result['bias'][key]:.4f}, 12-week out forecasts={forecasts[key]['forecasts']}")