"""Seasonal Decomposition by Moving Averages"""
import numpy as np
import pandas as pd
from pandas.core.nanops import nanmean as pd_nanmean
from statsmodels.tools.validation import PandasWrapper, array_like
from statsmodels.tsa.filters.filtertools import convolution_filter
from statsmodels.tsa.tsatools import freq_to_period
__all__ = [
"DecomposeResult",
"seasonal_decompose",
"seasonal_mean",
]
def _extrapolate_trend(trend, npoints):
"""
Replace nan values on trend's end-points with least-squares extrapolated
values with regression considering npoints closest defined points.
Parameters
----------
trend : ndarray
The trend series, which may contain NaN values at the start and/or
end that will be replaced.
npoints : int
The number of closest defined points on each end to use in the
extrapolating least-squares regression.
Returns
-------
ndarray
The trend array with NaN values at the endpoints replaced by
extrapolated values.
"""
front = next(i for i, vals in enumerate(trend) if not np.any(np.isnan(vals)))
back = (
trend.shape[0]
- 1
- next(i for i, vals in enumerate(trend[::-1]) if not np.any(np.isnan(vals)))
)
front_last = min(front + npoints, back)
back_first = max(front, back - npoints)
k, n = np.linalg.lstsq(
np.c_[np.arange(front, front_last), np.ones(front_last - front)],
trend[front:front_last],
rcond=-1,
)[0]
extra = (np.arange(0, front) * np.c_[k] + np.c_[n]).T
if trend.ndim == 1:
extra = extra.squeeze()
trend[:front] = extra
k, n = np.linalg.lstsq(
np.c_[np.arange(back_first, back), np.ones(back - back_first)],
trend[back_first:back],
rcond=-1,
)[0]
extra = (np.arange(back + 1, trend.shape[0]) * np.c_[k] + np.c_[n]).T
if trend.ndim == 1:
extra = extra.squeeze()
trend[back + 1 :] = extra
return trend
def seasonal_mean(x, period):
"""
Return means for each period in x
Parameters
----------
x : array_like
Time series values.
period : int
The number of periods per cycle, e.g., 12 for monthly data.
Returns
-------
ndarray
The mean value for each period, ignoring NaNs.
"""
return np.array([pd_nanmean(x[i::period], axis=0) for i in range(period)])
[docs]
def seasonal_decompose(
x,
model="additive",
filt=None,
period=None,
two_sided=True,
extrapolate_trend=0,
):
"""
Seasonal decomposition using moving averages.
Parameters
----------
x : array_like
Time series. If 2d, individual series are in columns. x must contain 2
complete cycles.
model : {"additive", "multiplicative"}, optional
Type of seasonal component. Abbreviations are accepted.
filt : array_like, optional
The filter coefficients for filtering out the seasonal component.
The concrete moving average method used in filtering is determined by
two_sided.
period : int, optional
Period of the series (eg, 1 for annual, 4 for quarterly, etc). Must be
used if x is not a pandas object or if the index of x does not have a
frequency. Overrides default periodicity of x if x is a pandas object
with a timeseries index.
two_sided : bool, optional
The moving average method used in filtering.
If True (default), a centered moving average is computed using the
filt. If False, the filter coefficients are for past values only.
extrapolate_trend : int or 'freq', optional
If set to > 0, the trend resulting from the convolution is
linear least-squares extrapolated on both ends (or the single one
if two_sided is False) considering this many (+1) closest points.
If set to 'freq', use `freq` closest points. Setting this parameter
results in no NaN values in trend or resid components.
Returns
-------
DecomposeResult
An object with seasonal, trend, and resid attributes.
See Also
--------
statsmodels.tsa.filters.bk_filter.bkfilter
Baxter-King filter.
statsmodels.tsa.filters.cf_filter.cffilter
Christiano-Fitzgerald asymmetric, random walk filter.
statsmodels.tsa.filters.hp_filter.hpfilter
Hodrick-Prescott filter.
statsmodels.tsa.filters.convolution_filter
Linear filtering via convolution.
statsmodels.tsa.seasonal.STL
Season-Trend decomposition using LOESS.
Notes
-----
This is a naive decomposition. More sophisticated methods should
be preferred.
The additive model is Y[t] = T[t] + S[t] + e[t]
The multiplicative model is Y[t] = T[t] * S[t] * e[t]
The results are obtained by first estimating the trend by applying
a convolution filter to the data. The trend is then removed from the
series and the average of this de-trended series for each period is
the returned seasonal component.
"""
pfreq = period
pw = PandasWrapper(x)
if period is None:
if isinstance(x, (pd.Series, pd.DataFrame)):
index = x.index
if isinstance(index, pd.PeriodIndex):
pfreq = index.freq
else:
pfreq = getattr(index, "freq", None) or getattr(
index, "inferred_freq", None
)
x = array_like(x, "x", maxdim=2)
nobs = len(x)
if not np.all(np.isfinite(x)):
raise ValueError("This function does not handle missing values")
if model.startswith("m"):
if np.any(x <= 0):
raise ValueError(
"Multiplicative seasonality is not appropriate "
"for zero and negative values"
)
if period is None:
if pfreq is not None:
pfreq = freq_to_period(pfreq)
period = pfreq
else:
raise ValueError(
"You must specify a period or x must be a pandas object with "
"a PeriodIndex or a DatetimeIndex with a freq not set to None"
)
if x.shape[0] < 2 * pfreq:
raise ValueError(
f"x must have 2 complete cycles requires {2 * pfreq} "
f"observations. x only has {x.shape[0]} observation(s)"
)
if filt is None:
if period % 2 == 0: # split weights at ends
filt = np.array([0.5] + [1] * (period - 1) + [0.5]) / period
else:
filt = np.repeat(1.0 / period, period)
nsides = int(two_sided) + 1
trend = convolution_filter(x, filt, nsides)
if extrapolate_trend == "freq":
extrapolate_trend = period - 1
if extrapolate_trend > 0:
trend = _extrapolate_trend(trend, extrapolate_trend + 1)
if model.startswith("m"):
detrended = x / trend
else:
detrended = x - trend
period_averages = seasonal_mean(detrended, period)
if model.startswith("m"):
period_averages /= np.mean(period_averages, axis=0)
else:
period_averages -= np.mean(period_averages, axis=0)
seasonal = np.tile(period_averages.T, nobs // period + 1).T[:nobs]
if model.startswith("m"):
resid = x / seasonal / trend
else:
resid = detrended - seasonal
results = []
for s, name in zip(
(seasonal, trend, resid, x), ("seasonal", "trend", "resid", None), strict=True,
):
results.append(pw.wrap(s.squeeze(), columns=name))
return DecomposeResult(
seasonal=results[0],
trend=results[1],
resid=results[2],
observed=results[3],
)
[docs]
class DecomposeResult:
"""
Results class for seasonal decompositions
Parameters
----------
observed : array_like
The data series that has been decomposed.
seasonal : array_like
The seasonal component of the data series.
trend : array_like
The trend component of the data series.
resid : array_like
The residual component of the data series.
weights : array_like, optional
The weights used to reduce outlier influence.
"""
def __init__(self, observed, seasonal, trend, resid, weights=None):
self._seasonal = seasonal
self._trend = trend
if weights is None:
weights = np.ones_like(observed)
if isinstance(observed, pd.Series):
weights = pd.Series(weights, index=observed.index, name="weights")
self._weights = weights
self._resid = resid
self._observed = observed
@property
def observed(self):
"""Observed data"""
return self._observed
@property
def seasonal(self):
"""The estimated seasonal component"""
return self._seasonal
@property
def trend(self):
"""The estimated trend component"""
return self._trend
@property
def resid(self):
"""The estimated residuals"""
return self._resid
@property
def weights(self):
"""The weights used in the robust estimation"""
return self._weights
@property
def nobs(self):
"""Number of observations"""
return self._observed.shape
[docs]
def plot(
self,
observed=True,
seasonal=True,
trend=True,
resid=True,
weights=False,
):
"""
Plot estimated components
Parameters
----------
observed : bool
Include the observed series in the plot
seasonal : bool
Include the seasonal component in the plot
trend : bool
Include the trend component in the plot
resid : bool
Include the residual in the plot
weights : bool
Include the weights in the plot (if any)
Returns
-------
matplotlib.figure.Figure
The figure instance containing the plot.
"""
from pandas.plotting import register_matplotlib_converters
from statsmodels.graphics.utils import _import_mpl
plt = _import_mpl()
register_matplotlib_converters()
series = [(self._observed, "Observed")] if observed else []
series += [(self.trend, "trend")] if trend else []
if self.seasonal.ndim == 1:
series += [(self.seasonal, "seasonal")] if seasonal else []
elif self.seasonal.ndim > 1:
if isinstance(self.seasonal, pd.DataFrame):
for col in self.seasonal.columns:
series += [(self.seasonal[col], "seasonal")] if seasonal else []
else:
for i in range(self.seasonal.shape[1]):
series += [(self.seasonal[:, i], "seasonal")] if seasonal else []
series += [(self.resid, "residual")] if resid else []
series += [(self.weights, "weights")] if weights else []
if isinstance(self._observed, (pd.DataFrame, pd.Series)):
nobs = self._observed.shape[0]
xlim = self._observed.index[0], self._observed.index[nobs - 1]
else:
xlim = (0, self._observed.shape[0] - 1)
fig, axs = plt.subplots(len(series), 1, sharex=True)
for i, (ax, (series, def_name)) in enumerate(zip(axs, series, strict=True)): # noqa: B020
if def_name != "residual":
ax.plot(series)
else:
ax.plot(series, marker="o", linestyle="none")
ax.plot(xlim, (0, 0), color="#000000", zorder=-3)
name = getattr(series, "name", def_name)
if def_name != "Observed":
name = name.capitalize()
title = ax.set_title if i == 0 and observed else ax.set_ylabel
title(name)
ax.set_xlim(xlim)
fig.tight_layout()
return fig
