INTRODUCTION AND FOUNDATIONAL CONCEPTS
Time series data represents a sequence of observations recorded at successive points in time. Unlike traditional datasets where each observation is independent, time series data exhibits temporal dependencies where the value at one time point often depends on previous values. This temporal structure is what makes time series analysis fundamentally different from standard statistical analysis and requires specialized techniques to extract meaningful insights.
The most common examples of time series data include sensor readings from industrial equipment that measure temperature, pressure, or vibration at regular intervals; sales numbers recorded daily, weekly, or monthly that show how a business's revenue fluctuates over time; and share values or stock prices that are recorded at specific trading intervals. Each of these examples contains valuable patterns that, when properly analyzed, can lead to accurate forecasts and actionable insights.
Time series data can be categorized into two primary types based on structure. Univariate time series contains a single variable measured over time, such as the daily closing price of a particular stock. Multivariate time series contains multiple variables measured simultaneously, such as the opening price, closing price, trading volume, and volatility all recorded for the same stock on the same day. Understanding which type of data you are working with is crucial because it determines which analytical techniques are most appropriate.
The fundamental characteristic that makes time series analysis challenging and interesting is autocorrelation, which means that observations are not independent but rather correlated with their own past values. A temperature reading at 2 PM is likely to be very close to the temperature reading at 1 PM because weather changes gradually. A stock price at 10:01 AM is likely to be similar to the price at 10:00 AM because markets do not typically jump dramatically in one minute. This autocorrelation structure is what allows us to make predictions about future values based on historical patterns.
KEY COMPONENTS OF TIME SERIES DATA
Before analyzing any time series, it is essential to understand the components that make up the data. Most real-world time series can be decomposed into distinct components, each representing different patterns in the data.
The trend component represents the long-term direction of the data. A sales number series might show an overall upward trend if a business is growing, indicating that sales are increasing over the long term despite short-term fluctuations. The trend can be linear, meaning it increases or decreases at a constant rate, or non-linear, meaning the rate of change itself changes over time. Identifying and understanding the trend is critical because it represents the underlying pattern that your forecasts should capture.
The seasonal component represents regular, repeating patterns that occur at fixed intervals. Many businesses experience seasonal patterns, such as increased sales during holiday periods or reduced sales during summer months. Temperature data often exhibits seasonal patterns where summer is consistently warmer than winter. Share values might show seasonal patterns related to quarterly earnings reports. Seasonality has a fixed period, meaning if data is collected monthly and shows a twelve-month seasonal pattern, that pattern will repeat every twelve months.
The noise or irregular component represents random fluctuations that cannot be explained by trend or seasonality. This is the unpredictable element in the data, such as a sudden equipment failure causing a spike in sensor readings, unexpected news causing a stock price jump, or an unplanned promotional event boosting sales. Distinguishing between true patterns and noise is one of the fundamental challenges in time series analysis.
Autocorrelation measures how strongly observations in the time series correlate with their own past values at different lags. If autocorrelation is high, it means that values separated by several time periods are still related to each other. Understanding autocorrelation helps in choosing appropriate forecasting models and in understanding how far back in history you need to look to predict future values.
Stationarity is a critical property that describes whether the statistical properties of a time series remain constant over time. A stationary time series has a constant mean, constant variance, and constant autocorrelation structure. Many traditional forecasting models, particularly ARIMA models, require the input data to be stationary. Most real-world time series are non-stationary because they have trends or changing means. Converting non-stationary data to stationary data through differencing or transformation is often a necessary preprocessing step.
DATA PREPARATION AND VALIDATION
The quality of your analysis depends entirely on the quality of your data. Time series data preparation is a complex process that involves multiple steps, each of which requires careful attention.
The first step in data preparation is understanding your data source and collection process. You must verify that the data was collected at regular intervals, understand what each variable represents, and confirm the units of measurement. It is also important to understand the business context, such as whether the data includes holidays, special events, or periods when the system was not operating normally.
The second step is handling missing values, which are common in real-world time series data. Missing values can occur due to sensor failures, data transmission errors, or systems being offline. Simple approaches like forward filling, where you copy the last known value forward, or linear interpolation, where you estimate the missing value based on surrounding values, can work for small gaps. For longer gaps, more sophisticated imputation methods may be necessary. The key principle is that the method you choose should preserve the underlying patterns in the data while not introducing artificial patterns.
Consider a simple temperature sensor that records readings every hour. If one reading is missing, linear interpolation would estimate the missing value as the average of the readings before and after the gap. However, if ten consecutive readings are missing, linear interpolation might not be appropriate because it assumes a linear relationship, when in reality temperature might change non-linearly during that period.
import numpy as np
def interpolate_missing_values(time_series, method='linear'):
"""
Handle missing values represented as NaN in the input array.
Parameters:
time_series: NumPy array with possible NaN values
method: 'linear' for linear interpolation or 'forward_fill' for
forward fill approach
Returns:
NumPy array with missing values filled
"""
if method == 'linear':
nans = np.isnan(time_series)
x = lambda z: z.nonzero()[0]
time_series[nans] = np.interp(
x(nans),
x(~nans),
time_series[~nans]
)
elif method == 'forward_fill':
mask = np.isnan(time_series)
idx = np.where(~mask, np.arange(len(mask)), 0)
np.maximum.accumulate(idx, axis=0, out=idx)
time_series[mask] = time_series[idx[mask]]
return time_series
Outlier detection is the next critical step. Outliers are values that deviate significantly from the expected pattern and can distort your analysis and forecasts. Outliers can be caused by measurement errors, sensor malfunctions, or genuine unusual events. A simple statistical approach identifies values that are more than three standard deviations from the mean as potential outliers, though this assumes a normal distribution. More sophisticated methods use the interquartile range, where values below the first quartile minus 1.5 times the interquartile range or above the third quartile plus 1.5 times the interquartile range are considered outliers. The appropriate method depends on your data distribution.
def detect_outliers_iqr(time_series):
"""
Detect outliers using the Interquartile Range method.
This method is robust to non-normal distributions and identifies
values that fall outside the typical range of the data.
Parameters:
time_series: NumPy array of values
Returns:
Boolean array where True indicates an outlier
"""
Q1 = np.percentile(time_series, 25)
Q3 = np.percentile(time_series, 75)
IQR = Q3 - Q1
lower_bound = Q1 - 1.5 * IQR
upper_bound = Q3 + 1.5 * IQR
return (time_series < lower_bound) | (time_series > upper_bound)
Normalization and scaling are important because different variables often have different units and ranges. A temperature sensor might record values between 0 and 100 degrees Celsius, while a sales number might be between 10,000 and 100,000 dollars. When these different variables are used together in machine learning models, the model might give inappropriate weight to the variable with the larger scale. Standardization, also known as z-score normalization, transforms each variable to have a mean of zero and a standard deviation of one. Min-max scaling transforms values to a range between zero and one. The choice depends on your specific use case, but standardization is more common for machine learning models.
def standardize_data(time_series):
"""
Apply z-score normalization to center data with unit variance.
This transformation is essential when combining multiple variables
with different scales in machine learning models.
Parameters:
time_series: NumPy array of values
Returns:
Tuple of (standardized array, mean, standard deviation)
"""
mean = np.mean(time_series)
std = np.std(time_series)
standardized = (time_series - mean) / std
return standardized, mean, std
Splitting your data into training and testing sets is crucial for developing and evaluating models. Unlike general machine learning where data points can be randomly shuffled, time series requires careful temporal splitting. You must ensure that training data comes before test data in time. A common approach is to use the first eighty percent of your historical data for training and the remaining twenty percent for testing. This preserves the temporal order and ensures that you are not accidentally using future information to predict the past, which would be unrealistic.
EXPLORATORY DATA ANALYSIS FOR TIME SERIES
Exploratory data analysis for time series goes beyond simple summary statistics and requires visualization and temporal analysis techniques.
Visualization is the most powerful tool for understanding time series data. A simple line plot showing how your variable changes over time immediately reveals trends and seasonality. If the plot shows an upward trend, downward trend, or no clear direction, this is valuable information. If the plot shows repeating patterns at regular intervals, seasonality is present. If there are sudden spikes or drops that break the pattern, these are outliers or special events.
import matplotlib.pyplot as plt
def plot_time_series(time_points, values, title, ylabel):
"""
Create a basic line plot for time series visualization.
This helps identify trends, seasonality, and anomalies at a glance.
Parameters:
time_points: Array of time indices or datetime objects
values: Array of observed values
title: String for plot title
ylabel: String for y-axis label
"""
plt.figure(figsize=(12, 6))
plt.plot(time_points, values, linewidth=1.5)
plt.title(title)
plt.ylabel(ylabel)
plt.xlabel('Time')
plt.grid(True, alpha=0.3)
plt.show()
Autocorrelation analysis measures how the time series correlates with lagged versions of itself. The autocorrelation function computes the correlation between observations at different time lags. If autocorrelation at lag one is high and close to one, it means that consecutive observations are very similar. If autocorrelation decreases quickly as lag increases, it means that only recent history matters for predicting the future. If autocorrelation has a pattern with peaks at regular intervals matching your suspected seasonal period, this confirms seasonality. The autocorrelation function plot, also called a correlogram, is an essential diagnostic tool.
def compute_autocorrelation(time_series, max_lag):
"""
Compute autocorrelation values for different lags.
This reveals how strongly the series correlates with its own
past values and helps identify seasonality patterns.
Parameters:
time_series: NumPy array of values
max_lag: Maximum lag to compute correlation for
Returns:
Array of autocorrelation values for each lag
"""
mean = np.mean(time_series)
variance = np.var(time_series)
autocorr = np.zeros(max_lag + 1)
for lag in range(max_lag + 1):
if lag == 0:
autocorr[lag] = 1.0
else:
covariance = np.mean(
(time_series[:-lag] - mean) *
(time_series[lag:] - mean)
)
autocorr[lag] = covariance / variance
return autocorr
Seasonal decomposition separates a time series into its trend, seasonal, and residual components. This is useful for understanding the relative strength of each component and for removing seasonality before fitting certain models. The classical decomposition method, also known as seasonal and trend decomposition using Loess (STL), fits a trend using a moving average and then removes the trend to identify the seasonal component. The remaining residual represents the noise or irregular component.
Statistical tests for stationarity help you determine whether your data has a constant mean and variance over time. The Augmented Dickey-Fuller test is a common statistical test where the null hypothesis is that the time series has a unit root, meaning it is non-stationary. A p-value less than 0.05 typically leads to rejecting this null hypothesis and concluding that the series is stationary. If your test shows non-stationarity, you can apply differencing, where you replace each value with the difference between it and the previous value, until the series becomes stationary.
def first_difference(time_series):
"""
Apply first-order differencing to make series more stationary.
Differencing is essential preprocessing for many time series models
that require stationary input data.
Parameters:
time_series: NumPy array of values
Returns:
Array of differenced values (one element shorter than input)
"""
return np.diff(time_series)
TRADITIONAL STATISTICAL APPROACHES TO TIME SERIES FORECASTING
Before the era of artificial intelligence and machine learning, statisticians developed sophisticated methods for time series analysis and forecasting. These methods remain valuable today and often serve as baseline comparisons for newer approaches.
The Moving Average method is one of the simplest approaches. It forecasts future values as the simple average of the most recent observations. A simple moving average with a window of seven days would average the values from the past seven days to predict the next day's value. This approach is extremely simple to understand and implement but cannot capture trends or seasonality. It is useful primarily for smoothing out short-term noise and is often used as a baseline for comparison with more sophisticated methods.
def simple_moving_average(time_series, window_size):
"""
Compute simple moving average for smoothing and denoising.
This basic approach is useful for removing noise while preserving
the general trend. Works well for stationary data without
significant seasonality.
Parameters:
time_series: NumPy array of values
window_size: Number of periods to average over
Returns:
Array of moving average values
"""
moving_avg = np.convolve(
time_series,
np.ones(window_size) / window_size,
mode='valid'
)
return moving_avg
Exponential Smoothing methods assign higher weights to more recent observations while giving less weight to older observations. The Simple Exponential Smoothing method uses a smoothing parameter, typically denoted alpha and set between zero and one, to control how quickly the method responds to new observations. A high alpha value means the forecast reacts quickly to new data, while a low alpha value means the forecast changes slowly. Exponential smoothing methods can be extended to handle trends with Double Exponential Smoothing and seasonality with Triple Exponential Smoothing, also known as Holt-Winters method. These methods work well when your data has clear trend or seasonal components that are relatively stable over time.
ARIMA, which stands for AutoRegressive Integrated Moving Average, is a powerful and flexible family of models that can handle a wide variety of time series patterns. The AutoRegressive component means the model uses past values to predict future values. The Integrated component means the model can handle non-stationary data by differencing. The Moving Average component means the model uses past forecast errors to improve predictions. An ARIMA model is specified with three parameters, ARIMA(p,d,q), where p is the number of autoregressive terms, d is the number of times differencing is applied, and q is the number of moving average terms. Finding the optimal values for these parameters typically involves examining the autocorrelation and partial autocorrelation functions and trying different combinations.
ARIMA models work well for univariate time series with clear patterns but have limitations for multivariate data and for very long time horizons. They also assume that the patterns in the historical data will continue into the future, which may not be true if the underlying process changes.
ARTIFICIAL INTELLIGENCE AND MACHINE LEARNING FOR TIME SERIES
Artificial intelligence and machine learning have revolutionized time series analysis by enabling models to learn complex, non-linear patterns from data. The advantage of AI approaches is that they require fewer assumptions about the underlying process and can adapt to changing patterns.
One of the most important concepts for applying machine learning to time series is the creation of supervised learning datasets from the raw sequential data. A time series is inherently a sequence, but most machine learning algorithms expect independent observations. To bridge this gap, you create overlapping windows or sequences from your time series where each window's inputs are historical values and the target is the next value or next several values. For example, if you have a series of daily sales numbers and you want to forecast sales three days into the future, you might create samples where each sample contains the past seven days of sales as input features and the sales on day eight as the target. Creating such datasets requires careful consideration of how much historical data is relevant, known as the look-back window, and how far into the future you want to predict, known as the forecasting horizon.
def create_supervised_dataset(time_series, look_back, forecast_horizon):
"""
Transform time series into supervised learning dataset.
This enables using machine learning algorithms on sequential data
by creating (input, target) pairs from the original series.
Parameters:
time_series: NumPy array of values
look_back: Number of past time steps to use as input
forecast_horizon: Number of future time steps to predict
Returns:
Tuple of (X, y) arrays where X contains input sequences
and y contains corresponding target values
"""
X = []
y = []
for i in range(len(time_series) - look_back - forecast_horizon + 1):
X.append(time_series[i:i + look_back])
y.append(time_series[i + look_back + forecast_horizon - 1])
return np.array(X), np.array(y)
Long Short-Term Memory networks, called LSTMs, are a specialized type of recurrent neural network designed to capture long-term dependencies in sequential data. Traditional neural networks struggle with sequences because information from distant past time steps becomes increasingly hard to propagate through the network. LSTMs address this through a sophisticated memory mechanism with gates that control what information is retained, forgotten, and used. An LSTM cell contains a memory cell that can store information over long periods, an input gate that controls whether new information should be added to memory, a forget gate that controls whether information should be discarded from memory, and an output gate that controls what information from memory should be used for prediction. This architecture makes LSTMs particularly powerful for time series forecasting where historical context from many time steps in the past might be relevant.
The training process for neural networks involves passing data through the network in small batches, computing predictions, comparing predictions to actual values using a loss function such as mean squared error, and using backpropagation to adjust the network's weights to reduce the error. This process is repeated for multiple epochs until the network learns the patterns in the data. Choosing the right batch size, learning rate, and number of epochs requires experimentation and monitoring of validation performance.
Transformer-based models represent the newest frontier in time series analysis. Transformers use an attention mechanism that allows the model to focus on different parts of the input sequence when making predictions. The attention mechanism computes weights for each position in the input sequence, indicating how relevant that position is for predicting the target. Unlike LSTMs that process sequences sequentially, transformers can process entire sequences in parallel, leading to faster training. Recent research has shown that transformer-based models can outperform LSTMs on many time series tasks, particularly when you have large amounts of historical data available.
Gradient boosting methods such as XGBoost and LightGBM have shown surprising effectiveness for time series forecasting. These methods build many decision trees sequentially, with each new tree learning to correct the errors of the previous trees. For time series, you create a supervised dataset as described earlier and then apply the gradient boosting algorithm. The advantage of these methods is that they handle non-linear relationships naturally, require less data than deep learning methods, and often provide good performance without extensive hyperparameter tuning. They are particularly effective when you have external features, such as day of week or holiday indicators, that affect your target variable.
Prophet is a forecasting tool developed by Facebook that combines classical statistical components with machine learning. Prophet automatically fits a model with trend, seasonality, and holiday components. It is designed to be robust to missing data and outliers and requires minimal data preprocessing. Prophet works well for business time series data with clear seasonal patterns and is particularly useful when you need a model that is easy to interpret and explain to business stakeholders.
STEP-BY-STEP FORECASTING WORKFLOW
A systematic approach to time series forecasting ensures that you develop robust, reliable models. The following steps should be followed in order, though you may need to iterate and return to earlier steps if results are unsatisfactory.
The first step is problem definition and data collection. You must clearly define what you are trying to forecast, the time horizon of interest, and how accurate the forecast needs to be. Different problems have different requirements. Forecasting sales for inventory planning might require high accuracy, while forecasting foot traffic for rough capacity planning might accept lower accuracy. You must also define your target variable clearly. For example, are you forecasting daily sales, weekly sales, or monthly sales? Different target definitions might require different approaches.
The second step is exploratory data analysis as discussed in the previous section. This includes visualization, checking for stationarity, identifying autocorrelation patterns, and detecting anomalies. This step builds your understanding of the data and informs your modeling choices.
The third step is data preparation including handling missing values, removing or adjusting outliers, and normalizing data. The specific techniques you use should be guided by your domain knowledge and your findings from exploratory analysis.
The fourth step is creating training and test sets with proper temporal ordering. For some methods, you might also create a validation set separate from the test set to tune hyperparameters during development.
The fifth step is selecting and implementing candidate models. It is good practice to start with simple models as baselines, such as simple moving average or exponential smoothing, and then move to more complex models. This helps you understand how much improvement more complex models actually provide.
The sixth step is training your models on the training data. For traditional statistical methods, this might mean fitting parameters. For machine learning methods, this means running optimization algorithms to learn weights or tree structures.
The seventh step is evaluating your models on the test set using appropriate metrics. For regression problems like time series forecasting, common metrics include mean absolute error, which is the average absolute difference between predicted and actual values; root mean squared error, which penalizes larger errors more heavily than smaller errors; and mean absolute percentage error, which gives error as a percentage of the actual value. The choice of metric should reflect what matters in your business context.
def evaluate_forecast(actual, predicted):
"""
Calculate standard metrics for evaluating forecast accuracy.
Parameters:
actual: Array of actual observed values
predicted: Array of predicted values
Returns:
Dictionary containing MAE, RMSE, and MAPE metrics
"""
mae = np.mean(np.abs(actual - predicted))
rmse = np.sqrt(np.mean((actual - predicted) ** 2))
mape = np.mean(
np.abs((actual - predicted) / actual)
) * 100
return {
'MAE': mae,
'RMSE': rmse,
'MAPE': mape
}
The eighth step is hyperparameter tuning. Most models have parameters that affect their behavior and performance. Grid search or random search can be used to try different parameter combinations and select the combination that performs best on validation data. Be careful not to overfit by tuning excessively on the test set.
The ninth step is creating prediction intervals or confidence intervals around your point forecasts. A point forecast gives a single predicted value, but in reality there is uncertainty in any forecast. Prediction intervals quantify this uncertainty by providing a range within which you expect the actual value to fall with a specified probability, such as ninety-five percent. Different methods have different approaches for computing prediction intervals.
The tenth step is monitoring and updating your model in production. As new data arrives, you should monitor whether the model's forecasts remain accurate. If accuracy degrades, you may need to retrain the model with more recent data or switch to a different approach.
FINDING CORRELATIONS BETWEEN TIME SERIES DATA SEQUENCES
Many real-world problems require understanding relationships between multiple time series. A retail company might want to know whether advertising spending influences sales numbers. A power utility might want to know whether temperature and time of day together drive electricity demand. A pharmaceutical company might want to know whether a particular drug's sales correlate with physician search interest. Identifying these correlations is essential for understanding causal relationships and building more accurate multivariate forecasting models.
The simplest approach is Pearson correlation, which measures the linear relationship between two variables. If you have two time series, both of which have been standardized to have mean zero and unit variance, the Pearson correlation coefficient ranges from negative one to positive one, where positive one indicates perfect positive correlation, negative one indicates perfect negative correlation, and zero indicates no linear relationship. Computing Pearson correlation is straightforward but has limitations. It only captures linear relationships, so two variables might have a strong non-linear relationship but show low Pearson correlation. Pearson correlation assumes both series are stationary; if both series have strong trends in the same direction, they might show high correlation even if changes in one do not actually cause changes in the other.
def compute_pearson_correlation(series1, series2):
"""
Compute Pearson correlation coefficient between two time series.
This measures the strength of linear relationship, with values
ranging from -1 to 1.
Parameters:
series1: First time series as NumPy array
series2: Second time series as NumPy array
Returns:
Correlation coefficient value
"""
mean1 = np.mean(series1)
mean2 = np.mean(series2)
numerator = np.sum((series1 - mean1) * (series2 - mean2))
denominator = np.sqrt(
np.sum((series1 - mean1) ** 2) *
np.sum((series2 - mean2) ** 2)
)
if denominator == 0:
return 0
return numerator / denominator
Lagged correlation analysis addresses one critical limitation of simple correlation: the two series might be related but with a time delay. For example, advertising spending might influence sales, but not immediately. There is typically a lag between when an advertisement is shown and when customers make purchases. To find lagged correlations, you compute the correlation between one series and the other series shifted by one, two, three, or more time steps. The lag at which the highest correlation occurs indicates the typical delay between cause and effect. This technique requires careful interpretation because just because correlation is highest at a particular lag does not necessarily mean causation. It could be that both series are responding to a third factor.
def compute_lagged_correlations(series1, series2, max_lag):
"""
Compute correlation between series1 and series2 at various lags.
This identifies delayed relationships between time series, such as
the delay between a cause and its effect.
Parameters:
series1: First time series as NumPy array
series2: Second time series as NumPy array
max_lag: Maximum lag to compute correlation for
Returns:
Array of correlation values for each lag
"""
correlations = []
for lag in range(max_lag + 1):
if lag == 0:
corr = compute_pearson_correlation(series1, series2)
else:
corr = compute_pearson_correlation(
series1[lag:],
series2[:-lag]
)
correlations.append(corr)
return np.array(correlations)
Cross-correlation is a more general concept that measures similarity between two series at different lags. While Pearson correlation measures only the strength of linear relationship, cross-correlation can reveal how similar the patterns in two series are, even if they have different magnitudes. Cross-correlation is computed by sliding one series across the other and computing the sum of products at each position. The position where the sum of products is largest indicates the lag at which the two series are most similar.
Dynamic Time Warping is a technique that measures similarity between two time series that may have different lengths or speeds. Unlike Pearson correlation which assumes both series are sampled at the same rate and have the same length, dynamic time warping can align sequences that are stretched or compressed in time. This is particularly useful for sensor data where the same event might take different amounts of time to complete. For example, a machine might perform a cycle in thirty seconds on one occasion and forty seconds on another, but the sequence of sensor readings should show the same pattern. Dynamic time warping computes a distance matrix that measures point-wise similarity between all pairs of points in the two series, then finds the warping path through this matrix that minimizes total distance while respecting the constraint that time can move forward but not backward.
def dynamic_time_warping(series1, series2):
"""
Compute Dynamic Time Warping distance between two series.
DTW measures similarity between series that may have different
lengths or temporal speeds, enabling pattern matching in complex
scenarios.
Parameters:
series1: First time series as NumPy array
series2: Second time series as NumPy array
Returns:
DTW distance (lower values indicate more similarity)
"""
n, m = len(series1), len(series2)
distance_matrix = np.full((n + 1, m + 1), np.inf)
distance_matrix[0, 0] = 0
for i in range(1, n + 1):
for j in range(1, m + 1):
cost = abs(series1[i-1] - series2[j-1])
distance_matrix[i, j] = cost + min(
distance_matrix[i-1, j],
distance_matrix[i, j-1],
distance_matrix[i-1, j-1]
)
return distance_matrix[n, m]
Granger Causality is a statistical hypothesis test that determines whether past values of one time series help predict another time series beyond what past values of that series alone can predict. The concept is based on the idea that if series X causes series Y, then past values of X should contain information useful for predicting Y. The test regresses Y on its own past values to get a baseline prediction, then regresses Y on its own past values plus past values of X and checks whether including X significantly reduces the prediction error. Granger causality requires that both series are stationary and has limitations; it can only identify temporal precedence, not true causation, and can be misleading if both series are driven by a third variable.
HANDLING SEASONALITY AND TREND IN MULTIVARIATE ANALYSIS
When analyzing multiple time series together, it becomes more complex to handle seasonality and trend. Different series might have different seasonal periods or different trend strengths, making it essential to address these components appropriately.
One approach is seasonal decomposition and deseasonalization before correlation analysis. If you decompose each series into trend, seasonal, and residual components, you can then compute correlations on just the residual components. This removes the confounding effects of seasonality and common trends and allows you to identify correlations based on shorter-term movements. This approach assumes that you want to understand correlations in detrended and deseasonalized data, which is appropriate for many applications.
Another approach is using differencing to remove both trend and seasonality. First differencing removes linear trends, while seasonal differencing, where you compute differences between observations separated by the seasonal period, removes seasonal effects. If your data has both trend and seasonality, you might apply both types of differencing. Once both series are differenced, you can compute correlations or use other analysis techniques.
BEST PRACTICES AND COMMON PITFALLS
Developing good time series forecasting practices helps avoid common mistakes that can lead to inaccurate models or overconfident predictions.
The first best practice is starting with simple models before moving to complex models. Simple moving average or exponential smoothing can provide a surprisingly good baseline. If you immediately jump to complex deep learning models without establishing a simple baseline, you will not understand how much improvement the complexity provides and whether it is justified.
The second best practice is always using temporal ordering when splitting data into training and test sets. Randomly shuffling time series data violates the principle that you cannot use future information to predict the past. Temporal splitting ensures that your model is evaluated on its ability to forecast into the future based only on historical information.
The third best practice is avoiding overfitting. Complex models with many parameters can memorize noise in the training data and perform poorly on new, unseen data. Monitoring validation performance during training and using regularization techniques can help prevent overfitting.
The fourth best practice is being careful about correlation and causation. Just because two time series are correlated does not mean one causes the other. You must think carefully about the underlying mechanisms and use additional analysis techniques like Granger causality to investigate potential causal relationships.
The fifth best practice is understanding and documenting your data. Knowing when data collection started, whether there are known anomalies or special events, whether the system underwent changes that might affect the data, and other contextual information is crucial for proper interpretation and modeling.
The sixth best practice is retesting and retraining models as new data arrives. Time series processes can change over time. A model that worked well historically might degrade as conditions change. Regular monitoring and updating ensures your forecasts remain accurate.
The seventh best practice is creating prediction intervals along with point forecasts. Stakeholders need to understand uncertainty in forecasts, not just a single predicted value.
The eighth best practice is avoiding data leakage in feature engineering. If you are creating features from the time series that should not be available at prediction time, your model will appear to perform well during development but fail in production. For example, if you include information about future values when creating features, your model cannot use those features in production.
PRODUCTION IMPLEMENTATION CONSIDERATIONS
Deploying time series models to production requires careful consideration of practical issues beyond model development.
You must implement data pipelines that collect, validate, and preprocess incoming data reliably. The preprocessing applied in production must exactly match the preprocessing applied during model training. If you standardized the training data, you must standardize production data using the same mean and standard deviation computed from training data, not computed from the current production data.
You must implement monitoring systems that track model performance over time. As new data arrives, you should compute forecast errors and detect if accuracy is degrading. When performance degrades beyond a threshold, you should trigger retraining or investigation into what has changed.
You must implement versioning and rollback procedures for models. When you train a new model version, you should carefully evaluate it against the current production model before switching. You should maintain previous versions in case you need to roll back.
You must consider computational requirements and latency. Some models like neural networks require substantial computation, while others like moving averages are instant. Your production system must meet both latency and throughput requirements.
You must implement logging and alerting to understand when and why your system fails. When a forecast is clearly wrong, logs should help you diagnose the cause.
RUNNING EXAMPLE: COMPLETE PRODUCTION-READY IMPLEMENTATION
This section provides a complete, production-ready implementation that demonstrates all key concepts discussed throughout this article. The implementation is designed to handle multiple use cases including sensor data, sales numbers, and share values.
import numpy as np
import json
from abc import ABC, abstractmethod
from dataclasses import dataclass
from datetime import datetime, timedelta
from typing import Tuple, List, Dict, Optional
from enum import Enum
class TimeSeriesModel(ABC):
"""
Abstract base class for all time series forecasting models.
This defines the interface that all models must implement,
ensuring consistent usage across different model types.
"""
@abstractmethod
def fit(self, X: np.ndarray, y: np.ndarray) -> None:
"""
Train the model on provided data.
Parameters:
X: Training features, shape (n_samples, n_features)
y: Training targets, shape (n_samples,)
"""
pass
@abstractmethod
def predict(self, X: np.ndarray) -> np.ndarray:
"""
Generate predictions for new data.
Parameters:
X: Features to predict on, shape (n_samples, n_features)
Returns:
Predicted values, shape (n_samples,)
"""
pass
@abstractmethod
def get_prediction_intervals(
self,
X: np.ndarray,
confidence: float = 0.95
) -> Tuple[np.ndarray, np.ndarray]:
"""
Provide upper and lower bounds for predictions.
Parameters:
X: Features to predict on
confidence: Confidence level for intervals (0-1)
Returns:
Tuple of (lower_bounds, upper_bounds)
"""
pass
class SimpleMovingAverageModel(TimeSeriesModel):
"""
Simple moving average forecasting model.
This baseline model computes the average of the most recent
observations to predict the next value. It serves as a useful
baseline for comparison with more complex models.
"""
def __init__(self, window_size: int = 7):
"""
Initialize the moving average model.
Parameters:
window_size: Number of recent observations to average
"""
self.window_size = window_size
self.mean = None
self.std = None
self.training_data = None
def fit(self, X: np.ndarray, y: np.ndarray) -> None:
"""
For moving average, fit just stores training data statistics.
The model computes mean and std for prediction intervals.
"""
self.training_data = y
self.mean = np.mean(y)
self.std = np.std(y)
def predict(self, X: np.ndarray) -> np.ndarray:
"""
For SMA, each prediction is the average of recent values.
Since X contains sequences of past values, we take the
average of each sequence as the prediction.
"""
if X.ndim == 1:
X = X.reshape(1, -1)
predictions = np.mean(X, axis=1)
return predictions
def get_prediction_intervals(
self,
X: np.ndarray,
confidence: float = 0.95
) -> Tuple[np.ndarray, np.ndarray]:
"""
Provide intervals based on standard deviation of training data.
"""
predictions = self.predict(X)
z_score = 1.96 if confidence == 0.95 else 2.576
margin = z_score * self.std
lower = predictions - margin
upper = predictions + margin
return lower, upper
class ExponentialSmoothingModel(TimeSeriesModel):
"""
Simple exponential smoothing model.
This model applies higher weight to recent observations through
an exponential decay factor. It adapts to gradual changes better
than simple moving average.
"""
def __init__(self, alpha: float = 0.3):
"""
Initialize the exponential smoothing model.
Parameters:
alpha: Smoothing factor between 0 and 1. Higher values
make the model more responsive to recent changes.
"""
if not 0 <= alpha <= 1:
raise ValueError('alpha must be between 0 and 1')
self.alpha = alpha
self.level = None
self.std = None
self.residuals = None
def fit(self, X: np.ndarray, y: np.ndarray) -> None:
"""
Fit exponential smoothing by computing smoothed level.
"""
self.level = y[0]
smoothed = np.zeros(len(y))
smoothed[0] = y[0]
for i in range(1, len(y)):
self.level = self.alpha * y[i] + (1 - self.alpha) * self.level
smoothed[i] = self.level
self.residuals = y - smoothed
self.std = np.std(self.residuals)
def predict(self, X: np.ndarray) -> np.ndarray:
"""
Exponential smoothing prediction is the current smoothed level.
"""
if X.ndim == 1:
X = X.reshape(1, -1)
batch_size = X.shape[0]
predictions = np.full(batch_size, self.level)
return predictions
def get_prediction_intervals(
self,
X: np.ndarray,
confidence: float = 0.95
) -> Tuple[np.ndarray, np.ndarray]:
"""
Provide intervals based on residual standard deviation.
"""
predictions = self.predict(X)
z_score = 1.96 if confidence == 0.95 else 2.576
margin = z_score * self.std
lower = predictions - margin
upper = predictions + margin
return lower, upper
class NeuralNetworkModel(TimeSeriesModel):
"""
Simple neural network for time series forecasting.
This implements a feedforward neural network with one hidden layer.
For production use cases, more sophisticated architectures like
LSTMs might be preferable, but this demonstrates the principle.
"""
def __init__(
self,
hidden_size: int = 64,
learning_rate: float = 0.01,
epochs: int = 100,
batch_size: int = 32
):
"""
Initialize neural network model.
Parameters:
hidden_size: Number of neurons in hidden layer
learning_rate: Step size for gradient updates
epochs: Number of training iterations
batch_size: Samples per training batch
"""
self.hidden_size = hidden_size
self.learning_rate = learning_rate
self.epochs = epochs
self.batch_size = batch_size
self.W1 = None
self.b1 = None
self.W2 = None
self.b2 = None
self.input_size = None
self.training_std = None
self.training_residuals = None
def _relu(self, x: np.ndarray) -> np.ndarray:
"""Rectified Linear Unit activation function."""
return np.maximum(0, x)
def _relu_derivative(self, x: np.ndarray) -> np.ndarray:
"""Derivative of ReLU for backpropagation."""
return (x > 0).astype(float)
def _forward(
self,
X: np.ndarray
) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
"""
Forward pass through the network.
Returns tuple of (hidden_activation, output, hidden_preactivation)
"""
hidden_preactivation = np.dot(X, self.W1) + self.b1
hidden_activation = self._relu(hidden_preactivation)
output = np.dot(hidden_activation, self.W2) + self.b2
return output, hidden_activation, hidden_preactivation
def _backward(
self,
X: np.ndarray,
y: np.ndarray,
output: np.ndarray,
hidden_activation: np.ndarray,
hidden_preactivation: np.ndarray
) -> None:
"""
Backward pass through the network for gradient computation.
"""
batch_size = X.shape[0]
output_error = output - y.reshape(-1, 1)
dW2 = np.dot(
hidden_activation.T,
output_error
) / batch_size
db2 = np.sum(output_error, axis=0, keepdims=True) / batch_size
hidden_error = np.dot(
output_error,
self.W2.T
) * self._relu_derivative(hidden_preactivation)
dW1 = np.dot(X.T, hidden_error) / batch_size
db1 = np.sum(hidden_error, axis=0, keepdims=True) / batch_size
self.W1 -= self.learning_rate * dW1
self.b1 -= self.learning_rate * db1
self.W2 -= self.learning_rate * dW2
self.b2 -= self.learning_rate * db2
def fit(self, X: np.ndarray, y: np.ndarray) -> None:
"""
Train the neural network model.
"""
self.input_size = X.shape[1]
self.W1 = np.random.randn(
self.input_size,
self.hidden_size
) * 0.01
self.b1 = np.zeros((1, self.hidden_size))
self.W2 = np.random.randn(
self.hidden_size,
1
) * 0.01
self.b2 = np.zeros((1, 1))
for epoch in range(self.epochs):
indices = np.random.permutation(len(X))
for i in range(0, len(X), self.batch_size):
batch_indices = indices[i:i + self.batch_size]
X_batch = X[batch_indices]
y_batch = y[batch_indices]
output, hidden_activation, hidden_preactivation = (
self._forward(X_batch)
)
self._backward(
X_batch,
y_batch,
output,
hidden_activation,
hidden_preactivation
)
output, _, _ = self._forward(X)
self.training_residuals = y - output.flatten()
self.training_std = np.std(self.training_residuals)
def predict(self, X: np.ndarray) -> np.ndarray:
"""
Generate predictions from the trained network.
"""
if X.ndim == 1:
X = X.reshape(1, -1)
output, _, _ = self._forward(X)
return output.flatten()
def get_prediction_intervals(
self,
X: np.ndarray,
confidence: float = 0.95
) -> Tuple[np.ndarray, np.ndarray]:
"""
Provide intervals based on training residuals.
"""
predictions = self.predict(X)
z_score = 1.96 if confidence == 0.95 else 2.576
margin = z_score * self.training_std
lower = predictions - margin
upper = predictions + margin
return lower, upper
class TimeSeriesProcessor:
"""
Comprehensive processor for time series data.
Handles data loading, preprocessing, dataset creation, model
training, and evaluation in a unified interface.
"""
def __init__(self, data: np.ndarray):
"""
Initialize the processor with raw time series data.
Parameters:
data: NumPy array of time series values
"""
self.raw_data = data.copy()
self.standardized_data = None
self.normalization_mean = None
self.normalization_std = None
self.training_data = None
self.test_data = None
self.train_size = None
def handle_missing_values(
self,
method: str = 'linear'
) -> np.ndarray:
"""
Handle missing values represented as NaN.
Parameters:
method: 'linear' for interpolation or 'forward_fill'
Returns:
Array with missing values filled
"""
processed = self.raw_data.copy()
if method == 'linear':
nans = np.isnan(processed)
x = lambda z: z.nonzero()[0]
processed[nans] = np.interp(
x(nans),
x(~nans),
processed[~nans]
)
elif method == 'forward_fill':
mask = np.isnan(processed)
idx = np.where(~mask, np.arange(len(mask)), 0)
np.maximum.accumulate(idx, axis=0, out=idx)
processed[mask] = processed[idx[mask]]
return processed
def detect_outliers(self, method: str = 'iqr') -> np.ndarray:
"""
Detect outliers in the time series.
Parameters:
method: 'iqr' for Interquartile Range method
Returns:
Boolean array indicating outlier positions
"""
if method == 'iqr':
Q1 = np.percentile(self.raw_data, 25)
Q3 = np.percentile(self.raw_data, 75)
IQR = Q3 - Q1
lower_bound = Q1 - 1.5 * IQR
upper_bound = Q3 + 1.5 * IQR
return (self.raw_data < lower_bound) | (
self.raw_data > upper_bound
)
return np.zeros(len(self.raw_data), dtype=bool)
def standardize_data(self) -> np.ndarray:
"""
Apply z-score standardization to the data.
Returns:
Standardized array
"""
self.normalization_mean = np.mean(self.raw_data)
self.normalization_std = np.std(self.raw_data)
self.standardized_data = (
(self.raw_data - self.normalization_mean) /
self.normalization_std
)
return self.standardized_data
def inverse_transform(self, data: np.ndarray) -> np.ndarray:
"""
Reverse standardization to original scale.
Parameters:
data: Standardized data array
Returns:
Data in original scale
"""
return (
data * self.normalization_std + self.normalization_mean
)
def compute_autocorrelation(
self,
max_lag: int
) -> np.ndarray:
"""
Compute autocorrelation function.
Parameters:
max_lag: Maximum lag to compute
Returns:
Array of autocorrelation values for each lag
"""
data = self.standardized_data
if data is None:
data = self.standardize_data()
mean = np.mean(data)
variance = np.var(data)
autocorr = np.zeros(max_lag + 1)
for lag in range(max_lag + 1):
if lag == 0:
autocorr[lag] = 1.0
else:
covariance = np.mean(
(data[:-lag] - mean) * (data[lag:] - mean)
)
autocorr[lag] = covariance / variance if variance > 0 else 0
return autocorr
def first_difference(self) -> np.ndarray:
"""
Apply first-order differencing.
Returns:
Differenced time series
"""
return np.diff(self.raw_data)
def split_train_test(
self,
train_ratio: float = 0.8
) -> Tuple[np.ndarray, np.ndarray]:
"""
Split data into training and test sets.
Parameters:
train_ratio: Proportion of data for training
Returns:
Tuple of (training_data, test_data)
"""
self.train_size = int(len(self.raw_data) * train_ratio)
self.training_data = self.raw_data[:self.train_size]
self.test_data = self.raw_data[self.train_size:]
return self.training_data, self.test_data
def create_supervised_dataset(
self,
data: np.ndarray,
look_back: int,
forecast_horizon: int
) -> Tuple[np.ndarray, np.ndarray]:
"""
Transform time series into supervised learning dataset.
Parameters:
data: Time series array
look_back: Number of past steps to use
forecast_horizon: Number of future steps to predict
Returns:
Tuple of (X, y) arrays for supervised learning
"""
X = []
y = []
for i in range(
len(data) - look_back - forecast_horizon + 1
):
X.append(data[i:i + look_back])
y.append(data[i + look_back + forecast_horizon - 1])
return np.array(X), np.array(y)
def compute_lagged_correlation(
self,
series1: np.ndarray,
series2: np.ndarray,
max_lag: int
) -> np.ndarray:
"""
Compute correlation at different lags.
Parameters:
series1: First time series
series2: Second time series
max_lag: Maximum lag to compute
Returns:
Array of correlation values at each lag
"""
correlations = []
for lag in range(max_lag + 1):
if lag == 0:
corr = self._compute_pearson_correlation(
series1,
series2
)
else:
corr = self._compute_pearson_correlation(
series1[lag:],
series2[:-lag]
)
correlations.append(corr)
return np.array(correlations)
def _compute_pearson_correlation(
self,
series1: np.ndarray,
series2: np.ndarray
) -> float:
"""
Compute Pearson correlation coefficient.
"""
mean1 = np.mean(series1)
mean2 = np.mean(series2)
numerator = np.sum(
(series1 - mean1) * (series2 - mean2)
)
denominator = np.sqrt(
np.sum((series1 - mean1) ** 2) *
np.sum((series2 - mean2) ** 2)
)
if denominator == 0:
return 0
return numerator / denominator
def compute_dynamic_time_warping(
self,
series1: np.ndarray,
series2: np.ndarray
) -> float:
"""
Compute Dynamic Time Warping distance.
Parameters:
series1: First time series
series2: Second time series
Returns:
DTW distance measure
"""
n, m = len(series1), len(series2)
distance_matrix = np.full((n + 1, m + 1), np.inf)
distance_matrix[0, 0] = 0
for i in range(1, n + 1):
for j in range(1, m + 1):
cost = abs(series1[i-1] - series2[j-1])
distance_matrix[i, j] = cost + min(
distance_matrix[i-1, j],
distance_matrix[i, j-1],
distance_matrix[i-1, j-1]
)
return distance_matrix[n, m]
def evaluate_model(
self,
model: TimeSeriesModel,
X_test: np.ndarray,
y_test: np.ndarray
) -> Dict[str, float]:
"""
Evaluate model performance on test data.
Parameters:
model: Trained TimeSeriesModel
X_test: Test features
y_test: Test targets
Returns:
Dictionary of evaluation metrics
"""
predictions = model.predict(X_test)
mae = np.mean(np.abs(predictions - y_test))
rmse = np.sqrt(np.mean((predictions - y_test) ** 2))
mape = np.mean(
np.abs((predictions - y_test) / y_test)
) * 100
return {
'MAE': mae,
'RMSE': rmse,
'MAPE': mape
}
class TimeSeriesForecaster:
"""
End-to-end forecasting system that combines all components.
This provides a high-level interface for loading data, training
models, making predictions, and monitoring performance.
"""
def __init__(self):
"""Initialize the forecaster system."""
self.processor = None
self.models = {}
self.best_model = None
self.best_model_name = None
self.predictions_cache = {}
def load_data(self, data: np.ndarray) -> None:
"""
Load time series data.
Parameters:
data: NumPy array of time series values
"""
self.processor = TimeSeriesProcessor(data)
def preprocess_data(
self,
handle_missing: bool = True,
remove_outliers: bool = False,
standardize: bool = True
) -> np.ndarray:
"""
Apply preprocessing pipeline.
Parameters:
handle_missing: Whether to handle NaN values
remove_outliers: Whether to detect outliers
standardize: Whether to standardize data
Returns:
Preprocessed data array
"""
if handle_missing:
self.processor.raw_data = (
self.processor.handle_missing_values()
)
if remove_outliers:
outliers = self.processor.detect_outliers()
self.processor.raw_data[outliers] = np.nan
self.processor.raw_data = (
self.processor.handle_missing_values()
)
if standardize:
return self.processor.standardize_data()
return self.processor.raw_data
def train_multiple_models(
self,
X_train: np.ndarray,
y_train: np.ndarray,
X_test: np.ndarray,
y_test: np.ndarray
) -> Dict[str, Dict[str, float]]:
"""
Train multiple model types and compare performance.
Parameters:
X_train: Training features
y_train: Training targets
X_test: Test features
y_test: Test targets
Returns:
Dictionary with metrics for each model
"""
results = {}
models_to_train = {
'SimpleMovingAverage': SimpleMovingAverageModel(
window_size=7
),
'ExponentialSmoothing': ExponentialSmoothingModel(
alpha=0.3
),
'NeuralNetwork': NeuralNetworkModel(
hidden_size=64,
learning_rate=0.01,
epochs=100
)
}
for name, model in models_to_train.items():
print(f'Training {name}...')
model.fit(X_train, y_train)
self.models[name] = model
metrics = self.processor.evaluate_model(
model,
X_test,
y_test
)
results[name] = metrics
print(f'{name} - MAE: {metrics["MAE"]:.4f}, '
f'RMSE: {metrics["RMSE"]:.4f}, '
f'MAPE: {metrics["MAPE"]:.2f}%')
best_model_name = min(
results,
key=lambda x: results[x]['MAE']
)
self.best_model_name = best_model_name
self.best_model = self.models[best_model_name]
return results
def generate_forecast(
self,
X_future: np.ndarray,
use_best: bool = True
) -> np.ndarray:
"""
Generate forecasts for future periods.
Parameters:
X_future: Features for future time periods
use_best: Whether to use best model or all models
Returns:
Forecast values
"""
if use_best:
return self.best_model.predict(X_future)
predictions = []
for model in self.models.values():
predictions.append(model.predict(X_future))
return np.mean(np.array(predictions), axis=0)
def get_forecast_intervals(
self,
X_future: np.ndarray,
confidence: float = 0.95
) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
"""
Get point forecasts with prediction intervals.
Parameters:
X_future: Features for future periods
confidence: Confidence level for intervals
Returns:
Tuple of (predictions, lower_bounds, upper_bounds)
"""
predictions = self.best_model.predict(X_future)
lower, upper = self.best_model.get_prediction_intervals(
X_future,
confidence
)
return predictions, lower, upper
def analyze_correlations(
self,
series1: np.ndarray,
series2: np.ndarray,
max_lag: int = 20
) -> Dict[str, any]:
"""
Analyze correlations between two time series.
Parameters:
series1: First time series
series2: Second time series
max_lag: Maximum lag for lagged correlations
Returns:
Dictionary containing correlation analysis results
"""
pearson_corr = self.processor._compute_pearson_correlation(
series1,
series2
)
lagged_corrs = self.processor.compute_lagged_correlation(
series1,
series2,
max_lag
)
best_lag = np.argmax(np.abs(lagged_corrs))
dtw_distance = self.processor.compute_dynamic_time_warping(
series1,
series2
)
return {
'pearson_correlation': pearson_corr,
'lagged_correlations': lagged_corrs,
'best_lag': best_lag,
'correlation_at_best_lag': lagged_corrs[best_lag],
'dtw_distance': dtw_distance
}
def run_sensor_data_example():
"""
Example use case: forecasting temperature sensor readings.
This demonstrates the complete workflow for a typical industrial
monitoring scenario.
"""
print('=' * 70)
print('EXAMPLE 1: Temperature Sensor Data Analysis')
print('=' * 70)
np.random.seed(42)
days = 365
trend = np.linspace(15, 25, days)
seasonality = 8 * np.sin(np.arange(days) * 2 * np.pi / 365)
noise = np.random.normal(0, 0.5, days)
temperature = trend + seasonality + noise
forecaster = TimeSeriesForecaster()
forecaster.load_data(temperature)
preprocessed = forecaster.preprocess_data(
handle_missing=True,
remove_outliers=False,
standardize=True
)
train_data, test_data = forecaster.processor.split_train_test(0.8)
look_back = 30
forecast_horizon = 1
X_train, y_train = forecaster.processor.create_supervised_dataset(
train_data,
look_back,
forecast_horizon
)
X_test, y_test = forecaster.processor.create_supervised_dataset(
test_data,
look_back,
forecast_horizon
)
results = forecaster.train_multiple_models(
X_train,
y_train,
X_test,
y_test
)
print(f'\nBest model: {forecaster.best_model_name}')
future_input = X_test[-1:, :]
forecast, lower, upper = forecaster.get_forecast_intervals(
future_input,
confidence=0.95
)
forecast_original = forecaster.processor.inverse_transform(forecast)
lower_original = forecaster.processor.inverse_transform(lower)
upper_original = forecaster.processor.inverse_transform(upper)
print(f'\nNext day temperature forecast:')
print(f' Point forecast: {forecast_original[0]:.2f} degrees')
print(f' 95% prediction interval: '
f'[{lower_original[0]:.2f}, {upper_original[0]:.2f}]')
def run_sales_data_example():
"""
Example use case: forecasting daily sales numbers.
This demonstrates the system applied to business metrics
with stronger seasonality.
"""
print('\n' + '=' * 70)
print('EXAMPLE 2: Daily Sales Data Analysis')
print('=' * 70)
np.random.seed(43)
days = 365
base_sales = 5000
trend = np.linspace(0, 2000, days)
seasonality = 1000 * np.sin(np.arange(days) * 2 * np.pi / 7)
monthly_seasonality = 500 * np.sin(
np.arange(days) * 2 * np.pi / 30
)
noise = np.random.normal(0, 200, days)
sales = (
base_sales + trend + seasonality +
monthly_seasonality + noise
)
sales = np.maximum(sales, 0)
forecaster = TimeSeriesForecaster()
forecaster.load_data(sales)
preprocessed = forecaster.preprocess_data(
handle_missing=False,
remove_outliers=False,
standardize=True
)
train_data, test_data = forecaster.processor.split_train_test(0.8)
look_back = 14
forecast_horizon = 1
X_train, y_train = forecaster.processor.create_supervised_dataset(
train_data,
look_back,
forecast_horizon
)
X_test, y_test = forecaster.processor.create_supervised_dataset(
test_data,
look_back,
forecast_horizon
)
results = forecaster.train_multiple_models(
X_train,
y_train,
X_test,
y_test
)
print(f'\nBest model: {forecaster.best_model_name}')
future_input = X_test[-1:, :]
forecast, lower, upper = forecaster.get_forecast_intervals(
future_input,
confidence=0.95
)
forecast_original = forecaster.processor.inverse_transform(forecast)
lower_original = forecaster.processor.inverse_transform(lower)
upper_original = forecaster.processor.inverse_transform(upper)
print(f'\nNext day sales forecast:')
print(f' Point forecast: ${forecast_original[0]:.2f}')
print(f' 95% prediction interval: '
f'[${lower_original[0]:.2f}, ${upper_original[0]:.2f}]')
def run_stock_price_example():
"""
Example use case: forecasting stock prices.
This demonstrates the system applied to financial data
with different characteristics.
"""
print('\n' + '=' * 70)
print('EXAMPLE 3: Stock Price Data Analysis')
print('=' * 70)
np.random.seed(44)
days = 252
returns = np.random.normal(0.0005, 0.02, days)
price = 100 * np.exp(np.cumsum(returns))
forecaster = TimeSeriesForecaster()
forecaster.load_data(price)
preprocessed = forecaster.preprocess_data(
handle_missing=False,
remove_outliers=False,
standardize=True
)
train_data, test_data = forecaster.processor.split_train_test(0.8)
look_back = 20
forecast_horizon = 1
X_train, y_train = forecaster.processor.create_supervised_dataset(
train_data,
look_back,
forecast_horizon
)
X_test, y_test = forecaster.processor.create_supervised_dataset(
test_data,
look_back,
forecast_horizon
)
results = forecaster.train_multiple_models(
X_train,
y_train,
X_test,
y_test
)
print(f'\nBest model: {forecaster.best_model_name}')
future_input = X_test[-1:, :]
forecast, lower, upper = forecaster.get_forecast_intervals(
future_input,
confidence=0.95
)
forecast_original = forecaster.processor.inverse_transform(forecast)
lower_original = forecaster.processor.inverse_transform(lower)
upper_original = forecaster.processor.inverse_transform(upper)
print(f'\nNext day stock price forecast:')
print(f' Point forecast: ${forecast_original[0]:.2f}')
print(f' 95% prediction interval: '
f'[${lower_original[0]:.2f}, ${upper_original[0]:.2f}]')
def run_correlation_example():
"""
Example use case: analyzing correlation between two time series.
This demonstrates how to identify relationships and potential
causal effects between different variables.
"""
print('\n' + '=' * 70)
print('EXAMPLE 4: Time Series Correlation Analysis')
print('=' * 70)
np.random.seed(45)
days = 365
advertising = 1000 + 500 * np.sin(
np.arange(days) * 2 * np.pi / 30
) + np.random.normal(0, 50, days)
advertising = np.maximum(advertising, 0)
lag_period = 7
sales_base = 5000 + 0.003 * np.concatenate(
[np.zeros(lag_period), advertising[:-lag_period]]
)
sales = (
sales_base + 1000 * np.sin(
np.arange(days) * 2 * np.pi / 7
) +
np.random.normal(0, 200, days)
)
sales = np.maximum(sales, 0)
forecaster = TimeSeriesForecaster()
forecaster.load_data(sales)
correlation_results = forecaster.analyze_correlations(
advertising,
sales,
max_lag=30
)
print(f'\nPearson correlation (lag 0): '
f'{correlation_results["pearson_correlation"]:.4f}')
print(f'Best lag: {correlation_results["best_lag"]} days')
print(f'Correlation at best lag: '
f'{correlation_results["correlation_at_best_lag"]:.4f}')
print(f'DTW distance: {correlation_results["dtw_distance"]:.2f}')
print(f'\nLagged correlations (sample lags):')
for lag in [0, 5, 7, 10, 15, 20]:
corr = correlation_results['lagged_correlations'][lag]
print(f' Lag {lag:2d}: {corr:.4f}')
if __name__ == '__main__':
run_sensor_data_example()
run_sales_data_example()
run_stock_price_example()
run_correlation_example()
print('\n' + '=' * 70)
print('ALL EXAMPLES COMPLETED SUCCESSFULLY')
print('=' * 70)
CONCLUSION
Time series analysis is a complex field that combines statistical theory, machine learning, and domain expertise. The techniques described in this article range from simple moving averages that anyone can understand and implement to sophisticated neural networks that can capture intricate patterns in data.
The key to successful time series forecasting is understanding your specific problem, preparing your data carefully, trying multiple approaches, and evaluating results honestly. Simple models often provide a good baseline, and more complex models should only be used if they demonstrably outperform simpler alternatives and justify their increased complexity and computational cost.
Artificial intelligence and machine learning have made time series analysis more accessible and powerful, but they do not eliminate the need for careful data preparation, understanding of underlying processes, and skeptical evaluation of results. The best practitioners combine deep understanding of both the time series techniques and the business domain they are working in.
As you apply these concepts to real-world problems such as sensor monitoring, sales forecasting, or financial analysis, remember that the most important aspects are getting high-quality data, asking clear questions, and iterating based on results. Start simple, validate assumptions, and only add complexity when justified by improved performance and business value.
