ChatWhole Learn — Free Programming & Data Science Tutorials

ChatWhole Learn

📈 Time Series Analysis in PySpark

Time Series Decomposition Architecture

Architecture Diagram

Detailed Explanation

What is Time Series Analysis?

Time series analysis involves analyzing time-ordered data points to extract meaningful statistics, identify patterns, and forecast future values. In PySpark, time series processing is parallelized across partitions, enabling analysis of petabyte-scale temporal datasets.

Key Takeaway: PySpark's distributed processing enables time series analysis at scales impossible with single-machine tools like pandas or R, while maintaining the same analytical rigor.

Temporal Windowing Operations

PySpark supports several windowing strategies for temporal aggregation:

Window Type	Description	Use Case
Tumbling	Fixed-size, non-overlapping	Hourly/daily aggregations
Hopping	Fixed-size, periodic start	Rolling averages, overlapping aggregations
Session	Dynamic, activity-based	User session analytics
Sliding	Row-count or time-based	Moving averages, lag features

Stationarity and Differencing

Most time series models require stationary data (constant mean and variance):

Test	Null Hypothesis	p-value Threshold
ADF (Augmented Dickey-Fuller)	Unit root exists	p < 0.05 = stationary
KPSS	Series is stationary	p > 0.05 = stationary
Phillips-Perron	Unit root exists	p < 0.05 = stationary

Differencing Techniques:

First-order differencing:
Seasonal differencing: where is the seasonal period
Log transformation: Stabilizes variance for exponential growth patterns

Time Series Decomposition

Decomposition separates a time series into components:

Component	Description	Detection Method
Trend	Long-term direction	Moving average, LOESS
Seasonal	Periodic patterns	Autocorrelation, spectral analysis
Residual	Unexplained noise	After removing trend and seasonality

Model Selection:

Additive: When seasonal variation is constant ()
Multiplicative: When seasonal variation increases with level ()

Feature Engineering for Time Series

Feature Type	Description	PySpark Function
Lag Features	Previous time step values	`lag(col, n)`
Rolling Statistics	Moving averages, std dev	`avg(col).over(Window)`
Calendar Features	Day of week, month, quarter	`dayofweek()`, `month()`
Exponential Weighted	Recency-weighted averages	`ewm()`
Fourier Features	Seasonal periodicity	`sin(2πft)`, `cos(2πft)`

Forecasting Methods

Method	Complexity	Accuracy	Parallelizable
ARIMA	O(n²)	High	Limited
Prophet	O(n log n)	High	Yes
Random Forest	O(n log n)	Medium-High	Yes
Gradient Boosting	O(n log n)	High	Yes
Deep Learning	O(n × d)	Variable	Yes

Best Practice: For distributed time series forecasting, prefer Prophet or ML-based approaches over ARIMA, which requires sequential processing.

Seasonal Patterns Detection

Autocorrelation Function (ACF):

Measures correlation between a time series and its lagged values
Significant spikes at lag indicate seasonality with period

Partial Autocorrelation Function (PACF):

Measures direct correlation at each lag, removing intermediate effects
Helps determine ARIMA order parameters (p, d, q)

Mathematical Foundations

Definition: Time Series

A time series is a sequence of observations indexed by time, where each represents the value at time .

ARIMA Model

The ARIMA(p, d, q) model combines autoregressive (AR), differencing (I), and moving average (MA) components:

where is the differenced series, are AR parameters, are MA parameters, and is white noise.

Stationarity Condition

A time series is weakly stationary if:

The Augmented Dickey-Fuller test statistic is:

Exponential Weighted Moving Average

EWMA gives exponentially decreasing weight to older observations:

where is the smoothing parameter.

Forecast Accuracy Metrics

Metric	Range	Best For
MAE	[0, ∞)	Interpretable scale
RMSE	[0, ∞)	Penalizes large errors
MAPE	[0, 100%]	Scale-independent
SMAPE	[0, 100%]	Symmetric errors

Key Insight

Time series decomposition assumes the observed series is a combination of interpretable components. STL (Seasonal and Trend decomposition using Loess) is robust to outliers and handles any seasonal period, making it preferred over classical decomposition for production forecasting.

Summary

Time series analysis in PySpark combines distributed windowing operations for feature engineering with statistical models for forecasting. Stationarity testing (ADF/KPSS) determines differencing requirements. Decomposition separates trend, seasonality, and residuals for interpretable analysis. Forecast accuracy metrics (MAE, RMSE, MAPE) guide model selection.

Key Concepts Table

Concept	Description	Performance Impact
Tumbling Window	Non-overlapping fixed-size intervals	O(n) processing per partition
Hopping Window	Periodic start with fixed size	O(n × window_size/step)
Lag Features	Previous time step values	O(n) with Window function
Rolling Aggregation	Moving statistics over window	O(n × window_size)
ADF Test	Augmented Dickey-Fuller stationarity test	O(n) with O(p) parameters
Seasonal Decomposition	Split into T+S+E components	O(n) with STL
Autocorrelation	Correlation at different lags	O(n × max_lag)
Prophet	Additive model with trend + seasonality	O(n log n) training
Calendar Features	Day/month/quarter encoding	O(n) computation
Fourier Features	Sinusoidal seasonality encoding	O(n × k) for k harmonics

Code Examples

Temporal Windowing and Feature Engineering

# Save as script.py, then run:
spark-submit \
  --master yarn \
  --deploy-mode client \
  --driver-memory 4g \
  --executor-memory 8g \
  --executor-cores 4 \
  script.py

python

from pyspark.sql import SparkSession
from pyspark.sql.functions import *
from pyspark.sql.window import Window
from pyspark.sql.types import *

spark = SparkSession.builder \
    .appName("TimeSeriesAnalysis") \
    .getOrCreate()

# Create sample time series data (sensor readings)
import random
random.seed(42)

timestamps = [f"2024-01-{d:02d} {h:02d}:00:00" 
              for d in range(1, 32) for h in range(24)]
n = len(timestamps)

data = [
    (timestamps[i], 
     20 + 5 * (i / 24) + 3 * ((i % 24) - 12) / 12 + random.gauss(0, 1),
     50 + 10 * (i / 24) + random.gauss(0, 2))
    for i in range(n)
]

df = spark.createDataFrame(data, ["timestamp", "temperature", "humidity"])
df = df.withColumn("ts", to_timestamp("timestamp"))

# Define windows for various temporal aggregations
# 24-hour moving average
window_24h = Window.orderBy("ts").rowsBetween(-23, 0)

# 7-day rolling statistics
window_7d = Window.orderBy("ts").rowsBetween(-167, 0)

# Add temporal features
ts_features = df \
    .withColumn("temp_24h_avg", avg("temperature").over(window_24h)) \
    .withColumn("temp_24h_std", stddev("temperature").over(window_24h)) \
    .withColumn("humidity_7d_avg", avg("humidity").over(window_7d)) \
    .withColumn("temp_lag_1h", lag("temperature", 1).over(Window.orderBy("ts"))) \
    .withColumn("temp_lag_24h", lag("temperature", 24).over(Window.orderBy("ts"))) \
    .withColumn("temp_diff_1h", col("temperature") - lag("temperature", 1).over(Window.orderBy("ts"))) \
    .withColumn("temp_ewm", avg("temperature").over(
        Window.orderBy("ts").rowsBetween(-23, 0))) \
    .withColumn("hour_of_day", hour("ts")) \
    .withColumn("day_of_week", dayofweek("ts")) \
    .withColumn("day_of_month", dayofmonth("ts")) \
    .withColumn("month", month("ts"))

ts_features.show(50, truncate=False)

Seasonal Decomposition and Differencing

# Save as script.py, then run:
spark-submit \
  --master yarn \
  --deploy-mode client \
  --driver-memory 4g \
  --executor-memory 8g \
  --executor-cores 4 \
  script.py

python

from pyspark.sql import SparkSession
from pyspark.sql.functions import *
from pyspark.sql.window import Window

spark = SparkSession.builder \
    .appName("SeasonalDecomposition") \
    .getOrCreate()

# Create daily sales data with trend and seasonality
import math
import random

random.seed(42)
dates = [f"2024-{m:02d}-{d:02d}" 
         for m in range(1, 13) for d in range(1, 29)]

data = []
for i, date in enumerate(dates):
    trend = 100 + 0.5 * i
    seasonality = 20 * math.sin(2 * math.pi * i / 7)  # Weekly pattern
    monthly_seasonality = 15 * math.sin(2 * math.pi * i / 30)  # Monthly pattern
    noise = random.gauss(0, 5)
    value = trend + seasonality + monthly_seasonality + noise
    data.append((date, value))

df = spark.createDataFrame(data, ["date_str", "value"])
df = df.withColumn("date", to_date("date_str"))

# Decomposition using moving averages
# Step 1: Trend (7-day moving average)
window_7 = Window.orderBy("date").rowsBetween(-3, 3)
df_with_trend = df.withColumn("trend", avg("value").over(window_7))

# Step 2: Deseasonalize (value - trend)
df_deseasonalized = df_with_trend.withColumn(
    "deseasonalized", col("value") - col("trend"))

# Step 3: Seasonal component (average of deseasonalized by day of week)
seasonal_pattern = df_deseasonalized \
    .withColumn("dow", dayofweek("date")) \
    .groupBy("dow") \
    .agg(avg("deseasonalized").alias("seasonal_factor"))

# Step 4: Join seasonal pattern back
df_seasonal = df_deseasonalized \
    .withColumn("dow", dayofweek("date")) \
    .join(seasonal_pattern, "dow") \
    .withColumn("residual", col("deseasonalized") - col("seasonal_factor"))

print("=== Decomposition Results ===")
df_seasonal.select("date", "value", "trend", "seasonal_factor", "residual") \
    .orderBy("date").show(20, truncate=False)

# Seasonal differencing
window_seasonal = Window.orderBy("date")
df_seasonal_diff = df.withColumn(
    "seasonal_diff", 
    col("value") - lag("value", 7).over(window_seasonal)
)

print("=== Seasonal Differencing (period=7) ===")
df_seasonal_diff.show(20, truncate=False)

ARIMA Model Fitting

# Save as script.py, then run:
spark-submit \
  --master yarn \
  --deploy-mode client \
  --driver-memory 4g \
  --executor-memory 8g \
  --executor-cores 4 \
  script.py

python

from pyspark.sql import SparkSession
from pyspark.sql.functions import *
from pyspark.sql.types import *
import numpy as np

spark = SparkSession.builder \
    .appName("ARIMA_Forecast") \
    .getOrCreate()

# Create monthly time series with trend and seasonality
import math
import random

random.seed(42)
months = 36
data = []
for i in range(months):
    trend = 100 + 2 * i
    season = 10 * math.sin(2 * math.pi * i / 12)
    noise = random.gauss(0, 3)
    data.append((f"2024-{(i%12)+1:02d}-01", trend + season + noise))

df = spark.createDataFrame(data, ["date_str", "value"])
df = df.withColumn("date", to_date("date_str"))

# Collect for local ARIMA fitting (ARIMA requires sequential processing)
values = [row.value for row in df.collect()]

# Simple AR(1) model estimation using least squares
def fit_ar1(series):
    """Fit AR(1) model: Y_t = c + phi * Y_{t-1} + epsilon"""
    n = len(series)
    y = np.array(series[1:])
    x = np.array(series[:-1])
    
    # Least squares estimation
    x_mean = np.mean(x)
    y_mean = np.mean(y)
    phi = np.sum((x - x_mean) * (y - y_mean)) / np.sum((x - x_mean) ** 2)
    c = y_mean - phi * x_mean
    
    # Residuals and forecast
    residuals = y - (c + phi * x)
    sigma = np.std(residuals)
    
    return c, phi, sigma, residuals

c, phi, sigma, residuals = fit_ar1(values)

print(f"AR(1) Model Parameters:")
print(f"  Intercept (c): {c:.4f}")
print(f"  AR coefficient (phi): {phi:.4f}")
print(f"  Residual std: {sigma:.4f}")

# Forecast next 6 months
forecast_values = []
last_value = values[-1]
for i in range(6):
    forecast = c + phi * last_value
    forecast_values.append(forecast)
    last_value = forecast

# Create forecast DataFrame
forecast_dates = [f"2025-{m:02d}-01" for m in range(1, 7)]
forecast_df = spark.createDataFrame(
    list(zip(forecast_dates, forecast_values)),
    ["date_str", "forecast_value"]
).withColumn("date", to_date("date_str"))

print("\n=== 6-Month Forecast ===")
forecast_df.show(truncate=False)

Prophet-like Distributed Forecasting

# Save as script.py, then run:
spark-submit \
  --master yarn \
  --deploy-mode client \
  --driver-memory 4g \
  --executor-memory 8g \
  --executor-cores 4 \
  script.py

python

from pyspark.sql import SparkSession
from pyspark.sql.functions import *
from pyspark.sql.window import Window
import numpy as np

spark = SparkSession.builder \
    .appName("ProphetStyleForecast") \
    .getOrCreate()

# Simulate Prophet's additive model: Y = T + S + E
# with distributed computation of components

# Create daily sales data
import math
import random

random.seed(42)
dates = [f"2024-{m:02d}-{d:02d}" 
         for m in range(1, 13) for d in range(1, 29)]

data = []
for i, date in enumerate(dates):
    # Trend: linear growth with change points
    if i < 100:
        trend = 100 + 0.5 * i
    else:
        trend = 150 + 0.8 * (i - 100)  # Change point at i=100
    
    # Weekly seasonality
    weekly = 15 * math.sin(2 * math.pi * i / 7)
    
    # Annual seasonality
    annual = 25 * math.sin(2 * math.pi * i / 365)
    
    noise = random.gauss(0, 5)
    value = trend + weekly + annual + noise
    data.append((date, value))

df = spark.createDataFrame(data, ["date_str", "y"])
df = df.withColumn("ds", to_date("date_str"))

# Step 1: Trend estimation using piecewise linear fit
# In production, this uses gradient descent optimization
window_trend = Window.orderBy("ds").rowsBetween(-14, 14)
df_trend = df.withColumn("trend", avg("y").over(window_trend))

# Step 2: Seasonal decomposition
# Weekly seasonality (period = 7)
df_with_dow = df_trend.withColumn("dow", dayofweek("ds"))

weekly_pattern = df_with_dow \
    .withColumn("detrended", col("y") - col("trend")) \
    .groupBy("dow") \
    .agg(avg("detrended").alias("weekly_effect"))

# Monthly seasonality (period = 30)
df_with_dom = df_trend.withColumn("dom", dayofmonth("ds"))

monthly_pattern = df_with_dom \
    .withColumn("detrended", col("y") - col("trend")) \
    .groupBy("dom") \
    .agg(avg("detrended").alias("monthly_effect"))

# Step 3: Join patterns and compute residuals
df_components = df_trend \
    .withColumn("dow", dayofweek("ds")) \
    .withColumn("dom", dayofmonth("ds")) \
    .join(weekly_pattern, "dow") \
    .join(monthly_pattern, "dom") \
    .withColumn("seasonal", col("weekly_effect") + col("monthly_effect")) \
    .withColumn("residual", col("y") - col("trend") - col("seasonal"))

print("=== Prophet-Style Decomposition ===")
df_components.select("ds", "y", "trend", "weekly_effect", "monthly_effect", "residual") \
    .show(30, truncate=False)

# Step 4: Residual-based uncertainty estimation
residual_stats = df_components.agg(
    stddev("residual").alias("residual_std"),
    avg("residual").alias("residual_mean")
).collect()[0]

print(f"\nResidual Statistics:")
print(f"  Mean: {residual_stats['residual_mean']:.4f}")
print(f"  Std: {residual_stats['residual_std']:.4f}")

# Step 5: Forecast with uncertainty intervals
forecast_days = 30
last_date = df.select(max("ds")).collect()[0][0]

future_dates = [f"2025-01-{d:02d}" for d in range(1, forecast_days + 1)]
future_df = spark.createDataFrame(
    [(d,) for d in future_dates], ["date_str"]
).withColumn("ds", to_date("date_str"))

# Simple linear trend extrapolation
trend_slope = 0.5  # From trend estimation
last_trend = df_components.orderBy(col("ds").desc()).select("trend").first()[0]

future_trend = future_df.withColumn(
    "trend", 
    lit(last_trend) + lit(trend_slope) * (datediff(col("ds"), lit(last_date))))
).withColumn("dow", dayofweek("ds")) \
 .join(weekly_pattern, "dow") \
 .withColumn("forecast", col("trend") + col("weekly_effect"))

# Add uncertainty intervals (95% CI = ±1.96σ)
z_score = 1.96
future_with_ci = future_trend \
    .withColumn("lower_ci", col("forecast") - lit(z_score * residual_stats['residual_std'])) \
    .withColumn("upper_ci", col("forecast") + lit(z_score * residual_stats['residual_std']))

print("\n=== 30-Day Forecast with 95% CI ===")
future_with_ci.select("ds", "forecast", "lower_ci", "upper_ci") \
    .orderBy("ds").show(30, truncate=False)

Anomaly Detection in Time Series

# Save as script.py, then run:
spark-submit \
  --master yarn \
  --deploy-mode client \
  --driver-memory 4g \
  --executor-memory 8g \
  --executor-cores 4 \
  script.py

python

from pyspark.sql import SparkSession
from pyspark.sql.functions import *
from pyspark.sql.window import Window

spark = SparkSession.builder \
    .appName("TimeSeriesAnomalyDetection") \
    .getOrCreate()

# Create time series with anomalies
import random
import math

random.seed(42)
n_points = 1000
data = []
for i in range(n_points):
    value = 50 + 10 * math.sin(2 * math.pi * i / 100) + random.gauss(0, 2)
    # Inject anomalies
    if i in [150, 300, 500, 750]:
        value += 30  # Point anomaly
    if 400 <= i <= 420:
        value += 15  # Contextual anomaly
    data.append((f"2024-01-01 {i//24:02d}:{i%24:02d}:00", value))

df = spark.createDataFrame(data, ["timestamp_str", "value"])
df = df.withColumn("ts", to_timestamp("timestamp_str"))

# Method 1: Z-Score based anomaly detection
window_stats = Window.orderBy("ts").rowsBetween(-50, 50)

df_zscore = df \
    .withColumn("rolling_mean", avg("value").over(window_stats)) \
    .withColumn("rolling_std", stddev("value").over(window_stats)) \
    .withColumn("z_score", 
        (col("value") - col("rolling_mean")) / col("rolling_std")) \
    .withColumn("is_anomaly_zscore", 
        when(abs(col("z_score")) > 3, True).otherwise(False))

# Method 2: IQR-based anomaly detection
window_iqr = Window.orderBy("ts").rowsBetween(-50, 50)

df_iqr = df \
    .withColumn("q1", expr("percentile_approx(value, 0.25)").over(window_iqr)) \
    .withColumn("q3", expr("percentile_approx(value, 0.75)").over(window_iqr)) \
    .withColumn("iqr", col("q3") - col("q1")) \
    .withColumn("lower_bound", col("q1") - 1.5 * col("iqr")) \
    .withColumn("upper_bound", col("q3") + 1.5 * col("iqr")) \
    .withColumn("is_anomaly_iqr",
        when((col("value") < col("lower_bound")) | 
             (col("value") > col("upper_bound")), True).otherwise(False))

# Method 3: Moving average deviation
window_ma = Window.orderBy("ts").rowsBetween(-24, 0)

df_ma = df \
    .withColumn("ma_24h", avg("value").over(window_ma)) \
    .withColumn("deviation", abs(col("value") - col("ma_24h"))) \
    .withColumn("dev_std", stddev(col("deviation")).over(window_ma)) \
    .withColumn("is_anomaly_ma",
        when(col("deviation") > 3 * col("dev_std"), True).otherwise(False))

# Combine methods (ensemble approach)
anomalies = df_zscore.select("ts", "value", "is_anomaly_zscore") \
    .join(df_iqr.select("ts", "is_anomaly_iqr"), "ts") \
    .join(df_ma.select("ts", "is_anomaly_ma"), "ts") \
    .withColumn("anomaly_score",
        (col("is_anomaly_zscore").cast("int") + 
         col("is_anomaly_iqr").cast("int") + 
         col("is_anomaly_ma").cast("int"))) \
    .withColumn("is_anomaly",
        when(col("anomaly_score") >= 2, True).otherwise(False))

print("=== Detected Anomalies ===")
anomalies.filter(col("is_anomaly") == True) \
    .select("ts", "value", "anomaly_score") \
    .orderBy("ts").show(20, truncate=False)

# Summary statistics
total_anomalies = anomalies.filter(col("is_anomaly") == True).count()
print(f"\nTotal anomalies detected: {total_anomalies}")
print(f"Anomaly rate: {total_anomalies / n_points * 100:.2f}%")

Performance Metrics

Operation	1K Points	100K Points	10M Points	100M Points
Moving Average (window=24)	< 1 sec	1-3 sec	10-20 sec	2-5 min
Lag Feature Generation	< 1 sec	1-2 sec	5-15 sec	1-3 min
Seasonal Decomposition	< 1 sec	2-5 sec	30-60 sec	5-10 min
ADF Stationarity Test	< 1 sec	< 1 sec	1-3 sec	5-10 sec
Autocorrelation (100 lags)	< 1 sec	1-2 sec	10-20 sec	2-5 min
Prophet-style Decomposition	1-3 sec	10-20 sec	2-5 min	30-60 min
ARIMA Fitting (local)	1-5 sec	5-15 sec	30-60 sec	5-10 min
Anomaly Detection (ensemble)	< 1 sec	2-5 sec	20-40 sec	3-8 min
Calendar Feature Extraction	< 1 sec	< 1 sec	1-3 sec	5-10 sec
Fourier Feature Generation	< 1 sec	1-2 sec	5-10 sec	1-2 min

Best Practices

Always sort by timestamp before applying window functions to ensure correct temporal ordering across partitions
Use repartition(num_partitions, timestamp_column) to co-locate temporal data and reduce shuffle during window operations
Pre-filter time ranges before feature engineering to avoid computing lag/rolling features on unnecessary data
Handle missing timestamps explicitly with fillna() or interpolation—gaps in time series break window calculations
Use pandas_udf for complex time series models (ARIMA, Prophet) that cannot be vectorized across partitions
Partition by time period (day, week, month) to enable partition pruning for time-range queries
Cache intermediate DataFrames when performing multiple temporal analyses on the same dataset
Avoid collecting large time series to the driver—use distributed algorithms for fitting models
Use approximate functions for real-time anomaly detection when exact quantiles are not required
Implement watermarking for streaming time series to handle late-arriving data gracefully
Validate temporal ordering with assert df.filter(col("ts") < lag("ts")).count() == 0 before analysis
Use calendar table joins for complex date logic instead of chaining multiple date functions

See also: ML Pipeline Integration (24), ML Feature Engineering (32), Advanced Aggregations (31), Geospatial Data (27)

Time Series Analysis in PySpark

📈 Time Series Analysis in PySpark

Time Series Decomposition Architecture

Architecture Diagram

Detailed Explanation

What is Time Series Analysis?

Temporal Windowing Operations

Stationarity and Differencing

Time Series Decomposition

Feature Engineering for Time Series

Forecasting Methods

Seasonal Patterns Detection

Mathematical Foundations

Definition: Time Series

ARIMA Model

Stationarity Condition

Exponential Weighted Moving Average

Forecast Accuracy Metrics

Key Insight

Summary

Key Concepts Table

Code Examples

Temporal Windowing and Feature Engineering

Seasonal Decomposition and Differencing

ARIMA Model Fitting

Prophet-like Distributed Forecasting

Anomaly Detection in Time Series

Performance Metrics

Best Practices

See Also

Rate & Feedback

Related Tutorials

Time Series Forecasting: Trend, Seasonality, ARIMA, and Prophet

Time Series Forecasting with LLMs

13. Window Operations in PySpark

Apache Spark Provider Integration with Airflow

PySpark at Scale: DataFrame API, Shuffles, and Adaptive Execution

Industrial IoT Platforms and Edge Computing

Continue Learning

Geospatial Data Processing in PySpark

Time Series Forecasting: Trend, Seasonality, ARIMA, and Prophet

Time Series Forecasting with LLMs

13. Window Operations in PySpark

Apache Spark Provider Integration with Airflow

Need Expert PySpark Help?