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

import math
from typing import Tuple, Dict, List
from dataclasses import dataclass
import numpy as np


x = np.array([0.3, 0.6, 1.1, 0.2, 0.8, 1.6, 1.5, 0.9, 0.9, 0.8], dtype=float)
y = np.array([2.8, 4.3, 6.3, 2.9, 5.3, 8.3, 8.0, 5.7, 5.4, 5.2], dtype=float)


ALPHA_SW = 0.05


def ols_fit_and_residuals(x: np.ndarray, y: np.ndarray) -> Tuple[float, float, np.ndarray, np.ndarray]:

    x = np.asarray(x, float)
    y = np.asarray(y, float)
    x_bar = float(np.mean(x))
    y_bar = float(np.mean(y))
    Sxx = float(np.sum((x - x_bar) ** 2))
    Sxy = float(np.sum((x - x_bar) * (y - y_bar)))
    a1 = Sxy / Sxx
    a0 = y_bar - a1 * x_bar
    y_hat = a0 + a1 * x
    u = y - y_hat
    return a0, a1, y_hat, u


def central_moments(u: np.ndarray) -> Tuple[float, float, float, float]:

    u = np.asarray(u, float)
    mu = float(np.mean(u))
    d = u - mu
    m2 = float(np.mean(d ** 2))
    m3 = float(np.mean(d ** 3))
    m4 = float(np.mean(d ** 4))
    return mu, m2, m3, m4


def skewness_kurtosis(u: np.ndarray) -> Tuple[float, float]:

    _, m2, m3, m4 = central_moments(u)
    if m2 <= 0:
        return float('nan'), float('nan')
    g1 = m3 / (m2 ** 1.5)
    g2 = m4 / (m2 ** 2) - 3.0
    return g1, g2


def adjusted_moment_tests(u: np.ndarray):

    n = float(len(u))
    g1, g2 = skewness_kurtosis(u)
    A_hat = math.sqrt(n * (n - 1.0)) / (n - 2.0) * g1
    E_hat = ((n - 1.0) / ((n - 2.0) * (n - 3.0))) * ((n + 1.0) * g2 + 6.0)
    var_A = 6.0 * n * (n - 1.0) / ((n + 1.0) * (n + 3.0) * (n - 2.0))
    var_E = 24.0 * n * (n - 1.0) ** 2 / ((n + 5.0) * (n + 3.0) * (n - 2.0) * (n - 3.0))
    two_sA = 2.0 * math.sqrt(var_A)
    two_sE = 2.0 * math.sqrt(var_E)
    decision = (abs(A_hat) <= two_sA) and (abs(E_hat) <= two_sE)
    return A_hat, E_hat, two_sA, two_sE, decision


_A_COEFFS_N10 = np.array([0.5739, 0.3291, 0.2141, 0.1224, 0.0399], dtype=float)
_W_CRIT_N10: Dict[float, float] = {0.10: 0.869, 0.05: 0.842, 0.02: 0.806, 0.01: 0.781}


def shapiro_wilk_w_n10(sample: np.ndarray) -> Tuple[float, float, List[Tuple[int, float, float, float, float, float]]]:

    x = np.sort(np.asarray(sample, float))
    n = x.size
    if n != 10:
        raise ValueError("Shapiro–Wilk (ця реалізація) підтримує лише n = 10.")
    S2 = float(np.sum((x - x.mean()) ** 2))
    rows: List[Tuple[int, float, float, float, float, float]] = []
    b = 0.0
    k = 5
    for i in range(k):
        left = float(x[i])
        right = float(x[-(i + 1)])
        diff = right - left
        ai = float(_A_COEFFS_N10[i])
        term = ai * diff
        b += term
        rows.append((i + 1, left, right, diff, ai, term))
    W = (b ** 2) / S2 if S2 > 0 else float('nan')
    return W, S2, rows


def shapiro_wilk_decision_n10(sample: np.ndarray, alpha: float = ALPHA_SW) -> Tuple[float, bool, float, float, List[Tuple[int, float, float, float, float, float]]]:
    if alpha not in _W_CRIT_N10:
        raise ValueError(f"Для n=10 доступні α: {sorted(_W_CRIT_N10.keys())}")
    W, S2, rows = shapiro_wilk_w_n10(sample)
    Wcrit = _W_CRIT_N10[alpha]
    return W, (W >= Wcrit), Wcrit, S2, rows


@dataclass
class NormalityResults:
    a0: float
    a1: float
    A_hat: float
    E_hat: float
    two_sA: float
    two_sE: float
    approx_moment_ok: bool
    W: float
    Wcrit: float
    W_ok: bool
    S2: float
    rows: List[Tuple[int, float, float, float, float, float]]
    u: np.ndarray
    u_sorted: np.ndarray
    u_mean: float


def run_full_check_and_print(x: np.ndarray, y: np.ndarray, alpha_sw: float = ALPHA_SW) -> NormalityResults:
    a0, a1, y_hat, u = ols_fit_and_residuals(x, y)

    A_hat, E_hat, two_sA, two_sE, approx_ok = adjusted_moment_tests(u)

    W, sw_ok, Wcrit, S2, rows = shapiro_wilk_decision_n10(u, alpha=alpha_sw)

    u_sorted = np.sort(u)
    u_mean = float(np.mean(u))

    gamma = 1.0 - alpha_sw

    print("\n=== OLS модель ===")
    print(f"a0 = {a0:.6f}, a1 = {a1:.6f}")

    print("\n=== Залишки U ===")
    for idx, val in enumerate(u, start=1):
        print(f"u[{idx:>2}] = {val: .6f}")
    print(f"Середнє ū = {u_mean:.6f}")

    print("\n=== Відсортовані залишки u_(i) (для SW) ===")
    for idx, val in enumerate(u_sorted, start=1):
        print(f"u_({idx:>2}) = {val: .6f}")

    print("\n=== Наближений критерій: Â*, Ê* ===")
    print(f"A^* = {A_hat:.4f},  2σ_A^* = {two_sA:.4f}  =>  {'OK' if abs(A_hat) <= two_sA else 'NOT OK'}")
    print(f"E^* = {E_hat:.4f},  2σ_E^* = {two_sE:.4f}  =>  {'OK' if abs(E_hat) <= two_sE else 'NOT OK'}")
    concl = "не відхиляємо H0 (нормальність)" if approx_ok else "відхиляємо H0 (є відхилення від нормальності)"
    print(f"Висновок (наближений критерій): {concl}")

    print("\n=== Шапіро–Вілк (n=10) — проміжні розрахунки ===")
    print(f"S^2 = sum (u_(i) − ū)^2 = {S2:.6f}")
    header = "  i    u_(i)        u_(n−i+1)    diff=(right−left)      a_i         a_i*diff"
    print(header)
    print("-" * len(header))
    b = 0.0
    for (i, left, right, diff, ai, term) in rows:
        b += term
        print(f"{i:3d}  {left: .6f}    {right: .6f}      {diff: .6f}        {ai: .4f}       {term: .6f}")
    print(f"b = sum(a_i * (u_(n−i+1) − u_(i))) = {b:.6f}")
    print(f"b^2 = {b*b:.6f}")

    print("\n=== Підсумок SW ===")
    print(f"W = b^2 / S^2 = {W:.6f}")
    print(f"Рівень довіри γ = {gamma:.2f} (α = {alpha_sw:.2f});  W_crit(γ) = {Wcrit:.3f}")
    print(f"Висновок: {'не відхиляємо H0' if sw_ok else 'відхиляємо H0'}")

    return NormalityResults(
        a0=a0, a1=a1,
        A_hat=A_hat, E_hat=E_hat, two_sA=two_sA, two_sE=two_sE, approx_moment_ok=approx_ok,
        W=W, Wcrit=Wcrit, W_ok=sw_ok, S2=S2, rows=rows,
        u=u, u_sorted=u_sorted, u_mean=u_mean,
    )


if __name__ == "__main__":
    _ = run_full_check_and_print(x, y, alpha_sw=ALPHA_SW)