import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.seasonal import seasonal_decompose

# Generate sample time series data
np.random.seed(42)
dates = pd.date_range('2020-01-01', periods=730, freq='D')

# Create components
trend = np.linspace(100, 200, 730)  # Upward trend
seasonal = 30 * np.sin(2 * np.pi * np.arange(730) / 365)  # Yearly seasonality
noise = np.random.normal(0, 10, 730)  # Random noise

sales = trend + seasonal + noise
ts = pd.Series(sales, index=dates)

# Decompose
decomposition = seasonal_decompose(ts, model='additive', period=365)

# Plot
fig, axes = plt.subplots(4, 1, figsize=(12, 10))

axes[0].plot(ts)
axes[0].set_title('Original Time Series')

axes[1].plot(decomposition.trend)
axes[1].set_title('Trend Component')

axes[2].plot(decomposition.seasonal)
axes[2].set_title('Seasonal Component')

axes[3].plot(decomposition.resid)
axes[3].set_title('Residual (Noise)')

plt.tight_layout()
plt.show()

from statsmodels.tsa.stattools import adfuller

def test_stationarity(series, name="Series"):
    """Test if series is stationary using Augmented Dickey-Fuller test."""
    result = adfuller(series.dropna())
    
    print(f"\n{name} - Stationarity Test:")
    print(f"  ADF Statistic: {result[0]:.4f}")
    print(f"  p-value: {result[1]:.4f}")
    print(f"  Stationary: {'Yes' if result[1] < 0.05 else 'No'}")
    
    return result[1] < 0.05

# Test our series
test_stationarity(ts, "Raw Sales")

# Make it stationary through differencing
ts_diff = ts.diff().dropna()
test_stationarity(ts_diff, "Differenced Sales")

# Differencing: subtract previous value to remove trend
# If today's value is 105 and yesterday's was 100, the differenced value is 5
ts_diff = ts.diff()

# Log transform: stabilize variance when fluctuations grow with the level
# (e.g., stock prices: a $10 swing means more at $50 than at $500)
ts_log = np.log(ts)

# Both combined: handles both growing variance AND trend
ts_log_diff = np.log(ts).diff()

# Plot transformations
fig, axes = plt.subplots(2, 2, figsize=(12, 8))

axes[0, 0].plot(ts)
axes[0, 0].set_title('Original')

axes[0, 1].plot(ts_diff)
axes[0, 1].set_title('Differenced')

axes[1, 0].plot(ts_log)
axes[1, 0].set_title('Log Transformed')

axes[1, 1].plot(ts_log_diff)
axes[1, 1].set_title('Log + Differenced')

plt.tight_layout()
plt.show()

def moving_average_forecast(series, window=7):
    """Forecast using simple moving average."""
    return series.rolling(window=window).mean()

# Calculate 7-day and 30-day moving averages
ma_7 = moving_average_forecast(ts, 7)
ma_30 = moving_average_forecast(ts, 30)

plt.figure(figsize=(12, 6))
plt.plot(ts.values[-100:], label='Actual', alpha=0.7)
plt.plot(ma_7.values[-100:], label='7-day MA', linewidth=2)
plt.plot(ma_30.values[-100:], label='30-day MA', linewidth=2)
plt.legend()
plt.title('Moving Average Smoothing')
plt.show()

from statsmodels.tsa.holtwinters import SimpleExpSmoothing, ExponentialSmoothing

# Simple Exponential Smoothing
ses = SimpleExpSmoothing(ts).fit(smoothing_level=0.3)
ses_forecast = ses.forecast(30)

# Holt-Winters (handles trend + seasonality)
hw = ExponentialSmoothing(
    ts, 
    trend='add', 
    seasonal='add', 
    seasonal_periods=365
).fit()
hw_forecast = hw.forecast(30)

# Plot
plt.figure(figsize=(12, 6))
plt.plot(ts.values[-100:], label='Historical')
plt.plot(range(100, 130), ses_forecast, label='Simple ES', linestyle='--')
plt.plot(range(100, 130), hw_forecast, label='Holt-Winters', linestyle='--')
plt.axvline(99, color='gray', linestyle=':')
plt.legend()
plt.title('Exponential Smoothing Forecasts')
plt.show()

from statsmodels.tsa.arima.model import ARIMA
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

# Plot ACF and PACF to determine p and q
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
plot_acf(ts_diff.dropna(), ax=axes[0], lags=40)
plot_pacf(ts_diff.dropna(), ax=axes[1], lags=40)
plt.tight_layout()
plt.show()

# Fit ARIMA model
# Train on first 700 days, test on last 30
train = ts[:700]
test = ts[700:]

model = ARIMA(train, order=(1, 1, 1))
fitted = model.fit()

print(fitted.summary())

# Forecast
forecast = fitted.forecast(len(test))

# Plot
plt.figure(figsize=(12, 6))
plt.plot(train.values[-100:], label='Training')
plt.plot(range(100, 100+len(test)), test.values, label='Actual')
plt.plot(range(100, 100+len(forecast)), forecast.values, label='Forecast', linestyle='--')
plt.axvline(99, color='gray', linestyle=':')
plt.legend()
plt.title('ARIMA Forecast')
plt.show()

# Using pmdarima for automatic order selection
# pip install pmdarima
from pmdarima import auto_arima

auto_model = auto_arima(
    train,
    start_p=0, start_q=0,
    max_p=5, max_q=5,
    d=None,  # Let it determine d
    seasonal=True,
    m=7,  # Weekly seasonality
    trace=True,
    error_action='ignore',
    suppress_warnings=True
)

print(auto_model.summary())

from statsmodels.tsa.statespace.sarimax import SARIMAX

# SARIMA with weekly seasonality
model = SARIMAX(
    train,
    order=(1, 1, 1),           # (p, d, q)
    seasonal_order=(1, 1, 1, 7) # (P, D, Q, m) m=7 for weekly
)

fitted = model.fit(disp=False)
forecast = fitted.forecast(len(test))

# Calculate RMSE
from sklearn.metrics import mean_squared_error
rmse = np.sqrt(mean_squared_error(test, forecast))
print(f"SARIMA RMSE: {rmse:.2f}")

# pip install prophet
from prophet import Prophet

# Prepare data in Prophet format
df = pd.DataFrame({
    'ds': dates,
    'y': sales
})

# Train/test split
train_df = df[:700]
test_df = df[700:]

# Fit Prophet
prophet = Prophet(
    yearly_seasonality=True,
    weekly_seasonality=True,
    daily_seasonality=False
)
prophet.fit(train_df)

# Make future dataframe
future = prophet.make_future_dataframe(periods=30)
forecast = prophet.predict(future)

# Plot
prophet.plot(forecast)
plt.title('Prophet Forecast')
plt.show()

# Plot components
prophet.plot_components(forecast)
plt.show()

# Define holidays
holidays = pd.DataFrame({
    'holiday': 'christmas',
    'ds': pd.to_datetime(['2020-12-25', '2021-12-25', '2022-12-25']),
    'lower_window': -2,  # 2 days before
    'upper_window': 1    # 1 day after
})

prophet = Prophet(holidays=holidays)

def create_features(df, target_col, n_lags=7):
    """Create features from time series for ML models."""
    df = df.copy()
    
    # Lag features
    for i in range(1, n_lags + 1):
        df[f'lag_{i}'] = df[target_col].shift(i)
    
    # Rolling statistics
    df['rolling_mean_7'] = df[target_col].rolling(7).mean()
    df['rolling_std_7'] = df[target_col].rolling(7).std()
    df['rolling_mean_30'] = df[target_col].rolling(30).mean()
    
    # Date features
    if isinstance(df.index, pd.DatetimeIndex):
        df['day_of_week'] = df.index.dayofweek
        df['month'] = df.index.month
        df['day_of_month'] = df.index.day
        df['is_weekend'] = (df.index.dayofweek >= 5).astype(int)
    
    return df.dropna()

# Prepare dataset
ts_df = pd.DataFrame({'sales': sales}, index=dates)
ts_ml = create_features(ts_df, 'sales', n_lags=7)

# Features and target
feature_cols = [col for col in ts_ml.columns if col != 'sales']
X = ts_ml[feature_cols]
y = ts_ml['sales']

# Time series split (don't shuffle!)
split_idx = int(len(X) * 0.8)
X_train, X_test = X[:split_idx], X[split_idx:]
y_train, y_test = y[:split_idx], y[split_idx:]

# Train models
from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor
from sklearn.metrics import mean_absolute_error, mean_squared_error

models = {
    'Random Forest': RandomForestRegressor(n_estimators=100, random_state=42),
    'Gradient Boosting': GradientBoostingRegressor(n_estimators=100, random_state=42)
}

for name, model in models.items():
    model.fit(X_train, y_train)
    predictions = model.predict(X_test)
    
    mae = mean_absolute_error(y_test, predictions)
    rmse = np.sqrt(mean_squared_error(y_test, predictions))
    
    print(f"{name}: MAE={mae:.2f}, RMSE={rmse:.2f}")

from sklearn.model_selection import TimeSeriesSplit

# Time series cross-validation
tscv = TimeSeriesSplit(n_splits=5)

# Visualize the splits
plt.figure(figsize=(12, 4))
for i, (train_idx, test_idx) in enumerate(tscv.split(X)):
    plt.plot(train_idx, [i] * len(train_idx), 'b-', linewidth=3, label='Train' if i==0 else '')
    plt.plot(test_idx, [i] * len(test_idx), 'r-', linewidth=3, label='Test' if i==0 else '')
plt.xlabel('Sample index')
plt.ylabel('CV iteration')
plt.legend()
plt.title('Time Series Cross-Validation')
plt.show()

# Use in model evaluation
from sklearn.model_selection import cross_val_score

rf = RandomForestRegressor(n_estimators=100, random_state=42)
scores = cross_val_score(rf, X, y, cv=tscv, scoring='neg_mean_squared_error')
rmse_scores = np.sqrt(-scores)
print(f"CV RMSE: {rmse_scores.mean():.2f} ± {rmse_scores.std():.2f}")

def recursive_forecast(model, X_last, n_steps, n_lags):
    """Recursively forecast multiple steps."""
    forecasts = []
    current_features = X_last.copy()
    
    for _ in range(n_steps):
        pred = model.predict(current_features.reshape(1, -1))[0]
        forecasts.append(pred)
        
        # Shift lag features
        for i in range(n_lags - 1, 0, -1):
            current_features[i] = current_features[i-1]
        current_features[0] = pred  # lag_1 = prediction
    
    return forecasts

def direct_forecast(df, target_col, horizons=[1, 7, 14]):
    """Train separate models for each forecast horizon."""
    models = {}
    
    for h in horizons:
        # Create target h steps ahead
        df_h = df.copy()
        df_h['target'] = df_h[target_col].shift(-h)
        df_h = df_h.dropna()
        
        feature_cols = [col for col in df_h.columns if col not in [target_col, 'target']]
        X = df_h[feature_cols]
        y = df_h['target']
        
        # Train model
        model = RandomForestRegressor(n_estimators=100, random_state=42)
        model.fit(X[:-h], y[:-h])  # Leave last h for testing
        
        models[h] = model
    
    return models

def evaluate_forecast(actual, predicted, name="Model"):
    """Calculate common time series metrics."""
    mae = mean_absolute_error(actual, predicted)
    mse = mean_squared_error(actual, predicted)
    rmse = np.sqrt(mse)
    mape = np.mean(np.abs((actual - predicted) / actual)) * 100
    
    print(f"\n{name} Metrics:")
    print(f"  MAE:  {mae:.2f}")
    print(f"  RMSE: {rmse:.2f}")
    print(f"  MAPE: {mape:.1f}%")
    
    return {'MAE': mae, 'RMSE': rmse, 'MAPE': mape}

# Compare models
evaluate_forecast(test, forecast, "ARIMA")

# WRONG: Using future data for scaling
from sklearn.preprocessing import StandardScaler

# Don't do this! The scaler learns mean/std from ALL data, including the future
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)  # Uses all data including test!

# CORRECT: Fit only on training (past) data
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)  # Learn stats from past only
X_test_scaled = scaler.transform(X_test)  # Apply past stats to future -- no leakage

# Test on multiple time periods, not just the most recent
# Use expanding or sliding window backtests:
# Window 1: Train on Jan-Jun, test on Jul
# Window 2: Train on Jan-Sep, test on Oct
# Window 3: Train on Jan-Dec, test on next Jan
# This reveals whether your model is robust across different market conditions