.. _replication-cmbsts:

CMBSTS — Multivariate Bayesian Structural Time Series (Menchetti & Bojinov 2022)
================================================================================

:Estimator: :doc:`../cmbsts` — :class:`mlsynth.CMBSTS`
:Source: Menchetti, Fiammetta, and Iavor Bojinov (2022), *"Estimating the
   Effectiveness of Permanent Price Reductions for Competing Products Using
   Multivariate Bayesian Structural Time Series Models,"* Annals of Applied
   Statistics 16(1): 414–435 [MenchettiBojinov2022]_.
:Reference implementation: the authors' ``CausalMBSTS`` R package
   (Bojinov & Menchetti 2020).
:Replication type: Cross-validation — the mlsynth port is checked cell-by-cell
   against the R package.
:Status: Verified — the port reproduces ``CausalMBSTS`` within Monte-Carlo
   error on the package's vignette and on the Florence supermarket study.

Validation strategy
-------------------

CMBSTS is a faithful NumPy/SciPy port of the ``CausalMBSTS`` R package: a
multivariate structural state space (trend, optional seasonal and cycle, and a
spike-and-slab regression), a Gibbs sampler with a Durbin–Koopman simulation
smoother for the latent states, and a posterior-predictive counterfactual. The
port is validated by running the genuine R package and comparing its output to
mlsynth's on identical inputs.

Seed-faithful agreement is not achievable: the R and Python samplers draw from
independent random-number streams, so individual draws differ. The honest target
is therefore agreement of the posterior summaries — the per-series temporal-
average effect and its credible bounds — within Monte-Carlo error.

The vignette (durable benchmark)
--------------------------------

The package vignette generates a bivariate weekly series,
:math:`y_{1t} = 3\sin(2t) + e`, :math:`y_{2t} = 2\cos(2t) + e`, with a fictional
:math:`+2` intervention, and fits a trend-plus-cycle model (cycle period 75).
The data ship as ``basedata/cmbsts_vignette.csv`` (generated by the R DGP with
``set.seed(1)``, so the noise is identical to the reference run). With the
estimator's default Inverse-Wishart prior
(:math:`0.01 \cdot \mathrm{diag}(\widehat{\mathrm{var}}(\mathbf{y}_{\text{pre}}))`)
and 2000 Gibbs iterations:

.. list-table::
   :header-rows: 1

   * - Quantity
     - R ``CausalMBSTS``
     - mlsynth ``CMBSTS``
   * - series 1 (treated) ATT
     - 2.4469
     - 2.457
   * - series 2 (group) ATT
     - 2.2756
     - 2.287
   * - series 1 95% CI lower
     - 0.352
     - 0.323
   * - series 1 95% CI upper
     - 4.567
     - 4.642

The posterior-mean effects agree to about :math:`0.01` and the credible bounds
to a few hundredths, and the injected :math:`+2` sits inside both bands. This is
captured durably in
`benchmarks/cases/cmbsts_vignette.py <https://github.com/jgreathouse9/mlsynth/blob/main/benchmarks/cases/cmbsts_vignette.py>`_.

The supermarket study
---------------------

The paper's empirical application is the Florence supermarket chain's permanent
discount on store-brand cookies, modelled pair-by-pair (a store brand and its
direct competitor) with a trend-plus-weekly-seasonal model and a regression block
of calendar dummies, a frozen store price, the competitor price, and ten
wine-control series. Table 3 reports the temporal-average effect per pair at the
one-month horizon, finding significant positive store-brand effects on pairs 4, 7
and 10 and no significant competitor effects. The authors' Supplement B data ship
in ``basedata/cmbsts_supermarket/``, and
`benchmarks/cases/cmbsts_supermarket.py <https://github.com/jgreathouse9/mlsynth/blob/main/benchmarks/cases/cmbsts_supermarket.py>`_
reproduces this at the one-month horizon.

Two things hold. The mlsynth store-brand effects cross-validate against the R
``CausalMBSTS`` package on identical controls and prior — pair 4 ``48.9`` vs
``47.4``, pair 7 ``80.0`` vs ``78.1``, pair 10 ``12.4`` vs ``12.3`` — within
Monte-Carlo error. And the substantive Table 3 result reproduces: large positive
store effects with no significant competitor effect. The controls here are
screened by dynamic time warping (the package's method) through the optional
``fastdtw`` package rather than the authors' ``MarketMatching``, so the control
set differs; the pair-10 effect is strictly significant in both implementations,
while pairs 4 and 7 sit on the zero boundary where Monte-Carlo noise alone moves
the lower credible bound across (the package and the port land on opposite sides).
The published strict significance on all three pairs uses the ``MarketMatching``
controls at ``niter = 2200``.

A reference-implementation difference worth noting: ``CausalMBSTS`` adds a single
observation-noise draw per posterior draw and recycles it across the forecast
horizon, so that noise does not average out of the temporal-average effect and
widens its credible interval. The mlsynth port reproduces this behaviour to match
the reference; a per-period-independent draw would give a narrower, arguably more
defensible interval.