chc.dlm

Stability tier: experimental.

Discount dynamic linear models: forward filtering with a learned observational variance, smoothing, backward sampling and the decomposition it gives, sampling with coefficients held to a sign, forecasting and monitoring. Experimental.

The model is West and Harrison's (1997): y_t = F_t' theta_t + nu_t with theta_t = G theta_(t-1) + omega_t, where the evolution variance is not a parameter but a discount of the information the last step carried. DynamicLinearModel is built from blocks, each with its own discount delta in (0, 1]:

• Polynomial: a level, or a level with a slope, and higher orders;
• Seasonal: Fourier harmonics of a period, which need not be an integer;
• Regression: coefficients on columns of x, which follow random walks.

With P_t = G C_(t-1) G' the evolved posterior covariance, the prior covariance is

• form="additive" (the default, West and Harrison's component discounting): R_t = P_t + blockdiag_i(P_t,ii (1 - delta_i) / delta_i), which inflates each block's own covariance and leaves the covariance between blocks as it was;
• form="multiplicative": R_t = D P_t D, D = diag(delta_i^(-1/2)), which also inflates the covariance between blocks.

With one discount for every block the multiplicative form is the single discount, R_t = P_t / delta; the additive form is not, since it still leaves the covariance between blocks uninflated. No third form is offered.

The observational variance is learned (Normal-Gamma): Prior.dof is n_0 and Prior.scale is S_0, the prior point estimate of the variance; an infinite dof means the variance is known and equal to scale. variance_discount (beta) lets it move: n_t = beta n_(t-1) + 1. The one-step forecast is Student t with beta n_(t-1) degrees of freedom, location f_t and squared scale Q_t, and DLMFit.log_likelihood sums its log density over the observed steps: that is the criterion for choosing discounts. The update is the Joseph form, so a covariance stays symmetric and positive semidefinite in floating point; a step whose Q_t is not positive and finite raises rather than being zeroed.

Missing observations are NaN in y. A missing step evolves and does not update: m_t = a_t, C_t = R_t, and the variance's degrees of freedom decay by beta.

What is exact and what is not:

• The filter and the forecasts are exact for the model. So is backward_sample (forward filtering, backward sampling) when beta = 1; below 1 it draws the variance path by the beta-gamma evolution run backwards and each state given the next at that step's variance (McAlinn and West 2019, appendix A.2).
• constrained_sample draws the posterior truncated to coefficients' signs by Gibbs sampling, which is exact only in the limit of its chains; its R-hat and effective sample size say whether a run can be read.
• smooth returns the exact means, covariances and lag-one covariances of those draws. Each covariance is in the units of E[V_t | D_T]: with beta = 1 that is the final estimate n_T S_T / (n_T - 2) at every step, not the S_t the filter had at t (a smoother that keeps each C_t in its own units misstates it by S_T / S_t); below 1 it is one quadrature a step. Only the marginals' shape is approximate below 1, where they are not Student t; SmoothedStates.exact says which.
• forecast holds the evolution variance of the first step ahead for every later one by default (evolution="constant"); "compounding" applies the discount at every step, which inflates the k-step state variance by delta^(-k). The two differ by a factor 3.9 in the state variance of a level 52 steps ahead at delta = 0.95.

What it does not do, and what to watch:

• Discounting is forgetting. With G = I and a single discount the posterior mean is weighted least squares with weights delta^(T - t) (Ameen and Harrison 1984), so the effective sample is (1 + delta) / (1 - delta) observations.
• A regressor that is zero carries no information, and its coefficient's variance still grows by 1 / delta a step. Regression(hold_when_idle=True) discounts only the coordinates whose regressor is non-zero at the step, a coordinate-wise form of restricted forgetting (Kulhavy 1987); a long idle run without it is logged as a warning (chc_event="dlm_windup").
• A coefficient that moves is not identified from a saturation that bends. On observational marketing data a time-varying coefficient and a nonlinear response can explain the same series (Dew et al. 2024, arXiv:2408.07678); this module fits the first and cannot tell them apart.
• A discount picked by one-step likelihood can be too high for the coefficients' intervals. On the random-walk world of scripts/bench_dlm.py, 500 series, the likelihood picked the coefficients' discount at 0.95 or above in 450 and below 0.9 in none, and the 90% intervals of a channel's contribution over 13 weeks covered 0.78 and 0.77; fixed at 0.85 or 0.9 they covered 0.88-0.90. An interval reported from the likelihood's pick alone is too narrow; the union over confidence_set's fits covered 0.94 there.
• A monitor has until the filter absorbs a change. Once the discount has taken a shift into the level, the one-step errors are white again. On scripts/bench_dlm_monitor.py's local level at delta = 0.9, 300 series, chc.gate.DriftAlarm at an average run length of 1000 caught a shift of two one-step deviations within 1000 steps in 74% of them, and one of three in all but one. West and Harrison's monitor at tau = 0.135 caught each within a median of one step, and on unchanged series alarmed within a median of 13; DriftAlarm alarmed by step 100 on 2.7% of them, under its bound of 10%. Those series are the model's own, so the bound holds there by construction.
• A known variance set too low reads as a change. With every deviation 1.25 times what the known variance allowed for, the alarm sounded by step 100 on 57% of the series; with the variance learned from a prior of five degrees of freedom, on 4.3%.

monitor_evalues turns the one-step errors into e-values for West and Harrison's (1986) alternatives -- a shift of the location, an inflation of the scale -- which chc.gate.DriftAlarm turns into an alarm with a guaranteed average run length. forward_filter(interventions=...) applies an extra discount at named steps, their feed-back intervention, and logs each (chc_event="dlm_intervention").

Over geos. GeoDLM stacks national blocks every geo reads and regional blocks each geo has its own copy of; forward_filter_geos filters a KPI a geo, y_(g,t) = F_(g,t)' theta_t + nu_(g,t) with nu_(g,t) ~ N(0, V w_g), w_g the geo's known variance relative to the others'. A geo's row of F reads the national blocks and its own regional ones, and a national and a regional regression block read the geo's columns of x from the first, so a geo's coefficient on a column both read is the national one plus its own deviation. The evolution and the discounts are this module's for the stacked model; the update is West and Harrison's for a vector with one variance scale (1997, section 16.4): with r_t geos observed, Q_t = F R_t F' + S_(t-1) W, A_t = R_t F' Q_t^-1, n_t = beta n_(t-1) + r_t and S_t = S_(t-1) (beta n_(t-1) + e' Q_t^-1 e) / n_t, and the one-step forecast of the observed geos is multivariate Student t. That is the update one geo at a time with no evolution between them, and its density the product of theirs (validation/geo_dlm.mac, steps 1 and 2).

• One geo is forward_filter's model of the national blocks and then the regional ones, and with every discount 1 the filter is the conjugate regression on the first state over every observation at once; the tests hold the first to 1e-12 and the second to 1e-10.
• The filter is dense: with p coordinates, the national ones and every geo's, a step costs O(p^3) and the fit holds 2 T covariances of p^2 entries. An observation of one geo adds nothing to the posterior precision between two others, and the multiplicative form, a diagonal scaling of the covariance, keeps that zero through the evolution while component discounting does not (step 3): a filter linear in the number of geos exists for the first form and is not built.
• smooth and backward_sample take a filter over geos and run back over its stacked state. They read the filter's moments, the evolution and the variance discount, and the variance's backward step reads n_t and S_t, not how many geos a step observed. Over three geos with missing observations the smoother is the joint Gaussian posterior of every state to 1e-8, the variance known or learned. There is no forecast or decomposition over geos yet. One variance scale serves every geo, each geo's share of it fixed.
• fit_geo_spread chooses the geos' spread, the regional prior's variance on named coordinates, by the marginal likelihood that the filter's log-likelihood is (type-II maximum likelihood), with each spread's likelihood-ratio interval. It draws the spread from its posterior under a prior flat on each standard deviation by importance sampling, so that GeoSpread.mixture carries a quantity's uncertainty about the spread as well as about the states. Where one spread is weakly identified the mixture is quadrature's within its Monte Carlo error, with the best inside the range, near its floor, and at its ceiling.

The implementation is NumPy in float64 and written from the published equations.

Polynomial dataclass

A polynomial trend: order=1 is a level, order=2 a level and a slope.

F = (1, 0, ..., 0) and G is the upper bidiagonal matrix of ones.

Raises:

Type	Description
ValueError	on an order below 1 or a discount outside (0, 1].

Seasonal dataclass

Fourier harmonics j of a period: a rotation by 2 pi j / period a step for each.

A harmonic takes two state coordinates, except j = period / 2, which takes one and flips sign every step. period need not be an integer (52.18 weeks to a year).

Raises:

Type	Description
ValueError	on a period not above 1, harmonics that are not distinct integers in [1, period / 2], or a discount outside (0, 1].

Regression dataclass

Coefficients on width columns of x, each a random walk (G = I).

The columns are taken from x in block order: the first regression block reads the first width columns, the next the columns after them.

hold_when_idle=True discounts, at each step, only the coordinates whose regressor is non-zero, so a coefficient's variance does not grow while its column is zero.

Raises:

Type	Description
ValueError	on a width below 1 or a discount outside (0, 1].

Prior dataclass

The state and the observational variance before the first observation.

theta_0 ~ T_dof(mean, covariance) and 1 / V ~ Gamma(dof / 2, dof * scale / 2): scale is the prior point estimate of the variance and covariance is stated at it. An infinite dof means the variance is known, equal to scale, and the state's prior is normal.

covariance must be positive definite. Every block's G is invertible, so every later prior and posterior covariance is positive definite too, which the smoother relies on.

Raises:

Type	Description
ValueError	on shapes that do not agree, a non-finite entry, a covariance that is not symmetric positive definite, or a scale or dof that is not positive.

DynamicLinearModel dataclass

Blocks, their prior, how their discounts combine, and the variance discount.

The state is the blocks' coordinates in order. See the module docstring for form.

Raises:

Type	Description
ValueError	on no blocks, a prior whose size is not the blocks' total, an unknown form, a variance_discount outside (0, 1], or a variance_discount below 1 with a known variance.

size property

The state's dimension.

regressors property

How many columns x must have.

DLMFit dataclass

The forward filter's output, one row per step, in West and Harrison's notation.

At step t: the prior theta_t | D_(t-1) has mean a and covariance R; the one-step forecast is Student t with one_step_dof degrees of freedom, location one_step_mean (f) and squared scale one_step_scale (Q), so its variance is Q dof / (dof - 2); the posterior has mean m and covariance C, the variance estimate S with n degrees of freedom. log_scores is the one-step forecast's log density at the observation, NaN on a missing step. C is stated at S_t, R at S_(t-1).

log_likelihood property

The sum of the one-step forecasts' log densities over the observed steps.

errors property

One-step errors y_t - f_t, NaN on a missing step.

SmoothedStates dataclass

theta_t | D_T for every step: means, covariances, and Cov(theta_t, theta_(t+1)).

A covariance is the variance, not the scale. With a variance discount of 1 (exact) the marginals are Student t with dof degrees of freedom, normal when infinite. Below 1 the variance moves, the marginals are scale mixtures of normals that are not Student t, and dof is that of a gamma matched to the smoothed precision's mean and variance, for an approximate t; the means and covariances stay exact.

PosteriorDraws dataclass

Joint draws of the states and the observational variance given all the data.

ConstrainedDraws dataclass

Draws of theta_(1:T) and V given D_T and the signs, from Gibbs chains after their warm-up. rhat is the largest rank-normalised split R-hat, bulk or tail, over every state at every step and the variance (Vehtari et al. 2021); ess the smallest bulk effective sample size over the same.

mixed property

Whether the draws can be read: rhat at most 1.01 and ess at least 100 a chain (Vehtari et al. 2021).

pooled property

The chains pooled into one set of draws, as decompose reads them.

DLMForecast dataclass

k-step-ahead marginal forecasts from the end of the data, k = 1, ..., horizon.

mean and scale are the Student t location and squared scale of y_(T+k), with dof degrees of freedom; state_mean and state_scale those of theta_(T+k).

GeoDLM dataclass

A discount DLM over geos geos: national blocks every geo reads and regional blocks each geo has its own copy of, filtered by forward_filter_geos.

prior is on the stacked state: the national blocks' coordinates, then each geo's copy of the regional blocks', in the geos' order; stacked_prior builds one whose geos' regional priors are alike and independent. relative_variance is each geo's observational variance relative to the others', 1 for every geo when omitted.

A national and a regional regression block read each geo's columns of x from the first, so where both read a column, a geo's coefficient on it is the national coefficient plus the geo's own deviation. The regional blocks' prior variance is then the spread of the geos about the national coefficient, and the regional discount how long that prior is remembered.

Raises:

Type	Description
ValueError	on a geos that is not a positive integer, a relative_variance that is not one positive finite number a geo, and what DynamicLinearModel refuses of the stacked model: no blocks, a prior of the wrong size, an unknown form, or a variance_discount it does not take.

stacked property

The stacked state's model: the national blocks, then the regional ones once a geo.

regressors property

How many columns each geo's x must have: the national regression blocks' widths summed, or the regional ones', whichever is more.

GeoDLMFit dataclass

forward_filter_geos's output over the stacked state, one row per step, in DLMFit's notation.

one_step_mean and one_step_scale are every geo's, observed or not: location F a and squared scale F R F' + S_(t-1) W, W the relative variances on the diagonal. log_scores is the observed geos' joint one-step log density, NaN on a step with none observed. n rises by the number of geos observed at a step.

log_likelihood property

The sum of the one-step forecasts' joint log densities over the steps with a geo observed.

errors property

One-step errors y - f, (T, geos), NaN where a geo is missing.

GeoSpread dataclass

How far the geos spread about the national coefficients, chosen by the marginal likelihood: fit_geo_spread's output.

variance[i] is the prior variance of every geo's regional coordinate pooled[i], stated at the prior's scale as Prior states a covariance: at a variance V the geos spread about the national coefficient with variance variance[i] V / scale. fit is the filter at variance, and lower and upper each spread's likelihood-ratio interval, the other spreads held at their best.

draws are spreads drawn from their posterior, (draws, pooled), under a prior flat on each standard deviation over the range searched, with self-normalised importance weights. mixture averages a quantity's posterior over them, so that it carries the spread's uncertainty as well as the states'.

at(variance)

The model with every geo's pooled coordinates at prior variance variance.

mixture(quantity)

A quantity's mean and variance over the spread's posterior.

quantity maps a filter to the quantity's posterior mean and variance at that filter's spread; the result is their mixture over the draws, sum_j w_j m_j and sum_j w_j (v_j + (m_j - mean)^2), one filter a draw of positive weight.

forward_filter(model, y, x=None, *, interventions=None)

Filter y through model, with x the regression blocks' columns.

interventions maps a step to an extra discount applied to the whole prior covariance at that step, R_t / delta: West and Harrison's feed-back intervention after a monitor's alarm. Each is logged.

Raises:

Type	Description
ValueError	on a y that is not a vector of finite values and NaN, an x whose shape is not (T, model.regressors) or that is not finite, or an intervention at a step outside y or with a discount outside (0, 1].
FloatingPointError	on a step whose one-step squared scale is not positive and finite, or whose update leaves a scale or a state that is not finite.

smooth(fit)

Retrospective moments of every state given all the data: Rauch-Tung-Striebel, with each step's covariance in the units E[V_t | D_T] of the variance the whole series supports. They are the moments of backward_sample's draws. A filter over geos is smoothed over its stacked state.

Raises:

Type	Description
ValueError	when the final degrees of freedom n_T is not above 2, so the covariances are infinite.

backward_sample(fit, draws, seed)

Forward filtering, backward sampling: joint draws of theta_(1:T) and V_(1:T) given D_T, with the volatility path when variance_discount < 1 (McAlinn and West 2019, appendix A.2). A filter over geos is sampled over its stacked state.

Raises:

Type	Description
ValueError	on draws below 1.

decompose(fit, draws)

The mean F_t' theta_t split into its parts, draw by draw: (draws, T, parts).

A Polynomial or Seasonal block is one part, its coordinates of F_t times its coordinates of the state: the level, or the seasonal effect. A Regression block is one part per column, the column times its coefficient. Parts follow the state's order, and a draw's parts add up to its F_t' theta_t to rounding.

Each part's interval is a quantile over the draws. A window's total, a difference between two windows, or a return per unit of a column is taken draw by draw first, so its interval carries the parts' dependence across steps and on each other.

Raises:

Type	Description
ValueError	on draws of another model's or another series' states.

constrained_sample(fit, signs, draws, seed, *, chains=4, warmup=None)

Draws of the states and the variance given the data with the coefficients of some columns of x held to a sign at every step: the model's posterior truncated to the signs, so the point estimate is the truncated posterior's mean, not a clipped one, and a channel's pull against its sign moves into the other states. signs maps a column of x to "positive" (at least 0) or "negative" (at most 0).

The discounts define the evolution variances W_t through the unconstrained filter, which depend on the design and not on y; with them the model is a Gaussian prior over paths, and the signs condition it. Gibbs sampling over three blocks, each chain started from an unconstrained backward draw moved into the signs (ADR 0042):

• the constrained coefficients' values with the other states integrated out, a Gaussian given V truncated to the signs, by one exact Hamiltonian trajectory of length pi / 2 (Pakman and Paninski 2014): the orbit is an ellipse, reflected where it meets a wall, with no step size and no rejection. The steps a zero evolution variance ties together (a discount of 1, or an idle held regressor) are one value;
• the other states given those paths, by forward filtering and backward sampling;
• V from its inverse gamma given the constrained values alone, the other states integrated out as well. Given every state, each evolution innovation would add a degree of freedom, which pins V to the path and slows its chain.

Truncating each backward-sampling step instead would not give this posterior, since the filter that feeds it never saw the signs. The values' covariance is held dense and factored once, so memory grows as the square of their number, and a trajectory meets more walls the harder the signs bind.

A constrained coefficient must be a Regression block of width one. rhat and ess say whether the chains mixed; the draws are not to be read when they are not ConstrainedDraws.mixed, and the run is then logged as a warning.

Raises:

Type	Description
ValueError	on no signs, a column that is not one of x's or not in a width-one regression block, a sign that is neither, form="multiplicative" (whose evolution couples the blocks), a variance discount below 1, a fit with interventions, draws below 1, chains below 2 or a negative warmup.

forecast(fit, horizon, x=None, *, evolution='constant')

Forecast horizon steps past the data, with x the regression columns over them.

evolution="constant" holds the first step's evolution variance W_(T+1) for every later step; "compounding" discounts again at every step. The variance's degrees of freedom decay by variance_discount a step and its estimate stays at S_T.

Under form="multiplicative" with blocks whose discounts differ, W = D P D - P is not a variance once the blocks are correlated enough, and holding it would drive the state variance indefinite; "constant" refuses such a W.

Raises:

Type	Description
ValueError	on a horizon below 1, an x whose shape is not (horizon, model.regressors), an unknown evolution, or, with evolution="constant", an evolution variance with a negative eigenvalue.

confidence_set(fits, parameters, level=0.95)

The fits a likelihood-ratio test at level does not reject, in the order given: those whose DLMFit.log_likelihood is within chi2_parameters(level) / 2 of the best. parameters is how many settings the fits vary over, two for a grid over two discounts.

An interval reported as the union of each member's interval, from its lowest lower end to its highest upper end, carries the uncertainty of the discounts as well as the states' (projection; Berger and Boos 1994). The best fit's interval alone does not, and a mixture weighted by the likelihood carries too little of it: the likelihood is flat where an evolution variance is weakly identified, and its best point there is a selection.

The union over-covers. On scripts/bench_dlm.py's 500 series a channel's 13-week contribution was covered 0.94 at a nominal 0.90, both on a random walk, where the pick covered 0.78, and on the discount model's own series, where it covered 0.85. The intervals there were 1.2 times as wide as at the true discounts.

Raises:

Type	Description
ValueError	on no fits, fits of other observations than the first's, a parameters that is not a positive integer, or a level outside (0, 1).

monitor_evalues(fit, *, shifts=(2.0,), inflations=(2.0,))

E-values for "the model's one-step forecasts hold", one row per observed step and one column per alternative, for chc.gate.DriftAlarm.

The standardised error u_t = e_t / sqrt(Q_t) is Student t under the model. Each alternative is a density for it -- the same t shifted by +h and by -h for each h in shifts, and scaled by k for each k in inflations -- and its e-value is the density ratio, whose expectation under the model is 1. West and Harrison's (1986) cumulative Bayes-factor monitor is e-CUSUM over the same ratios, and its usual threshold tau = 0.135 guarantees only that the average run length is at least 1 / tau, about 7.4 steps; chc.gate.DriftAlarm takes the average run length it guarantees as its threshold.

Raises:

Type	Description
ValueError	on no alternative, a shift that is not positive and finite, or an inflation that is not above 1 and finite.

stacked_prior(national, regional, geos)

The stacked state's prior: national on the national blocks and regional on each geo's copy of the regional blocks, every block independent of the others.

None stands for no blocks of that kind.

Raises:

Type	Description
ValueError	on two priors whose scale or dof differ, both None, or a geos that is not a positive integer.

forward_filter_geos(model, y, x=None)

Filter y, (T, geos), through model, with x, (T, geos, model.regressors), each geo's columns for the regression blocks. A missing observation is NaN.

Raises:

Type	Description
ValueError	on a y that is not (T, geos) of finite values and NaN, or an x whose shape is not (T, geos, model.regressors) or that is not finite.
FloatingPointError	on a step whose observed geos' one-step squared scale is not positive definite and finite, or whose update leaves a scale or a state that is not finite.

fit_geo_spread(model, y, x, pooled, draws, seed, *, level=0.9)

How far the geos spread about the national coefficients, by the marginal likelihood of y (type-II maximum likelihood, empirical Bayes), and draws of the spread from its posterior.

pooled are coordinates of a geo's regional blocks, counted from 0: each one's prior variance, the same in every geo and independent of every other coordinate, is the geos' spread on it, and the model's own is the vaguest considered. Each spread's log is searched from 1e-12 of the model's own up to it, from the model's own, by L-BFGS-B on central differences of GeoDLMFit.log_likelihood, which is exact, so the search is the only approximation. Each spread's interval at level is where the log-likelihood, the other spreads at their best, is within chi2_1(level) / 2 of the best, or the range's end where it stays so.

The posterior is under a prior flat on each standard deviation over the range (Gelman 2006), which unlike a flat prior on the log does not pile its mass at a spread of 0. draws spreads are drawn from it with seed by importance sampling in the range's logistic coordinates, log_spread = low + (high - low) expit(u): a Student t about the posterior's mode in u with its curvature there. Every draw is in the range, and the mode is inside it even where the best is at an end, as when the data would spread the geos wider than the model's own prior lets them. A quantity that carries the spread's uncertainty is the mixture of its posteriors over the draws, GeoSpread.mixture.

Raises:

Type	Description
ValueError	on pooled coordinates that are not distinct coordinates of the regional blocks, a pooled coordinate whose prior variance differs between geos or that has prior covariance with another, draws below 1, a level outside (0, 1), and what forward_filter_geos refuses.