# Testing for stability of a simulated dynamical system

### Background and question

I often work with simulations of dynamical systems and I usually track a single parameter $x$, such as the number of agents (for agents based models) or the error rate (for neural networks). This parameter usually has some interesting transient behavior. The type and length of this behavior depends on the random seed and simulation parameters. However, after the initial transient period, the system stabilizes and has small (in comparison to the size of the changes in the transient period) thermal fluctuations $\sigma$ around some mean value $x^*$. This mean value, and the size of thermal fluctuations depends on the simulation parameters, but not on the random seed, and are not known ahead of time.

The goal is to have a good way to automatically stop the simulation once it has transitioned from the initial transient period to the long-term behavior. We know that this transition will only happen once, and the transient period will have much more volatile behavior than the stable state.

### Naive approaches

The naive approach that first pops into my mind (which I have seen used as win conditions for some neural networks, for instance) is to pick to parameters $T$ and $E$, then if for the last $T$ timesteps there are not two points $x$ and $x'$ such that $x' - x > E$ then we conclude we have stabilized. This approach is easy, but not very rigorous. It also forces me to guess at what good values of $T$ and $E$ should be.

Alternatively we can use a time $T$ and confidence parameter $\alpha$ as in this answer and assume a normal distribution on the error. This will save us the effort of knowing the size of thermal fluctuations.

Let $y_t = x_{t + 1} - x_{t}$ be the change in the time series between timestep $t$ and $t + 1$. When the series is stable around $x^*$, $y$ will fluctuate around zero with some standard error. Take the last $T$, $y_t$'s and fit a Gaussian with confidence $\alpha$ using a function like Matlab's normfit. The fit will give us a mean $\mu$ with $\alpha$ confidence error on the mean $E_\mu$ and a standard deviation $\sigma$ with corresponding error $E_\sigma$. If $0 \in (\mu - E_\mu, \mu + E_\mu)$, then you can accept. If you want to be extra sure, then you can also renormalize the $y_t$s by the $\sigma$ you found (so that you now have standard deviation $1$) and test with the Kolmogorov-Smirnov test at the $\alpha$ confidence level, or one of the other strategies in this question.

### Notes

This question was originally posted to Cross Validated in a more time-series centered language. It generated some good discussion, but not an answer I was happy with. I tried MetaOptimize but the question generated zero interest there. I am not familiar with the tagging on Computational Science, so feel free to edit the tags.

-
What is it you're not happy with in your third "naïve approach"? – leftaroundabout Jan 9 '12 at 23:31
@leftaroundabout I would prefer if I could eliminate the $T$ parameter, since it often depends on the simulation parameters and thus I have to guess it. I would also prefer standard tools that I could reference (if such are available), instead of coming up with my own ad-hoc method and having to take space to explain and justify it. – Artem Kaznatcheev Jan 9 '12 at 23:36
If there is an "automated" method for this, it would be of enormous benefit in a number of different fields—including molecular simulations (my field). The fact that it would be so desirable and generally useful leads me to think that a rigorous method is either not available or very poorly publicized, if it exists. I'd be interested to see what the response is, if any. – aeismail Jan 10 '12 at 0:00
Any chance you could post a link to a typical trace for the parameter $x$? – Bill Barth Jan 10 '12 at 3:33

By your use of "thermal fluctuation", I assume you have some statistical mechanics framework. In that case you might use quantities that are known to converge monotonically, e.g., the entropy of the Boltzman-weighted probability distribution.

If you have a well-defined state space, you could look for an equivalent of the separatrix to determine whether you have crossed into an attractive region, e.g., compute the net flux from one region of state space to another (the 50:50 region could be the separatrix, which would be "transient"). This approach is sometimes taken to identify transition states of chemical reactions.

-
 I used "thermal fluctuation" in a loose sense, since I have more background in stat mech than in stats. Basically I meant random perturbations caused by ambient noise/randomness in the system. I will see if I can adapt your suggestion, thank you. – Artem Kaznatcheev Jan 10 '12 at 16:35 If you want to see an example of the idea for transition states check out this paper – Deathbreath Jan 15 '12 at 19:41

It seems to me that your problem is simply that you don't have a precise definition of "transient" and "stable" behavior. If you did, you could just check when the defiinition is satisfied. Until you use a precise definition, you don't even have a way of determining if one approach for distinguishing the two behaviors is better than another.

One possible source of a better definition: in the theory of continuous dynamical systems, there is a concept of center manifold; all solutions decay exponentially toward the center manifold but travel slowly along the center manifold. It might be possible to use this to come up with a quantitative definition of what you're calling transient behavior.

-
As far as I can tell, though, from a dynamics perspective, transience and stability are not numerically well-defined quantities. They're subject to interpretation (at least with respect to specifying values to test against). Perhaps that's part of what the OP is looking for here. – aeismail Jan 10 '12 at 11:00

As @aeismail comments, there doesn't seem to be a general method.

What I often do is take the autocorrelation and end the simulation when it is no longer changing. Advantages:

• Often part of a library (R, Fortran, Python, etc)
• Efficient
• Minimal parameter tweaking required
-
 how do you decide when the autocorrelation function is no longer changing? Wouldn't you require an arbitrary sensitivity there? (Since the autocorrelation of a stable system subject to noise will still fluctuate over time). Also, do you calculate the autocorrelation of your whole time series so far, or just the last $T$ time-steps? If you do the first, then your autocorrelation will for sure be constantly changing as you hit the stable region, since you will become more and more uncorrelated to the transient region. – Artem Kaznatcheev Jan 12 '12 at 18:14 If you need an error parameter for the sensitivity in deciding "no longer changing" and a window size $T$, then I don't understand the advantage over the naive approach. If I am misinterpreting, could you expand your answer to help me understand better? – Artem Kaznatcheev Jan 12 '12 at 18:15 I don't think that you are misinterpreting anything. If autocorrelation is on advantage vs. the naive approach depends on the problem. I usually find that autocorrelation method is less sensitive to parameter tweaking. – Eelvex Jan 12 '12 at 18:39