Time Series Analysis
& Forecasting

By the end of this chapter, you will be able to:

1. Decompose a time series into its core components: trend, seasonality, cyclical, and noise.
2. Test for stationarity using the ADF (Augmented Dickey-Fuller) and KPSS tests, and apply differencing/transformations to achieve stationarity.
3. Interpret ACF and PACF plots to identify appropriate model orders for AR, MA, and ARIMA models.
4. Build and tune ARIMA(p,d,q) models including seasonal variants (SARIMA) using the Box-Jenkins methodology.
5. Apply Exponential Smoothing methods: Simple, Holt's Linear Trend, and Holt-Winters Seasonal models.
6. Use Facebook Prophet for time series forecasting with trend changepoints, holiday effects, and custom seasonality.
7. Design LSTM networks using TensorFlow/Keras for sequence-to-one and sequence-to-sequence forecasting.
8. Understand multivariate time series analysis using Vector Autoregression (VAR).
9. Detect anomalies in time series data using statistical and ML-based approaches.
10. Evaluate forecasts with MAPE, RMSE, SMAPE, and other metrics, knowing their tradeoffs.
11. Apply time series methods to real-world Indian and global case studies, including Nifty50 forecasting, monsoon prediction, and energy demand.
12. Build end-to-end mini projects — Stock Price Forecaster and Weather Predictor.

Time series data is everywhere — from the minute-by-minute stock prices on the Bombay Stock Exchange to daily temperatures recorded by the India Meteorological Department, monthly GST revenues, and the second-by-second sensor readings in a smart factory. Unlike cross-sectional data where observations are independent, time series data has temporal ordering, and that ordering contains crucial information.

Time series analysis is the discipline of extracting meaningful statistics and characteristics from time-ordered data, while time series forecasting uses historical observations to predict future values. The applications span virtually every industry:

• Finance: Stock prices, exchange rates, portfolio risk (VaR)
• Weather: Temperature, rainfall, cyclone paths
• Retail: Sales demand, inventory optimization
• Healthcare: Patient vitals monitoring, epidemic forecasting
• Energy: Electricity demand, solar/wind output prediction
• Telecom: Network traffic, call volume
• Manufacturing: Sensor data for predictive maintenance

This chapter walks through the entire landscape — from classical decomposition and ARIMA to modern deep learning (LSTM) and Meta's Prophet. We adopt a practitioner's approach: understand the theory, implement from scratch, then use production-ready libraries.

The study of time series has deep roots in astronomy, economics, and meteorology.

Timeline of Key Developments

The Box-Jenkins methodology (1970) was the watershed moment. It provided a systematic approach — identify, estimate, diagnose — that remains the gold standard for classical time series modeling. George Box famously said, "All models are wrong, but some are useful," which perfectly captures the forecasting philosophy.

4.1 Time Series Components

Every time series Y(t) can be decomposed into four fundamental components:

• Trend (T): The long-term upward or downward movement. Example: India's GDP growth over decades.
• Seasonality (S): Regular, periodic patterns repeating at fixed intervals. Example: Ice cream sales peak every summer, Diwali shopping spikes every October/November.
• Cyclical (C): Longer-term fluctuations without fixed period, often linked to business cycles. Example: Economic booms and recessions lasting 5-10 years.
• Noise/Residual (ε): Random, unpredictable variation after removing other components.

Two decomposition models are common:

Use additive when seasonal fluctuations are roughly constant regardless of the level. Use multiplicative when seasonal variation grows proportionally with the level (e.g., airline passengers — the seasonal swing is larger when overall traffic is higher).

4.2 Stationarity

A time series is stationary when its statistical properties (mean, variance, autocorrelation) do not change over time. Most forecasting models require stationarity.

Strict vs. Weak Stationarity

• Strict stationarity: The joint distribution of (Y(t₁), Y(t₂), ..., Y(tₖ)) is identical to (Y(t₁+τ), Y(t₂+τ), ..., Y(tₖ+τ)) for all time shifts τ.
• Weak (covariance) stationarity: Mean is constant, variance is finite and constant, covariance depends only on the lag (not time). This is the practical requirement.

Testing for Stationarity

Augmented Dickey-Fuller (ADF) Test:

• H₀: Unit root exists (series is non-stationary)
• H₁: No unit root (series is stationary)
• If p-value < 0.05 → reject H₀ → series IS stationary

KPSS Test:

• H₀: Series is stationary
• H₁: Series has a unit root (non-stationary)
• If p-value < 0.05 → reject H₀ → series is NOT stationary

4.3 Differencing and Transformations

To make a non-stationary series stationary:

• First Differencing: Y'(t) = Y(t) - Y(t-1) removes linear trends.
• Second Differencing: Y''(t) = Y'(t) - Y'(t-1) removes quadratic trends.
• Seasonal Differencing: Y'(t) = Y(t) - Y(t-m) where m is the seasonal period.
• Log Transform: Stabilizes variance when it grows with the level.
• Box-Cox Transform: Generalized power transform Y(λ) = (Y^λ - 1)/λ.

4.4 ACF and PACF

Autocorrelation Function (ACF) measures the correlation between Y(t) and Y(t-k) for lag k. It includes indirect correlations through intermediate lags.

Partial Autocorrelation Function (PACF) measures the direct correlation between Y(t) and Y(t-k), removing the effect of all intermediate lags.

4.5 ARIMA Models

AR(p) — Autoregressive: Current value depends on p past values.

MA(q) — Moving Average: Current value depends on q past error terms.

ARIMA(p,d,q): Combines AR, differencing (d times), and MA. The I stands for "Integrated" — it reverses differencing.

4.6 SARIMA

SARIMA extends ARIMA for seasonal data: SARIMA(p,d,q)(P,D,Q,m) where uppercase letters are seasonal orders and m is the seasonal period (e.g., 12 for monthly data with annual seasonality).

4.7 Exponential Smoothing

• SES (Simple Exponential Smoothing): For level-only data — no trend, no seasonality.
• Holt's Linear: Adds a trend component.
• Holt-Winters: Adds both trend and seasonality. Comes in additive and multiplicative variants.

4.8 Prophet

Developed by Meta (Facebook), Prophet uses a decomposable additive model: y(t) = g(t) + s(t) + h(t) + ε(t) where g(t) is a piecewise-linear or logistic growth trend, s(t) is seasonality modeled with Fourier terms, and h(t) captures holiday effects. Prophet is designed to handle missing data, changepoints, and multiple seasonalities automatically.

4.9 LSTM for Time Series

Long Short-Term Memory networks are a type of recurrent neural network (RNN) that can learn long-range dependencies. They use a gating mechanism (forget, input, output gates) to selectively remember or forget information. For time series, the input is a sliding window of past observations, and the output is the forecast.

4.10 VAR — Vector Autoregression

When you have multiple interrelated time series (e.g., Nifty50 + USD/INR + gold prices), VAR models capture the linear interdependencies. Each variable is a linear function of past lags of itself AND all other variables.

4.11 Anomaly Detection in Time Series

Anomalies are observations that deviate significantly from expected patterns. Methods include: statistical thresholds (z-score, IQR), forecast-based (flag points with large residuals), isolation forests, autoencoders, and Prophet's built-in uncertainty intervals.

5.1 AR(p) Process

The process is stationary if all roots of the characteristic polynomial 1 - φ₁z - φ₂z² - ... - φₚzᵖ = 0 lie outside the unit circle.

5.2 MA(q) Process

An MA(q) process is always stationary (finite linear combination of white noise). It is invertible if all roots of 1 + θ₁z + θ₂z² + ... + θₑzᵍ = 0 lie outside the unit circle.

5.3 ARIMA(p,d,q)

5.4 SARIMA(p,d,q)(P,D,Q)ₘ

5.5 Exponential Smoothing — State Space Form

5.6 Evaluation Metrics

5.7 ADF Test Statistic

6.1 Deriving the AR(1) Autocorrelation Function

For AR(1): Y(t) = φY(t-1) + ε(t), with |φ| < 1 for stationarity.

Step 1: Variance — γ(0) = Var(Y(t)) = φ²Var(Y(t-1)) + σ² = φ²γ(0) + σ²

Therefore: γ(0) = σ² / (1 - φ²)

Step 2: Autocovariance at lag k — Multiply Y(t) = φY(t-1) + ε(t) by Y(t-k) and take expectations:

γ(k) = E[Y(t)·Y(t-k)] = φ·E[Y(t-1)·Y(t-k)] + E[ε(t)·Y(t-k)] = φ·γ(k-1) + 0 = φ·γ(k-1)

By recursion: γ(k) = φᵏ · γ(0)

Step 3: Autocorrelation: ρ(k) = γ(k)/γ(0) = φᵏ

6.2 Deriving MA(1) Autocorrelation

For MA(1): Y(t) = ε(t) + θε(t-1)

Step 1: γ(0) = Var(ε(t) + θε(t-1)) = (1 + θ²)σ²

Step 2: γ(1) = Cov(ε(t) + θε(t-1), ε(t-1) + θε(t-2)) = θσ²

Step 3: γ(k) = 0 for k ≥ 2 (non-overlapping noise terms)

6.3 Deriving Simple Exponential Smoothing

The SES forecast ŷ(t+1) = αy(t) + (1-α)ŷ(t) can be expanded recursively:

ŷ(t+1) = αy(t) + (1-α)[αy(t-1) + (1-α)ŷ(t-1)]

= αy(t) + α(1-α)y(t-1) + (1-α)²ŷ(t-1)

Continuing infinitely:

6.4 Box-Cox Transformation

Example 1: AR(1) Forecast by Hand

Given: AR(1) model Y(t) = 5 + 0.7·Y(t-1) + ε(t), with Y(100) = 20.

Task: Forecast Y(101), Y(102), Y(103).

Solution:

• ŷ(101) = 5 + 0.7 × 20 = 5 + 14 = 19.0
• ŷ(102) = 5 + 0.7 × 19.0 = 5 + 13.3 = 18.3
• ŷ(103) = 5 + 0.7 × 18.3 = 5 + 12.81 = 17.81

Note: Forecasts converge toward the long-run mean μ = c/(1-φ) = 5/(1-0.7) = 16.67.

Example 2: Simple Exponential Smoothing

Given: Sales data: [100, 120, 110, 130, 125]. α = 0.3. Initial forecast ŷ(1) = 100.

Solution:

• ŷ(2) = 0.3×100 + 0.7×100 = 30 + 70 = 100.0
• ŷ(3) = 0.3×120 + 0.7×100 = 36 + 70 = 106.0
• ŷ(4) = 0.3×110 + 0.7×106 = 33 + 74.2 = 107.2
• ŷ(5) = 0.3×130 + 0.7×107.2 = 39 + 75.04 = 114.04
• ŷ(6) = 0.3×125 + 0.7×114.04 = 37.5 + 79.83 = 117.33

Example 3: MAPE Calculation

Given: Actual = [100, 150, 200], Forecast = [110, 140, 180].

Solution:

• |100-110|/100 = 10/100 = 0.10
• |150-140|/150 = 10/150 = 0.0667
• |200-180|/200 = 20/200 = 0.10

MAPE = (100/3)(0.10 + 0.0667 + 0.10) = (100/3)(0.2667) = 8.89%

Example 4: Differencing by Hand

Given: Y = [10, 13, 18, 25, 34, 45]

First Difference: ΔY = [3, 5, 7, 9, 11] — still shows a trend

Second Difference: Δ²Y = [2, 2, 2, 2] — constant! Series is now stationary.

Since d=2 differences were needed, the original series has a quadratic trend (Y ≈ t²).

Example 5: ACF for MA(1)

Given: MA(1) model Y(t) = ε(t) + 0.6ε(t-1), σ² = 4.

Solution:

• γ(0) = (1 + 0.6²) × 4 = 1.36 × 4 = 5.44
• γ(1) = 0.6 × 4 = 2.4
• γ(k) = 0 for k ≥ 2
• ρ(1) = 2.4/5.44 = 0.441
• ρ(k) = 0 for k ≥ 2 — ACF cuts off after lag 1 ✓

10.1 Simple Exponential Smoothing from Scratch

PYTHON
import numpy as np

class SimpleExponentialSmoothing:
    """Simple Exponential Smoothing from scratch."""
    
    def __init__(self, alpha=0.3):
        self.alpha = alpha
        self.fitted_values = []
    
    def fit(self, series):
        """Fit the SES model to a time series."""
        self.series = np.array(series, dtype=float)
        n = len(self.series)
        self.fitted_values = np.zeros(n)
        
        # Initialize: first forecast = first observation
        self.fitted_values[0] = self.series[0]
        
        for t in range(1, n):
            self.fitted_values[t] = (
                self.alpha * self.series[t - 1] +
                (1 - self.alpha) * self.fitted_values[t - 1]
            )
        
        self.level = self.fitted_values[-1]
        return self
    
    def forecast(self, steps=1):
        """Forecast future values (SES gives flat forecast)."""
        last_smoothed = (
            self.alpha * self.series[-1] +
            (1 - self.alpha) * self.level
        )
        return np.full(steps, last_smoothed)
    
    def mape(self):
        """Calculate Mean Absolute Percentage Error."""
        actual = self.series[1:]
        forecast = self.fitted_values[1:]
        return np.mean(np.abs((actual - forecast) / actual)) * 100


# Example usage
data = [112, 118, 132, 129, 121, 135, 148, 148, 136, 119, 104, 118]
model = SimpleExponentialSmoothing(alpha=0.3)
model.fit(data)
print("Fitted:", np.round(model.fitted_values, 2))
print("Forecast (next 3):", np.round(model.forecast(3), 2))
print(f"MAPE: {model.mape():.2f}%")

10.2 AR(1) Model from Scratch

PYTHON
import numpy as np

class AR1Model:
    """AR(1) model: Y(t) = c + phi * Y(t-1) + epsilon."""
    
    def __init__(self):
        self.phi = None
        self.c = None
        self.sigma2 = None
    
    def fit(self, series):
        """Estimate AR(1) parameters using OLS."""
        y = np.array(series, dtype=float)
        n = len(y)
        
        # Set up regression: y(t) = c + phi * y(t-1)
        Y = y[1:]      # dependent variable
        X = y[:-1]      # lagged variable
        
        # Add intercept
        X_design = np.column_stack([np.ones(n - 1), X])
        
        # OLS: beta = (X'X)^{-1} X'Y
        beta = np.linalg.solve(X_design.T @ X_design, X_design.T @ Y)
        self.c = beta[0]
        self.phi = beta[1]
        
        # Residuals
        residuals = Y - X_design @ beta
        self.sigma2 = np.var(residuals, ddof=2)
        self.last_value = y[-1]
        
        return self
    
    def forecast(self, steps=5):
        """Generate multi-step ahead forecasts."""
        predictions = []
        current = self.last_value
        for _ in range(steps):
            next_val = self.c + self.phi * current
            predictions.append(next_val)
            current = next_val
        return np.array(predictions)
    
    def long_run_mean(self):
        """E[Y] = c / (1 - phi) for stationary process."""
        if abs(self.phi) >= 1:
            return float('inf')
        return self.c / (1 - self.phi)
    
    def is_stationary(self):
        """Check stationarity condition: |phi| < 1."""
        return abs(self.phi) < 1


# Example: Fit AR(1) to synthetic data
np.random.seed(42)
n = 200
y = np.zeros(n)
y[0] = 10
for t in range(1, n):
    y[t] = 3 + 0.7 * y[t-1] + np.random.normal(0, 1)

model = AR1Model()
model.fit(y)
print(f"Estimated: c={model.c:.4f}, phi={model.phi:.4f}")
print(f"True:      c=3.0000, phi=0.7000")
print(f"Stationary: {model.is_stationary()}")
print(f"Long-run mean: {model.long_run_mean():.4f} (true: {3/(1-0.7):.4f})")
print(f"Forecast: {np.round(model.forecast(5), 2)}")

10.3 Augmented Dickey-Fuller Test from Scratch

PYTHON
import numpy as np

def adf_test_manual(series, max_lags=1):
    """
    Manual ADF test implementation.
    Tests H0: unit root (non-stationary) vs H1: stationary.
    """
    y = np.array(series, dtype=float)
    n = len(y)
    
    # Compute first differences
    dy = np.diff(y)
    
    # Lagged level
    y_lag = y[:-1]
    
    # Trim to account for lags
    if max_lags > 0:
        # Include lagged differences
        dy_lags = np.column_stack([
            dy[max_lags - i - 1: n - 1 - i - 1]
            for i in range(max_lags)
        ])
        dy_trimmed = dy[max_lags:]
        y_lag_trimmed = y_lag[max_lags:]
        
        X = np.column_stack([
            np.ones(len(dy_trimmed)),   # intercept
            y_lag_trimmed,               # gamma * y(t-1)
            dy_lags                      # lagged differences
        ])
    else:
        dy_trimmed = dy
        y_lag_trimmed = y_lag
        X = np.column_stack([np.ones(len(dy_trimmed)), y_lag_trimmed])
    
    # OLS estimation
    beta = np.linalg.solve(X.T @ X, X.T @ dy_trimmed)
    residuals = dy_trimmed - X @ beta
    
    # Standard error of gamma
    sigma2 = np.sum(residuals**2) / (len(dy_trimmed) - len(beta))
    cov_matrix = sigma2 * np.linalg.inv(X.T @ X)
    se_gamma = np.sqrt(cov_matrix[1, 1])
    
    # Test statistic
    gamma = beta[1]
    t_stat = gamma / se_gamma
    
    # Critical values (approximate, for n > 100 with intercept)
    critical_values = {
        '1%':  -3.43,
        '5%':  -2.86,
        '10%': -2.57
    }
    
    print(f"ADF Test Statistic: {t_stat:.4f}")
    print(f"Estimated gamma:    {gamma:.6f}")
    for level, cv in critical_values.items():
        reject = "REJECT H0 (Stationary)" if t_stat < cv else "Fail to reject H0"
        print(f"  {level} critical value: {cv:.2f} → {reject}")
    
    return t_stat, gamma


# Test with a stationary series
np.random.seed(42)
stationary = np.random.normal(0, 1, 200).cumsum()  # Random walk (non-stationary)
print("=== Random Walk (Non-Stationary) ===")
adf_test_manual(stationary)

print("\n=== White Noise (Stationary) ===")
white_noise = np.random.normal(0, 1, 200)
adf_test_manual(white_noise)

10.4 ACF/PACF Computation

PYTHON
import numpy as np

def compute_acf(series, max_lag=20):
    """Compute the Autocorrelation Function."""
    y = np.array(series, dtype=float)
    n = len(y)
    mean = np.mean(y)
    var = np.sum((y - mean) ** 2) / n
    
    acf_values = []
    for k in range(max_lag + 1):
        if k == 0:
            acf_values.append(1.0)
        else:
            cov = np.sum((y[k:] - mean) * (y[:-k] - mean)) / n
            acf_values.append(cov / var)
    
    return np.array(acf_values)


def compute_pacf(series, max_lag=20):
    """Compute PACF using the Durbin-Levinson algorithm."""
    acf = compute_acf(series, max_lag)
    n_lags = max_lag
    pacf_values = np.zeros(n_lags + 1)
    pacf_values[0] = 1.0
    pacf_values[1] = acf[1]
    
    phi = np.zeros((n_lags + 1, n_lags + 1))
    phi[1, 1] = acf[1]
    
    for k in range(2, n_lags + 1):
        # Numerator
        num = acf[k] - sum(phi[k-1, j] * acf[k-j] for j in range(1, k))
        # Denominator
        den = 1.0 - sum(phi[k-1, j] * acf[j] for j in range(1, k))
        
        phi[k, k] = num / den
        
        for j in range(1, k):
            phi[k, j] = phi[k-1, j] - phi[k, k] * phi[k-1, k-j]
        
        pacf_values[k] = phi[k, k]
    
    return pacf_values


# Example
np.random.seed(42)
ar1_data = np.zeros(500)
for t in range(1, 500):
    ar1_data[t] = 0.8 * ar1_data[t-1] + np.random.normal(0, 1)

acf = compute_acf(ar1_data, 10)
pacf = compute_pacf(ar1_data, 10)

print("Lag  |  ACF    |  PACF")
print("-" * 30)
for k in range(11):
    print(f"  {k:2d}  | {acf[k]:6.3f}  | {pacf[k]:6.3f}")
print("\n→ ACF decays exponentially, PACF cuts off after lag 1 → AR(1) ✓")

11.1 LSTM Time Series Forecaster

TENSORFLOW / KERAS
import numpy as np
import tensorflow as tf
from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import LSTM, Dense, Dropout
from tensorflow.keras.callbacks import EarlyStopping
from sklearn.preprocessing import MinMaxScaler

# =========================================================
# 1. Generate / Load Data
# =========================================================
# Synthetic: trend + seasonality + noise
np.random.seed(42)
t = np.arange(0, 500)
series = 0.05 * t + 10 * np.sin(2 * np.pi * t / 50) + np.random.normal(0, 2, 500)
series = series.reshape(-1, 1)

# =========================================================
# 2. Scale Data
# =========================================================
scaler = MinMaxScaler(feature_range=(0, 1))
scaled = scaler.fit_transform(series)

# =========================================================
# 3. Create Sliding Window Dataset
# =========================================================
def create_sequences(data, window_size=30):
    X, y = [], []
    for i in range(window_size, len(data)):
        X.append(data[i - window_size:i, 0])
        y.append(data[i, 0])
    return np.array(X), np.array(y)

WINDOW = 30
X, y = create_sequences(scaled, WINDOW)

# Train-test split (80/20)
split = int(0.8 * len(X))
X_train, X_test = X[:split], X[split:]
y_train, y_test = y[:split], y[split:]

# Reshape for LSTM: (samples, timesteps, features)
X_train = X_train.reshape(-1, WINDOW, 1)
X_test = X_test.reshape(-1, WINDOW, 1)

# =========================================================
# 4. Build LSTM Model
# =========================================================
model = Sequential([
    LSTM(64, return_sequences=True, input_shape=(WINDOW, 1)),
    Dropout(0.2),
    LSTM(32, return_sequences=False),
    Dropout(0.2),
    Dense(16, activation='relu'),
    Dense(1)
])

model.compile(
    optimizer=tf.keras.optimizers.Adam(learning_rate=0.001),
    loss='mse',
    metrics=['mae']
)

model.summary()

# =========================================================
# 5. Train
# =========================================================
early_stop = EarlyStopping(
    monitor='val_loss', patience=10, restore_best_weights=True
)

history = model.fit(
    X_train, y_train,
    epochs=100,
    batch_size=32,
    validation_split=0.2,
    callbacks=[early_stop],
    verbose=1
)

# =========================================================
# 6. Evaluate
# =========================================================
y_pred_scaled = model.predict(X_test)
y_pred = scaler.inverse_transform(y_pred_scaled)
y_actual = scaler.inverse_transform(y_test.reshape(-1, 1))

rmse = np.sqrt(np.mean((y_actual - y_pred) ** 2))
mape = np.mean(np.abs((y_actual - y_pred) / y_actual)) * 100
print(f"\nRMSE: {rmse:.4f}")
print(f"MAPE: {mape:.2f}%")

# =========================================================
# 7. Multi-Step Forecast
# =========================================================
def forecast_multistep(model, last_window, steps, scaler):
    """Iteratively forecast multiple steps ahead."""
    predictions = []
    current = last_window.copy()
    
    for _ in range(steps):
        pred = model.predict(current.reshape(1, -1, 1), verbose=0)
        predictions.append(pred[0, 0])
        current = np.append(current[1:], pred[0, 0])
    
    predictions = np.array(predictions).reshape(-1, 1)
    return scaler.inverse_transform(predictions)

future = forecast_multistep(model, scaled[-WINDOW:, 0], 10, scaler)
print(f"\nNext 10 forecasts: {future.flatten().round(2)}")

11.2 Bidirectional LSTM with Attention

TENSORFLOW / KERAS
from tensorflow.keras.layers import (
    Bidirectional, LSTM, Dense, Dropout,
    Attention, Input, Concatenate, Flatten
)
from tensorflow.keras.models import Model

def build_attention_lstm(window_size, n_features=1):
    """LSTM with simple attention mechanism for time series."""
    
    inp = Input(shape=(window_size, n_features))
    
    # Bidirectional LSTM
    lstm_out = Bidirectional(
        LSTM(64, return_sequences=True)
    )(inp)
    lstm_out = Dropout(0.2)(lstm_out)
    
    # Second LSTM layer
    lstm_out2 = Bidirectional(
        LSTM(32, return_sequences=True)
    )(lstm_out)
    
    # Simple attention: learn which timesteps matter most
    attention = Dense(1, activation='tanh')(lstm_out2)
    attention = Flatten()(attention)
    attention = Dense(window_size, activation='softmax')(attention)
    
    # Apply attention weights
    # Reshape attention for element-wise multiplication
    from tensorflow.keras.layers import RepeatVector, Permute, Multiply
    attention = RepeatVector(64)(attention)
    attention = Permute([2, 1])(attention)
    
    context = Multiply()([lstm_out2, attention])
    context = tf.keras.layers.GlobalAveragePooling1D()(context)
    
    # Output
    out = Dense(32, activation='relu')(context)
    out = Dropout(0.2)(out)
    out = Dense(1)(out)
    
    model = Model(inputs=inp, outputs=out)
    model.compile(optimizer='adam', loss='mse', metrics=['mae'])
    
    return model

# Build and display
attn_model = build_attention_lstm(window_size=30)
attn_model.summary()

12.1 ARIMA with statsmodels + Prophet Pipeline

PYTHON
import numpy as np
import pandas as pd
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.stattools import adfuller, kpss
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
import warnings
warnings.filterwarnings('ignore')

# =========================================================
# 1. Load and Prepare Data
# =========================================================
# Using airline passengers dataset as example
from statsmodels.datasets import co2
data = co2.load().data
data = data.resample('M').mean().ffill()
data.columns = ['co2']
print(f"Data shape: {data.shape}")
print(data.head())

# =========================================================
# 2. Stationarity Tests
# =========================================================
def stationarity_report(series, name="Series"):
    """Run ADF and KPSS tests and print results."""
    print(f"\n=== Stationarity Report: {name} ===")
    
    # ADF Test
    adf_result = adfuller(series.dropna(), autolag='AIC')
    print(f"ADF Statistic: {adf_result[0]:.4f}")
    print(f"ADF p-value:   {adf_result[1]:.6f}")
    print(f"ADF → {'Stationary' if adf_result[1] < 0.05 else 'Non-Stationary'}")
    
    # KPSS Test
    kpss_result = kpss(series.dropna(), regression='ct')
    print(f"KPSS Statistic: {kpss_result[0]:.4f}")
    print(f"KPSS p-value:   {kpss_result[1]:.4f}")
    print(f"KPSS → {'Non-Stationary' if kpss_result[1] < 0.05 else 'Stationary'}")

stationarity_report(data['co2'], "CO2 Levels")
stationarity_report(data['co2'].diff().dropna(), "CO2 First Difference")

# =========================================================
# 3. Decomposition
# =========================================================
decomp = seasonal_decompose(data['co2'], model='additive', period=12)
# decomp.plot()  # Uncomment in Jupyter

# =========================================================
# 4. Auto ARIMA using pmdarima
# =========================================================
# pip install pmdarima
import pmdarima as pm

auto_model = pm.auto_arima(
    data['co2'],
    seasonal=True,
    m=12,
    stepwise=True,
    suppress_warnings=True,
    trace=True,
    error_action='ignore',
    information_criterion='aic'
)
print(f"\nBest model: {auto_model.summary()}")

# =========================================================
# 5. Manual ARIMA Fit
# =========================================================
train = data['co2'][:'2000']
test = data['co2']['2001':]

model = ARIMA(train, order=(1, 1, 1), seasonal_order=(1, 1, 1, 12))
fitted = model.fit()
print(fitted.summary())

# Forecast
forecast = fitted.forecast(steps=len(test))
rmse = np.sqrt(np.mean((test.values - forecast.values) ** 2))
print(f"\nSARIMA RMSE: {rmse:.4f}")

# =========================================================
# 6. Exponential Smoothing
# =========================================================
from statsmodels.tsa.holtwinters import ExponentialSmoothing

hw_model = ExponentialSmoothing(
    train,
    trend='add',
    seasonal='add',
    seasonal_periods=12
).fit()

hw_forecast = hw_model.forecast(steps=len(test))
hw_rmse = np.sqrt(np.mean((test.values - hw_forecast.values) ** 2))
print(f"Holt-Winters RMSE: {hw_rmse:.4f}")

12.2 Prophet Pipeline

PYTHON
# pip install prophet
from prophet import Prophet
import pandas as pd
import numpy as np

# =========================================================
# Prophet requires columns 'ds' (date) and 'y' (value)
# =========================================================
# Prepare data
df = data.reset_index()
df.columns = ['ds', 'y']

# Train-test split
train_df = df[df['ds'] < '2001-01-01']
test_df = df[df['ds'] >= '2001-01-01']

# =========================================================
# Build Prophet Model
# =========================================================
prophet_model = Prophet(
    yearly_seasonality=True,
    weekly_seasonality=False,
    daily_seasonality=False,
    changepoint_prior_scale=0.05,    # flexibility of trend
    seasonality_prior_scale=10.0,    # flexibility of seasonality
    interval_width=0.95              # 95% prediction interval
)

# Add Indian holidays (example)
# prophet_model.add_country_holidays(country_name='IN')

prophet_model.fit(train_df)

# =========================================================
# Forecast
# =========================================================
future = prophet_model.make_future_dataframe(periods=len(test_df), freq='M')
forecast = prophet_model.predict(future)

# Evaluate
pred = forecast[forecast['ds'] >= '2001-01-01']['yhat'].values
actual = test_df['y'].values
rmse = np.sqrt(np.mean((actual - pred[:len(actual)]) ** 2))
mape = np.mean(np.abs((actual - pred[:len(actual)]) / actual)) * 100
print(f"Prophet RMSE: {rmse:.4f}")
print(f"Prophet MAPE: {mape:.2f}%")

# =========================================================
# Components
# =========================================================
# prophet_model.plot_components(forecast)  # Uncomment in Jupyter

# =========================================================
# Anomaly Detection with Prophet
# =========================================================
anomalies = forecast[
    (df['y'].values < forecast['yhat_lower'].values) |
    (df['y'].values > forecast['yhat_upper'].values)
]
print(f"\nDetected {len(anomalies)} anomalies")

12.3 VAR Model for Multivariate Forecasting

PYTHON
from statsmodels.tsa.api import VAR
import pandas as pd
import numpy as np

# Simulated multivariate data: Nifty50 + USD/INR
np.random.seed(42)
n = 300
nifty = np.cumsum(np.random.normal(0.5, 50, n)) + 15000
usd_inr = np.cumsum(np.random.normal(0.01, 0.5, n)) + 75
gold = np.cumsum(np.random.normal(0.1, 20, n)) + 50000

df = pd.DataFrame({
    'Nifty50': nifty,
    'USD_INR': usd_inr,
    'Gold': gold
})

# Differencing for stationarity
df_diff = df.diff().dropna()

# Fit VAR
var_model = VAR(df_diff)
# Select optimal lag order
lag_order = var_model.select_order(maxlags=15)
print("Optimal lag orders:")
print(lag_order.summary())

# Fit with optimal lag
fitted_var = var_model.fit(maxlags=5, ic='aic')
print(fitted_var.summary())

# Forecast
forecast_diff = fitted_var.forecast(df_diff.values[-5:], steps=10)
forecast_df = pd.DataFrame(forecast_diff, columns=df.columns)

# Invert differencing
last_values = df.iloc[-1]
forecast_levels = forecast_df.cumsum() + last_values.values
print("\nVAR Forecast (levels):")
print(forecast_levels.round(2))

# Granger Causality Test
from statsmodels.tsa.stattools import grangercausalitytests
print("\n=== Granger Causality: USD/INR → Nifty50 ===")
gc_result = grangercausalitytests(
    df_diff[['Nifty50', 'USD_INR']].dropna(),
    maxlag=5, verbose=True
)

Case Study 1: Nifty50 Stock Index Forecasting

Challenge: Build a daily Nifty50 closing price forecaster for short-term (5-day) predictions.

Approach:

• Data: 10 years of daily Nifty50 closing prices from NSE India (nse-india.com).
• Stationarity: ADF test on log-returns (stationary). Raw prices are non-stationary (unit root).
• Models tested: ARIMA(2,1,2), GARCH(1,1) for volatility, LSTM with 60-day window.
• Feature engineering: Add moving averages (20, 50, 200-day), RSI, Bollinger Bands.

Results:

• ARIMA MAPE: 1.8% (5-day forecast)
• LSTM MAPE: 1.2% (same horizon)
• Key insight: Stock prices follow a random walk — even small improvements matter. LSTM captured regime changes (budget announcements, RBI rate decisions) better than ARIMA.

Lesson: Stock price prediction is extremely difficult (Efficient Market Hypothesis). Models are more useful for volatility estimation and risk management than price direction.

Case Study 2: IMD Monsoon Rainfall Prediction

Context: The Indian monsoon delivers ~70% of annual rainfall between June-September, critical for agriculture (which employs ~42% of the workforce). IMD issues Long Range Forecasts (LRF) in April.

Approach:

• Data: 150+ years of All-India monsoon rainfall (IITM Pune dataset).
• Predictors: ENSO (El Niño), IOD (Indian Ocean Dipole), snow cover, SST.
• Models: SARIMA(1,0,1)(1,1,1,12), Multiple regression with lagged climate indices, Neural network ensemble.
• Evaluation: Categorical accuracy (normal/excess/deficit).

Results: IMD's statistical model correctly categorizes monsoon in ~70% of years. The 2023 El Niño year was correctly predicted as "below normal" with a 3-month lead time.

Case Study 3: COVID-19 India Curve Forecasting

Context: During the 2020-2021 COVID-19 waves, ICMR and IIT research teams built forecasting models to predict daily cases and guide lockdown decisions.

Approach:

• Data: Daily confirmed cases, deaths, recoveries from covid19india.org API.
• Models: SIR/SEIR epidemiological models, ARIMA on log-transformed daily cases, LSTM with mobility data features, Prophet with custom changepoints for lockdown dates.

Key findings: Prophet with manual changepoints at lockdown start/end dates outperformed pure ARIMA. LSTM required careful windowing — the second wave (Delta variant, April 2021) had fundamentally different dynamics from the first wave, causing models trained on Wave 1 to severely underestimate Wave 2 peaks.

Lesson: Time series models assume the future resembles the past. Black swan events (new variants, policy changes) break this assumption. Ensemble methods with scenario analysis work best for crisis forecasting.

Case Study 1: Weather Forecasting (NOAA / ECMWF)

Context: Modern weather forecasting combines physics-based Numerical Weather Prediction (NWP) with statistical/ML post-processing. The ECMWF (European Centre for Medium-Range Weather Forecasts) produces the world's best 10-day forecasts.

ML Role:

• Post-processing NWP output using gradient boosting and neural networks.
• Google DeepMind's GraphCast (2023) uses GNNs to produce 10-day forecasts in under a minute (vs. hours for NWP), matching ECMWF accuracy.
• Time series methods (ARIMA, ETS) still used for localized short-term temperature/wind speed forecasting.

Metrics: RMSE for temperature (typically 1-2°C for 24-hour forecasts), skill scores relative to climatological baselines.

Case Study 2: Energy Demand Forecasting (Electricity Grids)

Context: Power grid operators (like ISO New England, or POSOCO in India) must forecast electricity demand hourly to dispatch generation units efficiently. Over-forecasting wastes fuel; under-forecasting causes blackouts.

Approach:

• Features: Historical load, temperature, humidity, day of week, holidays, special events.
• Models: SARIMA for baseline, gradient boosting for complex patterns, LSTM for intra-day load shape.
• Scale: Forecasts generated for 15-minute, hourly, daily, and weekly horizons.

Results: Modern ensembles achieve MAPE of 1-3% for day-ahead forecasting. The challenge grows with renewable energy integration — solar and wind output depends on weather, adding another layer of uncertainty.

Case Study 3: Retail Sales Forecasting (Walmart / Amazon)

Context: Walmart manages inventory for 4,700+ stores across the US. Accurate demand forecasting reduces waste, prevents stockouts, and saves billions in supply chain costs.

Approach:

• Walmart's M5 competition (2020) on Kaggle: 42,840 hierarchical time series for sales at item-store-state levels.
• Winning solutions used LightGBM with extensive feature engineering: lag features, rolling means, price changes, calendar events, SNAP benefits.
• Prophet used for interpretable decomposition of trends and holiday effects.

Results: Top solutions achieved WRMSSE (Weighted Root Mean Scaled Squared Error) of ~0.52. Key insight: simple lag features + gradient boosting beat deep learning in this competition.

15.1 FinTech — Cash Flow Forecasting

Startups: Razorpay, CashFlo, Recko (India)

Fintech startups help SMEs forecast cash flows using transaction history. Prophet models decompose revenue into trend + weekly + monthly seasonality. Alert systems flag anomalies (sudden drops in UPI collections) using residual-based detection.

15.2 AgriTech — Crop Yield Prediction

Startups: CropIn, Fasal, DeHaat (India)

These startups use satellite imagery time series (NDVI index over crop growth cycles) + weather time series to forecast crop yields. Models include: SARIMA for rainfall patterns, LSTM for multi-source sensor data fusion, and Prophet for price forecasting to advise farmers on optimal selling times.

15.3 HealthTech — Patient Volume Forecasting

Startups: Practo, 1mg

Hospital bed occupancy, emergency room visits, and appointment demand follow time series patterns with weekly seasonality (higher on Mondays), annual seasonality (flu season), and pandemic-driven anomalies. Accurate forecasting enables efficient staff scheduling.

15.4 E-Commerce — Demand Sensing

Startups: Meesho, Udaan, Flipkart

Indian e-commerce faces extreme demand spikes during sales events (Big Billion Days, Great Indian Festival). Time series models augmented with event indicators and real-time search trend data enable "demand sensing" — adjusting forecasts hours/days ahead of traditional weekly cycles.

16.1 RBI — Inflation Forecasting

The Reserve Bank of India uses time series models to forecast CPI inflation for its monetary policy decisions. The RBI's Quarterly Projection Model combines ARIMA with structural economic models. Accurate inflation forecasts determine repo rate changes that affect every Indian borrower.

16.2 NITI Aayog — GDP Forecasting

India's GDP forecasting uses VAR models with high-frequency indicators (IIP, PMI, GST collections, electricity consumption) as predictors. The "nowcasting" approach uses real-time data to estimate current-quarter GDP before official statistics are released.

16.3 Ministry of Health — Epidemic Surveillance

IDSP (Integrated Disease Surveillance Programme) monitors weekly disease incidence across India. Time series anomaly detection flags unusual spikes in dengue, malaria, and respiratory illness cases at the district level, triggering rapid response.

16.4 Smart Cities Mission — Traffic Forecasting

Cities like Pune, Hyderabad, and Bengaluru use traffic sensor time series data to forecast congestion, optimize signal timings, and plan infrastructure. SARIMA models capture daily and weekly patterns, while LSTM models incorporate weather and event data.

17.1 Manufacturing — Predictive Maintenance

Sensor time series (vibration, temperature, pressure) from machines are monitored using anomaly detection. LSTM autoencoders learn "normal" operating patterns and flag deviations. Companies like Tata Steel and Reliance Industries use this to prevent unplanned downtime worth crores per hour.

17.2 Telecom — Network Traffic Forecasting

Jio, Airtel, and Vi forecast network traffic to plan capacity. Time series at the cell tower level show strong daily patterns (peaks at 8-10 PM) and weekly patterns (weekends vs. weekdays). SARIMA and Prophet handle this well, while LSTM captures special events like cricket match streaming.

17.3 Banking — ATM Cash Demand

Banks forecast cash withdrawal patterns at each ATM to optimize replenishment schedules. The series shows strong patterns: salary days (1st and 15th), weekends, festivals, and even election days. Post-demonetization (Nov 2016), models had to rapidly adapt to fundamentally changed withdrawal patterns.

17.4 Pharma — Drug Demand Forecasting

Pharmaceutical companies forecast drug demand considering: seasonal illness patterns (flu season), patent expiry impacts (generic entry), and supply chain lead times. Cipla and Sun Pharma use hierarchical time series models (bottom-up from SKU level to national).

17.5 Aviation — Ticket Price Optimization

IndiGo, Air India, and SpiceJet use time series models for dynamic pricing. Historical booking curves (cumulative bookings vs. days before departure) are modeled per route-class combination. Forecasting residual demand enables revenue optimization (yield management).

Foundational Texts

1. Box, G.E.P., Jenkins, G.M., Reinsel, G.C., & Ljung, G.M. (2015). Time Series Analysis: Forecasting and Control, 5th Edition. Wiley. — The classic reference.
2. Hyndman, R.J. & Athanasopoulos, G. (2021). Forecasting: Principles and Practice, 3rd Edition. OTexts. Free online: otexts.com/fpp3 — The best modern textbook.
3. Hamilton, J.D. (1994). Time Series Analysis. Princeton University Press. — Graduate-level econometric treatment.
4. Shumway, R.H. & Stoffer, D.S. (2017). Time Series Analysis and Its Applications, 4th Edition. Springer.

Key Research Papers

5. Dickey, D.A. & Fuller, W.A. (1979). "Distribution of the Estimators for Autoregressive Time Series with a Unit Root." JASA, 74(366), 427-431.
6. Kwiatkowski, D. et al. (1992). "Testing the null hypothesis of stationarity against the alternative of a unit root." J. of Econometrics, 54(1-3), 159-178.
7. Hochreiter, S. & Schmidhuber, J. (1997). "Long Short-Term Memory." Neural Computation, 9(8), 1735-1780.
8. Taylor, S.J. & Letham, B. (2018). "Forecasting at Scale." The American Statistician, 72(1), 37-45. — The Prophet paper.
9. Salinas, D. et al. (2020). "DeepAR: Probabilistic Forecasting with Autoregressive Recurrent Networks." International J. of Forecasting, 36(3), 1181-1191.
10. Lim, B. et al. (2021). "Temporal Fusion Transformers for Interpretable Multi-horizon Time Series Forecasting." International J. of Forecasting, 37(4), 1748-1764.
11. Oreshkin, B.N. et al. (2020). "N-BEATS: Neural Basis Expansion Analysis for Interpretable Time Series Forecasting." ICLR 2020.
12. Lam, R. et al. (2023). "Learning Skillful Medium-Range Global Weather Forecasting." Science, 382(6677), 1416-1421. — GraphCast.
13. Ansari, A.F. et al. (2024). "Chronos: Learning the Language of Time Series." arXiv:2403.07815.

Indian Context References

14. Rajeevan, M. et al. (2012). "Analysis of variability and trends of extreme rainfall events over India using 104 years of gridded daily rainfall data." Geophysical Research Letters, 39(6).
15. Reserve Bank of India (2023). "Report on Currency and Finance 2022-23." RBI Publications. — Contains RBI's forecasting methodology.
16. Indian Meteorological Department (2024). "End of Season Report — Southwest Monsoon 2023." IMD, Ministry of Earth Sciences.
17. Sahoo, B.B. et al. (2019). "Long short-term memory (LSTM) recurrent neural network for low-flow hydrological time series forecasting." Acta Geophysica, 67, 1471-1481.

Software & Libraries

18. statsmodels: statsmodels.org — Python statistical models including ARIMA, VAR.
19. Prophet: facebook.github.io/prophet — Meta's forecasting library.
20. pmdarima: alkaline-ml.com/pmdarima — Auto-ARIMA for Python.
21. TensorFlow/Keras: tensorflow.org — LSTM implementation.
22. Nixtla: github.com/Nixtla — TimeGPT and StatsForecast libraries.
23. Darts: unit8co.github.io/darts — Unified time series library supporting many models.

Online Resources

24. Kaggle — M5 Forecasting Competition: kaggle.com/c/m5-forecasting-accuracy
25. IITM Pune — Indian Rainfall Data: tropmet.res.in
26. NSE India Historical Data: nseindia.com