import warnings; warnings.filterwarnings("ignore")  # Suppress warnings for clean output

# Let us check our Python and library versions
import sys
import pandas as pd
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt

print("Python version:", sys.version)
print("pandas version:", pd.__version__)
print("statsmodels version:", sm.__version__)
print("numpy version:", np.__version__)

# Data setup: let us load the Airline Passengers dataset
url = "https://raw.githubusercontent.com/jbrownlee/Datasets/master/airline-passengers.csv"
df = pd.read_csv(url, parse_dates=["Month"], index_col="Month")

print("Shape:", df.shape)
df.head()

# Let us plot the time series to see what we are working with
df.plot(y="Passengers", legend=False, figsize=(10,4), title="Monthly Airline Passengers")
plt.ylabel("Number of Passengers")
plt.xlabel("Date")
plt.show()

# Decompose time series into trend, seasonality, and residuals using statsmodels
from statsmodels.tsa.seasonal import seasonal_decompose

result = seasonal_decompose(df["Passengers"], model="multiplicative", period=12)
result.plot()
plt.show()

# Let us check for stationarity using rolling mean and rolling standard deviation
rollmean = df["Passengers"].rolling(window=12).mean()
rollstd = df["Passengers"].rolling(window=12).std()

plt.figure(figsize=(10,4))
plt.plot(df.index, df["Passengers"], label="Original")
plt.plot(df.index, rollmean, label="Rolling Mean")
plt.plot(df.index, rollstd, label="Rolling Std")
plt.legend()
plt.show()

# Test stationarity using Augmented Dickey-Fuller test
from statsmodels.tsa.stattools import adfuller

result = adfuller(df["Passengers"])
print("ADF Statistic:", result[0])
print("p-value:", result[1])
if result[1] < 0.05:
    print("This series is likely stationary.")
else:
    print("This series is likely not stationary.")
    

# Let us difference the data to make it stationary
df_diff = df["Passengers"].diff().dropna()

plt.figure(figsize=(10,4))
plt.plot(df_diff)
plt.title("First Difference of Monthly Passengers")
plt.show()

# Build a simple ARIMA model using statsmodels
from statsmodels.tsa.arima.model import ARIMA

train = df["Passengers"][:-12]
test = df["Passengers"][-12:]

model = ARIMA(train, order=(1,1,1))
fit = model.fit()
forecast = fit.forecast(steps=12)

plt.figure(figsize=(10,4))
plt.plot(train.index, train, label="Train", color="blue")
plt.plot(test.index, test, label="Test", color="green")
plt.plot(test.index, forecast, label="Forecast", color="red")
plt.legend()
plt.title("ARIMA Forecast vs Actuals")
plt.show()

# Fit SARIMA model, including yearly seasonality
from statsmodels.tsa.statespace.sarimax import SARIMAX

sarima_model = SARIMAX(train, order=(1,1,1), seasonal_order=(1,1,1,12))
sarima_fit = sarima_model.fit(disp=False)
sarima_forecast = sarima_fit.forecast(steps=12)

plt.figure(figsize=(10,4))
plt.plot(train.index, train, label="Train")
plt.plot(test.index, test, label="Test")
plt.plot(test.index, sarima_forecast, label="SARIMA Forecast", linestyle="dashed")
plt.legend()
plt.title("SARIMA with Seasonality")
plt.show()

# SARIMAX: Predict with an external variable (exogenous)just for demonstration, use a dummy one
exog = np.arange(len(train)).reshape(-1,1)  # Fake variable for demonstration
sarimax_model = SARIMAX(train, order=(1,1,1), seasonal_order=(1,1,1,12), exog=exog)
sarimax_fit = sarimax_model.fit(disp=False)
exog_forecast = np.arange(len(train), len(train)+12).reshape(-1,1)
sarimax_forecast = sarimax_fit.forecast(steps=12, exog=exog_forecast)

plt.figure(figsize=(10,4))
plt.plot(train.index, train, label="Train")
plt.plot(test.index, test, label="Test")
plt.plot(test.index, sarimax_forecast, label="SARIMAX Forecast", linestyle="dotted", color="magenta")
plt.legend()
plt.title("SARIMAX with Exogenous Variable")
plt.show()

# Evaluate error: Let us compare our forecast to real values using Mean Absolute Error
from sklearn.metrics import mean_absolute_error

arima_err = mean_absolute_error(test, forecast)
sarima_err = mean_absolute_error(test, sarima_forecast)
sarimax_err = mean_absolute_error(test, sarimax_forecast)

print("ARIMA MAE:", round(arima_err, 2))
print("SARIMA MAE:", round(sarima_err, 2))
print("SARIMAX MAE:", round(sarimax_err, 2))

# Load Shampoo Sales data
url2 = "https://raw.githubusercontent.com/jbrownlee/Datasets/master/shampoo.csv"
shampoo = pd.read_csv(url2, parse_dates=["Month"], index_col="Month")
shampoo.columns = ["Sales"]

print("Shape:", shampoo.shape)
shampoo.head()

# Plot Shampoo Sales time series
shampoo.plot(y="Sales", legend=False, figsize=(10,4), title="Monthly Shampoo Sales")
plt.ylabel("Sales")
plt.xlabel("Month")
plt.show()

# Your turn: Build and compare ARIMA and SARIMA on Shampoo Sales
train2 = shampoo["Sales"][:-6]
test2 = shampoo["Sales"][-6:]

arima2 = ARIMA(train2, order=(1,1,1)).fit()
sarima2 = SARIMAX(train2, order=(1,1,1), seasonal_order=(1,1,1,12)).fit(disp=False)

arima2_forecast = arima2.forecast(steps=6)
sarima2_forecast = sarima2.forecast(steps=6)

plt.plot(train2.index, train2, label="Train")
plt.plot(test2.index, test2, label="Test")
plt.plot(test2.index, arima2_forecast, label="ARIMA Forecast")
plt.plot(test2.index, sarima2_forecast, label="SARIMA Forecast")
plt.legend()
plt.title("Shampoo Sales Forecast Comparison")
plt.show()

# Extra tip: Avoid overfittingdo not make the model too complex for small data
order = (2,1,2)
seasonal_order = (1,1,1,12)

model_test = SARIMAX(train2, order=order, seasonal_order=seasonal_order)
fit_test = model_test.fit(disp=False)
forecast_test = fit_test.forecast(steps=6)

print("Test forecast values:")
print(forecast_test)

# Troubleshooting: What if you get an error like 'non-invertible starting MA parameters'? 
try:
    broken_model = SARIMAX(train2, order=(0,1,0), seasonal_order=(1,1,1,12)).fit(disp=False)
except Exception as e:
    print("Error message:", str(e))

# This happens if the model is too simple or does not fit the data well