I'm using tsfeatures from nixtla tsfeatures
Calculates various features from time series data. Python implementation of the R package tsfeatures.

Installation
You can install the released version of tsfeatures from the Python package index with:

pip install tsfeatures
Usage
The tsfeatures main function calculates by default the features used by Montero-Manso, Talagala, Hyndman and Athanasopoulos in their implementation of the FFORMA model.

from tsfeatures import tsfeatures
This function receives a panel pandas df with columns unique_id, ds, y and optionally the frequency of the data.



tsfeatures(panel, freq=7)
By default (freq=None) the function will try to infer the frequency of each time series (using infer_freq from pandas on the ds column) and assign a seasonal period according to the built-in dictionary FREQS:

FREQS = {'H': 24, 'D': 1,
         'M': 12, 'Q': 4,
         'W':1, 'Y': 1}
You can use your own dictionary using the dict_freqs argument:

tsfeatures(panel, dict_freqs={'D': 7, 'W': 52})
List of available features
Features		
acf_features	heterogeneity	series_length
arch_stat	holt_parameters	sparsity
count_entropy	hurst	stability
crossing_points	hw_parameters	stl_features
entropy	intervals	unitroot_kpss
flat_spots	lumpiness	unitroot_pp
frequency	nonlinearity	
guerrero	pacf_features	
See the docs for a description of the features. To use a particular feature included in the package you need to import it:

from tsfeatures import acf_features

tsfeatures(panel, freq=7, features=[acf_features])
 this is the output 	unique_id	nperiods	seasonal_period	trend	spike	linearity	curvature	e_acf1	e_acf10	seasonal_strength	peak	trough
0	10000495_US03_8_1131028	1	52	1	NaN	NaN	NaN	NaN	NaN	1	1	1
1	10000650_US03_10_1131023	1	52	1	NaN	NaN	NaN	NaN	NaN	1	1	1
2	10000687_US03_8_1131028	1	52	1	NaN	NaN	NaN	NaN	NaN	1	1	1
3	10000691_US03_8_1131035	1	52	1	NaN	NaN	NaN	NaN	NaN	1	1	1
4	10000696_US03_8_1131028	1	52	1	NaN	NaN	NaN	NaN	NaN	1	1	1
5	10000701_US03_8_1131028	1	52	1	NaN	NaN	NaN	NaN	NaN	1	1	1
6	10000701_US03_8_5643999	1	52	1	NaN	NaN	NaN	NaN	NaN	1	1	1
7	10000702_US03_8_1131028	1	52	1	NaN	NaN	NaN	NaN	NaN	1	1	1
8	11000263_US01_10_5076999	1	52	1	NaN	NaN	NaN	NaN	NaN	1	1	1
9	11000263_US01_10_6067493	1	52	0	0.000000e+00	0.000000e+00	-0.000000e+00	NaN	NaN	1	6	1
10	11000288_US01_10_5076999	1	52	0	NaN	0.000000e+00	NaN	NaN	NaN	1	2	1
11	11000288_US01_10_6665066	1	52	0	0.000000e+00	0.000000e+00	-0.000000e+00	NaN	NaN	1	2	3
12	11000326_US01_8_1131023	1	52	1	NaN	NaN	NaN	NaN	NaN	1	1	1
13	11000327_US01_8_1131037	1	52	1	NaN	NaN	NaN	NaN	NaN	1	1	1
14	11000329_US01_8_5072051	1	52	1	NaN	NaN	NaN	NaN	NaN	1	1	1
15	11000329_US01_8_5076999	1	52	0	1.665659e-69	1.025779e-16	2.134575e-18	-0.014356	0.510390	1	1	16
16	11000329_US01_8_5338738	1	52	0	4.789750e-70	2.192671e-16	1.150726e-18	-0.125017	0.551564	1	27	13
17	11000329_US01_8_5339183	1	52	0	2.456462e-69	4.526088e-17	-2.142536e-17	0.032159	0.062825	1	5	13
18	11000329_US01_8_5641466	1	52	0	7.145710e-71	1.276334e-17	3.771750e-19	-0.129093	0.149062	1	2	3
19	11000329_US01_8_6304687	1	52	1	NaN	NaN	NaN	NaN	NaN	1	1	1 which value would be the season_length? this is my ensemble model from statsforecast import StatsForecast
from statsforecast.models import AutoARIMA, AutoETS, 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(damped=True, 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', '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('DynamicOptimizedTheta') < 0).then(0).otherwise(pl.col('DynamicOptimizedTheta')).alias('DynamicOptimizedTheta'),
])

# Calculate the ensemble forecast
ensemble_forecast = forecasts_df[['AutoARIMA', 'AutoETS', '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('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', '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)
])

# 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",
    "EnsembleForecast", 
    "Ensemble-lo-95", 
    "Ensemble-hi-95",
    "AutoARIMA", 
    "AutoARIMA-lo-95", 
    "AutoARIMA-hi-95", 
    "AutoETS", 
    "AutoETS-lo-95", 
    "AutoETS-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) you see right here I put season_length= 12, I'm trying to make it automatic, so how do I do that from the output of stl_features models = [
    AutoARIMA(season_length=12),
    AutoETS(damped=True, season_length=12),
    DynamicOptimizedTheta(season_length=12)
]