Documentation Index

Fetch the complete documentation index at: https://nixtlaverse.nixtla.io/llms.txt

Use this file to discover all available pages before exploring further.

Step-by-step guide on using the CrostonOptimized Model with Statsforecast.

​

Table of Contents

• Introduction
• Croston Optimized Model
• Loading libraries and data
• Explore data with the plot method
• Split the data into training and testing
• Implementation of CrostonOptimized with StatsForecast
• Cross-validation
• Model evaluation
• References

​

Introduction

1. Exponential Smoothing: The Croston Optimized model uses exponential smoothing to capture the trend and seasonality in the intermittent demand data. This helps in identifying the underlying patterns and making more accurate forecasts.
2. Intermittent Demand Estimation: Since intermittent demand data often consists of long periods of zero demand, the Croston Optimized model employs a separate estimation process for the occurrence and size of non-zero demand values. It estimates the probability of occurrence and the average size of non-zero demand intervals, enabling better forecasting of intermittent demand.

​

Croston Optimized method

1. Initialization:
   • Let $(y_t)$ represent the intermittent demand time series data at time $t$.
   • Initialize two sets of variables: $(p_t)$ for the probability of occurrence and $(q_t)$ for the average size of non-zero demand intervals.
   • Initialize the forecast $(F_t)$ and forecast error $(E_t)$ variables as zero.
2. Calculation of $(p_t)$ and $(q_t)$:
   • Calculate the intermittent demand occurrence probability $(p_t)$ using exponential smoothing: $[p_t = \alpha + (1 - \alpha)(p_{t-1}),]$ where $(\alpha)$ is the smoothing parameter (typically set between 0.1 and 0.3).
   • Calculate the average size of non-zero demand intervals $(q_t)$ using exponential smoothing: $[q_t = \beta \cdot y_t + (1 - \beta)(q_{t-1}),]$ where $(\beta)$ is the smoothing parameter (typically set between 0.1 and 0.3).
3. Forecasting:
   • If $(y_t > 0)$ (non-zero demand occurrence):
     • Calculate the forecast $(F_t)$ as the previous forecast $(F_{t-1})$ divided by the average size of non-zero demand intervals $(q_{t-1})$: $[F_t = \frac{{F_{t-1}}}{{q_{t-1}}}]$
     • Calculate the forecast error $(E_t)$ as the difference between the actual demand $(y_t)$ and the forecast $(F_t)$: $[E_t = y_t - F_t]$
   • If $(y_t = 0)$ (zero demand occurrence):
     • Set the forecast $(F_t)$ and forecast error $(E_t)$ as zero.
4. Updating the model:
   • Update the intermittent demand occurrence probability $(p_t)$ and the average size of non-zero demand intervals $(q_t)$ using exponential smoothing as described in step 2.
5. Repeat steps 3 and 4 for each time point in the time series.

​

Some properties of the Optimized Croston Model

• It is accurate, even for time series with irregular demand.
• It is adaptable to different types of intermittent time series.
• It is easy to implement and understand.
• It is robust to outliers.

• Precision: The optimized Croston model has been shown to be more accurate than the classic Croston model for time series with irregular demand. This is because the optimized Croston model uses a probability function to forecast the average interval between orders, which is more accurate than the weighted average of historical orders.
• Adaptability: The optimized Croston model is also more adaptable to different types of intermittent time series. This is because the optimized Croston model uses a probability function to forecast the mean interval between orders, allowing it to accommodate different demand patterns.
• Ease of Implementation and Understanding: The optimized Croston model is easy to implement and understand. This is because the optimized Croston model is a modification of the classical Croston model, which is a well-known and well-understood model.
• Robustness: The optimized Croston model is also robust to outliers. This is because the optimized Croston model uses a probability function to forecast the mean interval between orders, which allows it to ignore outliers.

​

Loading libraries and data

Tip Statsforecast will be needed. To install, see instructions.

import matplotlib.pyplot as plt
import seaborn as sns
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
import plotly.graph_objects as go
plt.style.use('grayscale') # fivethirtyeight  grayscale  classic
plt.rcParams['lines.linewidth'] = 1.5
dark_style = {
    'figure.facecolor': '#008080',  # #212946
    'axes.facecolor': '#008080',
    'savefig.facecolor': '#008080',
    'axes.grid': True,
    'axes.grid.which': 'both',
    'axes.spines.left': False,
    'axes.spines.right': False,
    'axes.spines.top': False,
    'axes.spines.bottom': False,
    'grid.color': '#000000',  #2A3459
    'grid.linewidth': '1',
    'text.color': '0.9',
    'axes.labelcolor': '0.9',
    'xtick.color': '0.9',
    'ytick.color': '0.9',
    'font.size': 12 }
plt.rcParams.update(dark_style)


from pylab import rcParams
rcParams['figure.figsize'] = (18,7)

import pandas as pd
df=pd.read_csv("https://raw.githubusercontent.com/Naren8520/Serie-de-tiempo-con-Machine-Learning/main/Data/intermittend_demand2")

df.head()

• The unique_id (string, int or category) represents an identifier for the series.
• The ds (datestamp) column should be of a format expected by Pandas, ideally YYYY-MM-DD for a date or YYYY-MM-DD HH:MM:SS for a timestamp.
• The y (numeric) represents the measurement we wish to forecast.

df["unique_id"]="1"
df.columns=["ds", "y", "unique_id"]
df.head()

print(df.dtypes)

ds           object
y             int64
unique_id    object
dtype: object

df["ds"] = pd.to_datetime(df["ds"])

​

Explore Data with the plot method

from statsforecast import StatsForecast

StatsForecast.plot(df)

​

Autocorrelation plots

• Look for patterns in charts. Common patterns include repeating peaks and valleys, sawtooth patterns, and plateau patterns.
• Compare ACF and PACF charts. The PACF chart generally has fewer spikes than the ACF chart.
• Consider the length of the time series. ACF and PACF charts for longer time series will have more spikes.
• Use a confidence interval. The ACF and PACF plots also show confidence intervals for the autocorrelation values. If an autocorrelation value is outside the confidence interval, it is likely to be significant.

fig, axs = plt.subplots(nrows=1, ncols=2)

plot_acf(df["y"],  lags=30, ax=axs[0],color="fuchsia")
axs[0].set_title("Autocorrelation");

plot_pacf(df["y"],  lags=30, ax=axs[1],color="lime")
axs[1].set_title('Partial Autocorrelation')

plt.show();

​

Decomposition of the time series

• Level: This is the primary value that averages over time.
• Trend: The trend is the value that causes increasing or decreasing patterns in a time series.
• Seasonality: This is a cyclical event that occurs in a time series for a short time and causes short-term increasing or decreasing patterns in a time series.
• Residual/Noise: These are the random variations in the time series.

​

Additive time series

​

Multiplicative time series

from plotly.subplots import make_subplots

from statsmodels.tsa.seasonal import seasonal_decompose

def plotSeasonalDecompose(
    x,
    model='additive',
    filt=None,
    period=None,
    two_sided=True,
    extrapolate_trend=0,
    title="Seasonal Decomposition"):

    result = seasonal_decompose(
            x, model=model, filt=filt, period=period,
            two_sided=two_sided, extrapolate_trend=extrapolate_trend)
    fig = make_subplots(
            rows=4, cols=1,
            subplot_titles=["Observed", "Trend", "Seasonal", "Residuals"])
    for idx, col in enumerate(['observed', 'trend', 'seasonal', 'resid']):
        fig.add_trace(
            go.Scatter(x=result.observed.index, y=getattr(result, col), mode='lines'),
                row=idx+1, col=1,
            )
    return fig

plotSeasonalDecompose(
    df["y"],
    model="additive",
    period=24,
    title="Seasonal Decomposition")

​

Split the data into training and testing

train = df[df.ds<='2023-01-31 19:00:00']
test = df[df.ds>'2023-01-31 19:00:00']

train.shape, test.shape

((9500, 3), (500, 3))

​

Implementation of CrostonOptimized with StatsForecast

​

Load libraries

from statsforecast import StatsForecast
from statsforecast.models import CrostonOptimized

​

Instantiating Model

season_length = 24 # Hourly data
horizon = len(test) # number of predictions

# We call the model that we are going to use
models = [CrostonOptimized()]

• freq: a string indicating the frequency of the data. (See pandas’ available frequencies.)
• n_jobs: n_jobs: int, number of jobs used in the parallel processing, use -1 for all cores.
• fallback_model: a model to be used if a model fails.

sf = StatsForecast(models=models, freq='h')

​

Fit the Model

# fit the models
sf.fit(df=train)

StatsForecast(models=[CrostonOptimized])

result=sf.fitted_[0,0].model_
result

{'mean': array([27.41841685])}

​

Forecast Method

• h (int): represents the forecast h steps into the future. In this case, 500 hours ahead.

Y_hat = sf.forecast(df=train, h=horizon)
Y_hat

sf.plot(train, Y_hat, max_insample_length=500)

​

Predict method with confidence interval

• h (int): represents the forecast h steps into the future. In this case, 500 hours ahead.

forecast_df = sf.predict(h=horizon)
forecast_df

​

Cross-validation

​

Perform time series cross-validation

• df: training data frame
• h (int): represents h steps into the future that are being forecasted. In this case, 500 hours ahead.
• step_size (int): step size between each window. In other words: how often do you want to run the forecasting processes.
• n_windows(int): number of windows used for cross validation. In other words: what number of forecasting processes in the past do you want to evaluate.

crossvalidation_df = sf.cross_validation(df=df,
                                         h=horizon,
                                         step_size=50,
                                         n_windows=5)

• unique_id: series identifier
• ds: datestamp or temporal index
• cutoff: the last datestamp or temporal index for the n_windows.
• y: true value
• model: columns with the model’s name and fitted value.

crossvalidation_df

​

Model Evaluation

from functools import partial

import utilsforecast.losses as ufl
from utilsforecast.evaluation import evaluate

evaluate(
    test.merge(Y_hat),
    metrics=[ufl.mae, ufl.mape, partial(ufl.mase, seasonality=season_length), ufl.rmse, ufl.smape],
    train_df=train,
)

​

References

1. Changquan Huang • Alla Petukhina. Springer series (2022). Applied Time Series Analysis and Forecasting with Python.
2. Ivan Svetunkov. Forecasting and Analytics with the Augmented Dynamic Adaptive Model (ADAM)
3. James D. Hamilton. Time Series Analysis Princeton University Press, Princeton, New Jersey, 1st Edition, 1994.
4. Nixtla CrostonOptimized API
5. Pandas available frequencies.
6. Rob J. Hyndman and George Athanasopoulos (2018). “Forecasting Principles and Practice (3rd ed)”.
7. Seasonal periods- Rob J Hyndman.