Hameedi's occupation-number filter is fast. Exact-in-intent filtering of a 40-component binomial MSM over 2512 observations runs in 9.2598 seconds. The literature the paper cites stops at k of 8 to 10 and calls exact inference intractable past k = 13. The word "exact" is where it comes apart, and what breaks it is the paper's own orbit-invariant benchmark, the one built specifically to rule out labelling artefacts.
One disclosure before the argument. We could not reproduce the permutation-symmetry reduction as a drop-in exact filter. Each component carries its own switching probability gamma_l, so permuting components does not preserve the filtering dynamics. We ran the conventional full-state binomial MSM filter instead. Our own performance numbers further down test that substitute. They say nothing about the paper's exactness or complexity claims.
Why anyone runs this filter
The Markov-Switching Multifractal writes a return as drift plus sigma times the square root of a product of k latent multipliers, times a standard normal. Each multiplier switches independently, with switching probabilities that decay geometrically across scales through Eq. (10), gamma_l = 1 - (1 - gamma_1)^(b^(l-1)). The paper describes the result as "a long memory methodology for computing volatility of financial assets". Its motivation for caring cites evidence that MSM beats GARCH-type models in and out of sample on crude oil volatility.
Calibration goes through a forward filter. The joint state has D^k configurations. The naive predict step costs O(D^(2k)) per timestep, and Kronecker factorization of the transition kernel gets that to O(k D^(k+1)). Still exponential in k. The log-likelihood is only available pointwise: one full forward pass per candidate parameter vector. So the filter is the bottleneck in the whole estimation.
Hameedi's move is a symmetry argument. The Gaussian emission depends on the state only through the product of the k multipliers, and multiplication commutes, so the likelihood is invariant under the permutation group S_k acting on component labels (Proposition 2). Index states by occupation vectors instead: how many components sit in each of the D multiplier values. Burnside's Lemma gives the count. In the binomial case the orbit space has k + 1 elements, and an extra factor of k from a leave-one-out recursion leaves O(k^2). In the multinomial case it is (k + D - 1 choose D - 1), giving O(k^D). The filter carries component marginals and a sector distribution. Each step it rebuilds the sector distribution from the Poisson-multinomial generating function, then recovers marginals back out through a leave-one-out formula. Code is public. Validation is 20 trials per k for k in {2,3,4,5,6}, T = 2512, five multipliers {0.7, 0.9, 1.0, 1.1, 1.3}, one arbitrarily chosen non-ergodic initial distribution, seed 42, against the naive full-state filter.
Is the filter exact?
The abstract concedes the disagreement in the same breath as the claim: "while there are significant disagreements, the ground-truth recovery of the latter seems to improve on the former." The conclusion goes further, and states that "disagreement between the two posterior representations cannot by itself be interpreted as evidence that the sector filter is performing worse". Both sentences deserve a straight answer.
Start with what supports exactness. The relative log-likelihood discrepancy stays below 5e-4 in all 100 paired comparisons, with a maximum of 4.94e-4. At k = 2 the two filters are indistinguishable: mean absolute log-likelihood gap 0.013522, mean relative gap 0.000002, mean marginal L1 distance 0.003690, MAP agreement 1.00, mean Hamming 0.00.
Now the absolute numbers, which are the ones that matter for estimation. The mean absolute log-likelihood gap goes 0.680350 at k = 3, 0.768044 at k = 4, 0.899184 at k = 5, 1.550141 at k = 6. Mean L1 between recovered marginals climbs alongside it, from 0.003690 at k = 2 to 0.273904 at k = 6, with a worst single-trial L1 of 0.513606. Worst single-trial log-likelihood gap: 3.849454 nats. 4.94e-4 looks tiny because the log-likelihoods sit near 7900. By my own arithmetic rather than the paper's: double 3.85 nats into a likelihood-ratio statistic and you get 7.70. The critical value at one degree of freedom is 3.84. So 3.85 nats of filter error can flip a nested test. An MLE loop needs the naive filter as a check.
Hameedi's answer to the labelled disagreements is fair and I accept it. MAP agreement between two estimators is the wrong diagnostic when permutation-related configurations are identified by construction. The paper shows the argmax does not commute with the pushforward along the orbit map. So he rebuilds every metric to be orbit-invariant, and reports that the sector filter recovers the true configuration better: 0.6500 versus 0.5333 at k = 3, 0.6750 versus 0.5500 at k = 4, 0.6900 versus 0.5600 at k = 5, 0.6917 versus 0.4583 at k = 6. At k = 2 the naive filter is slightly ahead, 0.6750 against 0.6500. Labelled recovery is essentially tied throughout, 0.258 against 0.250 at k = 6.
The recovery defence cuts the wrong way. Two filters computing the same posterior would report identical recovery accuracy, so a gap of 0.6917 against 0.4583 says they are computing different objects.
And the metric designed to settle it settles it. Total variation between the sector posterior and the orbit projection of the naive posterior compares two distributions on the same quotient space. No relabelling can inflate it, and the paper itself calls it orbit-invariant by construction. It runs 0.0024, 0.0745, 0.1639, 0.1688, 0.2753 across k = 2 to 6. Orbit-level MAP agreement is 0.80 at k = 2 and exactly 0.00 for k = 4, 5 and 6. Mean orbit distance d_orbit is 1.5000 (sd 0.5130) at k = 4 and 2.5500 (sd 0.6048) at k = 6. The two filters' MAP configurations sit two and a half components apart on average. At k = 6 the two posteriors put 27.5% of their mass in different sectors, measured on the quotient space where labelling cannot be the culprit.
So the disagreement is real evidence, and it is evidence about the posterior itself.
Marginals in, joint out
The paper attributes the drift to posterior multiplicity and possibly to floating-point error from dividing by p-tilde_{l,D}, and states plainly that "a full numerical stability study is not carried out here." I would look somewhere else first.
The reduction rests on two stated properties: component independence and emission symmetry. The conclusion names them as "the permutation symmetry of its likelihood and the factorization of its transition kernel". Kernel symmetry is never claimed, and Eq. (10) makes gamma_l differ across l by construction, so the kernel fails to be S_k-symmetric. The algorithm propagates labelled component marginals for exactly that reason, each with its own D-by-D matrix. The lossy step is downstream. Algorithm 1 carries marginals from one timestep to the next, and rebuilds the sector distribution from the generating function of Eq. (24). That generating function is the one for independent components. But the Bayes update against a likelihood that depends only on the product induces dependence across components in the posterior. That dependence has nowhere to live in the state being carried forward.
My hypothesis, then, not the paper's: the error should be near zero at k = 2, where a pair of marginals plus sector weights nearly pins the joint, and should grow with k. The TV series is the test, and it goes 0.0024 at k = 2, 0.0745 at k = 3, 0.1639 at k = 4, 0.1688 at k = 5, 0.2753 at k = 6.
The experiment that separates the two explanations is cheap. Re-run the k = 6 comparison in extended precision. If TV falls toward zero, it is arithmetic. If it stays near 0.2753, the recursion is an approximation with an appealing regularising side effect, which would be a publishable result about regularisation.
Speed is polynomial in k only
O(k^D) is exponential in D, and the timing tables say so. At D = 2, filtering 2512 GSPC observations takes 0.3572 s at k = 5 and 9.2598 s at k = 40. At D = 3, k = 40 costs 255.0306 s. At D = 4, 2.3768 s at k = 5 becomes 5255.0459 s at k = 40. Arithmetic on those timings, mine and not the paper's. Take a derivative-free search over (m_0, b) needing a hundred likelihood evaluations. At D = 2 that is about 15 minutes (100 x 9.26 s). At D = 4 it is about 6.1 days (100 x 5255 s). The reduction buys deep cascades. D is where the cost still explodes.
Two further things about the timing. The tables report the sector filter's wall clock and the theoretical reference T k^D. We did not find a naive-filter runtime on the same k grid, so the headline comparison is against the literature rather than against a control the author ran. And Hameedi notes the algorithm over-performed the O(k^D) prediction at low D because NumPy dispatches to compiled C, BLAS and LAPACK. Honest disclosure, and also a reason to keep the ratio column out of any scaling-exponent argument. The GSPC sample is 2512 daily bars pulled with yfinance; we did not find a date range, so that particular series is not reproducible as stated.
Our own SPY overlay
We adapted the idea and ran a conventional full-state filter. Our run measures a full-state binomial MSM volatility overlay on SPY. Component-specific gamma_l breaks the permutation reduction, and the occupation-sector filter never generated a portfolio weight.
The strategy is a volatility-targeting overlay on SPY alone, long only, daily at the close, from 2020-01-01 to 2024-07-01. Each day we fit psi_k = (m_0, sigma, b) for k in {2,3,4,5,6} on an expanding history, with gamma_k fixed at 1/756. We select k by BIC with a penalty of 3 ln(T). We propagate the full-state posterior one day and annualize by sqrt(252). Exposure is min(1.5, 0.12 / v). Forecast volatility above 35% flattens the book. Position and gross caps are both 1.5. We charged 5 bps per side plus $0.004 a share; modelled slippage was zero.
The paper reports no strategy performance of its own, no forecast RMSE and no GARCH comparison, so there is nothing of theirs to set beside ours. Ours, over 2020-01-01 to 2024-07-01: total return 54.55%, Sharpe 0.90, Sortino 1.21, Calmar 0.50, max drawdown -20.21%, realized volatility 13.52%. Those figures speak to our SPY overlay implementation only and bear on none of the paper's exactness claims.
Three of our choices would explain a weak result before the model does. The window opens on the COVID crash, where a 35% gate spends its time in cash and the rest of the sample is largely one direction. Our MLE does not impose the binomial domain 1 < m_0 < 2, so a fit with m_0 at or above 2 makes the low multiplier non-positive and breaks the Gaussian variance the filter assumes. And the paper gives no initialisation for the real-data GSPC timing run. Algorithm 1 step 1 says only "Initialize component marginals p_lj,0", with no rule attached. Our initial component probabilities were therefore our own choice, which leaves the likelihood and our BIC selection of k not uniquely pinned down. Read our figures as a statement about that implementation.
We have argued before that volatility-model rankings move with the loss function you score them under (/articles/adding-egarch-asymmetry-buys-no-forecasting-edge-over-plain-garch). The desk-relevant gap here is the same one: a filter this fast is worth having, and nothing in the paper yet tells you what it forecasts.
One extended-precision re-run at k = 6 would tell us whether the 0.2753 is arithmetic or structure. Until someone does it, treat the sector filter as a fast approximation with a curious recovery advantage, and keep the naive filter for any k where you can still afford it.
Our backtest stops at 2024-07-01, and everything after that date is deliberately left untouched so the same strategy can be checked out of sample later.