Complexity, bifurcations, catastrophe

What makes a system complex?

It is a perplexing problem--both its description and its quantification. One might think that the description of a system as complex would suggest it has many subsystems each acting in accordance with its own rules, and interacting with each of the other subsystems in ways that we find difficult to describe. But there are systems involving very few "parts" which exhibit the kind of behaviour we call complex.

Another possible definition may stem from the notion of the compressibility of the system's information. Is it difficult to describe the sequential outputs in a manner that is simpler than merely listing all of our observations? A good random number generator would exhibit such behaviour, but we would not describe that as complex.

John Casti has proposed that the complexity of a system is at least partially dependent on the observer. He uses an example of a rock lying on the ground. To the layperson, there are only a limited number of ways to interact with the rock (kicking it, breaking it, throwing it, etc.) . To a trained geologist, there are more (as before, as well as mass-spectral geochem, x-radiography, electron probe, etc.). So the rock seems more complex to the geologist, but that additional complexity actually stems from the observer.

Another example can be share prices, where the complexity depends on the timing of your observations. If you look at the price of, say, Anadarko Petroleum over the past year, using closing prices only. (for disclosure--no position).

Then we can look at one-minute increments on a daily chart (Nov. 6, 2013).

Note that the variability within the two charts doesn't seem all that different despite the changes of scale. Lastly we could look at how trading in Anadarko looks over one second. One particular second, that is, between 3:59:59 and 4:00:00 on May 17 of this year.

Source

You probably wouldn't expect a lot of change over 1 s, but in this case you would be wrong: the price fell from about $90 to $0.01 in less than 50 ms. That's a loss of $1 billion in market capitalization per millisecond--keep losing money at that rate and before long you're talking real money!

This all seems to trigger a philosophical debate--is the complexity present when none of the observers are capable of seeing it? In the case of Anadarko, if you were a pension fund, your losses would have been real (although the trades were all cancelled and reversed after market close).

If the complexity of a system arises from within, then what characteristics do we ascribe to complexity. One characteristic is discontinuous behaviour, particularly when the inputs to the system are continuous. For instance, tectonic processes gradually cause stresses to accumulate in an area in a fairly uniform fashion, until a critical threshold is reached and the earthquake occurs.

The branch of mathematics that investigates the sudden onset of convulsions wrought by a slow change is called catastrophe theory. Catastrophe theory is generally considered to be a branch of bifurcation theory. By bifurcation we normally mean some change in the operations of a complex system. It could represent a transition from one stable state to another. It could also represent the development of new areas of stability in phase space (or their disappearance) or simply a change in the nature of a chaotic attractor.

In particular, sometimes the sudden appearance of a new mode of stability is brought about by the changing value of a slowly migrating parameter past a critical threshold. Such behaviour is called a catastrophe, in the mathematical sense.

In the past few years we have seen a major change in the mode of operations in the markets. In particular, the rapid growth of high-frequency trading has added complexity at timescales where such behaviour did not previously exist. This is another example of a catastrophe.

Inference of dynamics for complex systems, part 3

Equilibria in linear systems

It is common for linear systems to evolve towards a single fixed state.

The assumption that most dynamic systems have a single preferred fixed state for a given set of circumstances appears to be one of the drivers behind Central Bank policies. It is debatable whether the Central Banks recognize nonlinearities in socio-economic systems, or whether they do but are unable to express it it the typical 10-second soundbite that the average consumer of their "products" can absorb.

Certainly, the common example of using lower interest rates to combat unemployment sounds like the kind of thinking one expects about a system with a single equilibrium state. The notion that low interest rates ultimately lead to higher inflation reflects similar thinking.

Here are some examples in phase space of single equilibria in a simple linear system.

When all trajectories within a given region of phase space gradually and continually converge towards a single point, then the equilibrium can be described as being asymptotically stable, and is commonly referred to as a point attractor. Note that even though trajectories 1 and 3 appear to cross, they actually do not as one passes over the other in a three- (or higher-) dimensional projection.

Line crossings in a properly "unfolded" phase space portrait cannot cross, as each uniquely defined state represents a unique value and a unique sequence of values of the higher time derivatives of the time series from which they are projected. If two trajectories converged onto a single point, it would be impossible for them to diverge again--hence no crossings. If there appear to be crossings in a reconstructed phase space, then the phase space needs to be projected in more dimensions. Unfortunately, there are logistical problems with presenting even three (let alone more) dimensional figures, which is why I have limited the figures presented to two dimensions, with a caveat that apparent crossings mean we should really be looking at a projection in at least three dimensions.

Often we find that rather than continuous convergence, two states which originate close to one another as well as to a particular point in phase space tend to stay close together and near to the point. Such stability is called Lyapunov stability (sometimes a Lyapunov-stable area, or LSA). More formally, we might state that all states with a region of phase space (a) will tend to remain in another (possibly smaller) region of phase space (b). Once again, note that the apparent crossing of trajectories only suggests that we should construct this phase space in at least three dimensions.

Another form of equilibrium is the limit cycle. It represents a form of asymptotic stability whereby the final equilibrium is not a point but a continuous cycle. There is no reason that the cycle needs to form an ellipse--any closed shape is possible. In higher-dimensional projections, a limit cycle may be in the form of a shell, or a torus.

Equilibria in chaotic systems

A different form of equilibrium was discovered by Edward Lorenz in 1963. A system defined by three nonlinear differential equations reached a complex equilibrium in which any state evolved continually towards the "attractor"; however any two states starting arbitrarily close would diverge exponentially, even though at the same time they would both evolve in similar fashion through phase space.

This form of attractor was called a "strange attractor". Although it occupies a finite volume of phase space, the trajectory from any arbitrary point would evolve in unique fashion, so that all evolving trajectories from any arbitrary starting point would appear similar, yet would occupy subtly different regions of phase space.

Multiple equilibria (multistability)

In some complex systems, we observe any 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 (a bifurcation), resulting in chaotic changes in the evolution of the system [Parker and Chua, 1989]. Thus very complex behavior can arise in multistable systems.

Multistable behaviour generally arises from systems in which feedbacks, both negative and positive, impact a system which is perturbed by some sort of external forcing. Negative feedback tends to resist external forcing, resulting in stability in some regions of phase space. If the external forcing is sufficient to overcome the feedback, then positive feedback may actually accelerate the changes, causing the state to evolve rapidly towards another area of stability.

The smooth variation of one or more parameters in the system may result in a change in the type of or the number of attractors in the system, or even in the order in which the attractors are visited. Such a response is called a bifurcation. Bifurcations can represent a sudden transition within a system characterized by purely chaotic attractors to one with one or more LSA; between one LSA and multiple LSA; or between different configurations of LSA.

Bifurcations represent changes in the organization of the system, and their existence has been suggested by models [e.g., Ghil, 1994; Rahmstorf, 1995], and more rarely from observations [Livina et al., 2010; discussed here]  in future installments we will demonstrate such behaviour in natural systems. Initially, however, we will concentrate on interpretation of some of the phase space portraits presented last time.

References


Ghil, M. (1994), Cryothermodynamics: the chaotic dynamics of paleoclimate, Physica, 77D: 130-159.

Kauffman, S. (1993), The Origins of Order: Self-Organization and Selection in Evolution, Oxford Univ. Press, New York, 734 p.

Livina, V. N., Kwasniok, F., and Lenton, T. M. (2010), Potential analysis reveals changing number of climate states during the last 60 kyr. Climate of the Past, 6, 77-82, doi: 10.5194/cp-6-77-2010.

Parker, T. S. and L. O. Chua, (1989), Practical Numerical Algorithms for Chaotic Systems, Springer-Verlag, New York.

Rahmstorf, S. (1995), Bifurcation of the Atlantic thermohaline circulation in response to changes in the hydrologic cycle. Nature, 378, 145-149.

Recognizing change in complex systems: excursions vs. bifurcations

Continued from last time.

Once again, we have a fifteen-year plot of gold/copper ratio vs. silver/rough rice ratio. We are continuing our discussion of whether the event labelled C, which is still unfolding, is likely to be an excursion (which will then return to the region populated by most of the graph) or a bifurcation (which will lead us to a new area of Lyapunov stability somewhere new in phase space).

It appears to be at least a once-in-a-generation event. But how significant is it?

For this spike to represent a bifurcation as opposed to an excursion, what we would have to see the function settle in around a new area of phase space. Looking at the shape of the spike, the likely area for such a new area of stability will be centred near the "C", and could extend from the far right of the graph to the kink in the curve near (350, 1.5). For this to form with any degree of satisfaction is likely to take about two more years. If we don't see any sign of orbiting about the "C", then the most likely outcome is a return to the 15-year area of stability.

I found it interesting that on this graph, the silver spike of 1998 does not appear. It is lost in the middle of area of stability.

Is it possible that the late spike in silver is due to its comparison to a soft commodity, which perhaps haven't performed as well as metals? We can check this by changing our ratios slightly.

Same four commodities--but compared differently. This time we are looking at the silver($/oz) to copper ($/lb) ratio against gold to rough rice. There are three significant excursions, labelled A, B, and C.

Excursion A represents the spike in silver price in 1998 due to the Warren Buffet purchase. Excursion B is the rise in price in copper in 2006, which exceeded (in percentage terms) the gain in silver which occurred at the same time.

Excursion C has two phases. For the first six months or so, the excursion is dominated by increased gold prices (compared to rice), and for the last 13 months or so, the excursion records the outperformance of silver.

Once again, we won't be sure that this is a bifurcation as opposed to an excursion unless the phase space settles somewhere near or above the "C" for at least another two years.

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

It is obvious gold has not behaved well lately. Some of the stocks have been hit as well. We haven't looked at Detour Gold Corp. for awhile.

The recent price action has mirrored that of gold, with a sharp rise from the end of July and a sharper decline in the past few trading days.

The two-dimensional reconstructed phase space of the last eight months of the price of DGC-T.

The phase space portrait occupies a Lyapunov-stable area at around $30 for about six months, and from the beginning of August, the system evolves towards the $37 area. The system completes an orbit in the $37 range before leaving the area. This behaviour suggests that there may be another LSA in the $37 range, but I would feel more comfortable with it after a few more orbits.

What happens next? The two most likely scenarios would appear to be a return to the $30 LSA, or a return to the incipient LSA in the $37 range. Being long Detour, I would be encouraged by a reinforcement of the $37 LSA. Unfortunately, this outcome is nearly impossible. The reason is in the nature of the construction of the graph. Recall that the phase space reconstruction is generated by the price at closing of a particular day (on the vertical axis) against the closing price four trading days earlier (horizontal axis).

The coordinates of the last point are ($37.60, $31.00). In four days, the value of the x-coordinate is going to be $31 (today's closing price). In three days, the value of the x-coordinate is going to be $31.17 (yesterday's closing price). I don't know what the closing price of Detour will be on those days--but if it is at about the current price, then the state will lie within LSA30. The only possible hope for the Detour price state to reach the orbit at $37 would be for the price to recover to about $37 within the next three trading days--by Friday. Even then the state would not lie within LSA37, but at least the trajectory would indicated that that would be the likely outcome.

In fact, if we don't see some reassuring action tomorrow--when x falls to $33.60--we will be right on the outer edge of LSA30.

There is another, less desirable outcome. The price may continue falling to the previous LSA near the $23 level. Doh!

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

Time to liquidate more cocoa. Sadly, the good Ghanaian stuff is all gone. I am forced to take the South American stuff. But I can at least sweeten it a little with this.

Recognizing change in complex systems

"So I came downstairs and was surprised to find the dog reading the paper."

Such is the beginning of a typical shaggy dog story. If the story is true to form, the punchline would be something like the dog has been getting his news from the internet for years, or some such. If this were a true story, it would be that the dog is usually sleeping when I come downstairs. In such a case, I would presumably notice right away that the dog's behaviour was unusual.

One of the outstanding problems of complex dynamic systems is recognizing (hopefully in real time) when a change in behaviour is occurring.

Today we look at the long (ish)-term behaviour of some commodities with respect to one another. I haven't thought much about their economic significance. Today will just be arm-waving at charts, with a view to see if we can recognize the presence of change in the market fundamentals.

First up--some soft commodities. Below we see the ratio of cocoa prices to rough rice (contracts as defined in 1996).

The principal impact on the graph is from the first Ivoiran civil war. I find it interesting that the second war didn't really impact the price as much. Possibly this is because of rapidly increasing Ghanaian cocoa production. According to Voice of America, Ivoiran production is doomed due to years of under-investment.

The reconstructed phase space appears below. I have used the time delay method. To construct this figure I have smoothed the above data set through a 3-pt moving average filter, as the unfiltered data looked really noisy.

The fifteen years of data are confined to a comparatively small region of state space, with the only interesting feature the large excursion during the first civil war. This event lasted about two years. Note the magnitude, and more importantly, note the outcome--the system reverted back to the same area of phase space from which it began.

The second civil war by contrast is scarcely noticeable.

There are other ways to contruct phase space portraits. Instead of reconstructing them from a single time series we can build them from scatterplots of two (or more) time series. In studying geological systems I try to avoid this technique because it is difficult to build two time series of equivalent length and similar sample rates. For commodity time series this is less of a problem.

So let's go all out and plot the gold/copper ratio against the silver/rough rice ratio (all month-end closing prices).

IIRC, I have multiplied the silver price by 100 in the above chart. Silver and gold used $/oz, copper measured in $/lb.

Notice how for most of the fifteen year record, the states all plot within a relatively large oval (let's call this an LSA) within phase space. There have been three significant deviations from the LSA over the past fifteen years. The excursion marked A represents the strength in gold during the 2008-9 global meltdown. The excursion marked B is a rise in tandem of both silver and copper during the year 2006.

The last excursion (C) represents the recent rise in silver. This ongoing excursion appears to have lasted 17 months so far. The multi-billion dollar question is whether this is an excursion expected to revert to the  LSA, or has a bifurcation occurred, with the system evolving towards a new LSA somewhere else in phase space.

Refer once again to the phase space plot of cocoa/rough rice.

Bifurcations and excursions are both rare events. The excursions occur on at best a decadal scale. Bifurcations are more rare. Clearly they have happened in the past. For instance, the long-term gold/silver ratio was approximately 16 for centuries, but has not been near this ratio for several decades. A bifurcation occurred in the last century sometime. Perhaps we are in one now. Either way, I think it might be prudent to stock up on rice.

Disclosure: long gold, long silver, long copper, long cocoa, long rice. Sadly, I was forced to liquidate some of my cocoa holdings in the recent market turbulence.

But it was comforting.

Update (April, 2012) - I overlooked this story in the explanation for the rising price of cocoa in 2010. It seems to fit in with the first peak in the cocoa/rice ratio in 2010.

Bifurcation in gold-silver ratio! Where are we headed?

Is the recent exciting rise in the price of silver telling us something? Let's investigate.

The World Complex presents a two-dimensional reconstructed phase space portrait (using the time delay method) for monthly average gold and silver prices from January 1996 to April 2011. The time delay is twelve months, meaning that this plot shows the gold/silver ratio plotted against the ratio one year earlier.

Silver pricing data comes from here. Gold price monthly data from the World Gold Council website.

Don't get lost in the details. The important thing isn't every loop and swirl in the trajectory, nor is it even the trajectory itself. It is the space. For almost the entire data set, the system has cycled within a large area, which may have been an area of Lyapunov stability. The only exceptions are the excursion to the lower right area of the graph (an unusually high ratio coupled with an unusually high rate of increase in the gold/silver ratio, which happened in late 2008 to early 2009--I'm sure we all remember that one); and two excursions to the far left of the graph (in which the ratio is low and had also declined rapidly in the previous year)--these being the "Warren Buffett event" of 1998, and our present excursion.

Notice that the Warren Buffet excursion was short-lived. Indeed, the sudden reversal of the trajectory looks a little unnatural. We eagerly await to see whether our current excursion will be similarly short-lived, and marked by a similar reversal. Unfortunately, although we are eager, the time taken for such an event to play out is months to years, so patience is the order of the day.

The next several months will be crucial. If there are truly problems in silver supplies, it would be logical to expect to see the trajectory of our dynamic system leap into a new area of phase space--probably one beyond the boundaries of the graph.

One idea currently circulating is that bifurcations are preceded by a period of extreme stability (what we might term low volatility). It is not clear if the period of low volatility is dynamically necessary (clearing the road, as it were), or is merely an empirical observation. The volatility of the gold/silver ratio over the past couple of years has been breathtaking, so if a bifurcation is occurring, it would provide a counterargument to the extreme stability hypothesis.

Bifurcation to come in our economic decline?

I have been trying to finish this entry forever, but power has been out here for almost three days now due to an electrical storm a few nights ago. I have been running off portable batteries, charging once in awhile when we run the generator.

One debate I seem to have over and over again with one of my fellow expats here in Ghana is the success (or lack thereof) of central banks in controlling the unfolding economic decline.

This friend of mine insists that central banks will do what they have always done, which is oversee a gradual decline of economic power in the West, while ensuring enough volatility that it is hard for someone to sit in a long (or short) position and profit all the way down.

My thinking is informed by the behaviour I have observed in dynamic systems. Systems behave simply and change slowly over more than 99.9% of observable time, but behave unpredictably and change very rapidly during the remaining <0.1% of observable time.

The behaviour of dynamic systems, like climate, is nearly linear much of the time, but is characterized by episodic reorganizations. These reorganizations resemble car accidents, in that your life afterwards may have no similarity to your life beforehand. Consequently, trying to predict the behaviour of a system after a bifurcation is as difficult as predicting you will be in a car accident next week.

About all we can do is anticipate when the bifurcation will take place.

Central banks normally try to hide what they are doing, because they don't want the general population protecting themselves against the consequences of their actions. If you need deflation, you also need the general population to be hurt by it, otherwise there is no benefit. However this requires a prediction of what the general populace is going to do, and I think that this will prove the weak point of central bank operations.

For instance, the future operation of the system may well depend on people continuing to make their mortgage payments, even if the value of the house has fallen below the amount owing. Between the holy hell of MERS and all the recent foreclosure failings, it seems likely that more and more people are going to give up on those mortgage payments.

Whenever the central bank tries to interfere with the workings of the economy in order to bring about a desired change, the general population is going to have to act in the manner predicted by the genius economic modelers at these banks. But modeling is really hard, in this case because the behaviour of the population is unpredictable.

Here is an example. Some years ago the Canadian Mint introduced a one-dollar coin with a picture of a loon. The coin was made of nickel. The official opinion expressed was that Canadians would love these things and would call them “dollar coins” or “dollars”. But even before the coin was released, Canadians began calling them “loonies”, a name which has lasted to the present day.

A few years later, the Mint decided to introduce a two-dollar coin. This coin depicted a polar bear (now two polar bears). Once again, there was an official opinion—the public would call this coin the “Bear”. Of course, they were almost immediately known as “toonies”.

So much for the official prediction of public behaviour.

When the tipping point comes—and it likely won’t be recognized until afterwards—the changes will be broad and deep. I don’t know what our new world will be like. But I would guess it would have a lot less government and the economy would be more focussed on actual production and real wealth creation rather than the growth of paper assets.

Decline of variability when approaching bifurcation point: an alternative view of charting

One of the significant problems facing climate scientists is how to predict when the system is about to undergo a major change in organization. This point is important because a precise prediction of what is going to happen in response to continued alteration of the atmosphere is unlikely. What we may be able to do is look for signals of an impending global reorganization in the hopes that we will recognize same before the change suddenly occurs.

The basis of this discussion is that the climate system acts to resist external forcings, but if the forcing is strong enough, the system then evolves very rapidly towards some new equilibrium. We may or may not be able to predict what the new equilibrium will be.

It is generally viewed that a change in climate equilibrium would be bad for us as our entire system of global agriculture is geared towards our present equilibrium, which has lasted for about ten thousand years (which, come to think of it, is about as long as we have been practicing agriculture). We haven't put much effort into planning for agriculture in a new equilibrium state.

Once again, understanding the past is the key to predicting the future.

A number of papers have been presented recently, in which the authors have looked at climate signals prior to major climate shifts in the paleoclimate record. And they have noticed that prior to a major change in climatic regime, the short-term variability declines dramatically.

Carbonate content of marine sediments during a "climate tipping" event from greenhouse conditions to icehouse conditions ca. 34 million years ago. From Thompson and Sieber (in press). The lower figure tracks the rise of an autoregressive parameter--essentially an inverse measurement of variability.

The reason for this decline in variability is unclear, however it appears in both models and in observed time series. Given the detailed observations made of present-day climate, it is theoretically possible to anticipate when a sudden climate change is about to occur.

There are counter-arguments, however. Climate measurements are normally easier to make than ecological measurements.

What does this mean for the stock market?

Some well-known charting patterns show the same behaviour; most notably the triangle pattern.

This does seem to be at odds with some conclusions I've drawn earlier. I may have to reconsider.

Consider it evidence that dynamic systems are hard.






Reference

Thompson, J. M. T. and Sieber, J., in press. Predicting climate tipping as a noisy bifurcation: A review. International Journal of Bifurcation and Chaos.

Bifurcations in the paleoclimate system: an alternate view

In our last installment, we looked at a sequence of probability density plots of the 2D global ice volume state and inferred that the global climate has undergone changes in organization which we call bifucations.

A reviewer of a recently submitted paper asked why I use a 2D reconstruction when all of the metastable states appear to lie on a line, suggesting that one dimension might be enough.

Let's see what he means.

This is one frame of the animation in our last installment. (I suppose this should be Gipp, in rev.) The five Lyapunov-stable areas A2 to A6 all appear to lie on a line, approximately defined by y = x. So do we need two dimensions?

The above is a frame from a poster from the AGU fall meeting in 2008. The probability density is calculated from the last 750,000 years of the 2-d reconstructed phase space portrait from the ice volume proxy of Shackleton et al. (1990). The curved line is the 2-d phase space reconstruction over the last 125,000 years. At right we see the last 500,000 years of the proxy record.

By inspection we see that for the phase space portrait to plot in the A2 metastable state, the value of the proxy must be between 3.8 and 4.0. This is a necessary condition but is not sufficient for the state to plot on A2.

Similarly, for the A3 metastable state, the value of the proxy must lie between about 4.1 and 4.3, which is again necessary but not sufficient. To see why, consider the portions of the reconstructed state space near 128,000 years (128 ka), near 68 ka, and near 5 ka (the three blue circles on the reconstructed state space). At all three of these times, the proxy function has values between 4.1 and 4.3, yet none of these points correspond to the LSA--something which would not be obvious purely on the basis of the proxy function.

Such is the hazard of trying to determine the existence of metastable states on the basis of the time series without "unfolding" it into a more appropriate dimension (please note that 2-D is not an adequate projection either, but it does eliminate almost all of this sort of error).

The state space over the interval approximately from 58 ka until about 35 ka plots within the A4 LSA (note blue arrow connecting the appropriate regions on the two plots above). However, at the interval near 11 ka, the proxy function has the appropriate value for the A4 LSA (about 4.6)--nevertheless, the state space does not plot within the LSA, and again provides misleading information if one merely calculated probability density in one dimension.

The reason that the LSAs all line up is due to the topological equivalence between the two dimensional phase space plotted using the time lag method (as above) with one using the time derivative method. There are types of LSAs that would not plot along the y=x line--in the type of plot we have above we would see two corresponding patches representing the intersection of a toroidal limit cycle with the plane of the plot, each of which would be on opposite sides of the y=x line. The physical meaning of such an LSA would be of an oscillating glacial state--not impossible, but I can't interpret it from the record.

The LSAs A1 to A6 (in the animation) all represent relatively stable global ice volumes--achieved through a kind of dynamic stability that takes into consideration such factors as the slope of the continental shelves, anchoring points on the shelf or in coastal areas and lithospheric and asthenospheric reaction to ice distribution.

Livina et al. (2010) use a one-dimensional approach from a data set to infer the existence of multistability in the climate system over the past 60,000 years using proxies of Greenland temperatures over that interval.This is exactly the problem we are trying to solve (different data set though)!

The basic approach is to find the probability density of the data from the time series (what we have called one-dimensional here).

A brute force method would be as follows:

Just for fun I'm using a different data set here--this is believed to be a proxy for Himalayan monsoon strength over the past 2+ million years.

To find probability density by a brute force technique, you would draw numerous boxes (yellow) and calculate the length of time the proxy function plots within the box. As the box is slowly migrated from bottom to top, and these different lengths of time are plotted (green wiggle at right--which I have semi-arbitrarily sketched--sorry, I haven't calculated this), you have a probability density plot of the data. This answers the question of whether certain values are more common than others. The more "popular" values may represent metastable states (or Lyapunov-stable areas).

Now, to make things even more interesting you might do what I did in creating the animation in a previous post. Instead of calculating the probability density over the entire data set, you calculate it over a succession of "windows" of time (100,000 years duration, say--the choice depends on what you hope to show). Then you would have a probability density graph for each of these different windows (and there is no reason why the windows couldn't overlap).

If some of your windows had four probability peaks, and others had two or even one, then voila! you have demonstrated a bifurcation.

This is the essence of the argument of Livina et al. (2010). Now, dear reader, if you have been paying attention you will see the problem.

The technique just described overestimates the time spent actually occupying the stable states. To determine how much for the particular data set of Livina et al. (2010) would require some time with the data--which I'm not sure I can spare now.

Livina and her coauthors do not use the brute force method described above. They estimate the probability density using a Gaussian kernel estimator. Okay: I would do the same thing if I had a record as densely sampled as theirs. (There is a two-dimensional estimator that works on the same lines).

They calculate a potential based on the linear probability density estimate and a noise estimate, and approximate the shape of their (linear, remember) potential estimator by a polynomial of even order. The potential estimate [U(z)] has a well (local minimum) corresponding to peaks in probability.

People, people, please. If I've told you once, I've told you a thousand times*--use Chebyshev polynomials! Seriously. You'll thank me later.

The advantage of orthogonality is key. Chebyshev polynomials are orthogonal whereas standard polynomials are not. (For our purposes, the orthogonality means that the integral from -1 to 1 of the product of two Chebyshev polynomials (say Tn(x) and Tm(x) where m is not equal to n) is zero. Sorry, blogger is crappy at displaying this.

Because of the orthogonality of the individual functions, you can use a Fourier-type process to transform your original function into one composed of Chebyshev polynomials (technically this is a Fourier-Chebyshev transform). Such a transform has many similarities to a Fourier transform, such as the limitation on the highest-order Chebyshev polynomial (called the Nyquist Chebyshev number) being defined by your sample interval.

Codes for FC transform should be in Numerical Recipes. They were in the 1989 edition.

Data sets are noisy? Yup! (Paleoclimate sets are notoriously so). You can design filters in the Chebyshev number domain and eliminate some of that unpleasant high-Chebyshev-number noise (like high-frequency noise, only better).

You can also design filters in the Chebyshev number domain, which can be applied to your function without the need for a Fourier-Chebyshev transform by reverse transforming the filter into a series of convolution coefficients and convolving with the original data series. Just as in standard filter design, your filter has to be designed without sharp corners so that the convolution series remains small.

You can define the significance of the different probability peaks in terms of the relative significance of (in the two state case) the fourth order Chebyshev polynomial with respect to the first or second, or of the sixth order polynomial for a three state case.

Livina and her coauthors use a disingenous method of defining the number of wells in their potential function. Since their estimate for the potential is a polynomial of even order, their contention is that the number of wells is better calculated by counting the number of inflection points (where the second derivative is zero--or where the direction of curvature reverses) rather than the number of local minima.

It appears that they may be better off counting both. Here's why.

Here we see a potential function U(z) which is a polynomial of even order. The "+" represent regions of the curve where curvature is positive and the "-" represent segments of the curve with negative curvature.

There are four inflection points in the graph at left. Using equation 7 of Livina et al. (2010), we would infer the presence of three stable states.

By inspection we see only one real potential well, as well as one degenerate solution where both the first and second derivatives are zero. This type of point (just to the left of the origin) is of tremendous significance in dynamic systems analysis, but goes unnoticed by simply counting the inflections.

Compared to the Gaussian kernel estimator, looking at all critical points is a simple matter and is better not overlooked.

The Livina et al. (2010) methodology would be much improved by working in more than one dimension. I have worked in two, but for the data sets I am working on, going to three would be an extremely intense computational exercise. But for data sets with annual sampling over 60 ky or more, going to a three-dimensional probability estimator of the density of points in a three-dimensional phase space would be a huge step forward.

* I know, I know--I've never told you. But I'm telling you now.

References:

Gipp, M. R., 2008. Semiquantitave Characterization of Dynamic Response of Long Records Paleoclimate Records. [damn typos] American Geophysical Union, Fall Meeting 2008, abstract #PP13A-1428.

Gipp, M. R., in rev.  Imaging system change in reconstructed phase space: Evidence for multiple bifurcations in Quaternary climate. Submitted to Paleoceanography (in revision).

Livina, V. N., Kwasniok, F., and Lenton, T. M., 2010. Potential analysis reveals changing number of climate states during the last 60 kyr. Climate of the Past, 6: 77-82. doi:




Update November 14--see a new animation here

Probability density in phase space part 2: Detecting bifurcations

Continuing from last time.

We can use a succession of probability density plots (PDPs) in sequence to show whether the number of Lyapunov Stable Areas (LSAs) changes, or whether the locations of individual LSAs moves. If the PDPs are displayed quickly enough, they may be displayed as a movie.

The significance of this is that a change in the distribution of LSAs may be indicative of a bifurcation in the system. A bifurcation is suggestive of a reorganization of the system, often appearing as a sudden response to a gradual change in one or more system parameters which has crossed some critical threshold.

Let's look at at the bifurcations in the global ice volume suggested by a sequence of PDPs from the past two million years approximately (original data from Shackleton et al., 1990).



There are a number of book-keeping details in this video. The sliding scale bar at the right shows the interval of time represented by each frame. There is a little black bar representing 30,000 years of missing data (that section of the core was badly undersampled.

The big red interval on the time bar labelled "MPR" represents the time of the "Mid Pleistocene Revolution" (also called the Mid-Pleistocene Transition) which represents a major change in the operation of the climate system.

The probability density of every frame was calculated on the basis of 270,000 years of data, so the blocks that include the undersampled bit actually represent the available data over an interval of 300,000 years. There is a possibility that some bias may have resulted.

As is normal for geological projections, the earliest PDP is presented first. The movie therefore flows forwards in time, and we see that in general, over the past 2 million years, ice volumes have increased (low ice volume is at lower left and higher ice volume is at upper right). The system has become more complex, as there are generally more stable zones in the later part of the record as opposed to the former.

Much of the complexity can only be appreciated by going back and looking at the actual state space superimposed on top of the probability density plot. One way of establishing the complexity of the climate system is by computing the complexity of the simplest epsilon machine that reproduces the patterns observed within the data (see also here--will do an entry on this at some point).

There is another abstract here which expresses the same ideas with different terminology.

Next in the series we will compare this method with alternate approaches taken from the literature.