الوحدة 5 — ARIMA وSARIMA

الوحدة 3 عرّفت الاستقرارية والتفاضل وقراءة ACF/PACF؛ والوحدة 4 ثبّتت المرجع الذي يجب هزمه. الآن ندخل النموذج الإحصائيّ الكلاسيكيّ الأشهر: ARIMA وامتداده الموسميّ SARIMA. الهدف ليس استظهار الصيغ، بل بناء نموذج SARIMA على مبيعات صيدليّاتنا، وقراءة بواقيه، وتنبّؤ 28 يومًا مع فواصل ثقة قابلة للعرض على قسم الأعمال.

المكوّنات الثلاثة

اسم ARIMA(p, d, q) اختصار لثلاث فكرات مركّبة، كلّ منها بسيطة على حدة:

• AR(p) أو Autoregressive: يعتمد $y_t$ على قيمه السابقة $y_{t-1}, \dots, y_{t-p}$. الرتبة $p$ هي عدد التأخّرات المستعملة، وقيمتها تُقرأ من PACF (الوحدة 3).
• I(d) أو Integrated: تفاضل $d$ مرّة لتحقيق الاستقرارية. مبيعاتنا احتاجت $d = 1$ (اتّجاه).
• MA(q) أو Moving Average (لا يجب خلطه بمتوسّط متحرّك): يعتمد $y_t$ على أخطاء التنبّؤ السابقة $\varepsilon_{t-1}, \dots, \varepsilon_{t-q}$. الرتبة $q$ تُقرأ من ACF.

الصيغة العامّة بعد تفاضل $d$ مرّة، $z_t = \Delta^d y_t$، هي:

$$z_t = c + \phi_1 z_{t-1} + \dots + \phi_p z_{t-p} + \varepsilon_t + \theta_1 \varepsilon_{t-1} + \dots + \theta_q \varepsilon_{t-q}$$

المعاملات $\phi$ و$\theta$ تُقدَّر بالإمكان الأعظم (MLE). البرنامج يتكفّل بذلك، ونحن نختار الرتب.

SARIMA: إضافة الجزء الموسميّ

مبيعاتنا موسميّة أسبوعيًا وسنويًا؛ ARIMA البسيط لا يكفي. SARIMA(p, d, q)(P, D, Q)$_s$ يُضيف نفس الفكرات الثلاث لكن على تأخّرات موسميّة $s$: $P$ تأخّرات AR موسميّة، و$D$ تفاضل موسميّ، و$Q$ تأخّرات MA موسميّة. من ACF/PACF على السلسلة المُفاضلة مزدوجًا وجدنا في الوحدة 3:

• $p = 1$، $d = 1$، $q = 1$ للجزء غير الموسميّ
• $P = 0$، $D = 1$، $Q = 1$ للجزء الموسميّ، مع $s = 7$

فالمرشّح الأوّل هو SARIMA(1, 1, 1)(0, 1, 1)$_7$.

بناء النموذج بـstatsmodels

import pandas as pd
import numpy as np
from statsmodels.tsa.statespace.sarimax import SARIMAX

df = pd.read_csv("ventes_pharmacies.csv", parse_dates=["date"])
serie = df.set_index("date")["ventes_totales"].asfreq("D")

train = serie[:"2027-08-31"]
test = serie["2027-09-01":"2027-09-28"]

modele = SARIMAX(
    train,
    order=(1, 1, 1),
    seasonal_order=(0, 1, 1, 7),
    enforce_stationarity=False,
    enforce_invertibility=False,
)
resultat = modele.fit(disp=False)
print(resultat.summary())

النتيجة تعرض المعاملات المُقدَّرة، وخطأها المعياريّ، وقيم p الخاصّة بها (نتوقّع أن تكون p للمعامل الرئيسيّ < 0.05)، ومعياريّ AIC وBIC لاستعمالهما في المقارنة بين نماذج مرشّحة.

اختيار الرتب بـAIC، لا بالحدس

بدل تجريب رتبة واحدة، من الأفضل مسح شبكة صغيرة حول الرتب المقروءة من ACF/PACF، واختيار النموذج ذي أدنى AIC (Akaike Information Criterion). كلّ رتبة إضافيّة تخفض الخطأ لكن تعاقبها AIC بمعامل، فتحصل على مقايضة دقّة/تعقيد آليّة.

from itertools import product

resultats = []
for p, q, P, Q in product(range(0, 3), range(0, 3), range(0, 2), range(0, 2)):
    try:
        mod = SARIMAX(train, order=(p, 1, q), seasonal_order=(P, 1, Q, 7),
                      enforce_stationarity=False, enforce_invertibility=False)
        res = mod.fit(disp=False, maxiter=200)
        resultats.append({"p": p, "q": q, "P": P, "Q": Q, "aic": res.aic})
    except Exception as e:
        continue

tri = sorted(resultats, key=lambda r: r["aic"])[:5]
for r in tri:
    print(r)

قاعدة عمليّة: اختر النموذج ذا أدنى AIC، لكن إن كانت الفروق أقلّ من نقطتين، فضّل النموذج الأبسط. وتذكّر أنّ AIC مقياس تناسبيّ داخل مجموعة النماذج نفسها؛ لا يُقارَن AIC بين ARIMA وProphet أو LSTM. للمقارنة بينها نستعمل MAE على مجموعة اختبار (الوحدة 9).

تشخيص البواقي: خطوة لا تقفز عليها

نموذج SARIMA جيّد ينتج بواقٍ تشبه الضجيج الأبيض: بلا اتّجاه، بلا موسميّة متبقّية، ذات ارتباط ذاتيّ قريب من الصفر عند كلّ التأخّرات. أداة plot_diagnostics في statsmodels تعرض المطلوب في أربع لوحات:

resultat.plot_diagnostics(figsize=(12, 8))

نقرأها كالتالي: الرسم الأعلى-يسار (البواقي مع الزمن) يجب أن يبدو عشوائيًا؛ الأعلى-يمين (توزيع البواقي) يجب أن يقارب الطبيعيّ؛ الأسفل-يسار (Q-Q) يجب أن يمتدّ على القطر؛ الأسفل-يمين (ACF للبواقي) يجب ألّا يخرج عن حدود الثقة. إن رأينا قمّة بارزة في ACF عند تأخّر معيّن، فالنموذج فوّت بنية يجب إعادة النظر في الرتب.

اختبار Ljung-Box يُقنّن هذا التشخيص بقيمة p: نُريد p أكبر من 0.05 (فرضيّة صفريّة: البواقي لا ترتبط ببعضها، أي ضجيج أبيض).

التنبّؤ مع فواصل ثقة

هنا يكمن أحد أعظم ميّزات النماذج الإحصائيّة أمام الغابات والشبكات: فواصل الثقة مبنيّة نظريًا، لا بحيلة.

# التنبّؤ على 28 يومًا مع فواصل 95%.
prev = resultat.get_forecast(steps=28)
moyenne = prev.predicted_mean
inf, sup = prev.conf_int(alpha=0.05).iloc[:, 0], prev.conf_int(alpha=0.05).iloc[:, 1]

mae = np.mean(np.abs(moyenne.values - test.values))
print(f"MAE SARIMA(1,1,1)(0,1,1)_7 = {mae:.1f}")

# نُضاف الصفّ إلى جدول المقارنة (الوحدة 4).

في مثالنا نصل عادةً إلى MAE ≈ 118 على 28 يومًا، أي تحسّن 19% تقريبًا على naïf saisonnier ≈ 145 من الوحدة 4. جيّد لكن ليس مذهلًا؛ الوحدات 6 و7 و8 سيُحاولن تجاوزه.

فواصل تنبّؤ تتّسع بالأفق

فواصل SARIMA تتّسع مع الأفق: التنبّؤ باليوم 1 دقيق، والتنبّؤ باليوم 28 محاطٌ بشكّ أكبر بكثير. هذا صحيح رياضيًا ومفيد للأعمال: قسم الشراء يجب أن يعرف أنّ التوصية للأسبوع الرابع أقلّ ثقة، فيبني احتياطًا أعلى. إن رأيت فاصلًا ثابت العرض على أفق طويل، فأنت أمام خطأ في الحساب أو فرضيّة مقلوبة.

حدود SARIMA

SARIMA يفترض بنية خطّية ومكوّنات موسميّة ثابتة الطول ($s$ واحد). لا يُدير بشكل طبيعيّ الموسميّات المتعدّدة (أسبوعيّة + سنويّة معًا)، ولا يستوعب المتنبّئات الخارجيّة بلا امتداد SARIMAX. للحالات المعقّدة، انتقل إلى Prophet (الوحدة 6) أو التعزيز التدرّجي (الوحدة 7).

الخلاصة

• ARIMA(p, d, q) يُركّب AR وI وMA؛ SARIMA(p, d, q)(P, D, Q)$_s$ يُضيف الموسميّة.
• نبدأ من الرتب التي تُوحيها ACF/PACF (الوحدة 3)، ثمّ نمسح شبكة صغيرة ونختار أدنى AIC مع تفضيل الأبسط عند التقارب.
• التشخيص عبر plot_diagnostics واختبار Ljung-Box: البواقي يجب أن تبدو ضجيجًا أبيض؛ خلاف ذلك، النموذج فوّت بنية.
• SARIMA يُعطي فواصل تنبّؤ بنيويّة تتّسع مع الأفق، وهي ميزته الحقيقيّة أمام الشبكات؛ حدوده: بنية خطّيّة، وموسميّة وحيدة الطول.

الوحدة التالية: التمهيد الأسّي (Holt-Winters) وProphet، مع تركيز على نقاط قوّة كلّ منهما وحدوده الحقيقيّة.