Прикладна_економетрика_ЛР№11_Мірошниченко

import math
import matplotlib.pyplot as plt

DEFAULT_SEASON_LEN = 4
MAX_LAG_FOR_CORRELOGRAM = 12
DATASETS = {
    "demo": {
        "name": "Демонстраційний приклад (витрати на рекламу, квартали)",
        "y": [375, 371, 869, 1015, 357, 471, 992, 1020, 390, 355, 992, 905, 461, 454, 920, 927],
        "season_len": 4
    },
    "variant2": {
        "name": "Варіант 2",
        "y": [5.8, 4.5, 5.1, 9.1, 7.0, 5.0, 6.0, 10.1, 7.9, 5.5, 6.3, 10.8, 9.0, 6.5, 7.0, 11.1],
        "season_len": 4
    }
}

def title(s):
    print("\n" + "=" * 70)
    print(s)
    print("=" * 70)

def fmt(x, d=4):
    if x is None:
        return "-"
    try:
        return f"{float(x):.{d}f}"
    except Exception:
        return str(x)

def fmt_int_if_close(x, d=4):
    try:
        xf = float(x)
        if abs(xf - round(xf)) < 1e-12:
            return str(int(round(xf)))
        return f"{xf:.{d}f}"
    except Exception:
        return str(x)

def print_table(headers, rows):
    cols = len(headers)
    widths = [len(h) for h in headers]
    for r in rows:
        for j in range(cols):
            widths[j] = max(widths[j], len(str(r[j])))

    def line():
        s = "+"
        for w in widths:
            s += "-" * (w + 2) + "+"
        return s

    print(line())
    hrow = "|"
    for j, h in enumerate(headers):
        hrow += " " + h.ljust(widths[j]) + " |"
    print(hrow)
    print(line())

    for r in rows:
        row_s = "|"
        for j in range(cols):
            row_s += " " + str(r[j]).ljust(widths[j]) + " |"
        print(row_s)
    print(line())

def season_index(t1, p):
    return ((t1 - 1) % p) + 1

def mean(arr):
    return sum(arr) / len(arr) if arr else 0.0

def pearson_autocorr_shifted(y, lag):
    n = len(y)
    if lag <= 0 or lag >= n:
        return 0.0

    y1 = y[:-lag]
    y2 = y[lag:]
    m1 = mean(y1)
    m2 = mean(y2)

    num = 0.0
    s1 = 0.0
    s2 = 0.0
    for a, b in zip(y1, y2):
        da = a - m1
        db = b - m2
        num += da * db
        s1 += da * da
        s2 += db * db

    den = math.sqrt(s1 * s2)
    return num / den if den != 0 else 0.0

def moving_average(y, p):
    n = len(y)
    return [sum(y[i:i+p]) / p for i in range(n - p + 1)]

def centered_moving_average(ma):
    return [(ma[i] + ma[i+1]) / 2 for i in range(len(ma) - 1)]

def estimate_seasonals(y, p, kind):
    ma = moving_average(y, p)
    if p % 2 == 0:
        cma = centered_moving_average(ma)
        start_t = p // 2 + 1
    else:
        cma = ma[:]
        start_t = (p + 1) // 2

    t_cma = {start_t + i: cma[i] for i in range(len(cma))}

    raw_by_season = {i: [] for i in range(1, p + 1)}
    raw_by_t = {}

    for t, cm in t_cma.items():
        yt = y[t - 1]
        if kind == "additive":
            s = yt - cm
        else:
            s = yt / cm if cm != 0 else 0.0
        raw_by_t[t] = s
        raw_by_season[season_index(t, p)].append(s)

    avg = {}
    for i in range(1, p + 1):
        vals = raw_by_season[i]
        avg[i] = sum(vals) / len(vals) if vals else 0.0

    if kind == "additive":
        k = sum(avg.values()) / p
        seasonal = {i: avg[i] - k for i in avg}
    else:
        ssum = sum(avg.values())
        k = (p / ssum) if ssum != 0 else 1.0
        seasonal = {i: avg[i] * k for i in avg}

    return {
        "ma": ma,
        "cma": cma,
        "t_cma": t_cma,
        "raw_by_t": raw_by_t,
        "avg_seasonal": avg,
        "seasonal": seasonal
    }

def deseasonalize(y, seasonal, p, kind):
    tau = []
    for t in range(1, len(y) + 1):
        s = seasonal[season_index(t, p)]
        if kind == "additive":
            tau.append(y[t - 1] - s)
        else:
            tau.append(y[t - 1] / s if s != 0 else 0.0)
    return tau

def fit_linear_trend(tau):
    n = len(tau)
    t_vals = list(range(1, n + 1))
    sum_t = sum(t_vals)
    sum_t2 = sum(v * v for v in t_vals)
    sum_tau = sum(tau)
    sum_ttau = sum(t_vals[i] * tau[i] for i in range(n))

    den = n * sum_t2 - sum_t * sum_t
    a1 = (n * sum_ttau - sum_t * sum_tau) / den if den != 0 else 0.0
    a0 = (sum_tau / n) - a1 * (sum_t / n)
    return a0, a1

def build_yhat(a0, a1, seasonal, n, p, kind):
    yhat = []
    for t in range(1, n + 1):
        trend = a0 + a1 * t
        s = seasonal[season_index(t, p)]
        if kind == "additive":
            yhat.append(trend + s)
        else:
            yhat.append(trend * s)
    return yhat

def r_squared(y, yhat):
    ybar = mean(y)
    sst = sum((v - ybar) ** 2 for v in y)
    sse = sum((y[i] - yhat[i]) ** 2 for i in range(len(y)))
    r2 = 1.0 - (sse / sst) if sst != 0 else 0.0
    return r2, sse, sst

def forecast_2(a0, a1, seasonal, n, p, kind):
    out = []
    for h in (1, 2):
        t = n + h
        trend = a0 + a1 * t
        s = seasonal[season_index(t, p)]
        yh = (trend + s) if kind == "additive" else (trend * s)
        out.append((t, trend, s, yh))
    return out

def plot_series(y, name):
    t = list(range(1, len(y) + 1))
    plt.figure()
    plt.plot(t, y, marker="o")
    plt.title(f"Кореляційне поле (t–y): {name}")
    plt.xlabel("t")
    plt.ylabel("y")
    plt.grid(True)

def plot_correlogram(lags, rvals, name):
    plt.figure()
    plt.bar(lags, rvals)
    plt.title(f"Корелограма (autocorr): {name}")
    plt.xlabel("lag l")
    plt.ylabel("r_l")
    plt.grid(True, axis="y")

def plot_fact_vs_theor(y, yhat, name):
    t = list(range(1, len(y) + 1))
    plt.figure()
    plt.plot(t, y, marker="o", label="Факт y_t")
    plt.plot(t, yhat, marker="x", label="Теорія ŷ_t")
    plt.title(name)
    plt.xlabel("t")
    plt.ylabel("y")
    plt.grid(True)
    plt.legend()

def run_dataset(ds_key):
    ds = DATASETS[ds_key]
    name = ds["name"]
    y = ds["y"]
    p = ds.get("season_len", DEFAULT_SEASON_LEN)
    n = len(y)

    title(f"АНАЛІЗ ДАНИХ: {name}")
    print(f"n = {n}, сезон p = {p}")

    print("\n--- 1) Вихідні дані ---")
    rows = [[t, fmt_int_if_close(y[t-1], 4)] for t in range(1, n + 1)]
    print_table(["t", "y_t"], rows)

    print("\n--- 2) Кореляційне поле ---")
    plot_series(y, name)

    print("\n--- 3) Автокореляція ---")
    max_lag = min(MAX_LAG_FOR_CORRELOGRAM, n - 1)
    lags = list(range(1, max_lag + 1))
    rvals = [pearson_autocorr_shifted(y, lag) for lag in lags]

    rows = [[lag, fmt(r, 4), fmt(r, 6)] for lag, r in zip(lags, rvals)]
    print_table(["Lag", "r (4 знаки)", "r (6 знаків)"], rows)
    plot_correlogram(lags, rvals, name)

    title("АДИТИВНА МОДЕЛЬ")
    add = estimate_seasonals(y, p, kind="additive")

    print("\n--- 4.1) MA(p) та CMA (для p=4) ---")
    ma_rows = [[i, fmt(add["ma"][i-1], 2)] for i in range(1, len(add["ma"]) + 1)]
    print_table([f"MA індекс", f"MA({p})"], ma_rows)

    start_t = p // 2 + 1 if p % 2 == 0 else (p + 1) // 2
    cma_rows = [[start_t + i, fmt(add["cma"][i], 2)] for i in range(len(add["cma"]))]
    print_table(["t (для CMA)", "CMA_t"], cma_rows)

    print("\n--- 4.2) Оцінки сезонності: s_t = y_t - CMA_t ---")
    raw_rows = [[t, fmt(add["raw_by_t"][t], 2)] for t in sorted(add["raw_by_t"].keys())]
    print_table(["t", "s_t"], raw_rows)

    print("\n--- 4.3) Сезонні середні та скориговані S_i (сума = 0) ---")
    s_rows = []
    for i in range(1, p + 1):
        s_rows.append([i, fmt(add["avg_seasonal"][i], 2), fmt(add["seasonal"][i], 2)])
    print_table(["сезон i", "середнє (до корекц.)", "S_i (після корекц.)"], s_rows)

    print("S_i:")
    for i in range(1, p + 1):
        print(f"  S{i} = {fmt(add['seasonal'][i], 4)}")
    print("Контроль суми S_i =", fmt(sum(add["seasonal"].values()), 6))

    tau_add = deseasonalize(y, add["seasonal"], p, kind="additive")
    a0_add, a1_add = fit_linear_trend(tau_add)

    print("\n--- 4.4) Тренд τ̂(t)=a0+a1·t ---")
    print(f"a0 = {fmt(a0_add, 4)}")
    print(f"a1 = {fmt(a1_add, 7)}")
    print(f"Рівняння: τ̂(t) = {fmt(a0_add, 4)} + {fmt(a1_add, 7)}·t")

    yhat_add = build_yhat(a0_add, a1_add, add["seasonal"], n, p, kind="additive")
    r2_add, sse_add, sst_add = r_squared(y, yhat_add)

    print("\n--- 4.5) Якість (R²) ---")
    print(f"SSE = {fmt(sse_add, 4)}")
    print(f"SST = {fmt(sst_add, 4)}")
    print(f"R²  = {fmt(r2_add, 4)}  (тобто {fmt(r2_add * 100.0, 2)}%)")

    print("\n--- 4.6) Прогноз на t=[n+1, n+2] (адитивна) ---")
    fc = forecast_2(a0_add, a1_add, add["seasonal"], n, p, kind="additive")
    fc_rows = [[t, fmt(tr, 4), fmt(s, 4), fmt(yh, 4)] for (t, tr, s, yh) in fc]
    print_table(["t", "τ̂_t", "S_i", "ŷ_t прогноз"], fc_rows)

    plot_fact_vs_theor(y, yhat_add, f"Адитивна модель: факт vs теорія ({name})")

    title("МУЛЬТИПЛІКАТИВНА МОДЕЛЬ")
    mul = estimate_seasonals(y, p, kind="multiplicative")

    print("\n--- 5.1) Оцінки сезонності: s_t = y_t / CMA_t ---")
    raw_rows = [[t, fmt(mul["raw_by_t"][t], 4)] for t in sorted(mul["raw_by_t"].keys())]
    print_table(["t", "s_t"], raw_rows)

    print("\n--- 5.2) Сезонні середні та скориговані S_i (сума = p) ---")
    s_rows = []
    for i in range(1, p + 1):
        s_rows.append([i, fmt(mul["avg_seasonal"][i], 4), fmt(mul["seasonal"][i], 4)])
    print_table(["сезон i", "середнє (до корекц.)", "S_i (після корекц.)"], s_rows)
    print("Контроль суми S_i =", fmt(sum(mul["seasonal"].values()), 6))

    tau_mul = deseasonalize(y, mul["seasonal"], p, kind="multiplicative")
    a0_mul, a1_mul = fit_linear_trend(tau_mul)

    print("\n--- 5.3) Тренд τ̂(t)=a0+a1·t ---")
    print(f"a0 = {fmt(a0_mul, 4)}")
    print(f"a1 = {fmt(a1_mul, 7)}")
    print(f"Рівняння: τ̂(t) = {fmt(a0_mul, 4)} + {fmt(a1_mul, 7)}·t")

    yhat_mul = build_yhat(a0_mul, a1_mul, mul["seasonal"], n, p, kind="multiplicative")
    r2_mul, sse_mul, sst_mul = r_squared(y, yhat_mul)

    print("\n--- 5.4) Якість (R²) ---")
    print(f"SSE = {fmt(sse_mul, 4)}")
    print(f"SST = {fmt(sst_mul, 4)}")
    print(f"R²  = {fmt(r2_mul, 4)}  (тобто {fmt(r2_mul * 100.0, 2)}%)")

    print("\n--- 5.5) Прогноз на t=[n+1, n+2] (мультиплікативна) ---")
    fc = forecast_2(a0_mul, a1_mul, mul["seasonal"], n, p, kind="multiplicative")
    fc_rows = [[t, fmt(tr, 4), fmt(s, 4), fmt(yh, 4)] for (t, tr, s, yh) in fc]
    print_table(["t", "τ̂_t", "S_i", "ŷ_t прогноз"], fc_rows)

    plot_fact_vs_theor(y, yhat_mul, f"Мультиплікативна модель: факт vs теорія ({name})")

    title("ПОРІВНЯННЯ МОДЕЛЕЙ (за R²)")
    print(f"R² (адитивна)         = {fmt(r2_add, 4)}")
    print(f"R² (мультиплікативна) = {fmt(r2_mul, 4)}")


def main():
    run_dataset("demo")
    run_dataset("variant2")
    plt.show()

if __name__ == "__main__":
    main()