Catastrophe theory has been developed as a deterministic theory for systems that may respond to continuous changes in control variables by a discontinuous change from one equilibrium state to another. A key idea is that system under study is driven towards an equilibrium state. The behavior of the dynamical systems under study is completely determined by a so-called potential function, which depends on behavioral and control variables. The behavioral, or state variable describes the state of the system, while control variables determine the behavior of the system. The dynamics under catastrophe models can become extremely complex, and according to the classification theory of * Thom*, there are seven different families based on the number of control and dependent variables.

Let us suppose that the process y_{t} evolves over t = 1,…, T as

dy_{t} = -dV(y_{t}; α, β)dt/dy_{t} —– (1)

where V (y_{t}; α, β) is the potential function describing the dynamics of the state variable y_{t }controlled by parameters α and β determining the system. When the right-hand side of (1) equals zero, −dV (y_{t}; α, β)/dy_{t} = 0, the system is in equilibrium. If the system is at a non-equilibrium point, it will move back to its equilibrium where the potential function takes the minimum values with respect to y_{t}. While the concept of potential function is very general, i.e. it can be quadratic yielding equilibrium of a simple flat response surface, one of the most applied potential functions in behavioral sciences, a cusp potential function is defined as

−V(y_{t}; α, β) = −1/4y_{t}^{4} + 1/2βy_{t}^{2} + αy_{t} —– (2)

with equilibria at

-dV(y_{t}; α, β)dt/dy_{t} = -y_{t}^{3} + βy_{t} + α —– (3)

being equal to zero. The two dimensions of the control space, α and β, further depend on realizations from i = 1 . . . , n independent variables x_{i,t}. Thus it is convenient to think about them as functions

α_{x} = α_{0} +α_{1}x_{1,t} +…+ α_{n}x_{n,t} —– (4)

β_{x} = β_{0} + β_{1}x_{1,t} +…+ β_{n}x_{n,t} —– (5)

The control functions α_{x} and β_{x} are called normal and splitting factors, or asymmetry and bifurcation factors, respectively and they determine the predicted values of y_{t} given x_{i,t}. This means that for each combination of values of independent variables there might be up to three predicted values of the state variable given by roots of

-dV(y_{t}; α_{x}, β_{x})dt/dy_{t} = -y_{t}^{3} + βy_{t} + α = 0 —– (6)

This equation has one solution if

δ_{x} = 1/4α_{x}^{2} − 1/27β_{x}^{3} —– (7)

is greater than zero, δ_{x} > 0 and three solutions if δ_{x} < 0. This construction can serve as a statistic for bimodality, one of the catastrophe flags. The set of values for which the discriminant is equal to zero, δ_{x} = 0 is the bifurcation set which determines the set of singularity points in the system. In the case of three roots, the central root is called an “anti-prediction” and is least probable state of the system. Inside the bifurcation, when δ_{x} < 0, the surface predicts two possible values of the state variable which means that the state variable is bimodal in this case.

Most of the systems in behavioral sciences are subject to noise stemming from measurement errors or inherent stochastic nature of the system under study. Thus for a real-world applications, it is necessary to add non-deterministic behavior into the system. As catastrophe theory has primarily been developed to describe deterministic systems, it may not be obvious how to extend the theory to stochastic systems. An important bridge has been provided by the Itô stochastic differential equations to establish a link between the potential function of a deterministic catastrophe system and the stationary probability density function of the corresponding stochastic process. Adding a stochastic * Gaussian white noise* term to the system

dy_{t} = -dV(y_{t}; α_{x}, β_{x})dt/dy_{t} + σ_{yt}dW_{t} —– (8)

where -dV(y_{t}; α_{x}, β_{x})dt/dy_{t} is the deterministic term, or drift function representing the equilibrium state of the cusp catastrophe, σ_{yt} is the diffusion function and W_{t} is a * Wiener process*. When the diffusion function is constant, σ

_{yt}= σ, and the current measurement scale is not to be nonlinearly transformed, the stochastic potential function is proportional to deterministic potential function and probability distribution function corresponding to the solution from (8) converges to a probability distribution function of a limiting stationary stochastic process as dynamics of y

_{t}are assumed to be much faster than changes in x

_{i,t}. The probability density that describes the distribution of the system’s states at any t is then

f_{s}(y|x) = ψ exp((−1/4)y^{4} + (β_{x}/2)y^{2} + α_{x}y)/σ —– (9)

The constant ψ normalizes the probability distribution function so its integral over the entire range equals to one. As bifurcation factor β_{x} changes from negative to positive, the f_{s}(y|x) changes its shape from unimodal to bimodal. On the other hand, α_{x} causes asymmetry in f_{s}(y|x).