Time series data often exhibits patterns such as trends, seasonality, and cycles. Decomposing a time series into its components can help to better understand the underlying processes.
Time series are usually decomposed into three components:
Additive models assume that the components add together to form the time series: \[ Y_t = T_t + S_t + I_t \]
Multiplicative models assume that the components multiply together to form the time series: \[ Y_t = T_t \times S_t \times I_t \]
Sometimes, also a fourth component, Cyclic (\(C_t\)), is considered, which represents long-term oscillations that are not of fixed period.
This example uses the package statsmodels for time series decomposition.
import warnings
import numpy as np
import pandas as pd
import plotly.express as px
import plotly.graph_objs as go
import plotly.io as pio
import plotly.subplots as sp
from statsmodels.tsa.seasonal import seasonal_decompose, STL # used for decomposition
pio.renderers.default = (
"notebook" # set the default plotly renderer to "notebook" (necessary for quarto to render the plots)
)
warnings.filterwarnings("ignore")np.random.seed(42) # for reproducibility
date_range = pd.date_range(start="2010-01-01", end="2024-12-31", freq="MS")
trend = 100 + 0.9 * np.arange(len(date_range)) # nonlinear trend
seasonality = 121 * np.sin(np.pi * (date_range.month / 12)) # yearly seasonality
noise = np.random.normal(0, 21, len(date_range)) # random noise at each point
# Add outliers with a probability of 0.05
outlier_mask = np.random.rand(len(date_range)) < 0.05
outlier_values = np.random.normal(10, 100, len(date_range)) # large outliers
noise[outlier_mask] += outlier_values[outlier_mask]
sales = trend + seasonality + noise # final data as sum of all components
df = pd.DataFrame({"Date": date_range, "Sales": sales}).set_index("Date")fig = px.line(df, y="Sales", title="Sales Over Time")
fig.show()In the above plot, we can see a clear upward trend, but identifying seasonality is more difficult. To support this visual analysis, we add vertical lines at the beginning of each year.
fig = px.line(df, y="Sales", title="Sales Over Time")
for year in range(df.index.year.min(), df.index.year.max() + 1):
fig.add_vline(
x=pd.Timestamp(f"{year}-01-01"),
line_dash="dash",
line_color="gray",
opacity=0.5,
)
fig.show()Now it is apparent that there is a yearly seasonal pattern, with peaks around mid-year.
Algorithmic decomposition of a time series into its components is not a trivial task. Many algorithms exist to achieve this, each with its own advantages and disadvantages.
In the following, we use a very simple additive model based on moving averages to decompose our time series. First, the trend is estimated using a moving average with a smaller window size. Then, the detrended series is calculated by subtracting the trend from the original series. The seasonal component is estimated by averaging the values for each season (e.g., each month, year, ..) in the detrended time series. Finally, the residual component is calculated by subtracting both the trend and seasonal components from the original series.
Note that this is a very naive approach and according to the documentation of seasonal_decompose a more sophisticated method should be preferred. However, it is easy to understand and serves well for demonstration purposes.
decomposed_ts = seasonal_decompose(df, model="additive")decomposed_ts.plot(); # the object has a built-in plot method (based on matplotlib)
Helper function to create a plotly figure with four plots
def seasonal_decompose_plotly(decomposed_ts, title):
"""
Plots a time series decomposition using Plotly.
Parameters:
decomposed_ts (DecomposeResult): Decomposed time series result.
Returns:
plotly object
"""
# Extract components
observed = decomposed_ts.observed
trend = decomposed_ts.trend
seasonal = decomposed_ts.seasonal
resid = decomposed_ts.resid
# Create subplots
fig_decomp = sp.make_subplots(
rows=4,
cols=1,
shared_xaxes=True,
subplot_titles=["Observed", "Trend", "Seasonal", "Residual"],
)
fig_decomp.add_trace(
go.Scatter(x=observed.index, y=observed.values.flatten(), name="Observed"),
row=1,
col=1,
)
fig_decomp.add_trace(go.Scatter(x=trend.index, y=trend, name="Trend"), row=2, col=1)
fig_decomp.add_trace(go.Scatter(x=seasonal.index, y=seasonal, name="Seasonal"), row=3, col=1)
fig_decomp.add_trace(go.Scatter(x=resid.index, y=resid, name="Residual"), row=4, col=1)
fig_decomp.update_layout(height=900, title_text=title)
return fig_decompPlease note the different \(y\)-axes scales in the following plot.
seasonal_decompose_plotly(decomposed_ts, title="Seasonal Decomposition of Time Series (Additive)")A more robust method for time series decomposition is STL (Seasonal and Trend decomposition using Loess).
It not only allows for more flexibility in the decomposition process (e.g. look at seasonal changes over time), but is also more robust to outliers.
decomp_stl = STL(df, robust=True).fit() # robust=True to be less sensitive to expected outliersseasonal_decompose_plotly(decomp_stl, title="STL Decomposition")Having the decomposed components at hand, opens a variety of applications.
Using the trend and seasonal components from the STL decomposition, we can create a denoised version of the original time series.
decomp_stl_estimated = decomp_stl.seasonal + decomp_stl.trendfig = px.line(
decomp_stl_estimated,
title="Estimated versus Observed Sales",
labels={"value": "Sales", "index": "Date"},
)
fig.update_traces(name="Estimated", showlegend=True)
fig.add_trace(go.Scatter(x=df.index, y=df["Sales"], mode="lines", name="Observed"))
fig.show()In previous plot, we can see that the estimated values follow the observed values quite closely, indicating that the STL decomposition has effectively captured the underlying patterns in the data. But we can also see some points with larger deviations, which are likely outliers in the original data.
One way to quantify outliers is to apply a threshold to the residuals (e.g., 3 times the standard deviation of the residuals, also known as the 3-sigma rule).
resid = decomp_stl.resid
std_resid = resid.std()
fig = go.Figure()
fig.add_trace(go.Scatter(x=resid.index, y=resid, mode="lines", name="Residual"))
# Add horizontal lines for ±3 std. dev.
fig.add_hline(
y=resid.mean() + 3 * std_resid,
line_dash="dash",
line_color="blue",
opacity=0.7,
name="+3 Std. Dev.",
)
fig.add_hline(
y=resid.mean() - 3 * std_resid,
line_dash="dash",
line_color="blue",
opacity=0.7,
name="-3 Std. Dev.",
)
# Fill area between the lines
fig.add_traces(
[
go.Scatter(
x=resid.index,
y=[resid.mean() + 3 * std_resid] * len(resid),
mode="lines",
line=dict(width=0),
showlegend=False,
hoverinfo="skip",
),
go.Scatter(
x=resid.index,
y=[resid.mean() - 3 * std_resid] * len(resid),
mode="lines",
fill="tonexty",
fillcolor="rgba(0,0,255,0.08)",
line=dict(width=0),
showlegend=False,
hoverinfo="skip",
),
]
)
fig.update_layout(title="Residual with ±3 Std. Dev. Area", yaxis_title="Residual")
fig.show()So the outliers are:
outliers = df[(resid > resid.mean() + 3 * std_resid) | (resid < resid.mean() - 3 * std_resid)]
outliers| Sales | |
|---|---|
| Date | |
| 2014-09-01 | 447.316646 |
| 2016-07-01 | 139.561235 |
| 2018-03-01 | 489.942096 |
| 2021-06-01 | 227.748919 |
fig = px.line(
df,
y="Sales",
title="Sales with Outliers Highlighted",
labels={"value": "Sales", "index": "Date"},
)
fig.add_scatter(
x=outliers.index,
y=outliers["Sales"],
mode="markers",
marker=dict(color="red", size=10, symbol="x"),
name="Outliers",
)
fig.show()Using the decomposed components also allows for detrending the time series, which can be useful for further analysis.
# Detrend the observed series by subtracting the trend component
detrended = decomp_stl.observed["Sales"] - decomp_stl.trend
fig = px.line(
detrended,
title="Detrended Sales Over Time",
labels={"value": "Detrended Sales", "index": "Date"},
)
fig.show()Another important application is to deseasonalize the time series, which is often a prerequisite for forecasting models.
# Deseasonalize the observed series by subtracting the seasonal component
deseasonalized = decomp_stl.observed["Sales"] - decomp_stl.seasonal
fig = px.line(
deseasonalized,
title="Deseasonalized Sales Over Time",
labels={"value": "Deseasonalized Sales", "index": "Date"},
)
fig.show()