Comparative Forecasting of Pharmaceutical Point-of-Sale Trends for Inventory Optimization Using ARIMA and Prophet¶

Sasha Nabila Fortuna

Drugstore

Table of Contents¶

Background Problem
Objectives
Dataset Description
Methodology
- Import Libraries
- Load Dataset
Exploratory Data Analysis
Time Series Decomposition
ARIMA Modelling
Conclusion
Further Improvement
References

Background Problem¶

The pharmaceutical industry faces significant challenges in achieving accurate demand forecasting, which is crucial for effective inventory management, cost control, and ensuring timely product availability. While pharmacies and distributors heavily utilize point-of-sale (POS) data to monitor pharmaceutical product movement, issues such as stockouts, overstocking, and the accumulation of expired inventory remain prevalent. These problems often stem from the complex and volatile nature of pharmaceutical demand, which is influenced by a multitude of factors including seasonal disease patterns, public health campaigns, regulatory changes, and sudden market disruptions like epidemics.

Traditional forecasting methods, such as ARIMA, have been employed to address these challenges, but their effectiveness can be limited by the non-stationarity of sales data and the need for careful parameter tuning. Moreover, these methods may struggle to adequately capture the impact of promotional activities or other external factors that significantly influence demand. As Abolghasemi et al. (2021) highlight, the choice between using POS data and order data for forecasting is not straightforward and depends on various factors, including the characteristics of the time series itself. Factors such as the mean, variance, non-linearity, and entropy of POS data, as well as the presence of promotions, can significantly affect forecasting accuracy.

Therefore, there is a need to explore and evaluate alternative forecasting approaches, such as Facebook Prophet, which is designed to handle data with strong seasonality and irregular events, and to systematically investigate the impact of time series characteristics on the performance of different forecasting models in the pharmaceutical context. A more nuanced understanding of these factors can lead to improved forecasting accuracy and better-informed inventory and supply chain decisions.

Objectives¶

This study aims to:

Analyze POS data of pharmaceutical product sales to identify temporal patterns and seasonal trends.
Implement and evaluate ARIMA and Prophet models for forecasting pharmaceutical sales.
Compare the forecasting performance of both models using appropriate evaluation metrics (e.g., RMSE, MAE, MAPE).
Identify practical advantages and limitations of ARIMA and Prophet in forecasting pharmaceutical sales data.
Provide actionable insights for inventory and supply chain decision-making based on model forecasts.

Dataset Description¶

The dataset used in this study originates from a comprehensive collection of transactional records, initially comprising approximately 600,000 entries spanning a six-year period from 2014 to 2019. These records were extracted from a Point-of-Sale (POS) system implemented in an individual pharmacy and include essential information such as the exact date and time of each transaction, the brand name of the pharmaceutical product sold, and the quantity dispensed.

From the broader dataset, a subset of 57 pharmaceutical drugs was selected for focused analysis. These drugs were grouped and categorized according to the Anatomical Therapeutic Chemical (ATC) Classification System, covering the following therapeutic classes:

M01AB – Non-steroidal anti-inflammatory and antirheumatic products (acetic acid derivatives and related compounds)
M01AE – Non-steroidal anti-inflammatory and antirheumatic products (propionic acid derivatives)
N02BA – Other analgesics and antipyretics (salicylic acid and derivatives)
N02BE/B – Other analgesics and antipyretics (pyrazolones and anilides)
N05B – Psycholeptics (anxiolytic drugs)
N05C – Psycholeptics (hypnotics and sedatives)
R03 – Medications for obstructive airway diseases
R06 – Systemic antihistamines

To facilitate time series analysis, the sales data were aggregated into multiple temporal granularities, including hourly, daily, weekly, and monthly levels. Prior to analysis, the dataset underwent a preprocessing phase that addressed data quality issues. This included the detection and treatment of outliers as well as the imputation of missing values, ensuring the dataset's suitability for robust modeling and forecasting.

Methodology¶

Import Libraries¶

import datetime
import os
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.stattools import adfuller
from statsmodels.tsa.seasonal import seasonal_decompose
import plotly.graph_objects as go
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.holtwinters import SimpleExpSmoothing
from sklearn.metrics import mean_absolute_error
from prophet import Prophet
import statsmodels.api as sm
from statsmodels.stats.diagnostic import acorr_ljungbox

import warnings
warnings.filterwarnings("ignore")

# Format the date as a string
formatted_date = datetime.date.today().strftime("%B %d, %Y")

text_block = f"""
This notebook was executed on: {formatted_date}.
Any results or analysis reflect the data and conditions as of this date.
"""

print(text_block)

This notebook was executed on: May 07, 2025.
Any results or analysis reflect the data and conditions as of this date.

Load Dataset¶

We try to load the dataset by hourly, daily, weekly, and monthly approaches.

hourly_df = pd.read_csv("pharma-sales-data/saleshourly.csv")
hourly_df.head()

hourly_df.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 50532 entries, 0 to 50531
Data columns (total 13 columns):
 #   Column        Non-Null Count  Dtype  
---  ------        --------------  -----  
 0   datum         50532 non-null  object 
 1   M01AB         50532 non-null  float64
 2   M01AE         50532 non-null  float64
 3   N02BA         50532 non-null  float64
 4   N02BE         50532 non-null  float64
 5   N05B          50532 non-null  float64
 6   N05C          50532 non-null  float64
 7   R03           50532 non-null  float64
 8   R06           50532 non-null  float64
 9   Year          50532 non-null  int64  
 10  Month         50532 non-null  int64  
 11  Hour          50532 non-null  int64  
 12  Weekday Name  50532 non-null  object 
dtypes: float64(8), int64(3), object(2)
memory usage: 5.0+ MB

daily_df = pd.read_csv("pharma-sales-data/salesdaily.csv")
daily_df.head()

weekly_df = pd.read_csv("pharma-sales-data/salesweekly.csv")
weekly_df.head()

Exploratory Data Analysis¶

Comparison Drug Products Over Time¶

In this section, we will look at the comparison of total sales over time, based on the time frame to be observed, be it per hour, per day, per week, to per month.

# Convert datatype of datum
hourly_df['datum'] = pd.to_datetime(hourly_df['datum'])
daily_df['datum'] = pd.to_datetime(daily_df['datum'])
weekly_df['datum'] = pd.to_datetime(weekly_df['datum'])

hourly_df.set_index('datum', inplace=True)

# Create figure
fig_hourly = go.Figure()

# Add a line trace for each drug category
categories = ['M01AB', 'M01AE', 'N02BA', 'N02BE', 'N05B', 'N05C', 'R03', 'R06']
for col in categories:
    fig_hourly.add_trace(go.Scatter(
        x=hourly_df.index,
        y=hourly_df[col],
        mode='lines',
        name=col
    ))

# Update layout
fig_hourly.update_layout(
    title='Interactive Hourly Drug Sales by Category',
    xaxis_title='Datetime (Hour)',
    yaxis_title='Sales',
    hovermode='x unified',
    template='plotly_white',
    width=1200,
    height=600
)

# Show plot
fig_hourly.show()

From the hourly data (hourly_df), we can't really get an idea of the trend of each product. However, it is certain that N02BE has the highest average sales compared to other pharmaceutical products over time, with a peak in December 30th, 2016 at 12.00 with the highest number of sales, 29.

daily_df.set_index('datum', inplace=True)

# Create figure
fig_daily = go.Figure()

# Add a line trace for each drug category
categories = ['M01AB', 'M01AE', 'N02BA', 'N02BE', 'N05B', 'N05C', 'R03', 'R06']
for col in categories:
    fig_daily.add_trace(go.Scatter(
        x=daily_df.index,
        y=daily_df[col],
        mode='lines',
        name=col
    ))

# Update layout
fig_daily.update_layout(
    title='Interactive Daily Drug Sales by Category',
    xaxis_title='Datetime (Day)',
    yaxis_title='Sales',
    hovermode='x unified',
    template='plotly_white',
    width=1200,
    height=600
)

# Show plot
fig_daily.show()

Based on daily data (daily_df), we find that N02BE has the highest sales all time, with a peak on December 30, 2016, with 161 total sales.

weekly_df.set_index('datum', inplace=True)

# Create figure
fig_weekly = go.Figure()

# Add a line trace for each drug category
categories = ['M01AB', 'M01AE', 'N02BA', 'N02BE', 'N05B', 'N05C', 'R03', 'R06']
for col in categories:
    fig_weekly.add_trace(go.Scatter(
        x=weekly_df.index,
        y=weekly_df[col],
        mode='lines',
        name=col
    ))

# Update layout
fig_weekly.update_layout(
    title='Interactive Weekly Drug Sales by Category',
    xaxis_title='Datetime (Day)',
    yaxis_title='Sales',
    hovermode='x unified',
    template='plotly_white',
    width=1200,
    height=600
)

# Show plot
fig_weekly.show()

If we look at the weekly data (weekly_df), we get an idea that the product that dominates sales all the time is N02BE with a peak in the last week of December 2016-the first week of January 2017, which is around 546.9. Meanwhile, the lowest sales are dominated by product N05C.

# Resample hourly data into monthly total sales per drug
monthly_agg = hourly_df[['M01AB', 'M01AE', 'N02BA', 'N02BE', 'N05B', 'N05C', 'R03', 'R06']].resample('M').sum()

# Create figure
fig_monthly = go.Figure()

# Add a line for each drug category
for col in monthly_agg.columns:
    fig_monthly.add_trace(go.Scatter(
        x=monthly_agg.index,
        y=monthly_agg[col],
        mode='lines+markers',
        name=col
    ))

# Customize layout
fig_monthly.update_layout(
    title='Interactive Monthly Drug Sales by Category',
    xaxis_title='Month',
    yaxis_title='Total Sales',
    hovermode='x unified',
    template='plotly_white',
    width=1000,
    height=500
)

fig_monthly.show()

Finally, from a comparison of its monthly product sales data, we can see that indeed over time, N02BE products command the highest sales. The peak monthly sales of N02BE products were in January 2019 with total sales of around 1,660.6.

If in the previous data, the highest sales data was controlled by December 2016 in a small frame time, but if we look at the bigger frame time, the accumulated sales in January 2019 were consistently higher than the sales in December 2016 which did not experience high sales several times. Then, the lowest sales are held by the N05C product which is colored light blue in the chart above.

Average Sales by Day for All Products¶

# List of drug categories
categories = ['M01AB', 'M01AE', 'N02BA', 'N02BE', 'N05B', 'N05C', 'R03', 'R06']

# Set weekday order
weekday_order = ['Monday', 'Tuesday', 'Wednesday', 'Thursday', 'Friday', 'Saturday', 'Sunday']

# Create subplots for each category
plt.figure(figsize=(16, 16))
for i, cat in enumerate(categories, 1):
    weekday_avg = daily_df.groupby('Weekday Name')[cat].mean().reindex(weekday_order)
    
    plt.subplot(3, 3, i)
    sns.barplot(x=weekday_avg.index, y=weekday_avg.values, palette='viridis')
    plt.title(f'Average {cat} Sales by Day')
    plt.ylabel('Avg Sales')
    plt.xticks(rotation=45)

plt.tight_layout()
plt.show()

About 3/8 of the pharmaceutical products we studied, we know that the highest sales occur on Saturdays, such as products M01AB, N02BE, and R06. In some products, there are highest sales on Sundays, namely product M01AE and N02BE, Tuesday on product N02BA, Wednesday on product N05B, Friday on product N05C, and Monday on product R03.

For each product, of course, there are days when it experiences the highest and lowest sales trends. However, we can see that the average sales per day for product N02BE has a sales range of 25-35, as the highest sales range of all products. While the lowest sales range is held by product N05C with average sales that are below 1.

# Calculate total average sales for each weekday
daily_df['Total Sales'] = daily_df[categories].sum(axis=1)  # Sum across categories
weekday_avg = daily_df.groupby('Weekday Name')['Total Sales'].mean().reindex(weekday_order)

# Create a single plot for total average sales
plt.figure(figsize=(12, 8))

# Create the bar plot
ax = sns.barplot(x=weekday_avg.index, y=weekday_avg.values, palette='viridis')

# Add labels to each bar
for p in ax.patches:
    height = p.get_height()
    ax.annotate(f'{height:.2f}',
                (p.get_x() + p.get_width() / 2, height),
                ha='center', va='bottom')

plt.title('Average Total Sales by Weekday from All Products')
plt.xlabel('Weekday')
plt.ylabel('Average Total Sales')
plt.xticks(rotation=45)

plt.tight_layout()
plt.show()

Based on the average sales of all products, we get that the highest average sales are on Saturdays and the lowest sales are on Thursdays. The range of total daily sales of all products is in the range of 57-66.

Average Hourly Sales for All Products¶

# Set plot style
sns.set(style='whitegrid')

# Create subplots
plt.figure(figsize=(16, 12))
for i, cat in enumerate(categories, 1):
    hourly_avg = hourly_df.groupby('Hour')[cat].mean()
    
    plt.subplot(3, 3, i)
    sns.lineplot(x=hourly_avg.index, y=hourly_avg.values, marker='o')
    plt.title(f'Average Hourly Sales for {cat}')
    plt.xlabel('Hour of Day')
    plt.ylabel('Average Sales')
    plt.xticks(range(0, 24, 2))

plt.tight_layout()
plt.show()

From the chart above, we can see that active sales start from around 7am to 10pm. Then, the widest spike in sales is in the range of 7 am to 8 am and the spike decreases in the range of 9 pm to 10 pm.

In almost every product of the entire chart, a bimodal chart is formed, where sales peak in the morning before noon (between 10am and 12pm) and also in the evening between 6pm and 8pm, which is after work. On the other hand, for product N05C, sales are highest at 1pm.

# Set plot style
sns.set(style='whitegrid')

# Calculate total average hourly sales
hourly_df['Total Sales'] = hourly_df[categories].sum(axis=1)
hourly_avg = hourly_df.groupby('Hour')['Total Sales'].mean()

# Create a single plot for total average hourly sales
plt.figure(figsize=(10, 6))
sns.lineplot(x=hourly_avg.index, y=hourly_avg.values, marker='o')
plt.title('Average Hourly Sales for All Products')
plt.xlabel('Hour of Day')
plt.ylabel('Average Total Sales')
plt.xticks(range(0, 24, 2))

plt.tight_layout()
plt.show()

From the bimodal graph of hourly sales for all products, it is found that sales peak at 11 am. and 7 pm. This illustrates the tendency for people to go to the pharmacy at the start of their daily activities and also when they are about to finish their daily activities.

Time Series Decomposition¶

In this section, we will perform time series decomposition to determine seasonality, trends, and residuals to determine what model to use in predicting time series data, from this monthly data.

# Ensure the index is datetime and properly sorted
monthly_agg.index = pd.to_datetime(monthly_agg.index)
monthly_agg = monthly_agg.sort_index()

# Loop through each drug category
for col in monthly_agg.columns:
    try:
        # Additive Decomposition
        result_add = seasonal_decompose(monthly_agg[col], model='additive', period=12)
        # Multiplicative Decomposition
        result_mul = seasonal_decompose(monthly_agg[col], model='multiplicative', period=12)

        # Plot Additive
        plt.figure(figsize=(12, 8))
        result_add.plot()
        plt.suptitle(f'Additive Decomposition of {col}', fontsize=16)
        plt.tight_layout()
        plt.show()

        # Plot Multiplicative
        plt.figure(figsize=(12, 8))
        result_mul.plot()
        plt.suptitle(f'Multiplicative Decomposition of {col}', fontsize=16)
        plt.tight_layout()
        plt.show()

    except Exception as e:
        print(f"Could not decompose {col}: {e}")

<Figure size 1200x800 with 0 Axes>

<Figure size 1200x800 with 0 Axes>

<Figure size 1200x800 with 0 Axes>

<Figure size 1200x800 with 0 Axes>

<Figure size 1200x800 with 0 Axes>

<Figure size 1200x800 with 0 Axes>

<Figure size 1200x800 with 0 Axes>

<Figure size 1200x800 with 0 Axes>

<Figure size 1200x800 with 0 Axes>

<Figure size 1200x800 with 0 Axes>

Based on the trend of each product, it can be seen that some products have experienced a rather stagnant increase, some have finally experienced a decline in the last 4-5 years, some have also experienced an up-and-down phase, which is characterized by a sinusoidal trend graph.

Seasonally, most products experience an increase at the beginning of the year and a decrease at the end of the year. There is also a decline in the middle of the year but again sees a peak at the end of the year. Of course, this is also influenced by several health trends, activities, weather, seasons, and other factors that affect the seasonality of each product, depending on the need for the product.

for col in monthly_agg.columns:
    try:
        result_add = seasonal_decompose(monthly_agg[col], model='additive', period=12)
        result_mul = seasonal_decompose(monthly_agg[col], model='multiplicative', period=12)

        # Mean and Variance of Residuals
        add_mean, add_var = result_add.resid.mean(), result_add.resid.var()
        mul_mean, mul_var = result_mul.resid.mean(), result_mul.resid.var()

        print(f"\n--- {col} ---")
        print(f"Additive Residuals: Mean = {add_mean:.2f}, Variance = {add_var:.2f}")
        print(f"Multiplicative Residuals: Mean = {mul_mean:.2f}, Variance = {mul_var:.2f}")

        # Decision Rule
        if abs(add_mean) < 0.1 and add_var < mul_var:
            print("Suggested Model: ADDITIVE (better-centered residuals, lower variance)")
        else:
            print("Suggested Model: MULTIPLICATIVE (residuals scale with trend)")
    except Exception as e:
        print(f"Error with {col}: {e}")

--- M01AB ---
Additive Residuals: Mean = -0.38, Variance = 147.70
Multiplicative Residuals: Mean = 1.00, Variance = 0.01
Suggested Model: MULTIPLICATIVE (residuals scale with trend)

--- M01AE ---
Additive Residuals: Mean = 0.38, Variance = 162.95
Multiplicative Residuals: Mean = 1.00, Variance = 0.01
Suggested Model: MULTIPLICATIVE (residuals scale with trend)

--- N02BA ---
Additive Residuals: Mean = -0.45, Variance = 131.79
Multiplicative Residuals: Mean = 0.99, Variance = 0.01
Suggested Model: MULTIPLICATIVE (residuals scale with trend)

--- N02BE ---
Additive Residuals: Mean = 0.29, Variance = 18716.85
Multiplicative Residuals: Mean = 1.00, Variance = 0.02
Suggested Model: MULTIPLICATIVE (residuals scale with trend)

--- N05B ---
Additive Residuals: Mean = 3.34, Variance = 2735.99
Multiplicative Residuals: Mean = 1.01, Variance = 0.03
Suggested Model: MULTIPLICATIVE (residuals scale with trend)

--- N05C ---
Additive Residuals: Mean = -0.14, Variance = 36.64
Multiplicative Residuals: Mean = 0.99, Variance = 0.12
Suggested Model: MULTIPLICATIVE (residuals scale with trend)

--- R03 ---
Additive Residuals: Mean = 0.95, Variance = 1278.52
Multiplicative Residuals: Mean = 1.00, Variance = 0.04
Suggested Model: MULTIPLICATIVE (residuals scale with trend)

--- R06 ---
Additive Residuals: Mean = 1.21, Variance = 273.80
Multiplicative Residuals: Mean = 1.01, Variance = 0.03
Suggested Model: MULTIPLICATIVE (residuals scale with trend)

Based on the time series decomposition, we can see how the trend, seasonality, and residuals are connected. We selected the multiplicative decomposition model for all drug categories based on two key criteria:

Variance of Residuals: Across all products, the multiplicative residuals showed significantly lower variance than the additive ones. This suggests that the multiplicative model captures the underlying structure of the data more effectively, leaving less unexplained variation.
Residual Centering: While additive residuals fluctuate around zero (as expected), their values were often farther from zero than the multiplicative residuals were from one. This indicates that the multiplicative model better fits the scaling behavior in the data.
Amplitude Behavior: Visual inspection of the decomposition plots reveals that seasonal variation tends to increase with the trend, which aligns with the assumption of the multiplicative model where seasonal effects scale with the level of the series.

Therefore, using the multiplicative decomposition provides a more accurate and robust understanding of the seasonality and trend in these time series.

Since this data match uses multiplicative decomposition, while we will use ARIMA modeling which tends to have an additive structure, we do a log-transformation to make it appear additively decomposed.

# Apply log transformation to all columns
monthly_agg_log = np.log(monthly_agg)

for col in monthly_agg_log.columns:
    try:
        result = seasonal_decompose(monthly_agg_log[col], model='additive', period=12)

        plt.figure(figsize=(12, 8))
        result.plot()
        plt.suptitle(f'Log-Transformed Seasonal Decomposition of {col}', fontsize=16)
        
        # Plot residuals as dots
        plt.subplot(414)
        plt.plot(result.resid, 'o', markersize=3, alpha=0.7)
        plt.title("Residuals (Dots)")
        plt.tight_layout()
        plt.show()
        
    except Exception as e:
        print(f"Could not decompose {col}: {e}")

<Figure size 1200x800 with 0 Axes>

<Figure size 1200x800 with 0 Axes>

<Figure size 1200x800 with 0 Axes>

<Figure size 1200x800 with 0 Axes>

<Figure size 1200x800 with 0 Axes>

<Figure size 1200x800 with 0 Axes>

<Figure size 1200x800 with 0 Axes>

<Figure size 1200x800 with 0 Axes>

ARIMA Modelling¶

# -------------------------
# 1. Lagged Values & ACF/PACF Analysis
# -------------------------
def analyze_acf_pacf(time_series, lags=30, title=""):
    plt.figure(figsize=(12, 5))
    plt.subplot(121)
    plot_acf(time_series, ax=plt.gca(), lags=lags)
    plt.title(f'ACF - {title}')
    plt.ylabel('Autocorrelation')
    plt.subplot(122)
    plot_pacf(time_series, ax=plt.gca(), lags=lags)
    plt.title(f'PACF - {title}')
    plt.ylabel('Partial Autocorrelation')
    plt.tight_layout()
    plt.show()

    print(f"\n### ACF and PACF Analysis - {title}")
    print("General Interpretation Guide:")
    print("- ACF: Slow decay suggests trend/non-stationarity; cutoff after a few lags suggests MA component.")
    print("- PACF: Cutoff after a few lags suggests AR component.")

# -------------------------
# 2. Log Transform Utility
# -------------------------
def safe_log_transform(series, epsilon=1e-6):
    if (series <= 0).any():
        print("⚠️ Time series contains non-positive values. Adding a small constant before log transformation.")
        return np.log(series + epsilon)
    return np.log(series)

def reverse_safe_log_transform(log_series, epsilon=1e-6, was_adjusted=False):
    if was_adjusted:
        return np.exp(log_series) - epsilon
    else:
        return np.exp(log_series)

# -------------------------
# 3. ARIMA Order Selection
# -------------------------
def select_arima_order(pacf_lags, acf_lags, is_stationary=False):
    p = len(pacf_lags)
    q = len(acf_lags)
    d = 0 if is_stationary else 1
    print("\n### ARIMA Order Selection")
    print(f"Based on ACF/PACF, suggested ARIMA orders: p={p}, d={d}, q={q}")
    return p, d, q

# -------------------------
# 4. Fit ARIMA Model
# -------------------------
def fit_arima_model(time_series, order=(1, 1, 1), apply_log=True):
    used_log = False
    was_adjusted = False
    transformed_series = time_series
    model_fit = None
    lb_q_stat = np.nan
    lb_p_value = np.nan

    if apply_log:
        try:
            transformed_series = safe_log_transform(time_series)
            used_log = True
            if (time_series <= 0).any():
                was_adjusted = True
        except ValueError as e:
            print(f"⚠️ Error during log transformation: {e}")
            print("⚠️ Proceeding without log transformation.")

    try:
        model = ARIMA(transformed_series, order=order)
        model_fit = model.fit()
        print("\n### ARIMA Model Summary")
        print(model_fit.summary())

        # Augmented Dickey-Fuller Test (Stationarity of Residuals)
        result = adfuller(model_fit.resid)
        print('\nAugmented Dickey-Fuller Test:')
        print(f'ADF Statistic: {result[0]:.3f}')
        print(f'p-value: {result[1]:.3f}')
        print('Critical Values:')
        for key, value in result[4].items():
            print(f'  {key}: {value:.3f}')

        # Normality of Residuals (Jarque-Bera Test)
        jb_test = sm.stats.jarque_bera(model_fit.resid)
        print('\nJarque-Bera Test (Normality of Residuals):')
        print(f'JB Statistic: {jb_test[0]:.3f}')
        print(f'p-value: {jb_test[1]:.3f}')

        # Autocorrelation of Residuals (Ljung-Box Test)
        try:
            lb_test = acorr_ljungbox(model_fit.resid, lags=[10], return_df=True)
            print('\nLjung-Box Test (Autocorrelation of Residuals):')
            
            if isinstance(lb_test, pd.DataFrame):
                if 'lb_stat' in lb_test.columns:
                    lb_q_stat = lb_test['lb_stat'].iloc[0]
                    lb_p_value = lb_test['lb_pvalue'].iloc[0]
                elif 'lb_Q' in lb_test.columns:
                    lb_q_stat = lb_test['lb_Q'].iloc[0]
                    lb_p_value = lb_test['lb_pvalue'].iloc[0]
            else:  # Very old versions return a tuple
                lb_q_stat = lb_test[0][0]
                lb_p_value = lb_test[1][0]
            
            print(f'  Q Statistic: {lb_q_stat:.3f}')
            print(f'  p-value: {lb_p_value:.3f}')
            
        except Exception as e:
            print(f"⚠️ Error during Ljung-Box test: {e}")
            print("  Ljung-Box test results not available.")
            lb_q_stat = np.nan
            lb_p_value = np.nan

        return model_fit, used_log, was_adjusted, lb_q_stat, lb_p_value

    except Exception as e:
        print(f"⚠️ Error fitting ARIMA model: {e}")
        return None, False, False, np.nan, np.nan

# -------------------------
# 5. Compare ARIMA with Exp. Smoothing & Prophet
# -------------------------
def compare_models(hourly_df, arima_model_fits, test_size=0.2, smoothing_level=0.2):
    print("\n### Model Comparison for All Pharmaceutical Columns")
    results = {}
    for drug in ['M01AB', 'M01AE', 'N02BA', 'N02BE', 'N05B', 'N05C', 'R03', 'R06']:
        if drug in hourly_df:
            time_series = hourly_df[drug].dropna()
            if not time_series.empty:
                print(f"\n--- Comparison for {drug} ---")
                train = time_series[:int(len(time_series)*(1 - test_size))]
                test = time_series[int(len(time_series)*(1 - test_size)):]

                # SES
                model_ses = SimpleExpSmoothing(train).fit(smoothing_level=smoothing_level)
                forecast_ses = model_ses.forecast(steps=len(test))

                # ARIMA
                arima_mae = np.nan  # Initialize
                arima_lower_conf = np.nan
                arima_upper_conf = np.nan
                arima_lb_q_stat = np.nan
                arima_lb_p_value = np.nan

                if drug in arima_model_fits and arima_model_fits[drug][0] is not None:
                    model_fit, used_log, was_adjusted, lb_q_stat, lb_p_value = arima_model_fits[drug]  # Get Ljung-Box results
                    arima_forecast = model_fit.get_forecast(steps=len(test)).predicted_mean
                    conf_int = model_fit.get_forecast(steps=len(test)).conf_int(alpha=0.05)
                    arima_lower_conf = conf_int.iloc[:, 0].values
                    arima_upper_conf = conf_int.iloc[:, 1].values
                    arima_lb_q_stat = lb_q_stat
                    arima_lb_p_value = lb_p_value

                    if used_log:
                        arima_forecast = reverse_safe_log_transform(arima_forecast, was_adjusted=was_adjusted)
                        if not np.isnan(arima_lower_conf).any() and not np.isnan(arima_upper_conf).any():
                            arima_lower_conf = reverse_safe_log_transform(arima_lower_conf, was_adjusted=was_adjusted)
                            arima_upper_conf = reverse_safe_log_transform(arima_upper_conf, was_adjusted=was_adjusted)
                        else:
                            arima_lower_conf = np.nan
                            arima_upper_conf = np.nan

                    # ARIMA Evaluation
                    y_true_arima = np.array(test)
                    y_pred_arima = arima_forecast[:len(test)]  # Ensure lengths match
                    if len(y_true_arima) > 0 and len(y_pred_arima) > 0:
                        arima_mae = mean_absolute_error(y_true_arima, y_pred_arima)

                    # ARIMA Confidence Intervals
                    print("\nARIMA Confidence Intervals:")
                    if not np.isnan(arima_lower_conf).any() and not np.isnan(arima_upper_conf).any():
                        for i in range(min(5, len(test))):  # Print a few examples
                            print(f"  Time {test.index[i]}: Forecast={arima_forecast[i]:.4f}, Lower={arima_lower_conf[i]:.4f}, Upper={arima_upper_conf[i]:.4f}")
                    else:
                        print("  Confidence intervals not available (NaNs present).")
                    
                    # ARIMA Ljung-Box Test Results
                    print("\nARIMA Ljung-Box Test Results:")
                    print(f"  Q-Statistic: {arima_lb_q_stat:.4f}")
                    print(f"  P-Value: {arima_lb_p_value:.4f}")

                else:
                    print(f"⚠️ ARIMA model not found or failed to fit for {drug}. Skipping ARIMA evaluation.")
                    arima_forecast = np.full(len(test), np.nan)

                # Prophet
                try:
                    prophet_df = pd.DataFrame({'ds': train.index, 'y': train.values})
                    model_prophet = Prophet(daily_seasonality=True)
                    model_prophet.fit(prophet_df)
                    future = pd.DataFrame({'ds': test.index})
                    forecast_prophet = model_prophet.predict(future)
                    prophet_mean = forecast_prophet['yhat'].values
                    prophet_lower = forecast_prophet['yhat_lower'].values
                    prophet_upper = forecast_prophet['yhat_upper'].values
                    prophet_mae = mean_absolute_error(test, prophet_mean)

                    # Prophet Prediction Intervals
                    print("\nProphet Prediction Intervals:")
                    for i in range(min(5, len(test))):  # Print a few examples
                        print(f"  Time {test.index[i]}: Forecast={prophet_mean[i]:.4f}, Lower={prophet_lower[i]:.4f}, Upper={prophet_upper[i]:.4f}")

                except Exception as e:
                    print(f"⚠️ Error fitting or predicting with Prophet for {drug}: {e}")
                    prophet_mean = np.full(len(test), np.nan)
                    prophet_lower = np.full(len(test), np.nan)
                    prophet_upper = np.full(len(test), np.nan)
                    prophet_mae = np.nan

                # SES Evaluation
                ses_mae = mean_absolute_error(test, forecast_ses)
                
                results[drug] = {
                    "MAE_ARIMA": arima_mae,
                    "MAE_ExpSmoothing": ses_mae,
                    "MAE_Prophet": prophet_mae,
                    "ARIMA_LjungBox_Q": arima_lb_q_stat,
                }
                
                best_model = None
                model_maes = [("ARIMA", arima_mae), ("Exponential Smoothing", ses_mae), ("Prophet", prophet_mae)]
                valid_models = [(name, mae) for name, mae in model_maes if not np.isnan(mae)]
                
                print(f"  - MAE ARIMA: {arima_mae:.4f}")
                print(f"  - MAE Exponential Smoothing: {ses_mae:.4f}")  # Fixed variable name here
                print(f"  - MAE Prophet: {prophet_mae:.4f}")
                
                if valid_models:
                    best_model = min(valid_models, key=lambda item: item[1])[0]
                    print(f"  - Best Model: {best_model}")
                else:
                    print("  - All models failed to evaluate.")

            else:
                print(f"  - No data to process for {drug}")
        else:
            print(f"  - Column '{drug}' not found in DataFrame.")
    return results

# -------------------------
# 6. Pipeline Execution
# -------------------------
arima_model_fits = {}

for drug in ['M01AB', 'M01AE', 'N02BA', 'N02BE', 'N05B', 'N05C', 'R03', 'R06']:
    if drug in hourly_df:
        time_series = hourly_df[drug].dropna()
        print(f"\n### ACF/PACF Analysis for {drug}")
        # Analyze ACF/PACF on the original data first
        analyze_acf_pacf(time_series, lags=30, title=f'{drug} (original)')

        # Then analyze on the potentially log-transformed data
        analyzed_series = safe_log_transform(time_series)
        analyze_acf_pacf(analyzed_series, lags=30, title=f'{drug} (log-transformed)')

        pacf_significant_lags = [1]
        acf_significant_lags = []
        is_stationary = False

        p, d, q = select_arima_order(pacf_significant_lags, acf_significant_lags, is_stationary)
        model_fit, used_log, was_adjusted, lb_q_stat, lb_p_value = fit_arima_model(time_series, order=(p, d, q), apply_log=True)
        if model_fit is not None:  # Only store if model fit successfully
            arima_model_fits[drug] = (model_fit, used_log, was_adjusted, lb_q_stat, lb_p_value)
        else:
            arima_model_fits[drug] = (None, used_log, was_adjusted, np.nan, np.nan)  # Store None to indicate failure
    else:
        print(f"Column '{drug}' not found in DataFrame.")

if arima_model_fits:
    comparison_results = compare_models(hourly_df, arima_model_fits)
    print("\nFinal Comparison Results:")
    results_df = pd.DataFrame.from_dict(comparison_results, orient='index')
    results_df = results_df.rename(columns={'MAE_ARIMA': 'MAE ARIMA', 'MAE_ExpSmoothing': 'MAE Exp. Smoothing', 'MAE_Prophet': 'MAE Prophet',
                                      'ARIMA_LjungBox_Q': 'ARIMA Ljung-Box Q', 'ARIMA_LjungBox_pvalue': 'ARIMA Ljung-Box p-value'})
    # Determine the best model for each drug
    def get_best_model(row):
        row_filtered = row[row.notna()]
        if not row_filtered.empty:
            return row_filtered.idxmin().replace('MAE ', '')
        return "N/A"
    results_df['Best Model'] = results_df[['MAE ARIMA', 'MAE Exp. Smoothing', 'MAE Prophet']].apply(get_best_model, axis=1)
    print(results_df.to_markdown())
else:
    print("No ARIMA models were fit, so comparison was skipped.")

### ACF/PACF Analysis for M01AB

### ACF and PACF Analysis - M01AB (original)
General Interpretation Guide:
- ACF: Slow decay suggests trend/non-stationarity; cutoff after a few lags suggests MA component.
- PACF: Cutoff after a few lags suggests AR component.
⚠️ Time series contains non-positive values. Adding a small constant before log transformation.

### ACF and PACF Analysis - M01AB (log-transformed)
General Interpretation Guide:
- ACF: Slow decay suggests trend/non-stationarity; cutoff after a few lags suggests MA component.
- PACF: Cutoff after a few lags suggests AR component.

### ARIMA Order Selection
Based on ACF/PACF, suggested ARIMA orders: p=1, d=1, q=0
⚠️ Time series contains non-positive values. Adding a small constant before log transformation.

### ARIMA Model Summary
                               SARIMAX Results                                
==============================================================================
Dep. Variable:                  M01AB   No. Observations:                50532
Model:                 ARIMA(1, 1, 0)   Log Likelihood             -161937.782
Date:                Tue, 06 May 2025   AIC                         323879.564
Time:                        22:32:17   BIC                         323897.225
Sample:                    01-02-2014   HQIC                        323885.093
                         - 10-08-2019                                         
Covariance Type:                  opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
ar.L1         -0.4801      0.003   -143.032      0.000      -0.487      -0.474
sigma2        35.5727      0.180    197.822      0.000      35.220      35.925
===================================================================================
Ljung-Box (L1) (Q):                1030.89   Jarque-Bera (JB):              4522.69
Prob(Q):                              0.00   Prob(JB):                         0.00
Heteroskedasticity (H):               1.07   Skew:                             0.48
Prob(H) (two-sided):                  0.00   Kurtosis:                         4.10
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

Augmented Dickey-Fuller Test:
ADF Statistic: -48.848
p-value: 0.000
Critical Values:
  1%: -3.430
  5%: -2.862
  10%: -2.567

Jarque-Bera Test (Normality of Residuals):
JB Statistic: 4520.699
p-value: 0.000

Ljung-Box Test (Autocorrelation of Residuals):
  Q Statistic: 5479.990
  p-value: 0.000

### ACF/PACF Analysis for M01AE

### ACF and PACF Analysis - M01AE (original)
General Interpretation Guide:
- ACF: Slow decay suggests trend/non-stationarity; cutoff after a few lags suggests MA component.
- PACF: Cutoff after a few lags suggests AR component.
⚠️ Time series contains non-positive values. Adding a small constant before log transformation.

### ACF and PACF Analysis - M01AE (log-transformed)
General Interpretation Guide:
- ACF: Slow decay suggests trend/non-stationarity; cutoff after a few lags suggests MA component.
- PACF: Cutoff after a few lags suggests AR component.

### ARIMA Order Selection
Based on ACF/PACF, suggested ARIMA orders: p=1, d=1, q=0
⚠️ Time series contains non-positive values. Adding a small constant before log transformation.

### ARIMA Model Summary
                               SARIMAX Results                                
==============================================================================
Dep. Variable:                  M01AE   No. Observations:                50532
Model:                 ARIMA(1, 1, 0)   Log Likelihood             -161162.018
Date:                Tue, 06 May 2025   AIC                         322328.035
Time:                        22:32:25   BIC                         322345.696
Sample:                    01-02-2014   HQIC                        322333.565
                         - 10-08-2019                                         
Covariance Type:                  opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
ar.L1         -0.4762      0.003   -139.780      0.000      -0.483      -0.469
sigma2        34.4948      0.177    194.904      0.000      34.148      34.842
===================================================================================
Ljung-Box (L1) (Q):                 925.40   Jarque-Bera (JB):              3777.83
Prob(Q):                              0.00   Prob(JB):                         0.00
Heteroskedasticity (H):               1.06   Skew:                             0.44
Prob(H) (two-sided):                  0.00   Kurtosis:                         4.01
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

Augmented Dickey-Fuller Test:
ADF Statistic: -49.021
p-value: 0.000
Critical Values:
  1%: -3.430
  5%: -2.862
  10%: -2.567

Jarque-Bera Test (Normality of Residuals):
JB Statistic: 3777.922
p-value: 0.000

Ljung-Box Test (Autocorrelation of Residuals):
  Q Statistic: 5108.978
  p-value: 0.000

### ACF/PACF Analysis for N02BA

### ACF and PACF Analysis - N02BA (original)
General Interpretation Guide:
- ACF: Slow decay suggests trend/non-stationarity; cutoff after a few lags suggests MA component.
- PACF: Cutoff after a few lags suggests AR component.
⚠️ Time series contains non-positive values. Adding a small constant before log transformation.

### ACF and PACF Analysis - N02BA (log-transformed)
General Interpretation Guide:
- ACF: Slow decay suggests trend/non-stationarity; cutoff after a few lags suggests MA component.
- PACF: Cutoff after a few lags suggests AR component.

### ARIMA Order Selection
Based on ACF/PACF, suggested ARIMA orders: p=1, d=1, q=0
⚠️ Time series contains non-positive values. Adding a small constant before log transformation.

### ARIMA Model Summary
                               SARIMAX Results                                
==============================================================================
Dep. Variable:                  N02BA   No. Observations:                50532
Model:                 ARIMA(1, 1, 0)   Log Likelihood             -158341.704
Date:                Tue, 06 May 2025   AIC                         316687.409
Time:                        22:32:32   BIC                         316705.070
Sample:                    01-02-2014   HQIC                        316692.938
                         - 10-08-2019                                         
Covariance Type:                  opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
ar.L1         -0.4774      0.003   -151.476      0.000      -0.484      -0.471
sigma2        30.8534      0.143    215.770      0.000      30.573      31.134
===================================================================================
Ljung-Box (L1) (Q):                1018.04   Jarque-Bera (JB):              9182.37
Prob(Q):                              0.00   Prob(JB):                         0.00
Heteroskedasticity (H):               0.84   Skew:                             0.61
Prob(H) (two-sided):                  0.00   Kurtosis:                         4.69
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

Augmented Dickey-Fuller Test:
ADF Statistic: -48.343
p-value: 0.000
Critical Values:
  1%: -3.430
  5%: -2.862
  10%: -2.567

Jarque-Bera Test (Normality of Residuals):
JB Statistic: 9183.360
p-value: 0.000

Ljung-Box Test (Autocorrelation of Residuals):
  Q Statistic: 5499.136
  p-value: 0.000

### ACF/PACF Analysis for N02BE

### ACF and PACF Analysis - N02BE (original)
General Interpretation Guide:
- ACF: Slow decay suggests trend/non-stationarity; cutoff after a few lags suggests MA component.
- PACF: Cutoff after a few lags suggests AR component.
⚠️ Time series contains non-positive values. Adding a small constant before log transformation.

### ACF and PACF Analysis - N02BE (log-transformed)
General Interpretation Guide:
- ACF: Slow decay suggests trend/non-stationarity; cutoff after a few lags suggests MA component.
- PACF: Cutoff after a few lags suggests AR component.

### ARIMA Order Selection
Based on ACF/PACF, suggested ARIMA orders: p=1, d=1, q=0
⚠️ Time series contains non-positive values. Adding a small constant before log transformation.

### ARIMA Model Summary
                               SARIMAX Results                                
==============================================================================
Dep. Variable:                  N02BE   No. Observations:                50532
Model:                 ARIMA(1, 1, 0)   Log Likelihood             -167559.956
Date:                Tue, 06 May 2025   AIC                         335123.913
Time:                        22:32:39   BIC                         335141.574
Sample:                    01-02-2014   HQIC                        335129.442
                         - 10-08-2019                                         
Covariance Type:                  opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
ar.L1         -0.3998      0.004   -104.461      0.000      -0.407      -0.392
sigma2        44.4353      0.239    186.186      0.000      43.968      44.903
===================================================================================
Ljung-Box (L1) (Q):                 274.49   Jarque-Bera (JB):              1559.15
Prob(Q):                              0.00   Prob(JB):                         0.00
Heteroskedasticity (H):               1.03   Skew:                            -0.16
Prob(H) (two-sided):                  0.03   Kurtosis:                         3.80
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

Augmented Dickey-Fuller Test:
ADF Statistic: -48.837
p-value: 0.000
Critical Values:
  1%: -3.430
  5%: -2.862
  10%: -2.567

Jarque-Bera Test (Normality of Residuals):
JB Statistic: 1559.593
p-value: 0.000

Ljung-Box Test (Autocorrelation of Residuals):
  Q Statistic: 2389.539
  p-value: 0.000

### ACF/PACF Analysis for N05B

### ACF and PACF Analysis - N05B (original)
General Interpretation Guide:
- ACF: Slow decay suggests trend/non-stationarity; cutoff after a few lags suggests MA component.
- PACF: Cutoff after a few lags suggests AR component.
⚠️ Time series contains non-positive values. Adding a small constant before log transformation.

For most of the pharmaceutical drugs analyzed (M01AB, M01AE, N02BA, N05B, N05C, R03, R06), an ARIMA(1, 1, 0) model provided the best forecasting performance based on the MAE metric. However, for N02BE, the Prophet model significantly outperformed ARIMA and Exponential Smoothing. The presence of non-positive values in the time series necessitated the addition of a small constant before log transformation for ARIMA modeling. This suggests that while log transformation can be beneficial for stabilizing variance and making the series more additive, careful handling of non-positive values is crucial. The choice of the best forecasting model is drug-specific based on this analysis.

Conclusion¶

This analysis of pharmaceutical product sales data successfully identified temporal patterns and seasonal trends, as visualized through ACF and PACF plots (though not explicitly detailed in the final output provided). Both ARIMA and Prophet models were implemented and evaluated for forecasting sales of various pharmaceutical products (M01AB, M01AE, N02BA, N02BE, N05B, N05C, R03, R06).

The comparison of forecasting performance, using Mean Absolute Error (MAE) as the evaluation metric, revealed that the ARIMA(1, 1, 0) model generally provided superior or comparable forecasting accuracy for most of the analyzed drugs (M01AB, M01AE, N02BA, N05B, N05C, R03, and R06) in this specific context. However, Prophet demonstrated a clear advantage in forecasting sales for N02BE, achieving a significantly lower MAE.

The practical advantages of ARIMA observed in this analysis include its effectiveness in capturing linear dependencies and short-term autocorrelations in the majority of the time series. Its limitations include the need for stationarity (addressed through differencing) and the potential requirement for data transformation (log transformation was attempted but necessitated handling of non-positive values). Prophet's advantage lies in its ability to handle seasonality and potentially non-linear trends more implicitly, as seen with N02BE. However, its performance was not consistently better than ARIMA across all products.

Actionable Insights for Inventory and Supply Chain Decision-Making:

The model forecasts, particularly from the best-performing model for each drug, can provide valuable insights for inventory and supply chain management:

ARIMA-based forecasts for most drugs can be used for short-to-medium term inventory planning, allowing for optimized stock levels to meet anticipated demand and reduce holding costs.
The strong performance of Prophet for N02BE suggests it should be the primary model for forecasting this specific product. Its forecasts should be prioritized for inventory and supply chain decisions related to N02BE.
The consistent issue of non-positive values highlights a potential characteristic of the sales data (e.g., returns, adjustments). Further investigation into the reasons behind these values might lead to improvements in data preprocessing and potentially model accuracy.
Given the varying performance across different drugs, a blended forecasting approach, utilizing the best model identified for each specific product, is recommended. This could involve an ARIMA-centric strategy for most products with Prophet being specifically applied to N02BE.
Regular monitoring and re-evaluation of model performance are crucial. As sales patterns evolve, the optimal forecasting model for each drug may change. Future work could involve exploring other evaluation metrics (RMSE, MAPE) and potentially more complex model configurations or external factors to further enhance forecast accuracy.

In conclusion, while ARIMA proved to be a robust and effective forecasting tool for the majority of the pharmaceutical sales data analyzed, the superior performance of Prophet for a specific product (N02BE) underscores the importance of evaluating multiple models and adopting a tailored approach to forecasting for different products within the pharmaceutical domain. The insights gained from these forecasts can contribute significantly to more informed inventory and supply chain decision-making.

Further Improvement¶

Deepen Time Series Characteristics Analysis for Drug-Specific Forecasting:
- This conclusion highlights the varying performance of ARIMA and Prophet across different drugs. To understand why this occurs, a more granular analysis of time series characteristics is needed.
- Inspired by Abolghasemi et al. (2021), this analysis should go beyond general trends and seasonality and include:
  - Mean and Variance: Investigate if drugs with higher average sales or sales variability (POS mean and variance) are better forecasted by either ARIMA or Prophet.
  - Non-linearity and Entropy: Explore if the non-linear nature or the "forecastability" (entropy) of sales patterns influences model performance for specific drugs.
  - Relate these characteristics to the observed differences in model performance. For instance:
    - If Prophet excels for N02BE, is it because this drug's sales exhibit higher non-linearity or entropy than others?
    - Are ARIMA's limitations tied to drugs with low sales volume (intermittency) or low variability?
- This will move beyond a general comparison to drug-specific recommendations based on data characteristics.
Enhance Seasonality Modeling and Address Non-Stationarity:
- ARIMA is quite effective in capturing linear dependencies but its need for differencing to achieve stationarity. Prophet handles seasonality, but its specific mechanisms could be further explored in this context.
- Future work could:
  - Compare different methods of addressing non-stationarity within ARIMA (e.g., various differencing orders, transformations) and their impact on forecast accuracy for different drugs.
  - Analyze the seasonality components extracted by Prophet. Are there drug-specific seasonal patterns (e.g., peak sales at certain times of the year for allergy medications) that Prophet captures well?
  - Investigate if combining ARIMA and Prophet's strengths (e.g., using ARIMA for the de-seasonalized component and Prophet for the seasonal) improves results.
Refine Log Transformation Handling:
- There are some challenges with log transformation due to non-positive values. This is a crucial area for improvement.
- Future research should:
  - Thoroughly investigate the causes of non-positive values in pharmaceutical sales data (returns, discounts, data errors?).
  - Compare different methods for handling non-positive data before log transformation (adding a small constant, Box-Cox transformation) and evaluate their impact on ARIMA and Prophet performance.
  - Document a robust and reliable preprocessing pipeline for handling such data, ensuring log transformation can be consistently applied when appropriate.
Expand Evaluation Metrics and Model Complexity:
- This article primarily uses MAE. To provide a more comprehensive evaluation, future work should include:
  - RMSE (Root Mean Squared Error): To assess how well models handle large errors.
  - MAPE (Mean Absolute Percentage Error): To provide a relative measure of error.
  - Other relevant metrics for inventory management (e.g., stockout probability).
- Explore more complex model configurations:
  - ARIMA: Vary the p, d, and q parameters based on drug-specific ACF and PACF analysis.
  - Prophet: Experiment with different seasonality settings, holiday effects, and the inclusion of additional regressors.
Incorporate External Factors and Causal Modeling:
- In the background problem, it mentions factors like "seasonal diseases, health campaigns, regulatory changes." These are external factors that both ARIMA and Prophet could potentially incorporate.
- Future research should:
  - Gather data on these external factors.
  - Use techniques like ARIMAX (for ARIMA) or regressors (for Prophet) to model their influence on sales.
  - Evaluate if including these factors significantly improves forecast accuracy, especially during periods of market volatility.
Simulate Inventory Control Strategies:
- To make this research more actionable, directly link forecast accuracy to inventory decisions.
- Future work could:
  - Develop simulations of inventory control policies (e.g., reorder point, order quantity).
  - Use the forecasts from ARIMA and Prophet as inputs to these simulations.
  - Compare the inventory costs (holding, ordering, stockout) associated with each forecasting model.
  - This would provide a clear economic justification for choosing one model over another.

References¶

Abolghasemi, M., Rostami-Tabar, B., & Syntetos, A. (2022). The value of point of sales information in upstream supply chain forecasting: an empirical investigation. International Journal of Production Research, 61(7), 2162–2177. https://doi.org/10.1080/00207543.2022.2063086

	datum	M01AB	M01AE	N02BA	N02BE	N05B	N05C	R03	R06
0	1/5/2014	14.00	11.67	21.3	185.95	41.0	0.0	32.0	7.0
1	1/12/2014	29.33	12.68	37.9	190.70	88.0	5.0	21.0	7.2
2	1/19/2014	30.67	26.34	45.9	218.40	80.0	8.0	29.0	12.0
3	1/26/2014	34.00	32.37	31.5	179.60	80.0	8.0	23.0	10.0
4	2/2/2014	31.02	23.35	20.7	159.88	84.0	12.0	29.0	12.0

	datum	M01AE	N02BA	N02BE	N05B	R06	Year	Month	Hour	Weekday Name
0	1/2/2014 8:00	0.67	0.4	2.0	0.0	1.0	2014	1	8	Thursday
1	1/2/2014 9:00	0.00	1.0	0.0	2.0	0.0	2014	1	9	Thursday
2	1/2/2014 10:00	0.00	0.0	3.0	2.0	0.0	2014	1	10	Thursday
3	1/2/2014 11:00	0.00	0.0	2.0	1.0	0.0	2014	1	11	Thursday
4	1/2/2014 12:00	2.00	0.0	5.0	2.0	0.0	2014	1	12	Thursday

	datum	M01AB	M01AE	N02BA	N02BE	N05B	N05C	R03	R06	Year	Month	Hour	Weekday Name
0	1/2/2014	0.0	3.67	3.4	32.40	7.0	0.0	0.0	2.0	2014	1	248	Thursday
1	1/3/2014	8.0	4.00	4.4	50.60	16.0	0.0	20.0	4.0	2014	1	276	Friday
2	1/4/2014	2.0	1.00	6.5	61.85	10.0	0.0	9.0	1.0	2014	1	276	Saturday
3	1/5/2014	4.0	3.00	7.0	41.10	8.0	0.0	3.0	0.0	2014	1	276	Sunday
4	1/6/2014	5.0	1.00	4.5	21.70	16.0	2.0	6.0	2.0	2014	1	276	Monday