chc.dynamics_id¶
Stability tier: stable.
Causal identification of the residual's control channel, not just of a scalar effect.
chc.train fits a residual by prediction error. Under a confounded logging policy that is the
wrong object: if the historical action was chosen from a covariate that also drives the state,
then d(dx/dt)/du fitted on (x, u, x_next) is the observational response, and a
controller that optimises against it moves the state the wrong way. That is the same failure the
scalar headline benchmark shows -- and until now nothing protected the residual from it.
This module closes that gap with Robinson partialling-out lifted from a scalar effect to a
state-dependent matrix. With nuisances g(x,z) = E[y | x,z] and m(x,z) = E[u | x,z], where
y = (x_next - x)/dt - f_known(t, x, u) is the part the residual must explain::
y - g(x,z) = B_θ(x) (u - m(x,z)) + eps, E[eps | x, z, u - m] = 0
so the estimator regresses the residualised state rate on the residualised action, with
K-fold cross-fitting of both nuisances. Neyman orthogonality is what buys the guarantee: the
channel error is second order in nuisance error, which is exactly Results 18 (order transfer
p -> 2p) and 19 (debias every channel via cross-fit Robinson DML). This module is those
results' missing consumer, not a new one.
HONEST SCOPE, four limits worth stating before the code:
- Only the channel
B_θis interventional.a_θis fitted on the remainder and therefore absorbs whatever the omitted confounder contributes to the drift in-sample, so it is an observational-conditional drift. Planning is unbiased in the direction the optimiser moves; the predicted trajectory level still shifts if the confounder's distribution does. On real data this bites harder than it sounds: a confounder that trends with the state gets charged to positive feedback inx, so the fitted drift can come back unstable even where the channel is clean -- and an MPC uses the drift too. Check the drift's spectrum, and give measured exogenous drivers a place in the drift as regressors instead of leaving them only inadjust_for. - The residual must be control-affine (
chc.residual.ControlAffineResidual). A generalr_θ(x, u)has no partialling-out moment and gets no guarantee here. - With no adjustment set and no instrument nothing in the log identifies the channel, and an
instrument identifies it only along the directions it moves the action. The estimator reports
identified=Falsein both cases rather than returning a confident wrong answer -- that case belongs tochc.sensitivity, which prices the radius instead of pretending it away. - The premise is Markov: given
x,zandu, a transition's noise does not depend on the past. Where the noise persists within a unit and a lever persists too, the state is a common effect of the past lever and the past noise, soE[eps | x, z, u - m]is not 0 and the channel is biased by an amountchannel_errordoes not cover.persistence_checkreads the noise's persistence off the moment's residual; adjusting for the state, the levers andza period earlier removes the bias where the noise is AR(1) (ADR 0064).
Integrator = Literal['euler', 'rk4']
module-attribute
¶
Which one-step map the fit is asked to be consistent with. See fit_causal_residual.
ConfoundedControlAffineSystem
dataclass
¶
Euler-logged dx/dt = f_known(x,u) + A x + B_true u + C z under u = K z + eta.
The vector-state analogue of chc.causal.ConfoundedLinearSystem and the DGP
fit_causal_residual is scored against. z is the confounder: it drives the logged
action through K and the state rate through C, so with z withheld the response of
the rate to the action is B_true + C sigma_z^2 K^T / var(u) rather than B_true, where
var(u) = K K^T sigma_z^2 + G G^T + sigma_eta^2 collects every source of action variance.
A scalar channel with B=1, C=2, K=-1.5, sigma_z=1, sigma_eta=0.5 and no instrument
gives 1 - 2*1.5/2.5 = -0.2: the sign flip of the headline benchmark, reproduced per
channel. Turning the instrument on adds G G^T to that denominator and dilutes the bias
toward zero -- the exogenous shifter makes the log more informative, so the observational fit
looks less wrong while being no more identified. Worth knowing before reading a number off it.
The Euler discretisation is deliberately the one fit_causal_residual inverts, so a
failing recovery test means the identification is wrong, not that the integrator disagreed.
sample(n, key, known)
¶
Draw n transitions: x (n,d_x), z, u, x_next, w (instrument).
CausalDynamicsFit
dataclass
¶
A fitted causal residual next to what is and is not known about it.
identified is the load-bearing field: False means nothing the fit was given pins the
channel down from the log -- no adjustment set and no instrument, or an instrument whose moment
falls short of rank (instrument_rank) -- and the residual is then the observational fit,
kept only so the caller can compare.
error_cosines
property
¶
nu, the cosines of the periods the channel's error was read with where the fit was
given periods, counted over the distinct periods, and the degrees of freedom of the
t a test read off it is sized against; None where it was given none.
OmittedConfounderBound
dataclass
¶
How far a confounder the adjustment set left out could move a linear functional of the
channel, and how strong it would have to be to move it to null. Experimental.
bias_scale is identified; strength is the analyst's, as Gamma is
(docs/concepts/gamma.md), in partial R^2s rather than odds.
ChannelRegion
dataclass
¶
A region that holds the fitted channel's coefficients at a stated level, from one log.
Experimental: it may change or be withdrawn in any release. State s's coefficients
c_s (the channel's row s, raveled lever by lever, each lever's features in turn) lie in
{c : ||c_hat_s - c||_V <= radius[s]} for every state at once with probability at least
level. V = D'D + ridge P is the matrix the fit's least squares solved, D the
channel's design on the action residuals. channel_region says where it holds.
contains(channel)
¶
Whether every state's row of channel, shaped as channel, lies in its
ellipsoid.
operator_radius(state)
¶
A bound, at the region's level, on ||B_hat(x) - B(x)||_2 at the state x: the
largest move of the rate that the channel's error makes per unit action norm. It holds at
every state at once, since each is read off the same ellipsoids.
rate_error(states, reach)
¶
A bound, at the region's level, on ||(B_hat(x) - B(x)) a|| over the (N, n)
states and every action deviation a with |a_j| <= reach[j]: the rate error
that the plan's tube prices as one standard error (chc.decision's _model_error),
read here at a level. The largest of a convex form over a box is at a corner, so it reads
the 2^m corners.
stability_level(rho, state)
¶
The largest level at which operator_radius at state is within rho: the
probability the region gives that the channel's error there is at most rho, 0 where
no level reaches it.
A certainty-equivalent gain whose stability holds within a radius rho of the channel
(chc.regret.ce_explicit_constant_certificate) then stabilises the plant with at least
this probability. That certificate's rho is the true plant's, and its channel is the
discrete step's, dt times this one under Euler: pass its radius over dt. Read at
the estimate, rho is the estimate's, which the region does not check.
Solves sum_s (R_s A + b_s)^2 = (rho / g)^2 for A, g the operator radius per
unit of the radii, and A^2 = 2 log(n / delta) + log_det_ratio for delta
(validation/tube_levels.mac).
PersistenceCheck
dataclass
¶
Whether a fit's transition noise persists within units, read off its channel moment's
residual by persistence_check. Experimental.
The fit's premise is Markov: given the state, the levers and the adjustment set, a transition's noise does not depend on the past. Where the noise persists within a unit and a lever persists too, the state is a common effect of the past lever and the past noise, so given the state this period's lever and noise move together, and the channel is biased by an amount its error does not cover (ADR 0064).
ClosedLoopAttribution
dataclass
¶
Which coefficients of a control-affine fit the log identifies, on a tracked plant.
ClosedLoopAttributionCertificate
dataclass
¶
Two arms: an interaction the plant has, and one the tracking loop manufactures.
solve_channel_moment(state_residual, action_residual, states, *, instrument_action=None, degree=1, ridge=1e-06, weights=None)
¶
Solve E[(y_res - B_θ(x) u_res) (x) (z (x) phi(x))] = 0 for the channel.
The Robinson moment on its own, separated from the nuisance estimation
fit_causal_residual wraps around it. Public for two reasons: cross-fitted
ridge-polynomial nuisances are a default rather than a commitment, and the point of a debiased
score is that g and m may come from any learner -- gradient boosting, a neural net, a
model the caller already had -- so the caller needs a way in that does not go through ours.
It solves the moment it is handed and checks no rank. Along a direction of the channel that
instrument_action does not move, the moment holds at every channel, and the ridge picks
one. Whether an instrument identifies the channel is the caller's to settle before the solve,
which sees neither the raw shifter nor the covariates the nuisances took out, and a rank that
decides it reads both.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
state_residual
|
Array
|
|
required |
action_residual
|
Array
|
|
required |
states
|
Array
|
|
required |
instrument_action
|
Array | None
|
the action-shaped variable |
None
|
degree
|
int
|
monomial degree of |
1
|
ridge
|
float
|
Tikhonov term on the moment's Gram matrix, on each coefficient scaled by its regressor column's mean square, weighted as the rows are, so that it reads the same in any units of the actions and the state. Added as a constant, it outweighed actions logged in units a millionth of their own, and set their channel near zero. |
1e-06
|
weights
|
Array | None
|
one weight per transition, shape |
None
|
Returns:
| Type | Description |
|---|---|
Array
|
The channel coefficients, shape |
Array
|
attr: |
fit_causal_residual(known, data, dt, *, adjust_for=(), instrument=None, degree=1, channel_degree=None, nuisance_degree=2, folds=2, ridge=1e-06, seed=0, integrator='euler', drivers=(), weights=None, influence=False, clusters=None, periods=None)
¶
Fit a ControlAffineResidual whose channel is the interventional control response.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
known
|
Dynamics
|
the physics kept fixed; its own control dependence is known, not estimated. |
required |
data
|
dict[str, Array]
|
columns |
required |
adjust_for
|
tuple[str, ...]
|
observed confounders. Empty and no instrument => |
()
|
instrument
|
str | None
|
name of an exogenous action shifter, for when the confounder is latent. It
enters as the moment's instrument, not as the regressor, so this is real 2SLS. Expect a
variance premium, not a free lunch: identification here rides on however much of
the action the shifter explains, and on the reference DGP that is 18%, which costs
roughly an order of magnitude in channel error against adjusting for a logged
confounder (0.10 vs 0.002 at It identifies the channel only where its moment has rank: where the first stage's push
on the actions moves with the actions along every direction of the channel the actions
move (all but The rank does not grade an instrument: a column of noise drawn apart from the action
keeps it, and reads identified with a channel that misses. Its relevance grades it: on
the reference plant such columns read 0.002 to 0.035, the medians of eight at 500,
2000 and 8000 rows, of the order of |
None
|
degree
|
int
|
the feature degree of the drift, and of the channel unless |
1
|
channel_degree
|
int | None
|
the feature degree of the channel; |
None
|
nuisance_degree
|
int
|
flexibility of |
2
|
folds
|
int
|
cross-fitting folds. |
2
|
ridge
|
float
|
the Tikhonov term of the nuisances, the channel's moment, the drift regression and the instrument's first stage. The nuisances read standardised covariates. The drift regression and the first stage scale each slope's term by its column's variance, beside a free bias; the moment, which has no bias, scales each coefficient's by its column's mean square, on the log's raw actions. So the ridge reads the same in any units of the state, the actions, the drivers and the instrument. Added as a constant, it outweighed a column logged in units a millionth of its own: such actions read a channel near zero, and such a state a drift slope near zero. |
1e-06
|
integrator
|
Integrator
|
the one-step map the fitted field is made consistent with.
Up to 0.6.0 the fixed point was iterated instead -- add the leftover to the target,
refit -- and each state stopped once its defect's RMS stopped falling. That is not where
the fixed point is: the RMS flattens under the noise first, and the iteration contracts
the drift only by A log that falls by more in one step than RK4 can follow has no The default stays |
'euler'
|
drivers
|
tuple[str, ...]
|
exogenous columns that push the state and that nothing the plan does can move --
weather, a demand forecast. |
()
|
weights
|
Callable[[Array], Array] | None
|
a weight on each transition's channel moment, a function of its state alone.
Experimental: it may change or be withdrawn in any release. A callable is called with
Under a channel class that contains the truth every weight estimates the same channel,
and only the variance moves. The inverse of the rate's noise variance is the efficient
weight in the limit, but it loads the states where the noise is small, which can be
where the nuisances extrapolate: on a log whose rate noise ran Under a class that misses the truth, each weight estimates a different channel: the
projection of the true one under |
None
|
influence
|
bool
|
whether to keep each transition's influence on the fitted parameters,
|
False
|
clusters
|
ndarray | Array | None
|
each transition's cluster, |
None
|
periods
|
ArrayLike | None
|
each transition's period, a whole number of steps, on a log whose transitions
fall in one group, one unit's or the units' of one declared cluster, which leaves no
second cluster to sum within: one label a transition, two periods at least. The
channel's covariance then sums the scores of each period's transitions over them and
their states, and reads the sums' dependence across the periods by |
None
|
Returns:
| Name | Type | Description |
|---|---|---|
A |
CausalDynamicsFit
|
class: |
Raises:
| Type | Description |
|---|---|
ValueError
|
on a negative |
omitted_confounder_bound(fit, functional, *, cf_y, cf_d, rho=1.0, level=0.95, null=0.0)
¶
Bound a linear functional of the channel against a confounder adjust_for left out.
Experimental.
Chernozhukov, Cinelli, Newey, Sharma and Syrgkanis (2022) bound the bias of a linear
functional of a partially linear fit by |rho| sqrt(cf_y cf_d / (1 - cf_d)) sigma nu:
sigma^2 the rate's residual variance after the fit, nu^2 the second moment of the
functional's Riesz representer, both identified, and two shares that are not. cf_y is the
share of that residual variance the latent explains; cf_d the share of the long
representer's second moment the fit's short one misses, which for one lever is the latent's
partial R^2 with it. rho is the correlation between the two gaps, 1 at the worst. On a
linear Gaussian plant with one latent the bound is attained, for any functional of any number
of levers, and with two latents it is strict unless they move the lever and the rate in
proportion (validation/omitted_confounder_bound.mac). Over several states the bound sums
sigma_s nu_s, taking the shares as common to them; under Euler a functional of one
state's channel reads that state alone.
functional weighs the channel's coefficients, in its shape (n, m, n_features). One
lever's effect on one state is a single 1. Moving spend from lever a to lever b is
+1 at b and -1 at a, and since its value is linear in the channel, its
robustness_value is the confounding at which the move stops paying.
The sampling error follows DoubleML's: each confidence bound is one-sided at level, and
carries the estimate's influence and bias_scale's, whose nu^2 part is nu^2 -
alpha^2 for a representer alpha. Its sign convention for the estimate differs from this
one in the cross term of the two influences, which vanishes in expectation. A fit given
clusters has both influences summed within each cluster, as its channel's error is, and
one given periods reads them across its periods by cosines, each period's transitions
summed first. The bounds take a t's quantile with G - 1 degrees of freedom, G the
clusters, two-way the smaller dimension's count, or the rows where there are none; or nu,
the cosines', under periods (ADR 0075). At 5 units, a 5 % test of the channel read off its
error rejected 9.75 % and 14.75 % of panels against the normal's quantile, and 3.25 % and
5.75 % against t(4)'s (ADR 0062). DoubleML takes the normal's.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
fit
|
CausalDynamicsFit
|
an identified fit of |
required |
functional
|
Array
|
the weight on each channel coefficient. |
required |
cf_y
|
float
|
the share of the residual variance the latent explains, in |
required |
cf_d
|
float
|
the share of the long representer's second moment the short one misses, in
|
required |
rho
|
float
|
the correlation of the two gaps, in |
1.0
|
level
|
float
|
each confidence bound's one-sided level, in |
0.95
|
null
|
float
|
the value the robustness values measure the distance to. |
0.0
|
Raises:
| Type | Description |
|---|---|
ValueError
|
on a fit whose channel is not identified by adjustment -- an instrument's
representer is not the moment's -- or that is weighted, or kept no representer; on a
functional of another shape; on shares, |
channel_region(fit, level=0.95, *, bound=None)
¶
The self-normalised confidence ellipsoid of the fitted channel, at level.
Experimental: it may change or be withdrawn in any release. Abbasi-Yadkori, Pal and
Szepesvari (NeurIPS 2011, Theorems 1 and 2): for y_t = <X_t, c> + eta_t with X_t
predictable and eta_t conditionally R-sub-Gaussian, the ridge estimate c_hat
satisfies ||c_hat - c||_V <= R sqrt(2 log(det(V)^(1/2) det(V0)^(-1/2) / delta)) +
||c||_V0 with probability at least 1 - delta, V = V0 + X'X. The regressors need only
be predictable, so a log the plant was controlled in qualifies: a closed-loop policy's actions
move with the states and the noise before them. Here V0 = ridge P, the fit's own ridge,
one ellipsoid a state at delta / n.
It gives the plan's rate error a level (ChannelRegion.rate_error) and the
certainty-equivalent certificate a probability (ChannelRegion.stability_level),
beside the one-standard-error scales the tube carries, which it does not move.
Premises, and where they fail:
- No hidden confounder. The noise must be mean zero given the actions; a confounder the
fit does not adjust for pushes the rate with the action, and the ellipsoid misses the
channel, at any
level. - No persistent noise.
etamust be a martingale difference; AR noise is not, unless its lags are in the fit (persistence_check). - R is read off the log: each state's root mean square of the moment's residual, over
N - qrows, a plug-in for the theorem's knownR, which also counts the cross-fitted nuisances' error. The action residuals use nuisances fitted on other rows, so the regressors are predictable to the nuisance's first order only. - The prior term
||c||_V0isboundwhere given, else the estimate's own, a plug-in; at the default ridge it is about 1e-3 of the coefficients' effect on the rate.
Under a hidden confounder an instrument's sign-perturbed sums (Szentpeteri and Csaji, 2023) stay valid; they are not built here.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
fit
|
CausalDynamicsFit
|
an identified Euler fit by adjustment, unweighted, whose actions move every direction of the channel. |
required |
level
|
float
|
the joint level over the states, in |
0.95
|
bound
|
float | None
|
an upper bound, for every state, on the true coefficients' norm in the ridge's
metric, |
None
|
Raises:
| Type | Description |
|---|---|
ValueError
|
on a fit that is not identified, read through an instrument, fitted under
|
persistence_check(fit, units, periods)
¶
The lag-1 autocorrelation of fit's moment residual within units. Experimental.
units and periods label the fit's transitions, N each: the unit, of any one
sortable kind, and the period the transition starts in, a whole number of steps. Two
transitions of one unit whose periods differ by one are a pair. Each state's residual is
divided by its root mean square, so each state counts alike in any units. The pairs' products,
summed over the states, are summed within each unit before the test, since the transitions of
one unit need not be independent; the statistic is their total over its CR1 spread, each
unit's sum less its pairs' share of the total, against t(G - 1) (ADR 0062). Pairs of one
unit alone, or of a fit given periods, whose units share one group's noise, are read
across their periods instead: the products summed within each period, and the total over the
spread of their projections on nu cosines of the pairs' periods, against t(nu), as
the fit reads its error (ADR 0075). Pairs that all fall in one period leave no spread to read.
A rejection says the premise fails, not by how much the channel is off: that depends on how
far a lever persists as well (PersistenceCheck). Adjusting the fit for the state, the
levers and the adjustment set a period earlier removes the bias where the noise is AR(1): on
panels of 200 units over 100 periods whose noise and lever are AR(0.7), a channel of 0.8 read
0.0130 low, 0.0007 low with them, and 0.0250 low with the state's lag alone (ADR 0064).
Raises:
| Type | Description |
|---|---|
ValueError
|
the fit kept no moment residual, as where its channel is not identified; the labels are not one a transition; a period is not a whole number; or a unit has two transitions starting in one period. |
closed_loop_gain_attribution(states, actions, rates)
¶
Attribute a control-affine fit's coefficients to the controller and to the plant.
Result 41 left one item open: why the interaction coefficient b1 of
dx/dt = d + a*x + (b0 + b1*x)*u comes out large and negative on a setpoint-tracked zone.
The answer is a property of the log, not of the plant. A proportional loop
u = gain*(setpoint - x) + u0 puts every sample on an affine manifold x = c + m*u with
m = -1/gain, and restricted to that manifold the four-term class collapses to a quadratic in
the action::
d + a*x + (b0 + b1*x)*u = (d + a*c) + (a*m + b0 + b1*c)*u + (b1*m)*u^2
Only b1 reaches the u^2 term, so b1 is the identified coefficient and the pole is
not -- the inverse of the usual reading, and the reason Result 41 (a) found the reported pole
to be a units artefact while lambda(u) = a + b1*u held. Matching against an observed
response A + B*u + C*u^2 gives
``b1 = C/m = -gain*C``,
so the sign is decided by the curvature of the response in the action and the magnitude by the
controller: a tighter tracker reports a bigger interaction from identical physics.
proofs/closed_loop_attribution.v proves the matching identity and exhibits the explicit
one-parameter family that leaves d, a and b0 free.
This also refutes the standing guess that b1 = -1/gain. That is the manifold slope,
not the interaction; equating them forces C = 1/gain^2, a constraint on the plant rather
than an identity (guess_is_a_constraint_not_an_identity).
exploration_budget is what buys the rest back: variation off the manifold restores the
design's rank, and with it the separation of a from b0. Reading a pole off a fit whose
budget is ~0 is reading the regulariser.
closed_loop_attribution_certificate(gains=(0.5, 1.0, 2.6, 8.0), exploration=(0.0, 0.05, 0.25, 1.0), samples=4000, seed=0)
¶
Separate the interaction a plant HAS from the one a tracking loop MANUFACTURES.
Arm A: the plant really is d + a*x + (b0 + b1*x)*u. Under exact tracking the design is
singular, yet b1 comes back exactly while the drift a does not -- the non-identification
is real and it is the pole that suffers, not the interaction.
Arm B: the plant has no interaction, but its response carries a u^2 term the fitted
class cannot represent. The fit answers with b1 = -gain*curvature: a coefficient that is
entirely an artefact of misspecification amplified by the loop, growing linearly with the
controller's gain. That is the mechanism behind Result 41's large negative b1, and why
the number cannot be read as authority-falls-with-temperature without checking the budget.
The exploration sweep is the remedy: variation off the manifold restores the design's rank and
the drift with it. All three claims are measured here, and the middle one is the one that would
embarrass the entry if b1 turned out to be as unstable as a.