# Różnicowanie pierwszego rzędu (d=1):
# Y'_t = Y_t - Y_{t-1}

# Różnicowanie drugiego rzędu (d=2):
# Y''_t = Y'_t - Y'_{t-1} = (Y_t - Y_{t-1}) - (Y_{t-1} - Y_{t-2})

import pandas as pd
import numpy as np

# Przykład w Pythonie
dane = pd.Series([100, 102, 105, 107, 110, 115])

# Różnicowanie pierwszego rzędu
diff_1 = dane.diff().dropna()
print("Różnice rzędu 1:", diff_1.values)
# [2, 3, 2, 3, 5]

# Różnicowanie drugiego rzędu
diff_2 = dane.diff().diff().dropna()
print("Różnice rzędu 2:", diff_2.values)
# [1, -1, 1, 2]

# Zazwyczaj d=1 lub d=2 wystarcza do uzyskania stacjonarności

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

# Wczytanie danych
dane = pd.read_csv('dane_czasowe.csv', parse_dates=['data'])
szereg = dane['wartosc']

# Tworzenie wykresów ACF i PACF
fig, axes = plt.subplots(1, 2, figsize=(14, 4))

# ACF - pomaga określić q
plot_acf(szereg, ax=axes[0], lags=20)
axes[0].set_title('ACF - Autokorelacja')

# PACF - pomaga określić p
plot_pacf(szereg, ax=axes[1], lags=20)
axes[1].set_title('PACF - Częściowa Autokorelacja')

plt.tight_layout()
plt.show()

# Interpretacja:
# - Słupki poza niebieskim obszarem są statystycznie istotne
# - Szukamy, gdzie wykresy "ucinają się" lub zanikają

from statsmodels.tsa.stattools import adfuller

def test_stacjonarnosci(szereg, nazwa=''):
    """
    Przeprowadza test Augmented Dickey-Fuller.
    H0: Szereg jest niestacjonarny (ma pierwiastek jednostkowy)
    H1: Szereg jest stacjonarny
    """
    wynik = adfuller(szereg, autolag='AIC')

    print(f'Test ADF dla: {nazwa}')
    print(f'Statystyka ADF: {wynik[0]:.4f}')
    print(f'p-value: {wynik[1]:.4f}')
    print(f'Wartości krytyczne:')
    for klucz, wartosc in wynik[4].items():
        print(f'   {klucz}: {wartosc:.4f}')

    if wynik[1] <= 0.05:
        print('Wniosek: Szereg jest STACJONARNY (p <= 0.05)')
    else:
        print('Wniosek: Szereg jest NIESTACJONARNY (p > 0.05)')
        print('Zalecenie: Zastosuj różnicowanie (d >= 1)')

    return wynik[1]

# Użycie
p_value = test_stacjonarnosci(szereg, 'Sprzedaż')

import pandas as pd
import numpy as np
from statsmodels.tsa.arima.model import ARIMA
import matplotlib.pyplot as plt

# Wczytanie i przygotowanie danych
dane = pd.read_csv('sprzedaz.csv', parse_dates=['data'], index_col='data')
szereg = dane['wartosc']

# Podział na zbiór treningowy i testowy
train = szereg[:-12]  # Wszystko oprócz ostatnich 12 obserwacji
test = szereg[-12:]   # Ostatnie 12 obserwacji

# Budowa modelu ARIMA(1, 1, 1)
model = ARIMA(train, order=(1, 1, 1))
wynik = model.fit()

# Podsumowanie modelu
print(wynik.summary())

# Prognoza
prognoza = wynik.forecast(steps=12)
print("Prognoza:", prognoza.values)

# Prognoza z przedziałami ufności
prognoza_obj = wynik.get_forecast(steps=12)
prognoza = prognoza_obj.predicted_mean
przedzial = prognoza_obj.conf_int()

# Wizualizacja
plt.figure(figsize=(12, 6))

# Dane historyczne
plt.plot(train.index, train, label='Dane treningowe', color='blue')
plt.plot(test.index, test, label='Dane rzeczywiste', color='green')

# Prognoza
plt.plot(test.index, prognoza, label='Prognoza ARIMA',
         color='red', linestyle='--')

# Przedział ufności
plt.fill_between(test.index,
                 przedzial.iloc[:, 0],
                 przedzial.iloc[:, 1],
                 color='red', alpha=0.2,
                 label='95% przedział ufności')

plt.title('Prognoza ARIMA')
plt.xlabel('Data')
plt.ylabel('Wartość')
plt.legend()
plt.grid(True)
plt.show()

# pip install pmdarima

import pmdarima as pm
from pmdarima import auto_arima

# Automatyczne znajdowanie najlepszych parametrów
model_auto = auto_arima(
    train,
    start_p=0, max_p=5,      # Zakres dla p
    start_q=0, max_q=5,      # Zakres dla q
    d=None,                   # Automatyczne określenie d
    seasonal=False,           # Bez sezonowości (ARIMA)
    stepwise=True,            # Szybsze przeszukiwanie
    suppress_warnings=True,
    information_criterion='aic',  # Kryterium wyboru
    trace=True                # Pokazuj postęp
)

# Najlepszy model
print(f"Najlepszy model: ARIMA{model_auto.order}")
print(f"AIC: {model_auto.aic():.2f}")

# Prognoza
prognoza = model_auto.predict(n_periods=12)

# Diagnostyka
model_auto.plot_diagnostics(figsize=(12, 8))
plt.show()

# Reszty to różnice między wartościami rzeczywistymi a prognozami
reszty = wynik.resid

# Wykresy diagnostyczne
fig, axes = plt.subplots(2, 2, figsize=(12, 8))

# 1. Histogram reszt - powinny być normalnie rozłożone
axes[0, 0].hist(reszty, bins=25, density=True)
axes[0, 0].set_title('Histogram reszt')

# 2. Q-Q plot - punkty powinny leżeć na linii
from scipy import stats
stats.probplot(reszty, plot=axes[0, 1])
axes[0, 1].set_title('Q-Q Plot')

# 3. Reszty w czasie - powinny być losowe
axes[1, 0].plot(reszty)
axes[1, 0].axhline(y=0, color='r', linestyle='--')
axes[1, 0].set_title('Reszty w czasie')

# 4. ACF reszt - nie powinno być znaczących korelacji
plot_acf(reszty, ax=axes[1, 1], lags=20)
axes[1, 1].set_title('ACF reszt')

plt.tight_layout()
plt.show()

from statsmodels.stats.diagnostic import acorr_ljungbox
from scipy.stats import shapiro, jarque_bera

def diagnostyka_reszt(reszty):
    """Przeprowadza testy diagnostyczne na resztach modelu."""

    print("=== DIAGNOSTYKA RESZT ===\n")

    # 1. Test Ljung-Box - czy reszty są białym szumem
    lb_test = acorr_ljungbox(reszty, lags=[10, 20], return_df=True)
    print("Test Ljung-Box (H0: brak autokorelacji):")
    print(lb_test)
    print()

    # 2. Test Shapiro-Wilk - normalność
    stat, p = shapiro(reszty[:5000] if len(reszty) > 5000 else reszty)
    print(f"Test Shapiro-Wilk (H0: rozkład normalny):")
    print(f"Statystyka: {stat:.4f}, p-value: {p:.4f}")
    print(f"Wniosek: {'Normalny' if p > 0.05 else 'Nienormalny'}\n")

    # 3. Test Jarque-Bera - normalność
    stat, p = jarque_bera(reszty)
    print(f"Test Jarque-Bera:")
    print(f"Statystyka: {stat:.4f}, p-value: {p:.4f}")
    print(f"Wniosek: {'Normalny' if p > 0.05 else 'Nienormalny'}\n")

    # 4. Statystyki opisowe
    print(f"Średnia reszt: {reszty.mean():.4f} (powinna być ~0)")
    print(f"Odchylenie std: {reszty.std():.4f}")

diagnostyka_reszt(reszty)

from sklearn.metrics import mean_absolute_error, mean_squared_error
import numpy as np

def metryki_prognozy(rzeczywiste, prognoza):
    """Oblicza popularne metryki błędów prognozy."""

    # MAE - Mean Absolute Error (średni błąd bezwzględny)
    mae = mean_absolute_error(rzeczywiste, prognoza)

    # MSE - Mean Squared Error (średni błąd kwadratowy)
    mse = mean_squared_error(rzeczywiste, prognoza)

    # RMSE - Root Mean Squared Error
    rmse = np.sqrt(mse)

    # MAPE - Mean Absolute Percentage Error
    mape = np.mean(np.abs((rzeczywiste - prognoza) / rzeczywiste)) * 100

    print("=== METRYKI PROGNOZY ===")
    print(f"MAE:  {mae:.2f}")
    print(f"MSE:  {mse:.2f}")
    print(f"RMSE: {rmse:.2f}")
    print(f"MAPE: {mape:.2f}%")

    return {'MAE': mae, 'MSE': mse, 'RMSE': rmse, 'MAPE': mape}

# Użycie
metryki = metryki_prognozy(test.values, prognoza.values)

# Interpretacja MAPE:
# < 10% - Doskonała prognoza
# 10-20% - Dobra prognoza
# 20-50% - Akceptowalna prognoza
# > 50% - Słaba prognoza

from statsmodels.tsa.statespace.sarimax import SARIMAX

# Dane miesięczne z sezonowością roczną
# SARIMA(1, 1, 1)(1, 1, 1, 12)
model = SARIMAX(
    train,
    order=(1, 1, 1),              # (p, d, q)
    seasonal_order=(1, 1, 1, 12), # (P, D, Q, s)
    enforce_stationarity=False,
    enforce_invertibility=False
)

wynik = model.fit(disp=False)
print(wynik.summary())

# Auto SARIMA z pmdarima
model_auto = pm.auto_arima(
    train,
    seasonal=True,      # Włącz sezonowość
    m=12,               # Okres sezonowy
    start_P=0, max_P=2,
    start_Q=0, max_Q=2,
    D=1,                # Sezonowe różnicowanie
    trace=True,
    suppress_warnings=True
)

print(f"Najlepszy SARIMA: {model_auto.order}x{model_auto.seasonal_order}")

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.stattools import adfuller
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

# 1. Wczytanie danych
url = 'https://raw.githubusercontent.com/...'
dane = pd.read_csv(url, parse_dates=['Month'], index_col='Month')
szereg = dane['Passengers']

# 2. Wizualizacja
plt.figure(figsize=(12, 4))
plt.plot(szereg)
plt.title('Liczba pasażerów lotniczych (1949-1960)')
plt.xlabel('Data')
plt.ylabel('Tysiące pasażerów')
plt.show()

# 3. Statystyki opisowe
print(szereg.describe())
print(f"\nOkres: {szereg.index.min()} do {szereg.index.max()}")
print(f"Liczba obserwacji: {len(szereg)}")

# 4. Test stacjonarności
wynik_adf = adfuller(szereg)
print(f"ADF Statystyka: {wynik_adf[0]:.4f}")
print(f"p-value: {wynik_adf[1]:.4f}")
# Szereg niestacjonarny - potrzebne różnicowanie

# 5. Transformacja logarytmiczna (stabilizacja wariancji)
szereg_log = np.log(szereg)

# 6. Różnicowanie
szereg_diff = szereg_log.diff().dropna()

# 7. Ponowny test
wynik_adf2 = adfuller(szereg_diff)
print(f"\nPo różnicowaniu:")
print(f"p-value: {wynik_adf2[1]:.4f}")
# Teraz szereg jest stacjonarny

# 8. Wykresy ACF i PACF
fig, axes = plt.subplots(1, 2, figsize=(14, 4))
plot_acf(szereg_diff, ax=axes[0], lags=30)
plot_pacf(szereg_diff, ax=axes[1], lags=30)
plt.show()

# Na podstawie wykresów: p=2, q=2 (zanikanie w obu)

# 9. Podział danych
train = szereg_log[:-24]  # Treningowe
test = szereg_log[-24:]   # Testowe (2 lata)

# 10. Budowa modelu ARIMA(2, 1, 2)
model = ARIMA(train, order=(2, 1, 2))
wynik = model.fit()

# 11. Podsumowanie
print(wynik.summary())

# 12. Prognoza
prognoza_log = wynik.get_forecast(steps=24)
prognoza = np.exp(prognoza_log.predicted_mean)  # Odwrócenie log
rzeczywiste = np.exp(test)

# 13. Metryki
mape = np.mean(np.abs((rzeczywiste - prognoza) / rzeczywiste)) * 100
print(f"\nMAPE: {mape:.2f}%")

# 14. Wizualizacja
plt.figure(figsize=(12, 6))
plt.plot(np.exp(train), label='Treningowe')
plt.plot(rzeczywiste, label='Rzeczywiste')
plt.plot(prognoza, label='Prognoza', linestyle='--')
plt.legend()
plt.show()