The estimator that got worse the more data I gave it

I spent a while building a corporate treasury system whose operational state is a coupled bistable dynamical system. Five nodes — invoice matching, bulk transfer, credit risk, supplier disbursement, cash pooling — each with cubic dynamics and a threshold θ past which it flips state, and a coupling matrix α that lets one flipping pull the others over.

One agent in that system alerts on a Kramers escape rate. Another integrates the network forward to answer what-if questions. Both are arithmetic on θ and α. And θ and α were chosen, not measured — which makes every number downstream internally consistent and externally meaningless.

So I wrote the estimator. This is what it found, what I got wrong when I first wrote it up, and the two explanations I proposed for that and then had to throw away.

The estimator, and the result I published

The model is

du_k/dt = (1/τ_k)[ u_k(1−u_k)(u_k−θ_k) + Σ_j α_kj g_kj(u_j) + F_k ]

Two structural properties make it tractable. It is linear in α once θ is fixed, so the couplings profile out by least squares and the objective becomes a function of θ alone. And node k’s equation contains θ_k and no other threshold, so a five-dimensional search collapses into five independent one-dimensional ones.

Derivatives are never differenced — finite-differencing a noisy trajectory amplifies exactly the noise you are trying to see through. Instead the fit matches u(t+h) − u(t) against the trapezoid of the right-hand side between the two samples.

Given a trajectory driven through the interior with no noise, it recovers θ to a mean absolute error of 0.052 and α to 0.032, and reports all five thresholds identified. At the σ = 0.012 the simulator itself uses, it identifies nothing, however the system is driven.

I wrote that up as: observation noise is the binding constraint, not sample size and not excitation. And I named the fix — an estimator that separates process noise from observation noise, a state-space filter or multiple shooting over windows — as future work.

Then I went back and ran it properly. Most of that conclusion does not survive.

The correction I should have made first

Before anything else: none of what follows is a new discovery, and I want that established before the numbers rather than after them.

Estimating drift parameters of a diffusion from discretely sampled data is a well-worked field, and the bias of naive discretisation is its central known problem. The literature is blunt about it — a recent treatment notes that “most applications still use the Euler–Maruyama discretization, despite many proofs of its bias” — and the standard responses have names and decades behind them: Kessler’s Gaussian approximation, Ozaki’s local linearisation, Aït-Sahalia’s Hermite expansions, and simulated-likelihood approaches.

I did not know that when I wrote the estimator. I walked into a trap that a statistician working on discretely observed diffusions would have seen coming, and I walked in through a door marked “this is just least squares.”

What follows is therefore not a contribution to that literature. It is what the known problem looks like from inside an applied system, where the person fitting the model is not a diffusion statistician, the fit is one component among fifteen, and every diagnostic on the screen says the estimate is fine. That is a situation worth documenting precisely because the theory does not reach the people it needs to reach.

First: there is no observation noise

The simulator adds noise to the state and carries it forward. The recorded trajectory is the true state of a stochastic process, sampled exactly. There is no measurement error anywhere in the system.

So “separate process noise from observation noise” is a fix for a problem I did not have. There is nothing to separate. A state-space filter estimates a hidden state from corrupted observations; here the state is observed perfectly and the corruption is in the dynamics.

I want to be precise about how that error happened, because it is the kind that survives review. The phrase is standard. It names a real technique. It is the correct fix for a problem that looks exactly like this one from a distance. And it appears in the conclusion of a document whose earlier sections are careful, which is what makes it read as rigour. It is a sentence that describes a different system.

Second: the estimator is inconsistent, and cannot tell

This is the finding that changed how I think about the whole exercise.

I ran the fit across independent noise realisations at increasing sample sizes, on the same forcing, and looked at the sampling distribution of θ̂ rather than at any single fit.

As the data grows sixteen-fold, the standard deviation halves — exactly as it should — and the bias doubles. Node 3’s bias grows by a factor of eight while its scatter falls by more than half.

More data makes the estimate worse. Not noisier. Worse, and more confident about it.

That is the signature of an inconsistent estimator: it converges, but not to the truth, and the shrinking error bars are a faithful report of how tightly it is converging to the wrong number. Every internal diagnostic improves while the answer deteriorates. Residuals look unstructured. Condition numbers look fine.

Part of the mechanism is visible once you ask where the information lives. θ locates the unstable fixed point between the two attractors, so the only data that constrains it comes from excursions through the interior — transients, incidents, the moments a node is actually moving between states. A longer run at fixed forcing does not add more of those. It adds more time in the stationary phase, which contributes bias and no information, so the informative fraction of the sample shrinks as the sample grows.

Which means the instinct that “we need a longer window” is precisely backwards, and the diagnostics you would normally use to catch that are the ones that fail.

Third: two explanations, both wrong

Having established that the estimator is biased, I tried to say why. I proposed two mechanisms and measured both. Neither survived, and I am reporting them because the falsification is the useful part.

Hypothesis one: endogeneity. The trapezoid evaluates the right-hand side at both endpoints, so the regressor contains f(u(t+h)) — and u(t+h) contains the noise increment sitting in the residual. Regressor correlated with error term. That predicts a left-endpoint scheme, which evaluates f(u(t)) and is measurable before the increment happens, should do better.

It does, and the direction is right. Same estimator, same row filtering, only the quadrature changed:

Right-endpoint is uniformly worst, which is what endogeneity predicts. But look at the last block. That is noise-free data, where there is no increment to be correlated with, and the trapezoid still loses by a factor of six. Whatever is happening is mostly not about noise.

Hypothesis two: Jensen. If the endpoint evaluation is the problem, the mechanism should be curvature — the drift is a cubic, evaluating a convex function at a stochastically perturbed point biases it by roughly ½·f″·Var, and Var is σ²h. That predicts bias scaling as σ².

Fitting bias = C·σ^p gives p ≈ 0.2 to 0.3, nowhere near 2. And the bias is already large at σ = 0.003, eight times smaller noise than the system uses. Whatever this is, it is substantially present at nearly zero noise, and it reverses sign at high noise — so there are at least two effects fighting.

Jensen is not it either.

What the measurement does establish. Holding total simulated time constant and varying the step size, on noise-free data:

Both errors halve when the step halves. Both estimators are first-order in h — including the trapezoid, which is a second-order quadrature and is not delivering second-order accuracy here. And the ratio between them is flat at ≈ 5.7 across an eightfold range of step sizes, which is a stable enough constant to be a real property of the setup rather than an accident of one configuration.

So the honest position is: there is a robust, reproducible factor of roughly six between the two quadratures; it is a deterministic discretisation effect rather than a noise effect; it is first-order in the step size for both; and I do not know what produces it. Clipping is the live suspect — the simulator bounds the state into [0,1], the trajectory sits at the upper bound for a large fraction of the record, and the trapezoid straddles two rows where the left endpoint uses one — but driving amplitude does not move the gap monotonically, so that is a suspicion and not a result.

I would rather publish the falsifications than a third confident mechanism.

Fourth: it refused for the wrong reason, and was right anyway

The estimator declares a threshold identified only if the profile-residual interval is narrower than 0.25 and does not run to the edge of the search range. The interval walks outward from the optimum until the residual sum of squares has grown by 5%.

A flat 5% is not a confidence level. At the typical n ≈ 7,900 usable rows, inflating RSS by 5% corresponds to a χ² of 395 on one degree of freedom. A 95% profile interval needs a ratio of 1.000487, not 1.05. The shipped threshold is roughly ten times too wide in standard-error units — a twenty-sigma interval.

The refusal was not a measurement. It was a number I picked that happened to be enormous.

So I calibrated it properly — 1 + χ²₁(0.95)/(n−p) — and re-ran across seeds:

Calibrating the threshold converts refusals into answers, and the answers are wrong far more often than they claim. Node 2’s nominally-95% interval contains the true threshold half the time. That is the bias from the previous section becoming visible: once the interval is the right size, it is centred in the wrong place.

The over-wide interval was covering the truth 100% of the time for the same reason a net the size of the ocean catches every fish.

So the shipped system refused for a reason that does not survive inspection, and the refusal was correct. I do not think there is a version of this where that is a satisfying outcome. It is a control that worked by accident, and a control that works by accident stops working silently when the sample size changes.

That interaction — a miscalibrated identifiability test masking a biased estimator, each hiding the other — is the one part of this I have not seen described elsewhere, and it is the part I would most expect to be sitting in other people’s code.

What I think this generalises to

The specific finding is about one estimator on one system, and the underlying phenomenon has been in the statistics literature for decades. The shape of the failure is what travels.

Any method that estimates dynamical parameters by regressing observed increments on functions of the observed state is exposed to this, and a great deal of applied work does exactly that — Koopman operator methods, EDMD, SINDy, one-step regressions generally. The theory that covers it lives under “discretely observed diffusions,” which is not a phrase most people writing that code would think to search for.

And the failure is silent in every direction a practitioner normally looks. The fit converges. The residuals are unstructured. The condition number is fine. The error bars shrink as the dataset grows, which reads as the model getting better. To see it you have to already suspect it, because the experiment that reveals it is not a diagnostic on one fit — it is a Monte Carlo across realisations and sample sizes, which nobody runs by default.

If you take one operational thing from this: run your estimator at two sample sizes and compare the bias, not the residual. It costs one extra run and it is the only cheap test I know that distinguishes “noisy” from “wrong.”

Where this leaves the system

The headline stands, with different reasoning underneath it. Passive observation of a healthy treasury cannot calibrate its tipping points — not primarily because of noise, but because the estimator is biased under the noise process the system actually has, and because the data that would overcome that is the data a healthy system by definition does not produce. A treasury that never tips generates no evidence about where its tipping points are, and watching it for longer moves you further from the answer while making you more confident.

Three things change concretely. The refusal threshold has to be a function of sample size rather than a constant. The quadrature has to be the first-order one, which is the opposite of what numerical analysis would tell you if you asked it about integration rather than about estimation. And the estimator that would actually work is one of the ones the literature already named — Kessler, Ozaki, Aït-Sahalia, or simulated likelihood — rather than anything I was going to invent.

What has not changed is the conclusion I care most about. The escape rates and the forward integration in that system remain arithmetic on chosen numbers, and there is now a runnable check that says so rather than returning four decimal places it cannot support. The check was right. It was right for reasons I had not measured, two of the explanations I offered for it were wrong, and one of the three sentences I wrote about it described a system that does not exist.

I would rather publish that than the version where the estimator works.