Time series forecasting — decomposition, ARIMA আর exponential smoothing

সময়ের ডেটার হাতিয়ার আপনি আগেই চেনেন। AI-এর ভিত্তি কোর্সের পর্ব ২৬-এ DatetimeIndex, resample, rolling, shift আর diff দেখেছেন, আর পর্ব ৫৬-এ শিখেছেন, সময়ের ডেটা এলোমেলো করে ভাগ করা যায় না। এই কোর্সের পর্ব ২১-এ block bootstrap, পর্ব ২৪-এ prediction interval আর coverage, আর পর্ব ৪৬-এ periodic kernel দিয়ে সাপ্তাহিক অর্ডারের অনুমান এসেছে। এবার ধ্রুপদি পরিসংখ্যানের পথ: series-কে ভাগ করা, stationary বানানো, তারপর ARIMA আর exponential smoothing দিয়ে ভবিষ্যৎ বলা, আর সরল baseline-এর সামনে তা যাচাই করা। চারটি চালানো যায় এমন ব্লক পাতাতেই চলে (ব্রাউজারের Python, pandas, statsmodels আর matplotlib); প্রথমবার statsmodels নামতে প্রায় দশ সেকেন্ড লাগতে পারে।

১ · একটি series-এর চার উপাদান — trend, ঋতু, চক্র আর noise

একটি time series হলো সময়ের ক্রমে মাপা একটি মানের ধারা: y₁, y₂, …, y_T. সাধারণ tabular ডেটার সঙ্গে বড় পার্থক্য একটাই: সারিগুলো স্বাধীন নয়, আজকের মান গতকালের মানের সঙ্গে বাঁধা। Forecasting মানে এই বাঁধন কাজে লাগিয়ে y_T+1, y_T+2, …, y_T+h অনুমান করা; h-কে বলে horizon. ধ্রুপদি দৃষ্টিতে series চারটি উপাদানে গড়া:

• Trend: দীর্ঘমেয়াদি দিক, যেমন দোকানের ধীরে বাড়া বিক্রি।
• Seasonality: নির্দিষ্ট দৈর্ঘ্যের পুনরাবৃত্তি, যেমন সপ্তাহে (দৈর্ঘ্য m = ৭) বা বছরে। দৈর্ঘ্য জানা এবং স্থির।
• Cycle (চক্র): ওঠানামা, কিন্তু দৈর্ঘ্য অনির্দিষ্ট, যেমন অর্থনীতির মন্দা আর তেজি। দৈনিক ডেটার দুই বছরে সাধারণত trend-এর ভেতরেই মিশে থাকে।
• Remainder বা noise: বাকি সব, যার কিছু অংশ এলোমেলো আর কিছু অংশ (যেমন ঈদের ধাক্কা) আসলে মডেলে না-ধরা ঘটনা।

Additive না multiplicative? Additive মানে y = T + S + R: শুক্রবারে সবসময় প্রায় ৬০টি অর্ডার কম। Multiplicative মানে y = T × S × R: শুক্রবারে সবসময় প্রায় ১৫% কম, তাই দোকান বড় হলে খাদও গভীর হয়। ছবিতে ঋতুর ঢেউ যদি trend-এর সঙ্গে চওড়া হতে থাকে, তাহলে multiplicative. সহজ কৌশল: log নিলে log y = log T + log S + log R, অর্থাৎ multiplicative series additive হয়ে যায়। forecast করে শেষে exp নিয়ে ফেরত আসতে হয়, আর তখন গড়ের বদলে median পাওয়া যায়; এই সূক্ষ্ম তফাত বড় ভুলের মাপে গুরুত্বপূর্ণ হতে পারে।

২ · Decomposition — classical moving average আর STL

Classical decomposition (১৯২০-এর দশক থেকে প্রচলিত) তিন ধাপে চলে। (ক) m = ৭ দিনের কেন্দ্রিক চলমান গড় দিয়ে trend; এতে প্রথম আর শেষ তিনটি দিনের trend পাওয়া যায় না। (খ) y − T থেকে প্রতিটি বারের গড় নিয়ে ঋতু; ফলে সাত দিনের ছাঁদ দুই বছর ধরে হুবহু এক ধরা হয়। (গ) বাকিটা R. দুটি দুর্বলতা স্পষ্ট: ঋতু বদলাতে পারে না, আর একটি বড় outlier (ঈদের সপ্তাহ) চলমান গড়কে টেনে দেয়।

STL (Seasonal-Trend decomposition using Loess; Cleveland ও সহকর্মীরা, ১৯৯০) একই ভাগ করে loess নামের স্থানীয় মসৃণকরণ দিয়ে, বারবার পালা করে trend আর ঋতু হালনাগাদ করে। সুবিধা: ঋতু ধীরে বদলাতে পারে, ধারের দিনগুলোও বাদ যায় না, আর robust=True দিলে বড় residual-এর ওজন কমিয়ে outlier-এর টান কমায়। কত শক্ত ঋতু, তা মাপার একটি প্রচলিত সংখ্যা seasonal strength:

F_S = max(০, ১ − Var(R) / Var(S + R))

R-এর variance যদি S + R-এর তুলনায় নগণ্য হয়, F_S প্রায় ১: ঋতু প্রধান। আর R-ই সব হলে ০. একইভাবে trend strength F_T = max(০, ১ − Var(R)/Var(T + R)). নিচের ডেটা বাজারঘরের (কাল্পনিক দোকান) দুই বছরের দৈনিক অর্ডার, seeded কোডে তৈরি: ধীরে বাড়া এলোমেলো স্তর, শুক্রবারে খাদওয়ালা সাপ্তাহিক ছাঁদ যা একটু একটু বদলায়, দুটি কাল্পনিক ঈদ-ধাক্কা (তারিখগুলো আসল ঈদের নয়), আর আগের দিনের সঙ্গে বাঁধা noise. Index-এ freq="D" বসানো আছে; statsmodels এটি না পেলে frequency নিয়ে warning দেয়।

# বাজারঘর (কাল্পনিক দোকান): দুই বছরের দৈনিক অর্ডার, seeded
import numpy as np
import pandas as pd

def daily_orders(seed=51, n=728):
    rng = np.random.default_rng(seed)
    idx = pd.date_range("2024-01-01", periods=n, freq="D")  # freq = D
    level = 300 + np.cumsum(0.25 + rng.normal(0, 4, n))  # ধীরে বাড়া, এলোমেলো
    week = np.array([10, 15, 25, -5, -60, 5, 10.0])   # সোম..রবি; শুক্রবার কম
    dow, season = idx.day_of_week, np.empty(n)
    for i in range(n):                         # সাপ্তাহিক ছাঁদ ধীরে বদলায়
        week[dow[i]] += rng.normal(0, 1)
        season[i] = week[dow[i]]
    eid = np.zeros(n)                          # কাল্পনিক ঈদ-ধাক্কা, দুবার
    for start in (160, 510):
        eid[start:start + 10] = 150 * np.hanning(10)
    noise = np.zeros(n)
    for i in range(1, n):                      # AR(1) noise, phi = 0.5
        noise[i] = 0.5 * noise[i - 1] + rng.normal(0, 15)
    return pd.Series(np.round(level + season + eid + noise),
                     index=idx, name="orders")

শেষ ৫৬ দিন (৮ সপ্তাহ) শুরু থেকেই তালাবদ্ধ: decomposition, test আর model বাছাই সবই প্রথম ৬৭২ দিনে। নিচের ব্লক STL চালায়, দুটি strength আর প্রতিটি বারের গড় ঋতু-মান ছাপায়, classical পদ্ধতির trend-এ কয়টি দিন হারায় তা গোনে, আর চারটি অংশের ছবি আঁকে। ছবির লেখা ইংরেজিতে, কারণ পাতার matplotlib বাংলা অক্ষর ঠিকমতো আঁকে না।

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.seasonal import STL, seasonal_decompose
from bazar_ts import daily_orders

y = daily_orders().iloc[:-56]               # শুধু train অংশ
print(len(y), y.index.freqstr, y.index[0].date(), "->", y.index[-1].date())
stl = STL(y, period=7, robust=True).fit()
T, S, R = stl.trend, stl.seasonal, stl.resid
F_s = max(0, 1 - R.var() / (S + R).var())    # seasonal strength
F_t = max(0, 1 - R.var() / (T + R).var())    # trend strength
print(f"seasonal strength = {F_s:.3f}, trend strength = {F_t:.3f}")
prof = S.groupby(S.index.day_name().str[:3]).mean().round(1)
print(prof[["Mon", "Tue", "Wed", "Thu", "Fri", "Sat", "Sun"]].to_dict())
cl = seasonal_decompose(y, model="additive", period=7)
print("classical trend NaN at ends:", int(cl.trend.isna().sum()))
print("biggest STL residuals:")
print(R.abs().nlargest(3).round(1).to_string())
fig, ax = plt.subplots(4, 1, figsize=(9, 7), sharex=True)
fig.subplots_adjust(hspace=0.3)
for a, (s, name) in zip(ax, [(y, "orders"), (T, "trend"),
                             (S, "weekly seasonal"), (R, "remainder")]):
    a.plot(s.index, s.values, lw=0.8)
    a.set_ylabel(name)
ax[0].set_title("STL decomposition of daily orders (period = 7)")
plt.show()

672 D 2024-01-01 -> 2025-11-02
seasonal strength = 0.694, trend strength = 0.941
{'Mon': 8.2, 'Tue': 21.8, 'Wed': 27.3, 'Thu': -8.4, 'Fri': -53.0, 'Sat': -0.8, 'Sun': 4.8}
classical trend NaN at ends: 6
biggest STL residuals:
2024-06-14    136.4
2025-05-29    126.9
2025-05-30    116.7

৩ · Stationarity, differencing আর ADF/KPSS

ARIMA-জাতীয় model ধরে নেয় series stationary: গড়, variance আর “k দিন আগের সঙ্গে সম্পর্ক” সময়ের সঙ্গে বদলায় না। ওপরে-ওঠা অর্ডার stationary নয়: ২০২৪-এর গড় আর ২০২৫-এর গড় আলাদা। প্রতিকার differencing. Backshift operator B লিখলে By_t = y_t−1, তাই

(১ − B)y_t = y_t − y_t−1   আর   (১ − B^৭)y_t = y_t − y_t−7

প্রথমটি “গতকালের চেয়ে কত বেশি”, যা স্থির ঢাল বা এলোমেলো হাঁটা (random walk) স্তর সরায়। দ্বিতীয়টি seasonal differencing: “গত সপ্তাহের একই বারের চেয়ে কত বেশি”, যা সাপ্তাহিক ছাঁদ সরায়। দুটো একসঙ্গে (১ − B)(১ − B^৭)y_t. প্রতিটি differencing একটি করে মান খায় (৭-এরটি ৭টি), আর প্রয়োজনের বেশি করলে over-differencing হয়: noise বাড়ে আর ACF-এর lag ১-এ বড় ঋণাত্মক মান আসে।

Differencing লাগবে কি না, তা যাচাইয়ের দুটি প্রচলিত test আছে, আর তাদের null hypothesis উল্টো:

দুটি একসঙ্গে চালানোর যুক্তি: ADF বলে “unit root বাতিল করতে পারিনি” মানে unit root প্রমাণিত নয়; test-টির power কম হতে পারে। KPSS উল্টো দিক থেকে প্রশ্ন করে। দুটি একমত হলে আত্মবিশ্বাস বাড়ে। একটি ব্যবহারিক ফাঁদ: statsmodels-এর kpss p-value আনে একটি ছোট টেবিল থেকে, যা শুধু ০.০১ থেকে ০.১০ পর্যন্ত। আসল p এর বাইরে হলে এটি সীমার মানটি ফেরত দেয় আর একটি InterpolationWarning দেয়। সঠিক সামলানো মানে warning চাপা দিয়ে সংখ্যাটিকে হুবহু মান না ভাবা, বরং “p > ০.১০” বা “p < ০.০১” হিসেবে পড়া। নিচের ব্লক ঠিক তা করে: শুধু ওই একটি warning-শ্রেণি একটি সীমিত context-এ বন্ধ রাখে আর সীমা চিহ্নসহ ছাপায়।

import warnings
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import adfuller, kpss, acf, pacf
from statsmodels.tools.sm_exceptions import InterpolationWarning
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from bazar_ts import daily_orders

y = daily_orders().iloc[:-56]              # শেষ ৫৬ দিন তালাবদ্ধ
series = {"raw": y, "(1-B)": y.diff().dropna(),
          "(1-B^7)": y.diff(7).dropna(), "both": y.diff(7).diff().dropna()}

def kpss_p(s):                             # টেবিল শুধু 0.01 থেকে 0.10
    with warnings.catch_warnings():
        warnings.simplefilter("ignore", InterpolationWarning)
        stat, p, *_ = kpss(s, regression="c", nlags="auto")
    sign = ">" if p >= 0.10 else ("<" if p <= 0.01 else "=")
    return f"KPSS stat {stat:5.2f}, p {sign} {p:.2f}"

for name, s in series.items():
    print(f"{name:8s} n={len(s)}  ADF p = {adfuller(s, autolag='AIC')[1]:.3f}  {kpss_p(s)}")
np.set_printoptions(precision=2, suppress=True)
d = series["both"]
for name, s in [("raw", y), ("both", d)]:
    print(f"{name:5s} ACF  1-8:", acf(s, nlags=8)[1:])
    print(f"{name:5s} PACF 1-8:", pacf(s, nlags=8)[1:])
fig, ax = plt.subplots(1, 2, figsize=(9, 3))
fig.subplots_adjust(bottom=0.18, wspace=0.25)
plot_acf(d, lags=21, ax=ax[0], title="ACF of (1-B)(1-B^7) y")
plot_pacf(d, lags=21, ax=ax[1], title="PACF of (1-B)(1-B^7) y")
plt.show()

raw      n=672  ADF p = 0.921  KPSS stat  3.47, p < 0.01
(1-B)    n=671  ADF p = 0.000  KPSS stat  0.06, p > 0.10
(1-B^7)  n=665  ADF p = 0.000  KPSS stat  0.03, p > 0.10
both     n=664  ADF p = 0.000  KPSS stat  0.01, p > 0.10
raw   ACF  1-8: [0.88 0.79 0.78 0.77 0.75 0.82 0.89 0.81]
raw   PACF 1-8: [ 0.88  0.07  0.3   0.1   0.11  0.49  0.44 -0.45]
both  ACF  1-8: [-0.25  0.05 -0.08  0.03 -0.16  0.17 -0.44  0.07]
both  PACF 1-8: [-0.25 -0.02 -0.08 -0.01 -0.16  0.09 -0.42 -0.18]

৪ · ACF আর PACF পড়া — AR, MA আর তাদের ছাপ

ACF (autocorrelation function) lag k-তে y_t আর y_t−k-এর correlation দেয়, অর্থাৎ shift(k) করা কলামের সঙ্গে autocorrelation. PACF (partial ACF) একই প্রশ্ন করে, কিন্তু মাঝের ১ থেকে k − ১ lag-এর প্রভাব সরিয়ে: “গতকাল জানার পরেও পরশু কি নতুন কিছু বলে?” ছবিতে নীল ছায়া মোটামুটি ±১.৯৬/√n-এর সীমা; তার বাইরের দাগ “শূন্য নয়” বলে সন্দেহজনক। দুটি মৌলিক model:

• AR(p) (autoregressive): y_t = c + φ₁y_t−1 + … + φ_py_t−p + ε_t. আজকের মান আগের p দিনের মানের ওজনদার যোগফল আর নতুন ধাক্কা ε_t, যা white noise (গড় ০, স্থির variance, পরস্পর সম্পর্কহীন)।
• MA(q) (moving average model): y_t = c + ε_t + θ₁ε_t−1 + … + θ_qε_t−q. আজকের মান আগের q দিনের ধাক্কার ওপর নির্ভর করে। নামে মিল থাকলেও এটি চলমান গড় (rolling mean) নয়।

ব্লক ২-এর কাঁচা series-এ ACF ০.৮৮, ০.৭৯, … খুব ধীরে নামে: differencing দরকারের ক্লাসিক ছাপ। দুবার differencing-এর পরে ACF lag ১-এ −০.২৫ আর lag ৭-এ −০.৪৪, বাকিগুলো ছোট; PACF lag ৭, ৮-এ −০.৪২, −০.১৮ নিয়ে ধীরে ক্ষয়। টেবিল মেলালে: lag ১-এ একটি non-seasonal MA আর lag ৭-এ একটি seasonal MA. এটি বিখ্যাত “airline model” (০,১,১)(০,১,১)_৭-এর ছাপ। আমরা পরের অংশে তার চেয়ে একটু বড় (১,১,১)(১,১,১)_৭ চালাব আর দেখব বাড়তি AR পদ কিছু যোগ করে কি না।

এবার AR(1) দিয়ে forecast হাতে করুন। Mean-কেন্দ্রিক রূপে y_t − μ = φ(y_t−1 − μ) + ε_t, |φ| < ১ হলে stationary. ভবিষ্যতের ধাক্কার গড় ০, তাই h ধাপ পরের অনুমান ŷ_T+h = μ + φ^h(y_T − μ): বিচ্যুতি প্রতি ধাপে φ গুণ ছোট হয়ে গড়ে ফেরে। ভুলের variance σ²(১ + φ² + … + φ^২(h−১)), যা বাড়তে বাড়তে σ²/(১ − φ²)-এ থামে।

৫ · ARIMA আর SARIMA — Box–Jenkins-এর পথ

ARIMA(p, d, q) মানে: d বার differencing করার পরে series-টি ARMA(p, q). SARIMA (p, d, q)(P, D, Q)_m তার সঙ্গে ঋতুর স্তরে আরেক প্রস্থ যোগ করে: D বার seasonal differencing আর lag m, ২m, …-এ P টি AR ও Q টি MA পদ। Backshift-এ লিখলে (p = P = q = Q = ১, d = D = ১, m = ৭):

(১ − φB)(১ − ΦB^৭)(১ − B)(১ − B^৭)y_t = (১ + θB)(১ + ΘB^৭)ε_t

φ, Φ হলো non-seasonal আর seasonal AR coefficient, θ, Θ হলো MA coefficient, ε_t white noise যার variance σ². দুই পাশের গুণফল খুললে আজকের অর্ডার গতকাল, গত সপ্তাহ আর গত সপ্তাহের আগের দিনের মান ও ধাক্কার একটি রৈখিক যোগফল হয়। Fit হয় maximum likelihood দিয়ে (পর্ব ১৮): statsmodels-এর SARIMAX model-টিকে একটি state space রূপে লেখে, আর Kalman filter দিয়ে প্রতিটি দিনের এক-ধাপ অনুমানের ভুল থেকে likelihood গোনে। X অক্ষরটি exogenous (বাইরের) চলকের জন্য, যেমন ঈদের দিন চিহ্নিত একটি কলাম।

George Box আর Gwilym Jenkins-এর ১৯৭০ সালের বই Time Series Analysis: Forecasting and Control একটি চক্র প্রচলিত করে, যা আজও মূল কাঠামো: চেনা (plot, differencing, ACF/PACF) → অনুমান (MLE) → যাচাই (residual white noise কি না) → দরকারে আবার চেনা। কয়েকটি প্রার্থী model-এর মধ্যে AIC দিয়ে বাছাই করা যায়, তবে শুধু একই d আর D-এর model-এর মধ্যে, কারণ differencing বদলালে likelihood ভিন্ন ডেটার ওপর গোনা হয়। চূড়ান্ত বিচার holdout-এ (অংশ ৮)।

৬ · Exponential smoothing — simple, Holt আর Holt–Winters

Simple exponential smoothing (Robert Brown-এর ১৯৫০-এর দশকের কাজে প্রচলিত) একটি স্তর ℓ রাখে: ℓ_t = αy_t + (১ − α)ℓ_t−1, ০ < α ≤ ১. খুললে ওজন α, α(১ − α), α(১ − α)², …; α বড় মানে নতুন ডেটায় দ্রুত সাড়া। সব horizon-এর অনুমান একই ℓ_T. Holt-এর পদ্ধতি (Charles Holt, ১৯৫৭) একটি ঢাল b যোগ করে:

ℓ_t = αy_t + (১ − α)(ℓ_t−1 + b_t−1)
b_t = β(ℓ_t − ℓ_t−1) + (১ − β)b_t−1
ŷ_t+h = ℓ_t + h·b_t

নতুন স্তর হলো “আজ যা দেখলাম” আর “কাল যা ভেবেছিলাম” (ℓ_t−1 + b_t−1)-এর ওজনদার গড়। নতুন ঢাল হলো “স্তর আজ কতটা সরল” আর পুরোনো ঢালের গড়। Holt–Winters (Peter Winters, ১৯৬০) তৃতীয় সমীকরণে ঋতু-মান s যোগ করে, γ ওজন দিয়ে; additive রূপে স্তরের সমীকরণে y_t-এর জায়গায় y_t − s_t−m বসে।

ETS আর পরিসর। মূল Holt–Winters শুধু বিন্দু-অনুমান দেয়। Rob Hyndman ও সহকর্মীরা (২০০২; বই ২০০৮) একে ETS (Error, Trend, Seasonal) state space model হিসেবে লেখেন, যাতে likelihood, AIC আর prediction interval পাওয়া যায়। statsmodels-এ দুটি পথ: ExponentialSmoothing (ধ্রুপদি, parameter বাছে বর্গ-ভুল কমিয়ে) আর ETSModel (state space, likelihood দিয়ে)। আমরা প্রথমটি চালিয়ে তার simulate দিয়ে ১০০০টি সম্ভাব্য ভবিষ্যৎ টেনে পরিসর বের করব। একটি সুন্দর সম্পর্ক: simple exponential smoothing হুবহু ARIMA(০,১,১)-এর forecast দেয়, θ = α − ১ হলে; additive Holt একটি ARIMA(০,২,২)-এর সমান। দুই পরিবার আলাদা নয়, আংশিকভাবে একই model-এর দুই রূপ।

৭ · Baseline আর metric — হারানোর মতো কিছু চাই

একটি forecast ভালো কি না, তার একমাত্র সৎ প্রশ্ন: “সরল নিয়মের চেয়ে ভালো কি না?” দুটি baseline প্রায় বিনা খরচে পাওয়া যায়। Naive: সব ভবিষ্যৎ দিন = শেষ দেখা মান। Seasonal naive: প্রতিটি দিন = শেষ সপ্তাহের একই বার। যাচাই হয় সময় ধরে: শেষ অংশ holdout, আর model শুধু আগের অংশ দেখে (time-based split). একাধিক কাটার বিন্দুতে বারবার এই কাজ (rolling origin backtest) পর্ব ৫২-এর বিষয়।

ভুল e_t = y_t − ŷ_t ধরে চারটি মাপ। MAE = গড় |e|, অর্ডারের এককে। RMSE = √(গড় e²), বড় ভুলে বেশি শাস্তি। MAPE = গড় |e/y| × ১০০% (AI-এর ভিত্তি কোর্সের পর্ব ৬০): এককহীন, কিন্তু y শূন্যের কাছে হলে বিস্ফোরিত হয়, আর বেশি-অনুমান আর কম-অনুমানকে সমান শাস্তি দেয় না। sMAPE = গড় ২|e|/(|y| + |ŷ|) × ১০০% সেই সমস্যা কমাতে আনা, কিন্তু নামের “symmetric” পুরোপুরি সত্য নয়। Rob Hyndman আর Anne Koehler ২০০৬ সালে MASE প্রস্তাব করেন:

MASE = MAE_holdout / ( (১/(T − m)) ∑_t=m+1..T |y_t − y_t−m| )

হর হলো training ডেটায় এক-ধাপ (seasonal) naive-এর গড় পরম ভুল; m = ১ হলে সাধারণ naive. তাই MASE < ১ মানে আপনার forecast গড়ে training-এর ওই naive ভুলের চেয়ে ছোট ভুল করছে। এটি এককহীন, শূন্য-অর্ডারের দিনেও চলে, আর ভিন্ন series-এর মধ্যে তুলনা করা যায়। লেখকেরা একে সাধারণ ব্যবহারের জন্য সুপারিশ করেছিলেন; তবু কোনো একটি metric সব ব্যবসার প্রশ্নের উত্তর নয়, যেমন স্টক ফুরিয়ে যাওয়ার খরচ বেশি হলে quantile loss (পর্ব ২৪) বেশি মানানসই।

৮ · মুখোমুখি — SARIMA বনাম Holt–Winters বনাম baseline

এবার পুরো পরীক্ষা। প্রথম ৬৭২ দিনে SARIMA(১,১,১)(১,১,১)_৭ আর additive Holt–Winters (trend আর ঋতু দুটোই additive) fit করা হয়, তারপর ৫৬ দিন সামনে অনুমান, প্রতিটির সঙ্গে ৯৫% prediction interval. SARIMA-র পরিসর আসে model-এর সূত্র থেকে; Holt–Winters-এর পরিসর ১০০০টি simulated ভবিষ্যতের ২.৫ আর ৯৭.৫ percentile থেকে। চারটি metric নিজে লেখা, কোনো library ছাড়া। পরের ব্লক একই namespace-এ চলে (data-shared), তাই model আবার fit হয় না।

import numpy as np
import pandas as pd
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.tsa.holtwinters import ExponentialSmoothing
from bazar_ts import daily_orders

y = daily_orders()
train, test = y.iloc[:-56], y.iloc[-56:]         # সময় ধরে ভাগ, শেষ ৫৬ দিন
h, m, a = len(test), 7, test.values
sar = SARIMAX(train, order=(1, 1, 1), seasonal_order=(1, 1, 1, 7)).fit(disp=False)
sar_fc = sar.get_forecast(h)
sar_lo, sar_hi = sar_fc.conf_int(alpha=0.05).values.T
hw = ExponentialSmoothing(train, trend="add", seasonal="add",
                          seasonal_periods=m).fit()
sims = hw.simulate(h, repetitions=1000, error="add", anchor="end",
                   random_state=51)               # ১০০০টি সম্ভাব্য ভবিষ্যৎ
hw_lo, hw_hi = sims.quantile([0.025, 0.975], axis=1).values
print({k: round(float(v), 3) for k, v in sar.params.items()})
print({k: round(float(v), 3) for k, v in hw.params.items() if "smoothing" in k})
preds = {"naive": np.repeat(train.iloc[-1], h),
         "seasonal naive": np.tile(train.values[-m:], h // m),
         "SARIMA": sar_fc.predicted_mean.values, "Holt-Winters": hw.forecast(h).values}
scale = np.mean(np.abs(train.values[m:] - train.values[:-m]))   # MASE-এর হর
rows = {}
for k, p in preds.items():
    e = a - p
    rows[k] = {"MAE": np.mean(np.abs(e)), "RMSE": np.sqrt(np.mean(e ** 2)),
               "MAPE%": 100 * np.mean(np.abs(e) / a),
               "sMAPE%": 200 * np.mean(np.abs(e) / (np.abs(a) + np.abs(p))),
               "MASE": np.mean(np.abs(e)) / scale}
print(f"MASE scale (train seasonal-naive MAE) = {scale:.2f}")
print(pd.DataFrame(rows).T.round(2).to_string())
for k, lo, hi in [("SARIMA", sar_lo, sar_hi), ("HW", hw_lo, hw_hi)]:
    print(f"{k}: 95% coverage {np.mean((a >= lo) & (a <= hi)):.2f}, "
          f"width day 1 = {hi[0] - lo[0]:.0f}, day 56 = {hi[-1] - lo[-1]:.0f}")

{'ar.L1': 0.612, 'ma.L1': -0.924, 'ar.S.L7': 0.073, 'ma.S.L7': -0.943, 'sigma2': 351.142}
{'smoothing_level': 0.686, 'smoothing_trend': 0.0, 'smoothing_seasonal': 0.0}
MASE scale (train seasonal-naive MAE) = 22.84
                  MAE   RMSE  MAPE%  sMAPE%  MASE
naive           29.62  34.49   5.77    5.86  1.30
seasonal naive  24.52  29.02   4.75    4.84  1.07
SARIMA          13.37  16.35   2.61    2.61  0.59
Holt-Winters    17.54  21.78   3.37    3.45  0.77
SARIMA: 95% coverage 1.00, width day 1 = 73, day 56 = 179
HW: 95% coverage 1.00, width day 1 = 78, day 56 = 406

ছবিতে পরিসর দেখা আর residual যাচাই। Model ঠিক হলে residual white noise হওয়ার কথা। Ljung–Box test (Greta Ljung আর George Box, ১৯৭৮) প্রথম k টি residual autocorrelation একসঙ্গে পরীক্ষা করে: Q = n(n + ২)∑ r_j²/(n − j), H₀: সব শূন্য। ARMA-র residual-এ মুক্তির মাত্রা (degrees of freedom) থেকে আনুমানিক parameter-সংখ্যা বাদ যায়, তাই SARIMA-তে model_df=4. SARIMA-র প্রথম ৮টি residual (d + D·m = ১ + ৭) differencing-এর শুরুর অবস্থার কারণে অর্থহীনভাবে বড়, তাই বাদ।

import numpy as np
import matplotlib.pyplot as plt
from statsmodels.stats.diagnostic import acorr_ljungbox
from statsmodels.tsa.stattools import acf

fig, ax = plt.subplots(figsize=(9, 3.6))
fig.subplots_adjust(bottom=0.15)
ax.plot(train.index[-84:], train.values[-84:], color="gray", lw=1, label="train")
ax.plot(test.index, a, color="black", lw=1.3, label="actual (holdout)")
for k, lo, hi, c in [("SARIMA", sar_lo, sar_hi, "C0"),
                     ("Holt-Winters", hw_lo, hw_hi, "C1")]:
    ax.plot(test.index, preds[k], color=c, lw=1.3, label=k)
    ax.fill_between(test.index, lo, hi, color=c, alpha=0.15)
ax.set_title("56-day holdout: forecasts and 95% prediction intervals")
ax.set_ylabel("orders per day")
ax.legend(loc="upper left", fontsize=8)
plt.show()

r_sar = sar.resid.iloc[8:]                  # প্রথম d + D*s = 8টি বাদ
for name, r, df in [("SARIMA", r_sar, 4), ("Holt-Winters", hw.resid, 0)]:
    p = acorr_ljungbox(r, lags=[7, 14], model_df=df)["lb_pvalue"]
    print(f"{name:12s} Ljung-Box p (lag 7, 14) = {p.round(4).tolist()}")
    print(f"{'':12s} residual ACF 1-7 = {acf(r, nlags=7)[1:].round(2)}")
big = r_sar.abs().nlargest(4)
print("largest SARIMA residuals:", [f"{d.date()} {v:.0f}" for d, v in big.items()])

SARIMA       Ljung-Box p (lag 7, 14) = [0.0019, 0.0473]
             residual ACF 1-7 = [-0.05  0.08 -0.01  0.01 -0.11  0.03  0.  ]
Holt-Winters Ljung-Box p (lag 7, 14) = [0.0, 0.0]
             residual ACF 1-7 = [ 0.04 -0.05 -0.15 -0.13 -0.18 -0.01  0.1 ]
largest SARIMA residuals: ['2025-05-29 78', '2024-06-14 70', '2024-06-11 65', '2024-10-17 62']

৯ · AI-এর যুগে — forecasting কোড যাচাই

“আমার বিক্রির ডেটায় ARIMA চালিয়ে আগামী মাসের forecast দাও” বললে AI assistant কয়েক সেকেন্ডে কোড আর একটি সুন্দর ছবি দেয়। এই বিষয়ের ভুলগুলো প্রায়ই ছবিতে ধরা পড়ে না।

AI-এর দেওয়া forecasting কোড হাতে পেলে — যাচাইয়ের চেকলিস্ট

• ভাগ সময় ধরে, আর holdout কোনো ধাপে (scaling, differencing-এর পছন্দ, order বাছাই) ছোঁয়া হয়নি
• Naive আর seasonal naive-এর সঙ্গে একই holdout-এ তুলনা আছে, আর MASE ১-এর নিচে কি না দেখেছি
• একটি ছোট উদাহরণে (যেমন AR(1) বা Holt-এর এক ধাপ) সূত্র হাতে মিলিয়েছি
• Residual-এর ACF আর Ljung–Box দেখেছি, আর পরিসরের coverage আলাদা ডেটায় মেপেছি
• সব function আর parameter নিজের statsmodels সংস্করণের docs-এ খুঁজে পেয়েছি, আর কোনো warning চাপা দিইনি
• আসল বিক্রির গোপন সংখ্যা বা গ্রাহকের তথ্য AI-তে পেস্ট করিনি; দরকারে কাল্পনিক বা scaled ডেটা দিয়েছি