import math
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)


n = x.size
X = np.column_stack([np.ones(n), x]) 
beta = np.linalg.inv(X.T @ X) @ (X.T @ y)
y_hat = X @ beta
u = y - y_hat


den = float(np.sum(u ** 2))
dw = float(np.sum(np.diff(u) ** 2) / den) if den > 0 else float("nan")


def lag1_corr(u_vec: np.ndarray) -> float:
    u_vec = np.asarray(u_vec).ravel()
    if u_vec.size < 2:
        return float("nan")
    a = u_vec[1:]    
    b = u_vec[:-1]         
    da = a - np.mean(a)
    db = b - np.mean(b)
    num = float(np.sum(da * db))
    den = math.sqrt(float(np.sum(da * da)) * float(np.sum(db * db)))
    return num / den if den > 0 else float("nan")

r1 = lag1_corr(u)


d_L, d_U = 0.88, 1.32


def dw_status_label(d: float) -> str:
    if np.isnan(d):
        return "error"
    if d < d_L:
        return "reject_pos"  
    if d_U < d < (4 - d_U):
        return "accept"    
    if (4 - d_L) < d < 4:
        return "reject_neg"   
    if (d_L <= d <= d_U) or ((4 - d_U) <= d <= (4 - d_L)):
        return "inconclusive" 
    return "unknown"

dw_stat = dw_status_label(dw)

if dw_stat == "reject_pos":
    dw_verdict = "Відхиляємо H0: позитивна автокореляція (ρ>0)."
elif dw_stat == "accept":
    dw_verdict = "H0 не відхиляємо: автокореляції не виявлено."
elif dw_stat == "reject_neg":
    dw_verdict = "Відхиляємо H0: негативна автокореляція (ρ<0)."
elif dw_stat == "inconclusive":
    dw_verdict = "Сіра зона DW: остаточного рішення за DW немає; див. тест Бреуша–Годфрі."
elif dw_stat == "error":
    dw_verdict = "DW не обчислено (перевірте дані)."
else:
    dw_verdict = "Невідома ситуація для DW."


u_lag = u[:-1]
u_cur = u[1:]
rho_hat = float(np.sum(u_lag * u_cur) / np.sum(u_lag ** 2))
eps = u_cur - rho_hat * u_lag
df = n - 2
SSE = float(np.sum(eps ** 2))
S_eps2 = SSE / df
S_rho = math.sqrt(S_eps2 / np.sum(u_lag ** 2))
t_stat = rho_hat / S_rho if S_rho > 0 else float("inf")


t_crit = 2.306
bg_accept = (abs(t_stat) <= t_crit) 
bg_verdict = ("H0 не відхиляємо: автокореляції не виявлено."
              if bg_accept else
              "Відхиляємо H0: автокореляція є.")


if dw_stat == "inconclusive":
    final_concl = ("Остаточний висновок за Бреуша–Годфрі: "
                   f"{'автокореляції не виявлено' if bg_accept else 'виявлено автокореляцію'} "
                   f"(ρ̂ = {rho_hat:.6f}).")
elif dw_stat == "accept" and bg_accept:
    final_concl = "Автокореляції залишків не виявлено (обидва тести не відхилили H0)."
elif (dw_stat in ("reject_pos", "reject_neg")) and (not bg_accept):
    direction = "позитивну" if rho_hat > 0 else "негативну"
    final_concl = f"Виявлено {direction} автокореляцію AR(1) (обидва тести вказують на відхилення H0; ρ̂ = {rho_hat:.6f})."
elif (dw_stat in ("reject_pos", "reject_neg")) and bg_accept:
    direction_dw = "позитивну" if dw_stat == "reject_pos" else "негативну"
    final_concl = (f"Тести суперечливі: DW вказує на {direction_dw} автокореляцію, "
                   f"але Бреуша–Годфрі її не підтвердив (ρ̂ = {rho_hat:.6f}). "
                   f"Приймається обережний висновок за BG: на рівні α=0.05 "
                   f"автокореляцію не підтверджено.")
elif dw_stat == "accept" and (not bg_accept):
    direction = "позитивну" if rho_hat > 0 else "негативну"
    final_concl = (f"Тести суперечливі: DW не виявив автокореляцію, "
                   f"але Бреуша–Годфрі відхилив H0. Вважається, що є {direction} автокореляція "
                   f"(α=0.05; ρ̂ = {rho_hat:.6f}).")
else:
    final_concl = "Неможливо сформувати однозначний висновок — перевірте дані та критичні значення."


print("=== Перевірка передумови 3: відсутність автокореляції залишків ===")
print(f"n = {n}, модель: з константою (m = 2)")
print()
print("--- МНК-оцінки ---")
print(f"β̂0 = {beta[0]:.6f}, β̂1 = {beta[1]:.6f}")
print(f"Модель: ŷ = {beta[0]:.6f} + {beta[1]:.6f}·x")
print()
print("--- Тест Дарбіна–Уотсона ---")
print(f"d = {dw:.6f}")
print(f"Довідково: r(1) (ручний розрах.) ≈ {r1:.6f}, d ≈ 2(1 - r(1)) ≈ {2*(1-r1):.6f}")
print(f"Критичні межі (α=0.05, n=10, m=2): d_L = {d_L}, d_U = {d_U}, "
      f"інтервал прийняття H0: ({d_U:.2f}, {4-d_U:.2f})")
print(f"Вердикт DW: {dw_verdict}")
print()
print("--- Тест Бреуша–Годфрі (AR(1)) ---")
print(f"rho_hat = {rho_hat:.6f}")
print(f"t = {t_stat:.6f}, df = {df}, t_crit(α/2) ≈ {t_crit:.3f}")
print(f"Вердикт BG(AR1): {bg_verdict}")
print()
print("--- Підсумковий висновок ---")
print(final_concl)