28.2-dars: Klassik bashorat modellari

DARAJA YANGILANISHI:
  l_t = alfa * y_t + (1 - alfa) * l_{t-1}          0 < alfa < 1
      = l_{t-1} + alfa * (y_t - l_{t-1})           "xatoning alfa ulushicha tuzat"

BASHORAT:   f_{T+h} = l_T      (har qanday h uchun - TEKIS chiziq)

VAZNLAR (ochib yozsak):
  l_T = alfa*y_T + alfa(1-alfa)*y_{T-1} + alfa(1-alfa)^2*y_{T-2} + ...
  alfa = 0.2:  0.200  0.160  0.128  0.102  ...  (uzoq xotira, silliq)
  alfa = 0.6:  0.600  0.240  0.096  0.038  ...  (qisqa xotira, tez moslashadi)
  alfa -> 1: naive;   alfa -> 0: boshlang'ich daraja (o'rtachaga yaqin)

PARAMETRLAR: alfa va l_0 - bir qadamli xatolar kvadratlari yig'indisini
  (SSE) minimallashtirib topiladi

HOLT (chiziqli trend):
  l_t = alfa * y_t + (1 - alfa) * (l_{t-1} + b_{t-1})
  b_t = beta * (l_t - l_{t-1}) + (1 - beta) * b_{t-1}
  f_{T+h} = l_T + h * b_T                   <- trend CHEKSIZ davom etadi

DAMPED (so'nuvchi) TREND:  0 < phi < 1
  f_{T+h} = l_T + (phi + phi^2 + ... + phi^h) * b_T
  uzoq ufqda trend "tekislanadi" - amalda ko'pincha xavfsizroq

HOLT-WINTERS (mavsum m):
  qo'shiluvchi:     f_{T+h} = l_T + h*b_T + s_{T+h-m}
  ko'paytiruvchi:   f_{T+h} = (l_T + h*b_T) * s_{T+h-m}
  s_t = gamma * (y_t - l_t) + (1 - gamma) * s_{t-m}       (qo'shiluvchi)
  gamma = 0 -> mavsum profili O'ZGARMAS; gamma katta -> tez o'zgaradi

ETS(E, T, S):
  E - xato:    A (qo'shiluvchi) | M (ko'paytiruvchi)
  T - trend:   N (yo'q) | A (chiziqli) | Ad (so'nuvchi)
  S - mavsum:  N | A | M

  ETS(A,N,N) = SES          ETS(A,A,N) = Holt
  ETS(A,Ad,A) = so'nuvchi trend + qo'shiluvchi mavsum (ko'pincha kuchli sukut)

NIMA BERADI:
  holat fazosi (state space) ko'rinishi -> haqiqiy LIKELIHOOD
  -> AIC bilan model tanlash (bir xil ma'lumotda!)
  -> model asosidagi BASHORAT INTERVALLARI
  xato turi (A yoki M) nuqtali bashoratni deyarli o'zgartirmaydi,
  lekin intervallar va likelihood ni o'zgartiradi

AR(p) - AVTOREGRESSIYA: qator o'z o'tmishiga bog'liq
  y_t = c + phi_1*y_{t-1} + ... + phi_p*y_{t-p} + e_t
  AR(1), |phi| < 1: shokdan keyin o'rtachaga GEOMETRIK qaytadi

MA(q) - SIRPANUVCHI O'RTACHA (xatolar bo'yicha): o'tgan SHOKLARGA bog'liq
  y_t = c + e_t + theta_1*e_{t-1} + ... + theta_q*e_{t-q}
  shok faqat q qadam "eslanadi", keyin butunlay unutiladi

I(d) - INTEGRATSIYA: d marta farqlash
  d = 1:  y'_t = y_t - y_{t-1}  ga ARMA qo'llanadi

ARIMA(p,d,q):  (1 - phi_1 B - ... - phi_p B^p) (1 - B)^d y_t
                  = c + (1 + theta_1 B + ... + theta_q B^q) e_t

BASHORAT XULQI (uzoq ufqda):
  d = 0, c bor  -> qator o'rtachasiga qaytadi
  d = 1, c yo'q -> tekis chiziq (tasodifiy yurish kabi)
  d = 1, c bor  -> chiziqli trend (drift)
  d = 2         -> trend ham o'zgaradi (xavfli ekstrapolyatsiya)
  intervallar: d = 0 da cheklangan kenglikka yetadi, d >= 1 da cheksiz o'sadi

SARIMA(p,d,q)(P,D,Q)_m:
  mavsumiy AR:   (1 - Phi_1 B^m - ...)         "o'tgan yilning shu oyi"
  mavsumiy farq: (1 - B^m)^D                   y_t - y_{t-12}
  mavsumiy MA:   (1 + Theta_1 B^m + ...)       "o'tgan yilning shu oyidagi shok"
  hammasi oddiy qism bilan KO'PAYTIRILADI

"AIRLINE" MODELI:  SARIMA(0,1,1)(0,1,1)_12  (log qatorda)
  (1 - B)(1 - B^12) log y_t = (1 + theta B)(1 + Theta B^12) e_t
  faqat 2 parametr; oylik biznes qatorlari uchun klassik boshlang'ich nuqta

ESLATMA: m katta bo'lsa (kunlik ma'lumot + yillik mavsum, m = 365)
  SARIMA juda sekin va beqaror -> Fourier belgilar bilan ARIMAX yoki ML (28.3)

ACF 1-lagda cho'qqi, keyin uziladi        -> q = 1
PACF 1..p lagda cho'qqi, keyin uziladi    -> p
ACF 12-lagda manfiy cho'qqi               -> Q = 1 (mavsumiy MA)
PACF 12, 24-laglarda cho'qqi              -> P = 1

AIC = -2 * log L + 2 * k          (k - parametrlar soni)
BIC = -2 * log L + k * ln(n)      (kattaroq jarima -> soddaroq model)
dAIC < 2   - modellar deyarli teng -> soddasini oling

LJUNG-BOX testi (H0: 1..L laglarda avtokorrelyatsiya yo'q):
  Q = n(n+2) * sum_{k=1..L} r_k^2 / (n - k)   ~  chi^2 (L - model_df)
  model_df = p + q + P + Q     (baholangan ARMA parametrlari soni)
  L: oddiy qator uchun ~10, mavsumiy uchun 2m (24)

p < 0.05  -> qoldiqda tuzilish QOLGAN: tartibni oshiring, tashqi
             o'zgaruvchi qo'shing yoki almashtirishni o'zgartiring
p katta   -> "tuzilish topilmadi" (isbot emas - quvvat cheklangan)

QO'SHIMCHA: qoldiq ACF grafigi, qoldiq vs vaqt (dispersiya o'zgaradimi?),
qoldiqlarni tashqi omillar bo'yicha guruhlash (bayram oylari!)

1. MODEL TO'G'RI (tuzilish, tartib, mavsum turi)
2. PARAMETRLAR ANIQ MA'LUM (baholash xatosi hisobga olinmaydi)
3. XATOLAR NORMAL va dispersiyasi O'ZGARMAS

amalda uchalasi ham buziladi -> intervallar odatda TOR
"95%" interval backtestda 80-90% qamrashi odatiy hol

seasonal naive  <  ETS(A,Ad,A)  <  SARIMA  <  ARIMAX
(0 parametr)       (4 silliqlash    (3 ARMA     (+ exog koeffitsienti
                    parametri +      parametri   va tashqi ma'lumot
                    boshl. holat)    + d, D)     quvuri)

1. PARAMETRLARNI KAMROQ QAYTA BAHOLASH
   har N oynada fit, oraliqda: natija.apply(yangi_y)  (refit yo'q,
   faqat holat yangilanadi - juda tez)
2. ISSIQ START:  fit(start_params=oldingi.params) - iteratsiyalar kamayadi
3. SETKANI KICHRAYTIRISH: tartib bir marta (birinchi oynalarda) tanlanadi,
   backtest davomida qotiriladi
4. OYNALAR SONINI HISOBLAB TANLANG:
   10 oyna x 2 model x 1 s = 20 s;  60 oyna x setka 16 x 1 s = 16 daqiqa

REGRESSIYA + ARIMA XATOLAR:
  y_t = beta * x_t + u_t,     u_t ~ ARIMA(p,d,q)(P,D,Q)_m
  x_t: bayram kunlari soni, narx, aksiya, harorat, ish kunlari soni

SHART: x ning KELAJAKDAGI qiymatlari bashorat paytida MA'LUM bo'lishi kerak
  taqvim (bayram, Ramazon, ish kunlari)  -> aniq ma'lum        ✅
  rejalashtirilgan aksiya, narx           -> reja bo'yicha      ✅
  harorat, raqobatchi narxi                -> O'ZI bashorat      ⚠️
     (backtestda HAQIQIY harorat bilan sinasangiz - bu sizish:
      ishlab chiqarishda sizda faqat ob-havo BASHORATI bo'ladi)

y(t) = g(t) + s(t) + h(t) + e_t        (qo'shiluvchi regressiya, ARIMA emas)
  g(t) - bo'lakli chiziqli trend; o'zgarish nuqtalari (changepoints) avtomatik,
         ularning soni regulyarizatsiya bilan cheklanadi
  s(t) - mavsumlar Fourier qatorlari bilan: sum a_k cos(2 pi k t / P) + b_k sin(...)
         (haftalik P = 7, yillik P = 365.25 - bir vaqtda bir nechta)
  h(t) - bayramlar ro'yxati (sana + oldin/keyin oynasi)
  intervallar: trend o'zgarishlarini simulyatsiya qilib

KUCHLI: ko'p mavsum, bayramlar, bo'sh sanalar, tushunarli komponentlar
ZAIF:   avtokorrelyatsiyani modellamaydi (qisqa ufqda ARIMA/ETS yutishi mumkin),
        sozlamasiz trend ekstrapolyatsiyasi ba'zan g'alati;
        "avtomatik" = backtestsiz ishonish degani EMAS

# faqat g'oya - bu muhitda prophet o'rnatilmagan
from prophet import Prophet

df = pd.DataFrame({"ds": y.index, "y": np.log(y.to_numpy())})
bayramlar = pd.DataFrame({"holiday": "ramazon_hayiti",
                          "ds": pd.to_datetime(["2024-04-10", "2025-03-30"]),
                          "lower_window": -2, "upper_window": 1})
m = Prophet(holidays=bayramlar, yearly_seasonality=True, weekly_seasonality=True)
m.fit(df)
kelajak = m.make_future_dataframe(periods=90)
prognoz = m.predict(kelajak)[["ds", "yhat", "yhat_lower", "yhat_upper"]]

import warnings
import numpy as np
from statsmodels.tsa.holtwinters import ExponentialSmoothing, SimpleExpSmoothing
from statsmodels.tsa.exponential_smoothing.ets import ETSModel
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.stats.diagnostic import acorr_ljungbox

ly = np.log(y)                                       # y: freq="MS" li Series

with warnings.catch_warnings():
    warnings.simplefilter("ignore")                  # ConvergenceWarning va h.k.
    hw = ExponentialSmoothing(y, trend="add", damped_trend=True,
                              seasonal="mul", seasonal_periods=12).fit()
    ets = ETSModel(ly, error="add", trend="add", damped_trend=True,
                   seasonal="add", seasonal_periods=12).fit(disp=False)
    ar = ARIMA(x, order=(1, 0, 1), trend="c").fit()
    sar = SARIMAX(ly, exog=ramazon, order=(1, 1, 1),
                  seasonal_order=(0, 1, 1, 12)).fit(disp=False)

hw.forecast(12)
jad = ets.get_prediction(start=len(ly), end=len(ly) + 11).summary_frame(alpha=0.05)
f = sar.get_forecast(12, exog=ramazon_kelajak)       # exog kelajakda ma'lum!
f.predicted_mean, f.conf_int(alpha=0.05)

sar.aic, sar.bic                                     # faqat bir xil d, D ichida
acorr_ljungbox(sar.resid.iloc[13:], lags=[24], model_df=3)
sar.apply(ly_yangi, exog=ramazon_yangi)             # refit yo'q, holat yangilanadi
SARIMAX(...).fit(start_params=sar.params, disp=False)   # issiq start

qamrov = np.mean((haqiqiy >= past) & (haqiqiy <= yuqori))   # backtestda

SES: l_t = alfa*y_t + (1-alfa)*l_{t-1}; Holt +trend; HW +mavsum; ETS - likelihood
ARIMA(p,d,q): AR - o'tmish qiymatlar, MA - o'tmish shoklar, d - farqlash
SARIMA (P,D,Q)_m; airline = (0,1,1)(0,1,1)_12
AIC - faqat bir xil d, D va almashtirish ichida; Ljung-Box qoldiqda
interval qamrovi - backtestda o'lchanadi; ARIMAX - exog kelajakda ma'lum bo'lsin

"""Eksponensial silliqlash: SES noldan, Holt, damped trend, Holt-Winters, ETS."""

import warnings

import numpy as np
import pandas as pd
from scipy.optimize import minimize_scalar
from statsmodels.tsa.exponential_smoothing.ets import ETSModel
from statsmodels.tsa.holtwinters import ExponentialSmoothing, SimpleExpSmoothing

OY_MAVSUM = np.log([1.25, 1.18, 1.02, 0.90, 0.86, 0.95,
                    1.08, 1.06, 0.90, 0.88, 0.97, 1.15])


def ramazon_kunlari(sana):
    """Sintetik taqvim: har oyda nechta 'Ramazon kuni' bor (yiliga ~11 kun oldinga)."""
    kunlar = pd.date_range(sana[0], sana[-1] + pd.offsets.MonthEnd(0), freq="D")
    belgi = np.zeros(len(kunlar), dtype=int)
    for k in range(sana[0].year - 2011, sana[-1].year - 2011 + 1):
        bosh = pd.Timestamp(2011 + k, 8, 1) - pd.Timedelta(days=round(10.9 * k))
        belgi[(kunlar >= bosh) & (kunlar < bosh + pd.Timedelta(days=30))] = 1
    return pd.Series(belgi, index=kunlar).resample("MS").sum().to_numpy()


def elektr(seed=0, oylar=180):
    """Viloyatning oylik elektr iste'moli, GWh (sintetik)."""
    rng = np.random.default_rng(seed)
    sana = pd.date_range("2011-01-01", periods=oylar, freq="MS")
    t = np.arange(oylar)
    e = rng.normal(0, 0.02, oylar)
    shovqin = np.zeros(oylar)
    for i in range(1, oylar):
        shovqin[i] = 0.5 * shovqin[i - 1] + e[i]
    daraja = np.cumsum(rng.normal(0, 0.006, oylar))
    ly = (np.log(900) + 0.035 * t / 12 + daraja
          + (OY_MAVSUM - OY_MAVSUM.mean())[sana.month - 1]
          + 0.002 * ramazon_kunlari(sana) + shovqin)
    return pd.Series(np.exp(ly), index=sana, name="GWh")


def ses_noldan(y, alfa):
    """Bir qadamli prognozlar va SSE. l_0 = y_0."""
    daraja = y[0]
    prognoz = np.empty(len(y))
    for t, v in enumerate(y):
        prognoz[t] = daraja
        daraja = alfa * v + (1 - alfa) * daraja
    return prognoz, daraja, float(np.sum((y[1:] - prognoz[1:]) ** 2))


def main() -> None:
    rng = np.random.default_rng(1)
    print("=== 1. SES noldan va statsmodels bilan ===")
    talab = 40 + np.cumsum(rng.normal(0, 1.0, 120)) + rng.normal(0, 3.0, 120)
    q = pd.Series(talab, index=pd.date_range("2016-01-01", periods=120,
                                               freq="MS"))
    r = minimize_scalar(lambda a: ses_noldan(talab, a)[2],
                        bounds=(0.001, 0.999), method="bounded",
                        options={"xatol": 1e-6})
    _, oxirgi, sse = ses_noldan(talab, r.x)
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        sm = SimpleExpSmoothing(q, initialization_method="known",
                                initial_level=talab[0]).fit()
    print(f"  ehtiyot qism talabi, 120 oy: noldan alfa={r.x:.4f}, "
          f"prognoz {oxirgi:.3f}")
    print(f"  statsmodels:               alfa={sm.params['smoothing_level']:.4f}"
          f", prognoz {sm.forecast(1).iloc[0]:.3f}")
    for alfa in [0.2, 0.6]:
        w = alfa * (1 - alfa) ** np.arange(6)
        print(f"  alfa={alfa}: oxirgi 6 kuzatuv vaznlari "
              + " ".join(f"{v:.3f}" for v in w) + f"  (jami {w.sum():.3f})")

    print("\n=== 2. Oylik elektr: modellar va 24 oylik test ===")
    y = elektr()
    tr, te = y.iloc[:-24], y.iloc[-24:]
    ltr = np.log(tr)
    prognoz = {"seasonal naive": np.tile(tr.iloc[-12:].to_numpy(), 2)}
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        modellar = {
            "SES": ExponentialSmoothing(tr).fit(),
            "Holt (trend)": ExponentialSmoothing(tr, trend="add").fit(),
            "HW qo'shiluvchi": ExponentialSmoothing(
                tr, trend="add", seasonal="add", seasonal_periods=12).fit(),
            "HW ko'paytiruvchi": ExponentialSmoothing(
                tr, trend="add", seasonal="mul", seasonal_periods=12).fit(),
            "HW ko'p., damped": ExponentialSmoothing(
                tr, trend="add", damped_trend=True, seasonal="mul",
                seasonal_periods=12).fit(),
        }
        for nom, m in modellar.items():
            prognoz[nom] = m.forecast(24).to_numpy()
        ets = ETSModel(ltr, error="add", trend="add", damped_trend=True,
                       seasonal="add", seasonal_periods=12).fit(disp=False)
        prognoz["ETS(A,Ad,A) log"] = np.exp(ets.forecast(24).to_numpy())
    print(f"  {'model':<20} {'MAE, GWh':>9} {'1-yil':>7} {'2-yil':>7}")
    for nom, p in prognoz.items():
        x = np.abs(p - te.to_numpy())
        print(f"  {nom:<20} {x.mean():>9.2f} {x[:12].mean():>7.2f} "
              f"{x[12:].mean():>7.2f}")
    hw = modellar["HW ko'p., damped"].params
    print(f"  damped HW: alfa={hw['smoothing_level']:.3f}, "
          f"beta={hw['smoothing_trend']:.3f}, "
          f"gamma={hw['smoothing_seasonal']:.3f}, phi={hw['damping_trend']:.3f}")

    print("\n=== 3. ETS: AIC bilan tanlash (log qatorda, bir xil ma'lumot) ===")
    variantlar = [("A", None, False, "A"), ("A", "A", False, "A"),
                  ("A", "A", True, "A"), ("M", "A", True, "A")]
    natija = []
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        for xato, trend, damp, mavs in variantlar:
            m = ETSModel(ltr, error="add" if xato == "A" else "mul",
                         trend=None if trend is None else "add",
                         damped_trend=damp, seasonal="add",
                         seasonal_periods=12).fit(disp=False)
            nom = f"ETS({xato},{'N' if trend is None else trend}"
            nom += f"{'d' if damp else ''},{mavs})"
            p = np.exp(m.forecast(24).to_numpy())
            natija.append((m.aic, nom, np.mean(np.abs(p - te.to_numpy()))))
    eng = min(natija)[0]
    for aic, nom, mae in sorted(natija):
        print(f"  {nom:<14} AIC {aic:>9.2f}  dAIC {aic - eng:>6.2f}  "
              f"test MAE {mae:>6.2f}")

    print("\n=== 4. ETS(A,Ad,A) 95% bashorat intervali (24 oy) ===")
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        jad = ets.get_prediction(start=len(ltr), end=len(ltr) + 23).summary_frame(
            alpha=0.05)
    past, yuqori = np.exp(jad["pi_lower"].to_numpy()), np.exp(jad["pi_upper"].to_numpy())
    ichida = (te.to_numpy() >= past) & (te.to_numpy() <= yuqori)
    for h in [1, 6, 12, 24]:
        print(f"  h={h:>2}: [{past[h - 1]:7.1f}, {yuqori[h - 1]:7.1f}]  "
              f"kenglik {yuqori[h - 1] - past[h - 1]:6.1f}  "
              f"haqiqiy {te.iloc[h - 1]:7.1f}")
    print(f"  qamrov: {ichida.sum()}/24 = {ichida.mean():.3f} (nominal 0.95)")


if __name__ == "__main__":
    main()

=== 1. SES noldan va statsmodels bilan ===
  ehtiyot qism talabi, 120 oy: noldan alfa=0.1419, prognoz 32.512
  statsmodels:               alfa=0.1419, prognoz 32.512
  alfa=0.2: oxirgi 6 kuzatuv vaznlari 0.200 0.160 0.128 0.102 0.082 0.066  (jami 0.738)
  alfa=0.6: oxirgi 6 kuzatuv vaznlari 0.600 0.240 0.096 0.038 0.015 0.006  (jami 0.996)

=== 2. Oylik elektr: modellar va 24 oylik test ===
  model                 MAE, GWh   1-yil   2-yil
  seasonal naive           77.50   38.33  116.66
  SES                     156.22  159.23  153.21
  Holt (trend)            160.66  166.20  155.12
  HW qo'shiluvchi          49.14   30.84   67.44
  HW ko'paytiruvchi        51.59   26.65   76.54
  HW ko'p., damped         65.12   31.10   99.14
  ETS(A,Ad,A) log          73.88   33.60  114.16
  damped HW: alfa=0.620, beta=0.000, gamma=0.000, phi=0.992

=== 3. ETS: AIC bilan tanlash (log qatorda, bir xil ma'lumot) ===
  ETS(A,N,A)     AIC   -654.62  dAIC   0.00  test MAE  76.77
  ETS(A,A,A)     AIC   -653.03  dAIC   1.59  test MAE  48.24
  ETS(A,Ad,A)    AIC   -650.43  dAIC   4.19  test MAE  73.88
  ETS(M,Ad,A)    AIC   -649.77  dAIC   4.85  test MAE  73.88

=== 4. ETS(A,Ad,A) 95% bashorat intervali (24 oy) ===
  h= 1: [ 1439.0,  1597.4]  kenglik  158.4  haqiqiy  1546.6
  h= 6: [ 1079.8,  1283.7]  kenglik  203.9  haqiqiy  1169.7
  h=12: [ 1247.1,  1569.2]  kenglik  322.0  haqiqiy  1426.6
  h=24: [ 1197.7,  1639.3]  kenglik  441.6  haqiqiy  1524.2
  qamrov: 23/24 = 0.958 (nominal 0.95)

"""ARIMA: imzolar, parametrlarni tiklash, AIC setkasi, Ljung-Box va "AIC turli d da" tuzog'i."""

import warnings

import numpy as np
from statsmodels.stats.diagnostic import acorr_ljungbox
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.stattools import acf, adfuller, pacf


def simulyatsiya(rng, n, ar=(), ma=(), kuyish=200):
    """ARMA jarayoni: x_t = sum ar_i x_{t-i} + e_t + sum ma_j e_{t-j}."""
    e = rng.normal(0, 1, n + kuyish)
    x = np.zeros(n + kuyish)
    for t in range(n + kuyish):
        x[t] = e[t]
        for i, a in enumerate(ar, 1):
            if t - i >= 0:
                x[t] += a * x[t - i]
        for j, b in enumerate(ma, 1):
            if t - j >= 0:
                x[t] += b * e[t - j]
    return x[kuyish:]


def moslash(x, order, trend="n"):
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")          # optimizator ogohlantirishlari
        return ARIMA(x, order=order, trend=trend).fit()


def main() -> None:
    rng = np.random.default_rng(0)
    n = 400
    jarayonlar = {
        "AR(2) 0.6, 0.25": dict(ar=(0.6, 0.25)),
        "MA(1) 0.6": dict(ma=(0.6,)),
        "ARMA(1,1) 0.7, 0.4": dict(ar=(0.7,), ma=(0.4,)),
    }
    qatorlar = {k: simulyatsiya(rng, n, **v) for k, v in jarayonlar.items()}
    chegara = 1.96 / np.sqrt(n)
    print("=== 1. ACF va PACF imzolari (n=400, chegara +-0.098) ===")
    for nom, x in qatorlar.items():
        a, p = acf(x, nlags=5), pacf(x, nlags=5, method="ywm")
        yulduz = lambda v: "*" if abs(v) > chegara else " "
        print(f"  {nom}")
        print("    ACF  " + " ".join(f"{v:+.2f}{yulduz(v)}" for v in a[1:]))
        print("    PACF " + " ".join(f"{v:+.2f}{yulduz(v)}" for v in p[1:]))

    print("\n=== 2. Parametrlarni tiklash (baho +- SE, haqiqiy) ===")
    for (nom, x), order, haq in zip(qatorlar.items(),
                                    [(2, 0, 0), (0, 0, 1), (1, 0, 1)],
                                    [(0.6, 0.25), (0.6,), (0.7, 0.4)]):
        r = moslash(x, order)
        nomlar = [k for k in r.param_names if k != "sigma2"]
        qism = ", ".join(f"{k} {r.params[i]:+.3f}+-{r.bse[i]:.3f} ({h})"
                         for i, (k, h) in enumerate(zip(nomlar, haq)))
        print(f"  {nom:<20} {qism}")

    print("\n=== 3. AIC/BIC setkasi (ARMA(1,1) qatori, d=0) ===")
    x = qatorlar["ARMA(1,1) 0.7, 0.4"]
    setka = []
    for p in range(3):
        for q in range(3):
            r = moslash(x, (p, 0, q))
            setka.append((r.aic, r.bic, p, q))
    eng_aic = min(setka)[0]
    eng_bic = min(s[1] for s in setka)
    print(f"  {'(p,q)':<7} {'AIC':>8} {'dAIC':>6} {'BIC':>8} {'dBIC':>6}")
    for aic, bic, p, q in sorted(setka):
        print(f"  ({p},{q})   {aic:>8.1f} {aic - eng_aic:>6.1f} {bic:>8.1f} "
              f"{bic - eng_bic:>6.1f}")
    print(f"  AIC tanlovi: ({sorted(setka)[0][2]},{sorted(setka)[0][3]}), "
          f"BIC tanlovi: {min(setka, key=lambda s: s[1])[2:]}")

    print("\n=== 4. Qoldiqlar diagnostikasi: Ljung-Box (lag 10) ===")
    x = qatorlar["AR(2) 0.6, 0.25"]
    for order in [(1, 0, 0), (2, 0, 0)]:
        r = moslash(x, order)
        lb = acorr_ljungbox(r.resid, lags=[10], model_df=order[0] + order[2])
        p = float(lb["lb_pvalue"].iloc[0])
        a1 = acf(r.resid, nlags=2)
        print(f"  AR({order[0]}) qoldiq: ACF(1) {a1[1]:+.3f}, ACF(2) "
              f"{a1[2]:+.3f}, Ljung-Box p = {p:.4f} -> "
              f"{'qoldiqda tuzilish QOLGAN' if p < 0.05 else 'oq shovqinga o_xshaydi'}")

    print("\n=== 5. Tuzoq: AIC ni turli d orasida solishtirish ===")
    z = 50 + simulyatsiya(np.random.default_rng(7), 150, ar=(0.9,))
    print("  bitta STATSIONAR AR(1) phi=0.9 qatori (n=150), ikki birlikda:")
    print(f"  {'birlik':<8} {'AIC (1,0,0)':>12} {'AIC (2,0,0)':>12} "
          f"{'AIC (0,1,0)':>12}  AIC tanlovi")
    for birlik, k in [("GWh", 1.0), ("MWh", 1000.0)]:
        r10 = moslash(k * z, (1, 0, 0), trend="c")
        r20 = moslash(k * z, (2, 0, 0), trend="c")
        r01 = moslash(k * z, (0, 1, 0))
        aic = {"(1,0,0)": r10.aic, "(2,0,0)": r20.aic, "(0,1,0)": r01.aic}
        print(f"  {birlik:<8} {r10.aic:>12.1f} {r20.aic:>12.1f} "
              f"{r01.aic:>12.1f}  {min(aic, key=aic.get)}")
        if birlik == "GWh":
            farq_ichki, farq_d = r20.aic - r10.aic, r01.aic - r10.aic
        else:
            f1, f2 = r20.aic - r10.aic, r01.aic - r10.aic
            print(f"  d=0 ichida farq (2,0,0)-(1,0,0): {farq_ichki:+.2f} -> "
                  f"{f1:+.2f} (siljish {f1 - farq_ichki:+.2f})")
            print(f"  d lar orasida farq (0,1,0)-(1,0,0): {farq_d:+.2f} -> "
                  f"{f2:+.2f} (siljish {f2 - farq_d:+.2f}; "
                  f"2*ln(1000) = {2 * np.log(1000):.2f})")
    print(f"  sababi: log-likelihood dan chiqarilgan boshlang'ich qadamlar "
          f"d=0 -> {r10.loglikelihood_burn}, d=1 -> {r01.loglikelihood_burn}")

    print("\n=== 6. 40 ta shunday qator: d tanlovi va 12 qadamli xato ===")
    tanlov = {"AIC, GWh": [], "AIC, MWh": [], "ADF (p>=0.05 -> d=1)": []}
    xato = {"d=0 (haqiqiy)": [], "d=1": []}
    for s in range(40):
        z = 50 + simulyatsiya(np.random.default_rng(100 + s), 162, ar=(0.9,))
        tr, te = z[:150], z[150:]
        m0, m1 = moslash(tr, (1, 0, 0), trend="c"), moslash(tr, (0, 1, 0))
        xato["d=0 (haqiqiy)"].append(np.mean(np.abs(m0.forecast(12) - te)))
        xato["d=1"].append(np.mean(np.abs(m1.forecast(12) - te)))
        siljish = 2 * np.log(1000)            # MWh da (0,1,0)-(1,0,0) farqi shuncha kamayadi
        tanlov["AIC, GWh"].append(int(m1.aic < m0.aic))
        tanlov["AIC, MWh"].append(int(m1.aic - siljish < m0.aic))
        tanlov["ADF (p>=0.05 -> d=1)"].append(int(adfuller(
            tr, regression="c", autolag="AIC", result_object=False)[1] >= 0.05))
    for nom, v in tanlov.items():
        print(f"  {nom:<22} d=1 ni tanlagan ulush: {np.mean(v):.3f}")
    for nom, v in xato.items():
        print(f"  {nom:<22} MAE (12 qadam): {np.mean(v):.3f}")

if __name__ == "__main__":
    main()

=== 1. ACF va PACF imzolari (n=400, chegara +-0.098) ===
  AR(2) 0.6, 0.25
    ACF  +0.74* +0.66* +0.56* +0.51* +0.45*
    PACF +0.74* +0.25* +0.02  +0.05  +0.03
  MA(1) 0.6
    ACF  +0.49* +0.12* +0.15* +0.06  -0.01
    PACF +0.49* -0.15* +0.21* -0.13* +0.04
  ARMA(1,1) 0.7, 0.4
    ACF  +0.78* +0.49* +0.28* +0.14* +0.08
    PACF +0.78* -0.33* +0.08  -0.04  +0.04

=== 2. Parametrlarni tiklash (baho +- SE, haqiqiy) ===
  AR(2) 0.6, 0.25      ar.L1 +0.559+-0.053 0.6-bob, ar.L2 +0.264+-0.054 0.25-bob
  MA(1) 0.6            ma.L1 +0.630+-0.042 0.6-bob
  ARMA(1,1) 0.7, 0.4   ar.L1 +0.647+-0.048 0.7-bob, ma.L1 +0.414+-0.054 0.4-bob

=== 3. AIC/BIC setkasi (ARMA(1,1) qatori, d=0) ===
  (p,q)        AIC   dAIC      BIC   dBIC
  (1,1)     1153.2    0.0   1165.2    0.0
  (2,1)     1154.2    1.0   1170.2    5.0
  (1,2)     1154.3    1.1   1170.3    5.1
  (2,0)     1155.7    2.4   1167.6    2.4
  (2,2)     1156.1    2.9   1176.1   10.9
  (0,2)     1187.8   34.5   1199.7   34.5
  (1,0)     1198.6   45.3   1206.5   41.3
  (0,1)     1280.7  127.4   1288.7  123.4
  (0,0)     1580.6  427.4   1584.6  419.4
  AIC tanlovi: (1,1), BIC tanlovi: (1, 1)

=== 4. Qoldiqlar diagnostikasi: Ljung-Box (lag 10) ===
  AR(1) qoldiq: ACF(1) -0.199, ACF(2) +0.131, Ljung-Box p = 0.0000 -> qoldiqda tuzilish QOLGAN
  AR(2) qoldiq: ACF(1) -0.009, ACF(2) -0.035, Ljung-Box p = 0.2141 -> oq shovqinga o_xshaydi

=== 5. Tuzoq: AIC ni turli d orasida solishtirish ===
  bitta STATSIONAR AR(1) phi=0.9 qatori (n=150), ikki birlikda:
  birlik    AIC (1,0,0)  AIC (2,0,0)  AIC (0,1,0)  AIC tanlovi
  GWh             429.1        429.9        427.5  (0,1,0)
  MWh            2501.5       2502.3       2486.0  (0,1,0)
  d=0 ichida farq (2,0,0)-(1,0,0): +0.74 -> +0.72 (siljish -0.03)
  d lar orasida farq (0,1,0)-(1,0,0): -1.62 -> -15.53 (siljish -13.92; 2*ln(1000) = 13.82)
  sababi: log-likelihood dan chiqarilgan boshlang'ich qadamlar d=0 -> 0, d=1 -> 1

=== 6. 40 ta shunday qator: d tanlovi va 12 qadamli xato ===
  AIC, GWh               d=1 ni tanlagan ulush: 0.500
  AIC, MWh               d=1 ni tanlagan ulush: 1.000
  ADF (p>=0.05 -> d=1)   d=1 ni tanlagan ulush: 0.550
  d=0 (haqiqiy)          MAE (12 qadam): 1.574
  d=1                    MAE (12 qadam): 1.814

"""SARIMA: farqlangan qator ACF, kichik AIC setkasi, Ljung-Box va ARIMAX (Ramazon kunlari)."""

import warnings

import numpy as np
import pandas as pd
from statsmodels.stats.diagnostic import acorr_ljungbox
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.tsa.stattools import acf

OY_MAVSUM = np.log([1.25, 1.18, 1.02, 0.90, 0.86, 0.95,
                    1.08, 1.06, 0.90, 0.88, 0.97, 1.15])


def ramazon_kunlari(sana):
    """Sintetik taqvim: har oyda nechta 'Ramazon kuni' bor (yiliga ~11 kun oldinga)."""
    kunlar = pd.date_range(sana[0], sana[-1] + pd.offsets.MonthEnd(0), freq="D")
    belgi = np.zeros(len(kunlar), dtype=int)
    for k in range(sana[0].year - 2011, sana[-1].year - 2011 + 1):
        bosh = pd.Timestamp(2011 + k, 8, 1) - pd.Timedelta(days=round(10.9 * k))
        belgi[(kunlar >= bosh) & (kunlar < bosh + pd.Timedelta(days=30))] = 1
    return pd.Series(belgi, index=kunlar).resample("MS").sum().to_numpy()


def elektr(seed=0, oylar=180):
    """Viloyatning oylik elektr iste'moli, GWh (sintetik)."""
    rng = np.random.default_rng(seed)
    sana = pd.date_range("2011-01-01", periods=oylar, freq="MS")
    t = np.arange(oylar)
    e = rng.normal(0, 0.02, oylar)
    shovqin = np.zeros(oylar)
    for i in range(1, oylar):
        shovqin[i] = 0.5 * shovqin[i - 1] + e[i]
    daraja = np.cumsum(rng.normal(0, 0.006, oylar))
    ly = (np.log(900) + 0.035 * t / 12 + daraja
          + (OY_MAVSUM - OY_MAVSUM.mean())[sana.month - 1]
          + 0.002 * ramazon_kunlari(sana) + shovqin)
    return pd.Series(np.exp(ly), index=sana, name="GWh")


def sarima(y, order, exog=None):
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")          # ConvergenceWarning va h.k.
        return SARIMAX(y, exog=exog, order=order,
                       seasonal_order=(0, 1, 1, 12)).fit(disp=False)


def main() -> None:
    y = elektr()
    ly = np.log(y)
    ram = pd.Series(ramazon_kunlari(y.index), index=y.index, name="ramazon")
    n_te = 24
    ltr, te = ly.iloc[:-n_te], y.iloc[-n_te:]

    print("=== 1. Farqlangan qatorning ACF i: (1-B)(1-B^12) log y ===")
    dd = ltr.diff(12).diff().dropna()
    a = acf(dd, nlags=24)
    chegara = 1.96 / np.sqrt(len(dd))
    print(f"  n={len(dd)}, chegara +-{chegara:.3f}")
    print("  " + "  ".join(f"lag{k}:{a[k]:+.2f}{'*' if abs(a[k]) > chegara else ''}"
                           for k in [1, 2, 3, 11, 12, 13, 24]))
    print("  1-lag va 12-lagda manfiy cho'qqi -> MA(1) va mavsumiy MA(1)")

    print("\n=== 2. SARIMA(p,1,q)(0,1,1)_12: kichik AIC setkasi (d, D bir xil) ===")
    setka = []
    for p in [0, 1]:
        for q in [0, 1]:
            r = sarima(ltr, (p, 1, q))
            prog = np.exp(r.forecast(n_te).to_numpy())
            setka.append((r.aic, p, q, float(np.mean(np.abs(prog - te.to_numpy())))))
    eng = min(setka)[0]
    for aic, p, q, mae in sorted(setka):
        print(f"  ({p},1,{q})(0,1,1)12  AIC {aic:>9.2f}  dAIC {aic - eng:>5.2f}"
              f"  test MAE {mae:>6.2f}")
    _, p_eng, q_eng, _ = min(setka)
    order = (p_eng, 1, q_eng)

    print(f"\n=== 3. Qoldiqlar: SARIMA{order}(0,1,1)12 ===")
    r = sarima(ltr, order)
    qoldiq = r.resid.iloc[13:]                     # farqlash boshlanishi tashlanadi
    df = p_eng + q_eng + 1
    for lag in [12, 24]:
        lb = acorr_ljungbox(qoldiq, lags=[lag], model_df=df)
        print(f"  Ljung-Box lag {lag}: p = {float(lb['lb_pvalue'].iloc[0]):.4f}")
    qoldiq_ram = pd.Series(qoldiq.to_numpy(), index=qoldiq.index)
    ram_oy = ram.loc[qoldiq.index] >= 15
    print(f"  qoldiq o'rtachasi: Ramazonli oylar {qoldiq_ram[ram_oy].mean():+.4f},"
          f" boshqalar {qoldiq_ram[~ram_oy].mean():+.4f}")

    print("\n=== 4. ARIMAX: tashqi o'zgaruvchi - oydagi Ramazon kunlari ===")
    rx = sarima(ltr, order, exog=ram.iloc[:-n_te])
    b, se = rx.params["ramazon"], rx.bse["ramazon"]
    print(f"  koeffitsient: {b:.5f} +- {se:.5f} (haqiqiy 0.00200) -> "
          f"30 kunlik oy {100 * (np.exp(30 * b) - 1):.1f}% ko'proq")
    p_s = np.exp(r.forecast(n_te).to_numpy())
    p_x = np.exp(rx.forecast(n_te, exog=ram.iloc[-n_te:]).to_numpy())
    oy_ram = ram.iloc[-n_te:].to_numpy() >= 15
    print(f"  AIC: SARIMA {r.aic:.2f}, ARIMAX {rx.aic:.2f} (bir xil d, D, "
          f"bir xil ma'lumot)")
    lb = acorr_ljungbox(rx.resid.iloc[13:], lags=[24], model_df=df + 1)
    print(f"  ARIMAX Ljung-Box lag 24: p = {float(lb['lb_pvalue'].iloc[0]):.4f}")
    print(f"  {'model':<8} {'test MAE':>9} {'Ramazonli oylar':>16} {'boshqa':>7}")
    for nom, p in [("SARIMA", p_s), ("ARIMAX", p_x)]:
        x = np.abs(p - te.to_numpy())
        print(f"  {nom:<8} {x.mean():>9.2f} {x[oy_ram].mean():>16.2f} "
              f"{x[~oy_ram].mean():>7.2f}")
    print(f"  test davrida Ramazonli oylar: {int(oy_ram.sum())} ta; kelajak "
          f"exog qiymati oldindan ma'lum (taqvim)")


if __name__ == "__main__":
    main()

=== 1. Farqlangan qatorning ACF i: (1-B)(1-B^12) log y ===
  n=143, chegara +-0.164
  lag1:-0.20*  lag2:-0.05  lag3:-0.18*  lag11:-0.05  lag12:-0.39*  lag13:-0.04  lag24:+0.13
  1-lag va 12-lagda manfiy cho'qqi -> MA(1) va mavsumiy MA(1)

=== 2. SARIMA(p,1,q)(0,1,1)_12: kichik AIC setkasi (d, D bir xil) ===
  (1,1,1)(0,1,1)12  AIC   -591.62  dAIC  0.00  test MAE  45.03
  (0,1,1)(0,1,1)12  AIC   -588.52  dAIC  3.10  test MAE  55.95
  (1,1,0)(0,1,1)12  AIC   -584.19  dAIC  7.43  test MAE  60.54
  (0,1,0)(0,1,1)12  AIC   -571.70  dAIC 19.92  test MAE  71.00

=== 3. Qoldiqlar: SARIMA(1, 1, 1)(0,1,1)12 ===
  Ljung-Box lag 12: p = 0.6543
  Ljung-Box lag 24: p = 0.0490
  qoldiq o'rtachasi: Ramazonli oylar +0.0157, boshqalar -0.0013

=== 4. ARIMAX: tashqi o'zgaruvchi - oydagi Ramazon kunlari ===
  koeffitsient: 0.00209 +- 0.00026 (haqiqiy 0.00200) -> 30 kunlik oy 6.5% ko'proq
  AIC: SARIMA -591.62, ARIMAX -650.63 (bir xil d, D, bir xil ma'lumot)
  ARIMAX Ljung-Box lag 24: p = 0.2439
  model     test MAE  Ramazonli oylar  boshqa
  SARIMA       45.03           103.43   39.73
  ARIMAX       36.77            37.47   36.70
  test davrida Ramazonli oylar: 2 ta; kelajak exog qiymati oldindan ma'lum (taqvim)

"""Rolling-origin backtest: seasonal naive, ETS, SARIMA, ARIMAX; juftlashgan SE va interval qamrovi."""

import warnings

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

OY_MAVSUM = np.log([1.25, 1.18, 1.02, 0.90, 0.86, 0.95,
                    1.08, 1.06, 0.90, 0.88, 0.97, 1.15])


def ramazon_kunlari(sana):
    """Sintetik taqvim: har oyda nechta 'Ramazon kuni' bor (yiliga ~11 kun oldinga)."""
    kunlar = pd.date_range(sana[0], sana[-1] + pd.offsets.MonthEnd(0), freq="D")
    belgi = np.zeros(len(kunlar), dtype=int)
    for k in range(sana[0].year - 2011, sana[-1].year - 2011 + 1):
        bosh = pd.Timestamp(2011 + k, 8, 1) - pd.Timedelta(days=round(10.9 * k))
        belgi[(kunlar >= bosh) & (kunlar < bosh + pd.Timedelta(days=30))] = 1
    return pd.Series(belgi, index=kunlar).resample("MS").sum().to_numpy()


def elektr(seed=0, oylar=180):
    """Viloyatning oylik elektr iste'moli, GWh (sintetik)."""
    rng = np.random.default_rng(seed)
    sana = pd.date_range("2011-01-01", periods=oylar, freq="MS")
    t = np.arange(oylar)
    e = rng.normal(0, 0.02, oylar)
    shovqin = np.zeros(oylar)
    for i in range(1, oylar):
        shovqin[i] = 0.5 * shovqin[i - 1] + e[i]
    daraja = np.cumsum(rng.normal(0, 0.006, oylar))
    ly = (np.log(900) + 0.035 * t / 12 + daraja
          + (OY_MAVSUM - OY_MAVSUM.mean())[sana.month - 1]
          + 0.002 * ramazon_kunlari(sana) + shovqin)
    return pd.Series(np.exp(ly), index=sana, name="GWh")


class SarimaBashorat:
    """SARIMA(1,1,1)(0,1,1)12. Parametrlar har `har` oynada qayta baholanadi
    (oldingi parametrlardan boshlab), oraliqda faqat holat yangilanadi."""

    def __init__(self, exog=None, har=2):
        self.exog, self.har, self.natija, self.iter = exog, har, None, []

    def __call__(self, i, ly, h):
        ex = None if self.exog is None else self.exog.iloc[:len(ly)]
        ex_kel = None if self.exog is None else self.exog.iloc[len(ly):len(ly) + h]
        with warnings.catch_warnings():
            warnings.simplefilter("ignore")
            if self.natija is None or i % self.har == 0:
                m = SARIMAX(ly, exog=ex, order=(1, 1, 1),
                            seasonal_order=(0, 1, 1, 12))
                bosh = None if self.natija is None else self.natija.params
                self.natija = m.fit(start_params=bosh, disp=False)
                self.iter.append(self.natija.mle_retvals["iterations"])
                r = self.natija
            else:
                r = self.natija.apply(ly, exog=ex)   # refit yo'q
            f = r.get_forecast(h, exog=ex_kel)
        ci = f.conf_int(alpha=0.05).to_numpy()
        return f.predicted_mean.to_numpy(), ci[:, 0], ci[:, 1]


def snaive(i, ly, h):
    x = ly.to_numpy()
    f = np.array([x[len(x) - 12 + (k % 12)] for k in range(h)])
    s = np.std(x[12:] - x[:-12], ddof=1)          # mavsumiy farqlar std
    return f, f - 1.96 * s, f + 1.96 * s


def ets(i, ly, h):
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        m = ETSModel(ly, error="add", trend="add", damped_trend=True,
                     seasonal="add", seasonal_periods=12).fit(disp=False)
        jad = m.get_prediction(start=len(ly), end=len(ly) + h - 1).summary_frame(
            alpha=0.05)
    return jad["mean"].to_numpy(), jad["pi_lower"].to_numpy(), jad["pi_upper"].to_numpy()


def main() -> None:
    y = elektr()
    ly = np.log(y)
    ram = pd.Series(ramazon_kunlari(y.index), index=y.index)
    h, n_oyna = 6, 10
    boshlar = [len(y) - h * (i + 1) for i in range(n_oyna)][::-1]
    sarima, arimax = SarimaBashorat(), SarimaBashorat(exog=ram)
    usullar = {"seasonal naive": snaive, "ETS(A,Ad,A)": ets,
               "SARIMA": sarima, "ARIMAX (Ramazon)": arimax}
    mae = {k: [] for k in usullar}
    ichida = {k: [] for k in usullar}
    kenglik = {k: [] for k in usullar}
    for i, o in enumerate(boshlar):
        haq = y.iloc[o:o + h].to_numpy()
        for nom, fn in usullar.items():
            f, past, yuq = fn(i, ly.iloc[:o], h)
            mae[nom].append(np.mean(np.abs(np.exp(f) - haq)))
            ichida[nom].extend((np.log(haq) >= past) & (np.log(haq) <= yuq))
            kenglik[nom].append(np.mean(np.exp(yuq) - np.exp(past)))
    mae = {k: np.array(v) for k, v in mae.items()}

    print("=== 1. Backtest dizayni ===")
    print(f"  {n_oyna} oyna, ufq {h} oy, qadam {h} (qoplanmaydi); birinchi origin "
          f"{y.index[boshlar[0]].date()}, oxirgi {y.index[boshlar[-1]].date()}")
    print(f"  SARIMA parametrlari 12 oyda bir qayta baholandi: {len(sarima.iter)} "
          f"marta; optimizator iteratsiyalari {sarima.iter}")

    print("\n=== 2. Nuqtali bashorat: MAE (GWh), oynalar bo'yicha ===")
    print(f"  {'usul':<18} {'MAE':>7} {'min':>6} {'max':>6}")
    for nom, v in mae.items():
        print(f"  {nom:<18} {v.mean():>7.2f} {v.min():>6.1f} {v.max():>6.1f}")

    print("\n=== 3. Juftlashgan farq (qator - seasonal naive), SE oynalar bo'yicha ===")
    for nom in list(usullar)[1:]:
        d = mae[nom] - mae["seasonal naive"]
        se = d.std(ddof=1) / np.sqrt(len(d))
        hukm = ("sezilarli yaxshi" if d.mean() < -2 * se else
                "sezilarli yomon" if d.mean() > 2 * se else "farq sezilarli emas")
        print(f"  {nom:<18} {d.mean():+7.2f}  SE {se:5.2f}  -> {hukm}")

    print("\n=== 4. Qaror: eng yaxshisidan sezilarli yomon bo'lmagan eng sodda ===")
    eng = min(mae, key=lambda k: mae[k].mean())
    print(f"  eng past MAE: {eng}")
    munosib = []
    for nom in usullar:                            # soddalik tartibida
        d = mae[nom] - mae[eng]
        se = d.std(ddof=1) / np.sqrt(len(d)) if nom != eng else 0.0
        if nom == eng or d.mean() <= 2 * se:
            munosib.append(nom)
        if nom != eng:
            print(f"  {nom:<18} {d.mean():+7.2f}  SE {se:5.2f}")
    print(f"  munosiblar: {munosib} -> tanlov: {munosib[0]}")

    print("\n=== 5. 95% intervallar: haqiqiy qamrov (60 nuqta) ===")
    print(f"  {'usul':<18} {'qamrov':>7} {'o_rtacha kenglik':>17}")
    for nom in usullar:
        print(f"  {nom:<18} {np.mean(ichida[nom]):>7.3f} "
              f"{np.mean(kenglik[nom]):>17.1f}")
    n = len(ichida["SARIMA"])
    print(f"  nominal 0.95; {n} nuqtada binomial SE ~ "
          f"{np.sqrt(0.95 * 0.05 / n):.3f} (nuqtalar mustaqil emas - aslida kattaroq)")


if __name__ == "__main__":
    main()

=== 1. Backtest dizayni ===
  10 oyna, ufq 6 oy, qadam 6 (qoplanmaydi); birinchi origin 2021-01-01, oxirgi 2025-07-01
  SARIMA parametrlari 12 oyda bir qayta baholandi: 5 marta; optimizator iteratsiyalari [34, 14, 16, 26, 31]

=== 2. Nuqtali bashorat: MAE (GWh), oynalar bo'yicha ===
  usul                   MAE    min    max
  seasonal naive       53.35   23.9  114.8
  ETS(A,Ad,A)          37.27   12.7   84.1
  SARIMA               34.54   18.3   68.2
  ARIMAX (Ramazon)     34.51   19.4   75.4

=== 3. Juftlashgan farq (qator - seasonal naive), SE oynalar bo'yicha ===
  ETS(A,Ad,A)         -16.08  SE  7.00  -> sezilarli yaxshi
  SARIMA              -18.81  SE  7.07  -> sezilarli yaxshi
  ARIMAX (Ramazon)    -18.84  SE  6.18  -> sezilarli yaxshi

=== 4. Qaror: eng yaxshisidan sezilarli yomon bo'lmagan eng sodda ===
  eng past MAE: ARIMAX (Ramazon)
  seasonal naive      +18.84  SE  6.18
  ETS(A,Ad,A)          +2.76  SE  3.28
  SARIMA               +0.03  SE  3.53
  munosiblar: ['ETS(A,Ad,A)', 'SARIMA', 'ARIMAX (Ramazon)'] -> tanlov: ETS(A,Ad,A)

=== 5. 95% intervallar: haqiqiy qamrov (60 nuqta) ===
  usul                qamrov  o_rtacha kenglik
  seasonal naive       0.900             218.9
  ETS(A,Ad,A)          0.917             172.2
  SARIMA               0.967             174.9
  ARIMAX (Ramazon)     0.917             148.9
  nominal 0.95; 60 nuqtada binomial SE ~ 0.028 (nuqtalar mustaqil emas - aslida kattaroq)

eng = min([ARIMA(x, order=(1, 0, 0)).fit(), ARIMA(x, order=(0, 1, 0)).fit()],
          key=lambda r: r.aic)                                        # ⚠️
d = 1 if adfuller(x, result_object=False)[1] >= 0.05 else 0          # ✅ avval d,
eng = min([ARIMA(x, order=(p, d, q)).fit() for p in range(3)          #    keyin AIC
           for q in range(3)], key=lambda r: r.aic)

acorr_ljungbox(r.resid, lags=[24])                                    # ⚠️
acorr_ljungbox(r.resid.iloc[1 + 12:], lags=[24], model_df=p + q + P + Q)   # ✅

ExponentialSmoothing(y, trend="add").fit()                            # ⚠️
ExponentialSmoothing(y, trend="add", seasonal="mul", seasonal_periods=12).fit()  # ✅

past, yuqori = f.conf_int(alpha=0.05).T.to_numpy()                   # ⚠️ va'da
qamrov = np.mean((haq >= past) & (haq <= yuqori))                     # ✅ backtestda

SARIMAX(y, exog=haqiqiy_harorat).fit()      # backtestda haqiqiy harorat  ⚠️
SARIMAX(y, exog=harorat).fit()              # backtestda harorat BASHORATI ✅

for o in boshlar: eng = setka_qidir(y[:o])                           # ⚠️ soatlab
r = SARIMAX(y[:o], order=tanlangan, ...).fit(start_params=oldingi)   # ✅ har 12 oyda
r = r.apply(y[:o_keyingi])                                          #    oraliqda

prognoz = np.exp(f.predicted_mean)                   # ⚠️ bu mediana
prognoz = np.exp(f.predicted_mean + 0.5 * sigma2)    # ✅ taxminiy o'rtacha (lognormal)

1.  # SES da alfa = 1 bo'lsa, bashorat nimaga teng?
2.  # alfa = 0.3 bo'lsa, oxirgi kuzatuv va undan oldingisining vaznlari?
3.  # Holt ning 24 oylik bashorati qanday shaklda?
4.  # ARIMA(0,1,0) nimaning boshqa nomi?
5.  # ARIMA(1,0,0) + const ning uzoq ufqdagi bashorati?
6.  # farqlangan qatorning ACF i 1 va 12-lagda manfiy cho'qqi - qaysi model?
7.  # AIC (1,1,1) = -590, AIC (1,0,1) = -620 - (1,0,1) yaxshiroqmi?
8.  # Ljung-Box p = 0.001 nimani bildiradi?
9.  # SARIMA(1,1,1)(0,1,1)12 da Ljung-Box uchun model_df?
10. # 95% interval backtestda 0.81 qamradi - nima deysiz?
11. # ARIMAX ga backtestda haqiqiy haroratni berish - muammo nima?
12. # 60 oyna x 16 model x 1.5 s - qancha vaqt?

1.  model = ExponentialSmoothing(oylik_savdo).fit()     # mavsumli qator
    prognoz = model.forecast(12)

2.  a = ARIMA(y, order=(1, 0, 1)).fit()
    b = ARIMA(y, order=(1, 1, 1)).fit()
    eng = a if a.aic < b.aic else b

3.  lb = acorr_ljungbox(r.resid, lags=[24])            # SARIMA(1,1,1)(0,1,1)12

4.  print("95% interval:", f.conf_int(alpha=0.05))     # hisobotda kafolat sifatida

5.  for o in boshlar:                                    # 60 oyna
        eng = min(setka_16_model(ly[:o]), key=lambda r: r.aic)