Testing for stability of a simulated dynamical system

Question

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.

Answer 1

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.

Answer 2

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.

Answer 3

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

Answer 4

I like the autocorrelation idea, but there you will also need a parameter to decide when the autocorrelation changes are "small enough".

Idea: running averages:

Calculate a couple running averages: $$ m_{5} = \sum_{i}^{i-5} x_i \\ m_{20} = \sum_{i}^{i-20} x_i \\ m_{100} = \sum_{i}^{i-100} x_i \\ ... $$

If your system is at total equilibrium, then all of those should be the same, if you have only "thermal" fluctuations, then the changes will be in the shorter running averages but not in the longer ones. Depending on how different your simulations really are, you might be able to decide your system status on these running averages. You actually have a tradeoff at hand between numerical cost and how robust your criterion is..

This is a primitive version of the autocorrelation approach, but might just do the trick.

You will need a parameter though no matter what.