Probability density in phase space
Dynamic systems can be interpreted qualitatively by describing their phase spaces in terms of the type, order, and number of attractors (or areas of stability) that are traced out as the system evolves through time. Some systems may be described by a number of disjoint Lyapunov-stable areas (LSA), each separated in phase space by a separatrix (Kauffman, 1993). At any given time, the state of the system occupies only one such LSA, so that their number therefore constitutes the total number of alternative long-term behaviors, or equilibrium states, of the system.
Since an LSA is likely to be smaller than the total allowable range of states, the system tends to become boxed into an LSA unless it is subjected to external forcing. When the state approaches a separatrix, small perturbations can trigger a change to a nearby state, which can result in chaotic changes in the evolution of the system.
An important goal, then, is showing the presence of LSAs. But how do we show them?
Here is the original data--one million years of a proxy record for global ice volume (from Shackleton et al., 1990).
And here is a 2-D reconstructed state space of the last 500 thousand years of the above data set.
This is a mess. There is too much detail.
An approach I tried in Gipp (2001) (also appearing in Gipp, 2010) is to take a phase space reconstruction and smooth it out by calculating the probability density over some window of time (most recently I have been using windows of 270,000 years, but used windows of about 750,000 years in the 2001 paper).
There are two advantages to this approach: 1) simplification, making the portraits easier to characterize; and 2) reducing the impact of noise and dating errors, as the nature of the plot smoothes each point evenly over a region in phase space.
Regions of high probability may result from multiple visits to the same region of phase space, or by a drop the rate of evolution of the system (i.e., several consecutive states occupy the same general region of phase space). In either case, the region of higher probability can be interpreted as representing a region of stability, although whether it is asymptotically stable or Lyapunov-stable cannot be determined.
The method is straightforward, but computationally intensive.
Consider the following plot.
The plot represents the 2-D phase space from 130,000 years before present to 9,000 years before present. Yellow dots are plotted at 10,000-year intervals. The drop from 130 (thousand years ago) to 120 (thousand years ago) represents rapid deglaciation, to a global condition similar to that of the present. The interglacial lasted until about 70 (thousand years ago) whereupon global ice volume increased in two stages, reaching a maximum at about 20 (thousand years ago). It has been downhill ever since (the cores from which this data set are taken did not record the last 9,000 years of data--but we are reasonably sure from other sources that global ice volume has declined since then.
We take an area the size of the green square, and superimpose it over each of the evenly spaced red dots, and for each dot calculate the length of time the state space occupies the green square that is centred there. For instance, in the image above, the time spent in the green box is zero (at least over the interval observed above.
If we were to place the green box on the dot with coordinates (3.9, 3.9), we would find that the state space occupied the box for nearly 40,000 years. By contrast, the box placed at (4.0, 4.4) would be occupied for only about 3,000 years.
Here is what the intermediate plot might look like:
This data is actually from Raymo et al. (2004) and is a 2-D reconstruction of a 30 thousand year interval of the 13-C isotopic record of a deep-sea core.
The length of time between each blot is 3,000 years. The large numbers are the length of time that the state (the plot) spent within the 2x2 box centred on the number.
The 16.5 near the bottom left tells us that the function spent 16.5 thousand years in the 2x2 square centred on that number 16.5.
In order to come up with one set of contourable probability densities, I would string together nine of these plots (totalling 270 thousand years), sum them up, and contour the results, which would look something like this:
The regions in phase space labelled A2, A3, etc. all represent regions in phase space with a high probability of being occupied by a system state at a randomly selected time, over the interval from 729 thousand years ago until 459 thousand years ago. In practical terms it means that most of the time, the ice volume was either near A2 (interglacial), A3 (weak glacial), A4, A5, or A6 (all of which are fairly strong glacial conditions). Each of these labelled regions is a region of stability.
The stability most likely occurs because the system dynamics lead to negative feedbacks dominating while within the areas of stability, and positive feedbacks dominate while the system is evolving from one region of stability to another.
Therefore we would state that in the interval (or window) from 729-459 thousand years ago, there were five areas of stability in the global ice volume 2-D state. In particular, there were five particular global ice volumes (likely there were five global configurations for the ice volume) which were metastable.
But the best is still to come . . .
Gipp, M. R., 2010. Epsilon machine computation in paleoclimate proxy data: Implications for testing climate models. Geocanada 2010 Proceedings.
Gipp, M. R., 2001. Interpretation of climate dynamics from phase space portraits: Is the climate system strange or just different? Paleoceanography, 16: 335-351.
Kauffman, S., 1993. The Origins of Order: Self-Organization and Selection in Evolution, Oxford Univ. Press, New York, 734 p.
Raymo, M.E., D.W. Oppo, B.P. Flower, D.A. Hodell, J. McManus, K.A. Venz, K.F. Kleiven, and K. McIntyre, 2004. Stability of North Atlantic water masses in face of pronounced climate variability during the Pleistocene. Paleoceanography, 19, PA2008, doi:10.1029/2003PA000921.
Shackleton, N. J., A. Berger and W. R. Peltier, 1990. An alternative astronomical calibration of the Lower Pleistocene timescale based on ODP site 677. Trans. Royal Soc. Edinb. Earth Sci., 81, 251-261.