Trend models

Module for trend models in prophetverse.

FlatTrend

Bases: TrendEffectMixin, BaseEffect

Flat trend model.

The mean of the target variable is used as the prior location for the trend.

Parameters:

Name	Type	Description	Default
changepoint_prior_scale	float	The scale of the prior distribution on the trend changepoints. Defaults to 0.1.	0.1

Source code in src/prophetverse/effects/trend/flat.py

class FlatTrend(TrendEffectMixin, BaseEffect):
    """Flat trend model.

    The mean of the target variable is used as the prior location for the trend.

    Parameters
    ----------
    changepoint_prior_scale : float, optional
        The scale of the prior distribution on the trend changepoints. Defaults to 0.1.
    """

    def __init__(self, changepoint_prior_scale: float = 0.1) -> None:
        self.changepoint_prior_scale = changepoint_prior_scale
        super().__init__()

    def _fit(self, y: pd.DataFrame, X: pd.DataFrame, scale: float = 1):
        """Initialize the effect.

        Set the prior location for the trend.

        Parameters
        ----------
        y : pd.DataFrame
            The timeseries dataframe

        X : pd.DataFrame
            The DataFrame to initialize the effect.
        """
        self.changepoint_prior_loc = y.mean().values

    def _transform(self, X: pd.DataFrame, fh: pd.Index) -> dict:
        """Prepare input data (a constant factor in this case).

        Parameters
        ----------
        idx : pd.PeriodIndex
            the timeseries time indexes

        Returns
        -------
        dict
            dictionary containing the input data for the trend model
        """
        idx = X.index
        return jnp.ones((len(idx), 1))

    def _sample_params(self, data: Any, predicted_effects: Dict[str, jnp.ndarray]):
        """Sample parameters from the prior distribution.

        Parameters
        ----------
        data : jnp.ndarray
            A constant vector with the size of the series time indexes

        Returns
        -------
        dict
            dictionary containing the sampled parameters for the trend model
        """
        return {
            "trend_flat_coefficient": numpyro.sample(
                "trend_flat_coefficient",
                dist.Gamma(
                    rate=self.changepoint_prior_loc / self.changepoint_prior_scale**2,
                    concentration=self.changepoint_prior_loc,
                ),
            ),
        }

    def _predict(  # type: ignore[override]
        self,
        data: jnp.ndarray,
        predicted_effects: dict,
        params: dict,
    ) -> jnp.ndarray:
        """Apply the trend.

        Parameters
        ----------
        constant_vector : jnp.ndarray
            A constant vector with the size of the series time indexes

        Returns
        -------
        jnp.ndarray
            The forecasted trend
        """
        # Alias for clarity
        constant_vector = data

        coefficient = params["trend_flat_coefficient"]

        return constant_vector * coefficient

PiecewiseLinearTrend

Bases: TrendEffectMixin, BaseEffect

Piecewise Linear Trend model.

This model assumes that the trend is piecewise linear, with changepoints at regular intervals. The number of changepoints is determined by the changepoint_interval and changepoint_range parameters. The changepoint_interval parameter specifies the interval between changepoints, while the changepoint_range parameter specifies the range of the changepoints.

This implementation is based on the Prophet_ library. The initial values (global rate and global offset) are suggested using the maximum and minimum values of the time series data.

Parameters:

Name	Type	Description	Default
changepoint_interval	int	The interval between changepoints.	required
changepoint_range	int	The range of the changepoints.	required
changepoint_prior_scale	Distribution	The prior scale for the changepoints.	required
offset_prior_scale	float	The prior scale for the offset. Default is 0.1.	0.1
squeeze_if_single_series	bool	If True, squeeze the output if there is only one series. Default is True.	True
remove_seasonality_before_suggesting_initial_vals	bool	If True, remove seasonality before suggesting initial values, using sktime's detrender. Default is True.	True

Source code in src/prophetverse/effects/trend/piecewise.py

class PiecewiseLinearTrend(TrendEffectMixin, BaseEffect):
    """Piecewise Linear Trend model.

    This model assumes that the trend is piecewise linear, with changepoints
    at regular intervals. The number of changepoints is determined by the
    `changepoint_interval` and `changepoint_range` parameters. The
    `changepoint_interval` parameter specifies the interval between changepoints,
    while the `changepoint_range` parameter specifies the range of the changepoints.

    This implementation is based on the `Prophet`_ library. The initial values (global
    rate and global offset) are suggested using the maximum and minimum values of the
    time series data.


    Parameters
    ----------
    changepoint_interval : int
        The interval between changepoints.
    changepoint_range : int
        The range of the changepoints.
    changepoint_prior_scale : dist.Distribution
        The prior scale for the changepoints.
    offset_prior_scale : float, optional
        The prior scale for the offset. Default is 0.1.
    squeeze_if_single_series : bool, optional
        If True, squeeze the output if there is only one series. Default is True.
    remove_seasonality_before_suggesting_initial_vals : bool, optional
        If True, remove seasonality before suggesting initial values, using sktime's
        detrender. Default is True.


    """

    def __init__(
        self,
        changepoint_interval: int,
        changepoint_range: int,
        changepoint_prior_scale: dist.Distribution,
        offset_prior_scale=0.1,
        squeeze_if_single_series: bool = True,
        remove_seasonality_before_suggesting_initial_vals: bool = True,
        **kwargs,
    ):
        self.changepoint_interval = changepoint_interval
        self.changepoint_range = changepoint_range
        self.changepoint_prior_scale = changepoint_prior_scale
        self.offset_prior_scale = offset_prior_scale
        self.squeeze_if_single_series = squeeze_if_single_series
        self.remove_seasonality_before_suggesting_initial_vals = (
            remove_seasonality_before_suggesting_initial_vals
        )
        super().__init__()

    def _fit(self, y: pd.DataFrame, X: pd.DataFrame, scale: float = 1):
        """Initialize the effect.

        Set the prior location for the trend.

        Parameters
        ----------
        y : pd.DataFrame
            The timeseries dataframe

        X : pd.DataFrame
            The DataFrame to initialize the effect.

        scale : float, optional
            The scale of the timeseries. For multivariate timeseries, this is
            a dataframe. For univariate, it is a simple float.
        """
        super()._fit(X=X, y=y, scale=scale)

        t_scaled = self._index_to_scaled_timearray(
            y.index.get_level_values(-1).unique()
        )
        self._setup_changepoints(t_scaled)
        self._setup_changepoint_prior_vectors(y)
        self._index_names = y.index.names
        self._series_idx = None
        if y.index.nlevels > 1:
            self._series_idx = y.index.droplevel(-1).unique()

    def _fh_to_index(self, fh: pd.Index) -> Union[pd.Index, pd.MultiIndex]:
        """Convert an index representing the fcst horizon to multiindex if needed.

        If there's a single timeseries, just returns the fh.

        Parameters
        ----------
        fh : pd.Index
            The timeindex representing the forecasting horizon.

        Returns
        -------
        Union[pd.Index, pd.MultiIndex]
            The fh for all time series passed during fit
        """
        if self._series_idx is None:
            return fh

        idx_list = self._series_idx.to_list()
        idx_list = [x if isinstance(x, tuple) else (x,) for x in idx_list]
        # Create a new multi-index combining the existing levels with the new time index
        new_idx_tuples = list(
            map(
                lambda x: (
                    *x[0],
                    x[1],
                ),
                # Create a cross product of current indexes
                # and dates in fh
                itertools.product(idx_list, fh.to_list()),
            )
        )
        return pd.MultiIndex.from_tuples(new_idx_tuples, names=self._index_names)

    def _transform(self, X: pd.DataFrame, fh: pd.Index) -> dict:
        """
        Prepare the input data for the piecewise trend model.

        Parameters
        ----------
        X: pd.DataFrame
            The exogenous variables DataFrame.
        fh: pd.Index
            The forecasting horizon as a pandas Index.

        Returns
        -------
        jnp.ndarray
            An array containing the prepared input data.
        """
        idx = self._fh_to_index(fh)
        return self.get_changepoint_matrix(idx)

    def _sample_params(self, data: Any, predicted_effects: Dict[str, jnp.ndarray]):

        changepoint_matrix = data

        offset = numpyro.sample(
            "offset",
            dist.Normal(self._offset_prior_loc, self._offset_prior_scale),
        )
        changepoint_coefficients = numpyro.sample(
            "changepoint_coefficients",
            dist.Laplace(self._changepoint_prior_loc, self._changepoint_prior_scale),
        )

        if changepoint_matrix.ndim == 3:
            changepoint_coefficients = changepoint_coefficients.reshape((1, -1, 1))
            offset = offset.reshape((-1, 1, 1))

        return {
            "changepoint_coefficients": changepoint_coefficients,
            "offset": offset,
        }

    def _predict(
        self,
        data: jnp.ndarray,
        predicted_effects: Dict[str, jnp.ndarray],
        params: dict,
    ) -> jnp.ndarray:
        """
        Compute the trend based on the given changepoint matrix.

        Parameters
        ----------
        data: jnp.ndarray
            The changepoint matrix.
        predicted_effects: Dict[str, jnp.ndarray]
            Dictionary of previously computed effects. For the trend, it is an empty
            dict.

        Returns
        -------
        jnp.ndarray
            The computed trend.
        """
        # alias for clarity
        changepoint_matrix = data
        changepoint_coefficients = params["changepoint_coefficients"]
        offset = params["offset"]

        trend = (changepoint_matrix) @ changepoint_coefficients + offset

        if trend.ndim == 1 or (
            trend.ndim == 3 and self.n_series == 1 and self.squeeze_if_single_series
        ):
            trend = trend.reshape((-1, 1))

        return trend

    def get_changepoint_matrix(self, idx: pd.PeriodIndex) -> jnp.ndarray:
        """
        Return the changepoint matrix for the given index.

        Parameters
        ----------
        idx: pd.PeriodIndex
            The index for which to compute the changepoint matrix.

        Returns
        -------
            jnp.ndarray: The changepoint matrix.
        """
        t_scaled = self._index_to_scaled_timearray(idx)
        changepoint_matrix = self._get_multivariate_changepoint_matrix(t_scaled)

        # If only one series, remove the first dimension
        if changepoint_matrix.shape[0] == 1:
            if self.squeeze_if_single_series:
                changepoint_matrix = changepoint_matrix[0]

        return changepoint_matrix

    @property
    def n_changepoint_per_series(self):
        """Get the number of changepoints per series.

        Returns
        -------
        int
            Number of changepoints per series.
        """
        return [len(cp) for cp in self._changepoint_ts]

    @property
    def n_changepoints(self):
        """Get the total number of changepoints.

        Returns
        -------
        int
            Total number of changepoints.
        """
        return sum(self.n_changepoint_per_series)

    def _get_multivariate_changepoint_matrix(self, t_scaled) -> jnp.ndarray:
        """
        Get the changepoint matrix.

        The changepoint matrix has shape (n_series, n_timepoints, total number of
        changepoints for all series). A mask is applied so that for index i at dim 0,
        only the changepoints for series i are non-zero at dim -1.

        Parameters
        ----------
        t_scaled: jnp.ndarray
            Transformed time index.

        Returns
        -------
        jnp.ndarray
            The changepoint matrix.
        """
        changepoint_ts = np.concatenate(self._changepoint_ts)
        changepoint_design_tensor_list = []
        changepoint_mask_tensor_list = []
        for i, n_changepoints in enumerate(self.n_changepoint_per_series):
            A = _get_changepoint_matrix(t_scaled, changepoint_ts)

            start_idx = sum(self.n_changepoint_per_series[:i])
            end_idx = start_idx + n_changepoints
            mask = np.zeros_like(A)
            mask[:, start_idx:end_idx] = 1

            changepoint_design_tensor_list.append(A)
            changepoint_mask_tensor_list.append(mask)

        changepoint_design_tensor: np.ndarray = np.stack(
            changepoint_design_tensor_list, axis=0
        )
        changepoint_mask_tensor: np.ndarray = np.stack(
            changepoint_mask_tensor_list, axis=0
        )
        return changepoint_design_tensor * changepoint_mask_tensor

    def _setup_changepoints(self, t_scaled) -> None:
        """
        Set changepoint variables.

        This function has the collateral effect of setting the following attributes:
        - self._changepoint_ts.

        Parameters
        ----------
        t_scaled: jnp.ndarray
            Transformed time index.

        Returns
        -------
            None
        """
        changepoint_intervals = _to_list_if_scalar(
            self.changepoint_interval, self.n_series
        )
        changepoint_ranges = _to_list_if_scalar(self.changepoint_range, self.n_series)

        changepoint_ts = []
        for changepoint_interval, changepoint_range in zip(
            changepoint_intervals, changepoint_ranges
        ):
            changepoint_ts.append(
                _get_changepoint_timeindexes(
                    t_scaled,
                    changepoint_interval=changepoint_interval,
                    changepoint_range=changepoint_range,
                )
            )

            if len(changepoint_ts[-1]) == 0:
                raise ValueError(
                    "No changepoints were generated. Try increasing the changing"
                    + f" the changepoint_range. There are {len(t_scaled)} timepoints "
                    + f" in the series, changepoint_range is {changepoint_range} and "
                    + f"changepoint_interval is {changepoint_interval}."
                )

        self._changepoint_ts = changepoint_ts

    def _setup_changepoint_prior_vectors(self, y: pd.DataFrame) -> None:
        """
        Set up the changepoint prior vectors for the model.

        Parameters
        ----------
        y: pd.DataFrame
            The input DataFrame containing the time series data.

        Returns
        -------
            None
        """
        if self.remove_seasonality_before_suggesting_initial_vals:
            detrender = Detrender()
            y = y - detrender.fit_transform(y)

        self._global_rates, self._offset_prior_loc = (
            self._suggest_global_trend_and_offset(y)
        )
        self._changepoint_prior_loc, self._changepoint_prior_scale = (
            self._get_changepoint_prior_vectors(global_rates=self._global_rates)
        )
        self._offset_prior_scale = self.offset_prior_scale

    def _get_changepoint_prior_vectors(
        self,
        global_rates: jnp.array,
    ) -> Tuple[jnp.array, jnp.array]:
        """
        Return the prior vectors for the changepoint coefficients.

        Parameters
        ----------
        global_rates: jnp.array
            The global rates for each series.

        Returns
        -------
        Tuple[jnp.array, jnp.array]
            A tuple containing the changepoint prior location vector and the
            changepoint prior scale vector.
        """
        n_series = len(self.n_changepoint_per_series)

        def zeros_with_first_value(size, first_value):
            x = jnp.zeros(size)
            x.at[0].set(first_value)
            return x

        changepoint_prior_scale_vector = np.concatenate(
            [
                np.ones(n_changepoint) * cur_changepoint_prior_scale
                for n_changepoint, cur_changepoint_prior_scale in zip(
                    self.n_changepoint_per_series,
                    _to_list_if_scalar(self.changepoint_prior_scale, n_series),
                )
            ]
        )

        changepoint_prior_loc_vector = np.concatenate(
            [
                zeros_with_first_value(n_changepoint, estimated_global_rate)
                for n_changepoint, estimated_global_rate in zip(
                    self.n_changepoint_per_series, global_rates
                )
            ]
        )

        return jnp.array(changepoint_prior_loc_vector), jnp.array(
            changepoint_prior_scale_vector
        )

    def _suggest_global_trend_and_offset(
        self, y: pd.DataFrame
    ) -> Tuple[jnp.array, jnp.array]:
        """
        Suggest the global trend and offset for the given time series data.

        Parameters
        ----------
        y: pd.DataFrame
            The time series data.

        Returns
        -------
        Tuple[jnp.array, jnp.array]
            A tuple containing the global trend and offset.
        """
        t = self._index_to_scaled_timearray(y.index.get_level_values(-1).unique())

        y_array = series_to_tensor(y)

        global_rate = _enforce_array_if_zero_dim(
            (y_array[:, -1].squeeze() - y_array[:, 0].squeeze())
            / (t[0].squeeze() - t[-1].squeeze())
        )
        offset_loc = y_array[:, 0].squeeze() - global_rate * t[0].squeeze()

        return global_rate, offset_loc

n_changepoint_per_series property

Get the number of changepoints per series.

Returns:

Type	Description
int	Number of changepoints per series.

n_changepoints property

Get the total number of changepoints.

Returns:

Type	Description
int	Total number of changepoints.

get_changepoint_matrix(idx)

Return the changepoint matrix for the given index.

Parameters:

Name	Type	Description	Default
idx	PeriodIndex	The index for which to compute the changepoint matrix.	required

Returns:

Type	Description
jnp.ndarray: The changepoint matrix.

Source code in src/prophetverse/effects/trend/piecewise.py

def get_changepoint_matrix(self, idx: pd.PeriodIndex) -> jnp.ndarray:
    """
    Return the changepoint matrix for the given index.

    Parameters
    ----------
    idx: pd.PeriodIndex
        The index for which to compute the changepoint matrix.

    Returns
    -------
        jnp.ndarray: The changepoint matrix.
    """
    t_scaled = self._index_to_scaled_timearray(idx)
    changepoint_matrix = self._get_multivariate_changepoint_matrix(t_scaled)

    # If only one series, remove the first dimension
    if changepoint_matrix.shape[0] == 1:
        if self.squeeze_if_single_series:
            changepoint_matrix = changepoint_matrix[0]

    return changepoint_matrix

PiecewiseLogisticTrend

Bases: PiecewiseLinearTrend

Piecewise logistic trend model.

This logistic trend differs from the original Prophet logistic trend in that it considers a capacity prior distribution. The capacity prior distribution is used to estimate the maximum value that the time series trend can reach.

It uses internally the piecewise linear trend model, and then applies a logistic function to the output of the linear trend model.

The initial values (global rate and global offset) are suggested using the maximum and minimum values of the time series data.

Parameters:

Name	Type	Description	Default
changepoint_interval	int	The interval between changepoints.	required
changepoint_range	int	The range of the changepoints.	required
changepoint_prior_scale	Distribution	The prior scale for the changepoints.	required
offset_prior_scale	float	The prior scale for the offset. Default is 0.1.	10
squeeze_if_single_series	bool	If True, squeeze the output if there is only one series. Default is True.	required
remove_seasonality_before_suggesting_initial_vals	bool	If True, remove seasonality before suggesting initial values, using sktime's detrender. Default is True.	required
capacity_prior	Distribution	The prior distribution for the capacity. Default is a HalfNormal distribution with loc=1.05 and scale=1.	None

Source code in src/prophetverse/effects/trend/piecewise.py

class PiecewiseLogisticTrend(PiecewiseLinearTrend):
    """
    Piecewise logistic trend model.

    This logistic trend differs from the original Prophet logistic trend in that it
    considers a capacity prior distribution. The capacity prior distribution is used
    to estimate the maximum value that the time series trend can reach.

    It uses internally the piecewise linear trend model, and then applies a logistic
    function to the output of the linear trend model.


    The initial values (global rate and global offset) are suggested using the maximum
    and minimum values of the time series data.


    Parameters
    ----------
    changepoint_interval : int
        The interval between changepoints.
    changepoint_range : int
        The range of the changepoints.
    changepoint_prior_scale : dist.Distribution
        The prior scale for the changepoints.
    offset_prior_scale : float, optional
        The prior scale for the offset. Default is 0.1.
    squeeze_if_single_series : bool, optional
        If True, squeeze the output if there is only one series. Default is True.
    remove_seasonality_before_suggesting_initial_vals : bool, optional
        If True, remove seasonality before suggesting initial values, using sktime's
        detrender. Default is True.
    capacity_prior : dist.Distribution, optional
        The prior distribution for the capacity. Default is a HalfNormal distribution
        with loc=1.05 and scale=1.
    """

    def __init__(
        self,
        changepoint_interval: int,
        changepoint_range: int,
        changepoint_prior_scale: float,
        offset_prior_scale=10,
        capacity_prior: dist.Distribution = None,
        **kwargs,
    ):
        if capacity_prior is None:
            capacity_prior = dist.TransformedDistribution(
                dist.HalfNormal(
                    0.2,
                ),
                dist.transforms.AffineTransform(loc=1.05, scale=1),
            )

        self.capacity_prior = capacity_prior

        super().__init__(
            changepoint_interval,
            changepoint_range,
            changepoint_prior_scale,
            offset_prior_scale=offset_prior_scale,
            **kwargs,
        )

    def _suggest_global_trend_and_offset(
        self, y: pd.DataFrame
    ) -> Tuple[jnp.array, jnp.array]:
        """
        Suggest the global trend and offset for the given time series data.

        This implementation considers the max and min values of the timeseries
        to anallytically estimate the global rate and offset.

        Parameters
        ----------
        y: pd.DataFrame
            The input time series data.

        Returns
        -------
        Tuple[jnp.array, jnp.array]
            A tuple containing the suggested global rates
            and offset.
        """
        t_arrays = self._index_to_scaled_timearray(
            y.index.get_level_values(-1).unique()
        )
        y_arrays = series_to_tensor(y)

        if hasattr(self.capacity_prior, "loc"):
            capacity_prior_loc = self.capacity_prior.loc
        else:
            capacity_prior_loc = y_arrays.max() * 1.05

        global_rates, offset = _suggest_logistic_rate_and_offset(
            t=t_arrays.squeeze(),
            y=y_arrays.squeeze(),
            capacities=capacity_prior_loc,
        )

        return global_rates, offset

    def _sample_params(self, data, predicted_effects):
        """
        Sample params for the effect.

        Use super to sample the changepoint coefficients and offset, and then sample
        the capacity using the capacity prior.

        Parameters
        ----------
        data : Any
            The input data.
        predicted_effects : Dict[str, jnp.ndarray]
            The predicted effects

        Returns
        -------
        dict
            The sampled parameters.
        """
        with numpyro.plate("series", self.n_series, dim=-3):
            capacity = numpyro.sample("capacity", self.capacity_prior)

        return {
            "capacity": capacity,
            **super()._sample_params(data=data, predicted_effects=predicted_effects),
        }

    def _predict(  # type: ignore[override]
        self,
        data: Any,
        predicted_effects: Dict[str, jnp.ndarray],
        params: Dict[str, jnp.ndarray],
    ) -> jnp.ndarray:
        """
        Compute the trend for the given changepoint matrix.

        Parameters
        ----------
        changepoint_matrix: jnp.ndarray
            The changepoint matrix.

        Returns
        -------
        jnp.ndarray
            The computed trend.
        """
        trend = super()._predict(
            data=data, predicted_effects=predicted_effects, params=params
        )

        capacity = params["capacity"]
        if self.n_series == 1:
            capacity = capacity.squeeze()

        trend = capacity / (1 + jnp.exp(-trend))
        return trend

TrendEffectMixin

Mixin class for trend models.

Trend models are effects applied to the trend component of a time series.

Attributes:

Name	Type	Description
t_scale	float	The time scale of the trend model.
t_start	float	The starting time of the trend model.
n_series	int	The number of series in the time series data.

Source code in src/prophetverse/effects/trend/base.py

class TrendEffectMixin:
    """
    Mixin class for trend models.

    Trend models are effects applied to the trend component of a time series.

    Attributes
    ----------
    t_scale: float
        The time scale of the trend model.
    t_start: float
        The starting time of the trend model.
    n_series: int
        The number of series in the time series data.
    """

    _tags = {"skip_predict_if_no_match": False, "supports_multivariate": True}

    def _fit(self, y: pd.DataFrame, X: pd.DataFrame, scale: float = 1) -> None:
        """Initialize the effect.

        Set the time scale, starting time, and number of series attributes.

        Parameters
        ----------
        y : pd.DataFrame
            The timeseries dataframe

        X : pd.DataFrame
            The DataFrame to initialize the effect.
        """
        # Set time scale
        t_days = convert_index_to_days_since_epoch(
            y.index.get_level_values(-1).unique()
        )
        self.t_scale = (t_days[1:] - t_days[:-1]).mean()
        self.t_start = t_days.min() / self.t_scale
        if y.index.nlevels > 1:
            self.n_series = y.index.droplevel(-1).nunique()
        else:
            self.n_series = 1

    def _index_to_scaled_timearray(self, idx):
        """
        Convert the index to a scaled time array.

        Parameters
        ----------
        idx: int
            The index to be converted.

        Returns
        -------
        float
            The scaled time array value.
        """
        if idx.nlevels > 1:
            idx = idx.get_level_values(-1).unique()

        t_days = convert_index_to_days_since_epoch(idx)
        return (t_days) / self.t_scale - self.t_start