Ecological Classical Political Economy
Workshop Notes. Otto-Friedrich-Universität Bamberg, April 4-6 2022
Oriol Vallès Codina, Research Fellow, Leeds University Business School, JUST2CE Project.
April 4 - 6 2022
Schedule
| Monday, Apr 4 | Tuesday, April 5 | Wednesday, April 6 |
|---|---|---|
| Basics and Varieties of Dynamical Systems | Lotka-Volterra Prey-Predator Oscillations | Early Agent-Based Models in Complexity Science: Cellular Automata and the Ising Model |
| Discrete Models: Logistic Growth, Stability and Chaos | Economic Application 1: The Goodwin Model of Capital Accumulation and Class Struggle | Economic Application of the van der Pol Oscillator: Kaldor’s Business Cycle |
| Economic Application: Ricardian Growth in a World of Limited Resources | Economic Application 2: The Flaschel-Semmler Model of Multi-Sector Growth |
Complex Dynamical Systems
Belousov-Zhabotinsky Chemical Reactions
Phase Transition in the Ising Model
Complexity theory studies highly organized but decentralized systems composed of very large numbers of individual components, with common features that transcend its differences in scale, material components, and organizing laws. They are:
- complex, showing the potential to configure their component parts in an astronomically large number of ways,
- adaptive, undergoing constant change in response to environmental stimulus and their own development,
- self-organizing, exhibiting a strong tendency to achieve recognizable, stable patterns in their configuration, and
- non-equilibrium systems, avoiding stable, strictly self-reproducing states
Such models first emerged as applications of corporate-managerial problems in the mid-1950s out of systems thinking, a causality-driven, holistic approach that describes the interactions between the components of a system and their interaction with its environment.
For Brian Arthur (2014, 92), complex adaptive systems are characterized by six features that posit difficulties for the traditional mathematics used in economics:
- dispersed interaction among (heterogeneous) agents,
- existence of no global controller,
- cross-cutting hierarchical organization,
- continual adaptation,
- perpetual novelty, and
- out-of-equilibrium dynamics
Five “Big Ideas” that distinguish complexity economics from traditional (neoclassical) economics (table 4.1, Beinhocker 2006, 97)
| Complexity Economics | Neoclassical Economics | |
|---|---|---|
| Dynamics | Open, dynamic, nonlinear systems, far from equilibrium | Closed, static, linear systems in equilibrium |
| Agents | Modeled individually; use inductive rules of thumb to make decisions; have incomplete information; are subject to errors and biases; learn and adapt over time | Modeled collectively; use complex deductive calculations to make decisions; have complete information; make no errors and have no biases; have no need for learning or adaptation (are already perfect) |
| Networks | Explicitly model interactions between individual agents; networks of relationships change over time | Assume agents only interact indirectly through market mechanisms (e.g., auctions) |
| Emergence | No distinction between micro- and macroeconomics; macro patterns are emergent result of micro-level behaviors and interactions | Micro- and macroeconomics remain separate disciplines |
| Evolution | The evolutionary process of differentiation, selection, and amplification provides the system with novelty and is responsible for its growth in order and complexity | No mechanism for endogenously creating novelty, or growth in order and complexity |
Table 2.1, Mirowski (1996), p.15
- Romanticism wholist, historicist, intedeterminist, emotivist, experimental, diversity, aesthetic detachment, organicism, perspectival
- Classicism atomist, atemporal, determinist, stoic, formal, uniformity, interventionist control, mechanicism, aperspectival
Causal loop diagram of a model examining the growth or decline of a life insurance company
Basics of Dynamical Systems
Definitions
In many sciences dynamical systems play a central role in the description of reality. The basic idea is to represent the state of a system at a particular time \(t\) in terms of measurements, i.e. different variables or dimensions, which is called the state vector \(x(t)\). The space in which the state vector is defined (usually a finite-dimensional Euclidean space) is called the state space. The basic idea of a dynamical system is that certain laws of motion govern the evolution of the state vector.Hence we can define a dynamical system as a state space \(S\), a set of times \(T\), and a rule \(R\) that regulates the temporal evolution of the system: \[ R: S \times T \to S \]
- State space \(S\) can be continuous or discrete (e.g. cellular automata or Markov chains)
- Rule \(R\) can be deterministic (including deterministic chaos) or stochastic (e.g. Markov chains, but also the FitzHugh-Nagumo model of a spiking neuron). Rules \(R\) are often characterized by a set of ordinary differential equations that deterministically relate the current state of a system to a future state of the system.
- Discrete Maps
- Dynamical systems that operate over discrete time and a continuous state space, of the form \[f: S \to S\]. If \(x_t\) describes the system at time \(t\), a discrete map is characterized by \[x_{t+1} = f(x_t)\]
The most famous example is the logistic map: \[ x_{t+1}= rx (1 - \frac{x}{K}) \]
- Continuous-Time Flows
- Deterministic dynamical systems on a manifold \(M\) (subset of Euclidean space) \[\phi: R \times M \to M\], with an orbit \(x(t) = \phi_t (x(0))\). They are generally written as a set of differential equations, \[\dot{x} = \phi(x)\]. Flows can generate cycles in two dimensions and chaos in three dimensions. Turbulence can also be described by a dynamical system operating on a complex manifold. ODEs are deterministic, which implies that orbits cannot intersect themselves.
- Equilibrium
- \(x^*\) is the equilibrium or fixed point of a dynamical system if: \[x^* = f(x^*) \quad \text{in the discrete case}\] \[f(x^*) = 0 \quad \text{in the continuous case}\]
- Kinds of Fixed Points
- Fixed points can be:
- a sink if orbits lead to it
- a source if orbits run away from it,
- a center if they are center of stable cycles, or
- a saddle if orbits in one direction lead to it while they run away from it in the other.
- Attractor
- Roughly speaking, a subset of the state space to which orbits originating from typical initial conditions tend as time increases: e.g. “butterfly-effect” Lorenz Attractor, Rossler Attractor.
- Basin of Attraction
- Set of initial conditions leading to long-time behavior that approaches that attractor. The qualitative behavior of the long-time motion of a given system can be fundamentally different depending on which basin of attraction the initial condition lies in (e.g. attractors can correspond to periodic, quasiperiodic or chaotic behaviors of different types).
Dynamical Behavior
Linear Systems
In vectorial form, \[ \begin{pmatrix} x_1(t+1) \\ \ldots \\ x_N(t+1) \\ \end{pmatrix} = A \begin{pmatrix} x_1(t) \\ \ldots \\ x_N(t) \\ \end{pmatrix} \] In matricial form, \[ x_{t+1} = A x_t \] In ‘most’ cases, the matrix \(A\) has an eigenspace, characterized by \(N\) eigenvalues \(\lambda_1,\ldots,\lambda_N\) with corresponding eigenvectors \(w_1,\ldots,w_N\): \[ A w_i = \lambda_i w_i, \qquad i = 1,\ldots,N \] Computing \[ \det (A - \lambda I) =0 \] produces the characteristic equation with eigenvalues as their roots.
A general analytical solution for the linear discrete-time system is: \[ \begin{pmatrix} x_1(t+1) \\ \ldots \\ x_N(t+1) \\ \end{pmatrix} = b_1 \lambda^t_1 w_1 + \ldots + b_N \lambda^t_N w_N \]
In continuous time, nontrivial solutions will be of the form \[ x(t) = e^{\lambda t} u \] where \(u\) is the (eigen)vector of real or complex constants and \(\lambda\) is the real or complex eigenvalue.
Stability
The fixed point is:
- a source, i.e. is unstable if the eigenvalues are real and positive (above one module in the discrete-time case),
- a sink, i.e. is stable if the eigenvalues are real and negative (below one module in the discrete-time case),
- a saddle-point if one eigenvalue is real positive and the other is real negative,
- a center producing oscillations if the eigenvalues are complex conjugates with zero real component. The non-negative real component will add amplification or dampening to the oscillations
In the long term, i.e. for \(t \to \infty\), there are two qualitatively distinct cases: convergence to or oscillations around the fixed point.
- Convergence
- Let \(\lambda\) be the eigenvalue of \(A\) with the largest absolute value, i.e. the dominant eigenvalue, and let \(w\) be the corresponding eigenvector. Then, if \(\lambda\) is a real number, the long-term behaviour is given by: \[ \begin{pmatrix} x_1(t) \\ \ldots \\ x_N(t) \\ \end{pmatrix} = b \lambda^t w \qquad \text{for } t \to \infty \] where \(bw\) is the component of the initial condition in the direction of the eigenvector \(w\). If \(\lambda\) is real positive (negative), the dynamic variables will eventually grow (shrink) at a rate \(\lambda\), and the vector of variables will be a multiple of \(w\), the dominant eigenvector indicating the balanced-growth path of the system. In economics, the balanced-growth path (Von Neumann 1945) is called equilibrium and the dominant eigenvector \(w\) plays the role to the Sraffian `standard commodity’ (Sraffa 1960), while the dominant eigenvalue corresponds to the maximum rate of expansion (Shaikh 2016).
- Oscillations
- If \(\lambda\) is complex, \(\lambda_+ = r + is = c(\cos \phi + i \sin \phi)\) with \(\phi \neq 0\), then \(w = u + iv\) with \(v \neq 0\), and the complex conjugate \(\lambda_- r - is\) is also an eigenvalue of \(A\) with corresponding eigenvector \(w = u - iv\). The real-valued solution will eventually converge to oscillating behavior with exponentially increasing (or decreasing) amplitude: \[ \begin{pmatrix} x_1(t) \\ \ldots \\ x_N(t) \\ \end{pmatrix} = c^t \left[ b_1 \left(\cos \phi t \cdot u - \sin \phi t \cdot v \right) + b_2 \left( b_1 \cos \phi t \cdot u - \sin \phi t \cdot v \right) \right] \quad \text{for } t \to \infty \] where \(c = \sqrt{r^2 + s^2}\) is the absolute value of \(\lambda\) and \(b_1\) and \(b_2\) are determined from the initial condition.
Nonlinear Systems
This procedure is very helpful to compute the stability of fixed points of nonlinear systems. In order to do so, we compute the Jacobian matrix of the system evaluated at its fixed points. For instance, in two dimensions we have, for a continuous-time system: \[ \begin{cases} \dot{x} = f_1(x,y) \\ \dot{y} = f_2(x,y) \end{cases} \] The Jacobian is \[ J_{(x^*,y^*)} (f) = \begin{pmatrix} \frac{\partial f_1}{\partial x} & \frac{\partial f_1}{\partial y} \\ \frac{\partial f_2}{\partial x} & \frac{\partial f_2}{\partial y} \end{pmatrix} \] Computing \[ \det (J - \lambda I) =0 \] produces the characteristic equation with eigenvalues \(\lambda\) as their roots.
Here you have a practical example including a phase portrait.
Trace-Determinant Method
We can write the determinant of the two-dimensional system as \[ \det(A - \lambda I) = \lambda^2 - \underbrace{(a+d)}_{\text{trace}} \lambda + \underbrace{(ad-bc)}_{\text{determinant}} = \lambda^2 + T \lambda + D \] The eigenvalues of \(A\) are given by \[ \lambda = \frac{T \pm \sqrt{T^2 - 4D}}{2} \] - If \(T^2 - 4D > 0\) we have two distinct real eigenvalues - If \(T^2 - 4D < 0\) we have two complex eigenvalues that are conjugate - If \(T^2 - 4D = 0\) we have repeated eigenvalues
Classification of Phase Portraits in Two Dimensions in terms of the trace and the determinant of A
Bifurcation Theory
- Definition
- A bifurcation occurs when a small smooth quantitative change made to the parameter values (the bifurcation parameters) of a system causes a sudden ``qualitative” or topological change in its behavior. Bifurcations occur in both continuous systems (described by ODEs, DDEs or PDEs) and discrete systems (described by maps). Saddle-node bifurcations may be associated with hysteresis loops and catastrophes Rosser (2000).
- Saddle-Node
- In a saddle-node bifurcation, two fixed points (or equilibria) of a dynamical system collide and annihilate each other. An example is \[\dot{x} = r + x^2,\] with fixed point \(x^* = \sqrt{-r}\). If \(r<0\), there are two equilibrium points, a stable point at \(-\sqrt{-r}\) and an unstable point at \(\sqrt{-r}\). At \(r=0\) there is the bifurcation with one equilibrium point with saddle-node stability. If \(r>0\) there are no equilibrium points.
- Pitchfork
- A pitchfork bifurcation is a particular type of local bifurcation where the system transitions from one fixed point to three fixed points. Pitchfork bifurcations can be supercritical, where a stable fixed point becomes unstable while ``giving birth” to two stable fixed points. Example: \[\dot{x} = rx - x^3; \quad x^* = \begin{cases} 0 \\ \pm \sqrt{r} \end{cases}\] In the subcritical case, an unstable fixed point becomes stable while giving birth to two unstable fixed points. Example: \[\dot{x} = rx + x^3; \quad x^* = \begin{cases} 0 \\ \pm \sqrt{-r} \end{cases}\]
- Transcritical
- A transcritical bifurcation is one in which a fixed point exists for all values of a parameter and is never destroyed. However, such a fixed point interchanges its stability with another fixed point as the parameter is varied. In other words, both before and after the bifurcation, there is one unstable and one stable fixed point. However, their stability is exchanged when they collide. So the unstable fixed point becomes stable and vice versa. The normal form is \[\dot{x} = r x - x^2\]
- Hopf
- A Hopf bifurcation is a critical point where a system’s stability switches and a periodic solution (i.e. oscillations) arises.
Hopf bifurcations occur in the Lotka–Volterra model of predator–prey interaction (known as paradox of enrichment), the Hodgkin–Huxley model for nerve membrane, the Selkov model of glycolysis, the Belousov–Zhabotinsky reaction, the Lorenz attractor, and the Brusselator.
- Period-Doubling
- A period doubling bifurcation in a discrete dynamical system is a bifurcation in which a slight change in a parameter value in the system’s equations leads to the system switching to a new behavior with twice the period of the original system. With the doubled period, it takes twice as many iterations as before for the numerical values visited by the system to repeat themselves. A period doubling cascade is a sequence of doublings and further doublings of the repeating period, as the parameter is adjusted further and further. The most famous example is the logistic map, where period doubling leads to chaos in the form of aperiodic oscillations.
Exponential Growth in a Limited-Resource World
| Model | Dimensions | Topic |
|---|---|---|
| Bhaduri-Harris | 1 | Complex Dynamics of the Simple Ricardian System |
Predator-Prey Oscillations with 2+ Dimensions
| Model | Dimensions | Topic | Prey | Predator |
|---|---|---|---|---|
| Goodwin | 2 | Distribution | employment rate | labor share of income |
| Flaschel-Semmler | \(2N\) | Growth | prices/profits | quantities/capital |
R Setup
Check first if our version of R is the latest, otherwise
download it from the Comprehensive R Archive Network (CRAN) in https://www.r-project.org/ >> Download
>> CRAN >> Select your mirror.
Check for useful R Cheatsheets in https://www.rstudio.com/resources/cheatsheets/.
# Download and install RStudio
# Also latest version of R from https://ftp.fau.de/cran/
version## _
## platform x86_64-apple-darwin17.0
## arch x86_64
## os darwin17.0
## system x86_64, darwin17.0
## status
## major 4
## minor 1.3
## year 2022
## month 03
## day 10
## svn rev 81868
## language R
## version.string R version 4.1.3 (2022-03-10)
## nickname One Push-UpWe are using Markdown for our notes: check https://www.markdownguide.org/ for its basic and
extended syntax and https://bookdown.org/yihui/rmarkdown/, the R Markdown:
The Definitive Guide by Yihui Xie.
When writing code, remember to be nice to your future self (and others)!
Logistic Growth
\[ x_{t+1} = r x_t (1 - x_t) \] where \(x_t\) is a number between zero and one that represents the ratio of existing population to the maximum possible population \(K\). The values of interest for the parameter \(r\) are those in the interval [0,4]. This nonlinear difference equation is intended to capture two effects:
- reproduction where the population will increase at a rate proportional to the current population when the population size is small.
- density-dependent starvation where the growth rate will decrease at a rate proportional to the value obtained by taking the theoretical “carrying capacity” \(K\) of the environment less the current population.
Fixed Points
The fixed points are found by computing the solutions of \(f(x)=x\): \[x = rx - rx^2\] \[r x^2 + (1-r)x = x(rx+1-r) = 0\] There are two fixed points, but only one is non-trivial:
\[ x^* = \begin{cases} 0 \\ \frac{r-1}{r} = 1 - \frac{1}{r} \end{cases} \]
Stability Analysis
In order to derive the stability of the fixed points, we take the derivative of \(f(x) = rx(1-x)\) and evaluate it at the fixed points. The derivative of \(f(x)\) is \[f'(x) = r (1-2x)\] evaluated at \(x^*\): \[f'(0)=r\] Hence, \(x^*_1 = 0\) will be stable if \(r < 1\) and unstable if \(r > 1\). \[f'(\frac{r-1}{r}) = 2-r\] The second fixed point \(x^*_2 = 1 - 1/r\) will be unstable if \(r<1\), stable if \(1<r<3\), and unstable again if \(r>3\).
# LOGISTIC MAP
x <- c()
x[1] <- 0.05 # initial condition
r <- 3.5 # parameter
for (t in 1 : 100) x[t+1] <- r * x[t] * (1 - x[t])
plot(x, type='l', lwd=2, xlab = 'time', ylab = 'population')
# We can write a function logistic
# Three local variables: initial condition; reproduction rate (r), and number of iterations
logistic <- function (init, par, time.length)
{
x <- c()
x[1] <- init
for (t in 1 : time.length) x[t+1] <- par * x[t] * (1 - x[t])
return(x)
}
plot(logistic(init = 0.25,
par = 3.56995,
time.length = 20),
type='l', lwd=2, xlab='time', ylab = 'population')
grid()
Dynamical Behavior
By varying the parameter \(r\), the following behavior is observed.
Starvation
With \(r\) between 0 and 1, the population will eventually die, independent of the initial population.
library(colorspace)
tf <- 15
plot(NA, xlim=c(0,tf), ylim=c(0,1), xlab = 'time (t)', ylab = 'Population/Capacity Ratio (x)')
grid()
rval <- seq(0.1, 1, by = 0.25)
t_0 <- c(0.5, .8, .25, .6)
for(k in 1:length(rval))
{
lines(1:(tf+1), logistic(t_0[k], rval[k], tf), col = rainbow_hcl(n = length(rval))[k], lwd=2)
}
legend('topright', paste('r:', rval), pch = 16, col = rainbow_hcl(n = length(rval)))
Logistic Growth
With \(r\) between 1 and 2, the population will quickly approach the value \(\frac{r-1}{r}\), independent of the initial population.
plot(NA, xlim = c(0, tf), ylim = c(0,1), xlab='time (t)', ylab='Population/Capacity Ratio (x)')
grid()
rval <- seq(1.1, 2, by = 0.25)
t_0 <- c(0.01, .8, .25, .01)
for(k in 1 : length(rval))
{
lines(1:(tf+1), logistic(t_0[k], rval[k], tf), col = rainbow_hcl(n = length(rval))[k], lwd=2)
}
mtext(paste(' r*:',rval), side = 4, at = (rval - 1)/rval, col = rainbow_hcl(n = length(rval)), las=1, cex=0.8)
legend('topleft', paste('r:', rval), pch=16, col = rainbow_hcl(n = length(rval)), bty='n')
Convergence to Stability with Dampening Oscillations
With \(r\) between 2 and 3, the population will also eventually approach the same value \(\frac{r-1}{r}\), but first will fluctuate around that value for some time.
plot(NA, xlim = c(0, tf), ylim = c(0,1), xlab='time (t)', ylab='Population/Capacity Ratio (x)')
grid()
rval<-seq(2.1,3,by=0.1)
for(k in 1 : length(rval))
{
lines(1:(tf+1), logistic(0.8, rval[k], tf), col = rainbow_hcl(n = length(rval))[k], lwd=2)
}
mtext(paste(' r*:',rval), side = 4, at = (rval - 1)/rval, col = rainbow_hcl(n = length(rval)), las=1, cex=0.8)
legend('topleft', paste('r:', rval), pch=16, col = rainbow_hcl(n = length(rval)), bty='n')
Periodic Oscillations
With \(r\) between 3 and \(1 + \sqrt{6} \sim 3.44949\), from almost all initial conditions the population will approach permanent oscillations between two values. These two values are dependent on \(r\).
plot(NA, xlim = c(0, tf), ylim = c(.4 , .9), xlab='time (t)', ylab='Population/Capacity Ratio (x)')
grid()
rval <- seq(3, 3.45, by = 0.05)
for(k in 1:length(rval))
{
lines(1:(tf+1), logistic(0.5, rval[k], tf), col = rainbow_hcl(n = length(rval))[k], lwd=2)
}
legend('topleft',paste('r:',rval), pch = 16, col = rainbow_hcl(n = length(rval)), bty = 'n')
With \(r\) between 3.44949 and 3.54409 (approximately), from almost all initial conditions the population will approach permanent oscillations among four values.
plot(NA, xlim = c(0, tf), ylim = c(0, .9), xlab='time (t)', ylab='Population/Capacity Ratio (x)')
grid()
rval<-seq(3.45, 3.65, by = 0.025)
for(k in 1:length(rval))
{
lines(1:(tf+1), logistic(0.5,rval[k],tf), col = rainbow_hcl(n = length(rval))[k], lwd=2)
}
legend('topleft', paste('r:',round(rval,2)),pch = 16,col = rainbow_hcl(n = length(rval)),bty = 'n')
Onset of Chaos: the Period-Doubling Cascade
With \(r\) increasing beyond 3.54409, from almost all initial conditions the population will approach oscillations among 8 values, then 16, 32, etc. This behavior is an example of a period-doubling cascade.
tf <- 50
plot(NA, xlim = c(0,tf), ylim=c(0,1), xlab = 'time (t)',ylab = 'Population/Capacity Ratio (x)')
grid()
# rval <- seq(3.65 , 3.8, by = 0.025)
rval <- rep(3.9, 4)
lines(1:(tf+1), logistic(0.499, rval[1], tf), col = rainbow_hcl(n=length(rval))[1], lwd = 2)
lines(1:(tf+1), logistic(0.501, rval[2], tf), col = rainbow_hcl(n=length(rval))[2], lwd = 2)
lines(1:(tf+1), logistic(0.51, rval[3], tf), col = rainbow_hcl(n=length(rval))[3], lwd = 2)
lines(1:(tf+1), logistic(0.4995, rval[4], tf), col = rainbow_hcl(n=length(rval))[4], lwd = 2)
legend('topleft', paste('r:',round(rval,2)), pch=16, col = rainbow_hcl(n=length(rval)), bty='n')
At \(r \sim 3.56995\) is the onset of chaos, at the end of the period-doubling cascade. From almost all initial conditions, we no longer see oscillations of finite period.
plot(NA, xlim = c(0,4), ylim = c(0,1), xlab = 'Parameter r', ylab = 'Stable Population Values', main = 'Bifurcation Graph')
grid(lwd = 1.5, nx=7)
for(r in seq(0.1,4,by=.005))
{
yy <- as.numeric(names(table(logistic(0.5, r, 100)[50:100])))
points(rep(r, times = length(yy)), yy, pch = 16, cex = .1)
}
Iteration Maps
plot(NA, xlim = c(0, 1), ylim = c(0, 1), xlab = 'x[t]', ylab = 'x[t+1]', main = 'Iteration Map')
grid(lwd = 1.5)
library(scales)
# lines(rep(tmp,each=2)[-c(1,202)],rep(tmp,each=2)[-c(1,2)],col=rainbow_hcl(r0[k]))
r0 <- seq(3, .5, by = -.5)
for(k in 1 : length(r0))
{
tmp <- logistic(runif(1), r0[k], 100)
lines(rep(tmp, each = 2)[-c(1,202)], rep(tmp,each=2)[-c(1,2)], col = alpha(rainbow_hcl(length(r0))[k],0.8), lwd = 3)
lines(seq(0,1,by=.01),r0[k]*seq(0,1,by=.01)*(1-seq(0,1,by=.01)),col=rainbow_hcl(length(r0))[k],lwd=2)
points(tmp[1], tmp[2])
points(tmp[100], tmp[101], col = rainbow_hcl(length(r0))[k])
}
abline(a=0,b=1,lty=2)
legend('topleft',paste('r:',round(r0,2)),pch=16,col=rainbow_hcl(n=length(r0)),bty='n')
plot(NA, xlim = c(0, 1), ylim = c(0, 1), xlab = 'x[t]', ylab = 'x[t+1]', main = 'Iteration Map')
grid(lwd = 1.5)
# lines(rep(tmp,each=2)[-c(1,202)],rep(tmp,each=2)[-c(1,2)],col=rainbow_hcl(r0[k]))
r0 <- seq(3, 3.5, by = .1)
for(k in 1:length(r0))
{
tmp <- logistic( runif(1), r0[k], 100)
lines(rep(tmp, each = 2)[-c(1,202)], rep(tmp, each=2)[-c(1,2)], col = alpha(rainbow_hcl(length(r0))[k],0.6), lwd=3)
lines(seq(0,1,by=.01), r0[k]*seq(0,1,by=.01)*(1-seq(0,1,by=.01)), col = rainbow_hcl(length(r0))[k], lwd=2)
}
abline(a=0,b=1,lty=2)
legend('topleft',paste('r:',round(r0,2)),pch=16,col=rainbow_hcl(n=length(r0)),bty='n')
plot(NA,xlim=c(0,1),ylim=c(0,1),xlab='x[t]',ylab='x[t+1]',main='Iteration Map')
grid(lwd=1.5)
# lines(rep(tmp,each=2)[-c(1,202)],rep(tmp,each=2)[-c(1,2)],col=rainbow_hcl(r0[k]))
r0<-seq(3.5,4,by=.1)
for(k in 1:length(r0))
{
tmp<-logistic(runif(1),r0[k],100)
lines(rep(tmp,each=2)[-c(1,202)],rep(tmp,each=2)[-c(1,2)],col=alpha(rainbow_hcl(length(r0))[k],0.6),lwd=3)
lines(seq(0,1,by=.01),r0[k]*seq(0,1,by=.01)*(1-seq(0,1,by=.01)),col=rainbow_hcl(length(r0))[k],lwd=2)
}
abline(a=0,b=1,lty=2)
legend('topleft',paste('r:',round(r0,2)),pch=16,col=rainbow_hcl(n=length(r0)),bty='n')
The Complex Dynamics of the Simple Ricardian System
Visual depiction of the Ricardian one-sector model following Pasinetti
The idea that the process of capitalist accumulation ultimately terminates in a stationary state is deeply rooted in the classical tradition of political economy. Ricardo, in particular, visualized the process of accumulation that is almost inexorably driven toward its end because the capitalists gradually lose command over the investible surplus as their profit dwindles due to increasing pressure of rent on limited natural resources like land. Bhaduri and Harris (1987)
Consider an economy where labor \(N\) is applied in fixed proportion to less and less fertile land as the margin of cultivation is extended in the process of capital accumulation. Accordingly, there exists a falling marginal product curve for labor with fixed ``doses” of land. This curve is assumed to be linear of the form: \[ MP_t = \frac{dY_t}{dN_t} = a - bN_t \quad a>0, \quad b>0 \] By integration, the total product curve is: \[ Y_t = aN_t - bN^2/2 \] The average product is: \[ AP_t = \frac{Y_t}{N_t} = a - bN_t/2 \] Given the margin of cultivation reach at time \(t\), the level of land rent \(R_t\) emerges as the difference between the average and the marginal product of labor, \[ R_t = \left( \frac{Y_t}{N_t} - \frac{dY}{dN_t} \right) N_t = bN^2_t/2 \] Profit \(P_t\) is the residual after payment of rent and replacement of the wage fund \(W_t\) advanced to hire labor. Thus \[ P_t = Y_t - R_t - W_t \] It is the dynamics of the wage fund that represents the process of accumulation in this simple economy. At a given wage rate \(w\) we have \[ W_t = wN_t \] Classical condition imposes the dynamical behavior Accumulation of the wage fund comes entirely from reinvestment of profits accruing to capitalists: \[ W_{t+1}-W_t = P_t \] which yields \[ N_{t+1}-N_t = \frac{P_t}{w} = \frac{a}{w} N_t - \frac{b}{w} N^2_t - N_t \] Finally we obtain the discrete dynamic equation: \[ N_{t+1}= \frac{a}{w} N_t - \frac{b}{w} N^2_t \] There exists an equilibrium at \[ N_{t+1} = N_t = N^* = \frac{a-w}{b} \] In such an equilibrium, the marginal product of labor \(MP^*\) is equal to the given real wage rate \(w\) so that all profit disappears \(P^*=0\) and the whole surplus over the wage bill accrues to the landlords as rent – the Ricardian stationary state.
Since the positivity of the marginal product requires \(a > bN_t\), corresponding to all economically meaningful employment levels \(N_t\), we can define a new variable \(n_t\) such that \[ n_t \equiv \frac{bN_t}{a} \] which lies between 0 and 1. Consequently, the dynamic equation can be rewritten as: \[ n_{t+1} = A n_t (1-n_t) \] where \(A = a/w\). This is the equation for the logistic map. It has two fixed points, \(n^*=0\) and \(n^* = \frac{a-w}{a}\).Economic Interpretation
The degree of complexity of the dynamic behavior of the simple Ricardian system depends critically on the value of the ratio \(A\). If the system is to be capable of generating any surplus over wages, then the production coefficient \(a\) must exceed the wage rate \(w\) and thus \(A>1\) for the economy to be viable.Rewriting \(a/w\) as \[ \frac{a-w+w}{w}=1+ \frac{a-w}{w}= 1 + \epsilon_0 \] where \(\epsilon_0 \equiv (a-w)/w\).
When this rate of exploitation exceeds a relatively high critical value, that is \(\epsilon_0>3\), that the simple Ricardian system generates the dynamics resembling chaos.\(\epsilon_0\) corresponds to the rate of exploitation measured in corn at the very initial state of the system, i.e. at ``primitive accumulation”. It is also the maximum rate of exploitation that the system supports.
The dynamic behavior of the system is governed in a remarkably simple way by the distributional relations.
The value of the crucial ratio \(A\) can be expressed in terms of the distribution of income between rent and wages in the Ricardian stationary states (when profit is zero): \[ \frac{R^*}{W^*} = \frac{1}{2} \left(\frac{a-w}{w}\right) = \frac{1}{2} (A-1) \]
Predator-Prey Oscillations with 2+ Dimensions
| Model | Dimensions | Topic | Prey | Predator |
|---|---|---|---|---|
| Goodwin | 2 | Distribution | employment rate | labor share of income |
| Flaschel-Semmler | \(2N\) | Growth | prices/profits | quantities/capital |
Lotka-Volterra / Predator-Prey Oscillations
The Lotka-Volterra equations are frequently used to describe the dynamics of biological systems in which two species interact.
The Lotka-Volterra equations, also known as the predator-prey equations, are a pair of first-order nonlinear differential equations, frequently used to describe the dynamics of biological systems in which two species interact, one as a predator and the other as prey. Originally proposed in the theory of autocatalytic chemical reactions.
Brusselator
The Brusselator is a theoretical model for a type of autocatalytic reaction. The Brusselator model was proposed by Ilya Prigogine and collaborators at the Université Libre de Bruxelles.
It is characterized by the reactions: \[ A \to X \] \[ 2 X + Y \to 3 X\] \[ B + X \to Y + D \] \[ X \to E \] Under conditions where \(A\) and \(B\) are in vast excess and can thus be modeled at constant concentration, the rate equations become \[ \frac{d}{dt} \{X\} = \{A\} + \{X\}^2 \{Y\} - \{B\}\{X\} - \{X\} \] \[ \frac{d}{dt} \{Y\} = \{B\} \{X\} - \{X\}^2 \{Y\} \] \[ \dot{x} = a + x^2 y - (1+b) \] \[ \dot{y} = bx - x^2 y \] where, for convenience, the rate constants have been set to 1.
The Brusselator has a fixed point at \[ \{ X \} = A; \qquad x^* = a \] \[ \{ Y \} = \frac{B}{A} \qquad y^* = \frac{b}{a}\]
The fixed point becomes unstable when \[ B > 1 + A^2 \] leading to an oscillation of the system. Unlike the Lotka–Volterra equation, the oscillations of the Brusselator do not depend on the amount of reactant present initially. Instead, after sufficient time, the oscillations approach a limit cycle.
Spatial Brusselator chemicall reactions
Model
The populations change through time according to the pair of equations: \[ \begin{cases} \dot{x} = \alpha x - \beta xy \\ \dot{y} = \delta x y - \gamma y \end{cases} \] \[ \frac{\dot{x}}{x} = \alpha - \beta y \] \[ \frac{\dot{y}}{y} = \delta x - \gamma \]
where \(x\) and \(y\) are the prey and predator populations, respectively; \(\dot{x}\) and \(\dot{y}\) are their growth rates; and \(\alpha\), \(\beta\), \(\delta\), and \(\gamma\) are parameters describing the interactions between the species.Assumptions
- The prey population finds ample food at all times.
- The food supply of the predator population depends entirely on the size of the prey population.
- The rate of change of population is proportional to its size.
- During the process, the environment does not change in favour of one species, and genetic adaptation is inconsequential.
- Predators have limitless appetite.
Dynamical Behavior
The Lotka-Volterra model has two fixed points: an extinction equilibrium and a characteristic limit cycle (with a shape that depends on its four parameters and oscillation amplitude that depends on the initial values for \(x\) and \(y\)).
Population equilibrium occurs in the model when neither of the population levels is changing, i.e. when both of the derivatives are equal to 0: \[x(\alpha - \beta y) = 0\] \[-y (\gamma - \delta x) =0\] The fixed points \((x^*,y^*)\) are thus \[x^*_1 = (0,0) \quad \text{extinction}\] and \[x^*_2=(\frac{\gamma}{\delta},\frac{\alpha}{\beta}) \quad \text{oscillations}\]
The stability of the fixed points can be determined by performing a linearization using partial derivatives. The Jacobian matrix of the predator–prey model is \[J(x,y) = \begin{pmatrix} \alpha - \beta y & - \beta x \\ \delta y & \delta x - \gamma \end{pmatrix}\]
Thus, for the first fixed point at the origin (extinction): \[J(0,0) = \begin{pmatrix} \alpha & 0 \\ 0 & \gamma \end{pmatrix}\] with eigenvalues \(\lambda_1=\alpha\) and \(\lambda_2=-\gamma\).
In the model \(\alpha\) and \(\gamma\) are always greater than zero, and as such the sign of the eigenvalues above will always differ. Hence the fixed point at the origin is a saddle point.
The stability of this fixed point is of significance. If it were stable, non-zero populations might be attracted towards it, and as such the dynamics of the system might lead towards the extinction of both species for many cases of initial population levels. However, as the fixed point at the origin is a saddle point, and hence unstable, it follows that the extinction of both species is difficult in the model. (In fact, this could only occur if the prey were artificially completely eradicated, causing the predators to die of starvation. If the predators were eradicated, the prey population would grow without bound in this simple model.) The populations of prey and predator can get infinitesimally close to zero and still recover.
For the second fixed point (oscillations), the Jacobian yields \[J(\frac{\gamma}{\delta},\frac{\alpha}{\beta}) = \begin{pmatrix} 0 & -\frac{\beta \gamma}{\delta} \\ \frac{\alpha \delta}{\beta} & 0 \end{pmatrix}\] with conjugate pure imaginary eigenvalues \[\lambda = \pm i \sqrt{\alpha \gamma}\]. As the eigenvalues are both purely imaginary and conjugate to each others, this fixed point is elliptic, so the solutions are periodic, oscillating on a small ellipse around the fixed point, with frequency \[\omega = \sqrt{\lambda_1 \lambda_2} = \sqrt{\alpha \gamma}\]. As illustrated in the circulating oscillations in the figure above, the level curves are closed orbits surrounding the fixed point: the levels of the predator and prey populations cycle and oscillate without damping around the fixed point.
Lotka-Volterra oscillations for two sets of parameters
Phase-space plot for the predator prey problem for various initial conditions of the predator population
Computation
The Lotka-Volterra model is in continuous time. First, we attempt to compute the problem as if it was in discrete time. In order to do so, we write the equations as: \[x_{t+1} = x_t + \Delta t f_x(x,y)\] \[y_{t+1} = y_t + \Delta t f_y(x,y)\]
alpha <- 1
beta <- 1
gamma <- 1
delta <- 1
x <- y <- c() # initialization
x[1] <- 1 # initial conditions
y[1] <- 1.2
orbit <- data.frame(x = rep(x, 10), y = rep(y, 10))
foo <- list(x, y, orbit)
interval <- .01
timelength <- 50
for(t in 1 : (timelength/interval))
{
x[t+1] <- x[t] + interval * (alpha * x[t] - beta * x[t] * y[t])
y[t+1] <- y[t] + interval * (delta * x[t] * y[t] - gamma * y[t])
}
lotka.volterra <- function (pars, time.length, interval, init)
{
x <- y <- c() # initialization
x[1] <- init[1]
y[1] <- init[2]
for(t in 1 : (time.length/interval))
{
x[t+1] <- x[t] + interval * (pars[1] * x[t] - pars[2] * x[t] * y[t])
y[t+1] <- y[t] + interval * (pars[3] * x[t] * y[t] - pars[4] * y[t])
}
orbit <- data.frame(x = x, y = y)
return(orbit)
}
parameters <- c(1,1,3,1)
orbit <- lotka.volterra(pars = parameters, time.length = 10, interval = 0.0001, init = c(1, 1.2))
timelength <- nrow(orbit)
interval <- 0.01
plot(orbit$x, type='l' , xlim = c(0, timelength), ylim=c(0, max(orbit)), xlab='time', ylab='population', col='steelblue', lwd=2)
lines(orbit$y, col='tomato', lwd=2)
alpha <- parameters[1]
beta <- parameters[2]
gamma <- parameters[3]
delta <- parameters[4]
abline(h = delta/gamma, lty=2, col = 'steelblue', lwd = 2) # fixed point for prey
abline(h = alpha/beta, lty=2, col = 'tomato', lwd = 2) # fixed point for predator
legend('topleft',c('prey','predator'), col = c('steelblue','tomato'), pch=16, bty='n')
# Phase Space
plot(orbit$x, orbit$y, type='l')
points(delta/gamma,alpha/beta)
Let’s try with another set of parameters and initial conditions:
alpha<-1.5
beta<-1
gamma<-3
delta<-1
x<-y<-c()
x[1]<-0.5 # initial conditions
y[1]<-1.2
interval<-.01
timelength<-50
for(t in 1:(timelength/interval))
{
x[t+1] <- x[t]+interval*(alpha*x[t] - beta*x[t]*y[t])
y[t+1] <- y[t]+interval*(delta*x[t]*y[t] - gamma*y[t])
}
plot(seq(0,timelength,by=interval),x,type='l',xlim=c(0,timelength),ylim=c(0,20),xlab='time',ylab='population',col='steelblue',lwd=2)
lines(seq(0,timelength,by=interval),y,col='tomato',lwd=2)
abline(h=delta/gamma,lty=2,col='steelblue',lwd=2) # fixed point x
abline(h=alpha/beta,lty=2,col='tomato',lwd=2) # fixed point y
legend('topleft',c('prey','predator'),col=c('steelblue','tomato'),pch=16,bty='n')
# Phase Space
plot(x,y,type='l')
points(delta/gamma,alpha/beta)
We are not observing closed orbits and stable oscillations. How so?
The problem is discretization, which amplifies the oscillations
indefinitely by adding ‘energy’ to the system at each step. In order to
address this issue, we need to compute the solution of the system of
continuous differential equations. Fortunately, there is already a
R package for that called deSolve, which has
the function ode.
library(deSolve)
# We define the function to give to ode
LotVmod <- function (Time, State, Pars) {
with(as.list(c(State, Pars)), {
dx = x*(alpha - beta*y)
dy = -y*(gamma - delta*x)
return(list(c(dx, dy)))
})
}
# Parameters of the Model
Pars <- c(alpha = 2, beta = .5, gamma = .2, delta = .6)
State <- c(x = 10, y = 10)
Time <- seq(0, 100, by = .01)
# Compute the Solution and Plot it
solution <- as.data.frame(ode(func = LotVmod, y = State, parms = Pars, times = Time))
plot(solution$time,solution$x,xlab = "time", ylab = "population",col='steelblue',lwd=2,type='l',ylim=c(0,max(solution[,c('x','y')])))
lines(solution$time,solution$y,col='tomato',lwd=2)
abline(h=Pars['delta']/Pars['gamma'],lty=2,col='steelblue',lwd=2) # fixed point x
abline(h=Pars['alpha']/Pars['beta'],lty=2,col='tomato',lwd=2) # fixed point y
legend("topright", c("prey","predator"),col=c('steelblue','tomato'),pch=16,bty='n')
# Phase Space
plot(solution$x,solution$y,type='l')
points(Pars['delta']/Pars['gamma'],Pars['alpha']/Pars['beta']) # fixed point
Let’s try the former set of parameters now.
Pars <- c(alpha = 1, beta = 1, gamma = 1, delta = 1)
State <- c(x = 1, y = 1.2)
Time <- seq(0, timelength, by = interval)
# Compute the Solution and Plot it
solution <- as.data.frame(ode(func = LotVmod, y = State, parms = Pars, times = Time))
plot(solution$time,solution$x,xlab = "time", ylab = "population",col='steelblue',lwd=2,type='l',ylim=c(0,max(solution[,c('x','y')])))
lines(solution$time,solution$y,col='tomato',lwd=2)
abline(h=Pars['delta']/Pars['gamma'],lty=2,col='steelblue',lwd=2) # fixed point x
abline(h=Pars['alpha']/Pars['beta'],lty=2,col='tomato',lwd=2) # fixed point y
abline(v=seq(0,timelength,by=2*pi),lty=2) # the period is 2*pi/frequency, where frequency is 1
legend("bottomright", c("prey","predator"),col=c('steelblue','tomato'),pch=16,bty='n')
# Phase Space
plot(solution$x,solution$y,type='l')
points(Pars['delta']/Pars['gamma'],Pars['alpha']/Pars['beta']) # fixed point
Let’s try random parameters for 4 random initial conditions:
Pars <- c(alpha = 3, beta = 1.2, gamma = 1.3, delta = 0.67)
timelength<-20
Time <- seq(0, timelength, by = interval)
# Compute the Solution and Plot it
solution <- list()
for(i in 1 : 4)
{
State <- c(x = runif(1,0,10), y = runif(1,0,10))
solution[[i]]<-as.data.frame(ode(func = LotVmod, y = State, parms = Pars, times = Time))
}
ymax <- max(sapply(solution,function(x){max(x[,c('x','y')])}))
par(mfrow=c(2,2))
for(i in 1 : 4)
{
plot(NA,xlab = "time", ylab = "population",ylim=c(0,ymax),xlim=c(0,timelength))
grid()
lines(solution[[i]]$time,solution[[i]]$x,col=rainbow_hcl(n=4)[i],lwd=2)
lines(solution[[i]]$time,solution[[i]]$y,col=rainbow_hcl(n=4)[i],lwd=2,lty=2)
abline(h=Pars['delta']/Pars['gamma'],lty=2,col='steelblue',lwd=2) # fixed point x
abline(h=Pars['alpha']/Pars['beta'],lty=2,col='tomato',lwd=2) # fixed point y
}
# Phase Space
plot(NA,xlim=c(0,ymax),ylim=c(0,ymax),xlab='prey',ylab='predator',main=paste0('Phase Space for (A,B,G,D)=(',paste(as.vector(Pars),collapse=","),')'))
grid()
for(i in 1:4){lines(solution[[i]]$x,solution[[i]]$y,col=rainbow_hcl(n=5)[i],lwd=2)}
for(i in 1:4){points(solution[[i]]$x[1],solution[[i]]$y[1],col=rainbow_hcl(n=5)[i],lwd=2)}
The Goodwin Model of Class Struggle
The main thrust of the Goodwin model is that distributional conflict produces endogenous cycles in employment and the labor share.
The key element through which this happens is the fact that out of a balanced growth path –a path where all variables grow at the same rate, namely the growth rate of labor productivity in the conventional wage share model– wages and labor productivity might grow at different rates. In particular, upward or downward pressure on wages may result from tightness in the labor market. The employment rate can be used as a measure of labor market tightness.
The economic intuition is the following. Suppose the economy is expanding, and employment increases. Higher labor demand generates wage inflation which, as long as real wages increase more than labor productivity, increases the wage share in output. If, consistently with the Classical view, workers do not save, the resulting decrease in the profit share will act in reducing future investment and output – i.e. a profit squeeze. But then the economy is down, and lower labor demand will then correspond to lower output, leading the way to lower wage inflation or even wage deflation. The labor share will decrease. But a higher profit share will produce a surge in investment, which will generate higher employment, thus improving workers’ bargaining power and consequently wages. At this point, the wage share has increased, and the cycle can repeat itself [Tavani, Notes on Goodwin].
Assumptions
- steady technical progress (disembodied)
- steady growth in the labour force
- only two factors of production, labour and capital, both homogeneous and non-specific
- all quantities real and net
- all wages consumed, all profits saved and invested (classical condition; constant proportional savings)
- a constant capital-output ratio
- a real wage rate which rises in the neighbourhood of full employment, i.e. \(\frac{w}{w} = -\gamma + \rho v\)
Model
\(q\) is output, \(k\) is capital, \(w\) is the wage rate, \(a=a_0 e^{\alpha t}\) is labour productivity, where \(\alpha\) is constant, \(\sigma\) is the capital-output ratio (inverse of capital productivity), \(u=w/a\) is the labour share of income, \(n=n_0 e^{\beta t}\) is labour supply with constant \(\beta\), \(l=q/a\) is employment, and \(v=l/n\) is the employment rate.
Investment \(\dot{k}\) comes from profits, which are the capital share of income \(l-w/a\) times output \(q\), \[ \dot{k} = (l - w/a)q \] The profit rate is \[ \frac{\dot{k}}{k} = \frac{(l - w/a)q}{k} = \frac{l - w/a}{\sigma} \] \[ \frac{\dot{q/l}}{q/l} = \frac{\dot{q}}{q}- \frac{\dot{l}}{l} = \alpha \] so that \[ \frac{\dot{l}}{l} = \frac{l - w/a}{\sigma} - \alpha \] Then, we compute for the labour share of income \[\frac{\dot{u}}{u} = \frac{\dot{w}}{w} - \alpha = - (\alpha + \gamma) + \rho v\] We finally obtain the equivalent of the Lotka-Volterra oscillations for class struggle, where the wage share of income \(u\) plays the role of predator and the employment rate \(v\) plays the role of prey: \[ \dot{v} = v \left( \underbrace{[\frac{1}{\sigma}- (\alpha+\beta)]}_{\alpha^*} - \underbrace{\frac{1}{\sigma}}_{\beta^*} u \right) = v (\alpha^* - \beta^* u) \] \[ \dot{u} = u \left( -\underbrace{[\alpha + \gamma]}_{\gamma^*} + \underbrace{\sigma}_{\delta^*} v \right) = u (-\gamma^* + \delta^* v) \]
``In this form we recognize the Volterra case of prey and predator. To some extent, the similarity is purely formal, but not entirely so. It has long seemed to me that Volterra’s problem of the symbiosis of two populations –partly complementary, partly hostile– is helpful in the understanding of the dynamical contradictions of capitalism, especially when stated in a more or less Marxian form.” Goodwin (1982)
Clockwise Cycles in the Employment-Distribution Space Garcı́a Molina and Herrera Medina (2010)
alpha <- 0.25 # labor productivity growth rate
sigma <- 1 # capital/output ratio
beta <- 0.25 # population growth rate
gamma <- 0.25 # intercept of the wage curve close to full employment
parameters <- c(1/sigma - alpha - beta, # alpha*
1/sigma, # beta*
sigma, # delta*
alpha + gamma) # gamma*
orbit <- lotka.volterra(pars = parameters, time.length = 30, interval = 0.0001, init = c(0.5, 0.75))
plot(orbit$x, type='l' , xlim = c(0, nrow(orbit)), ylim=c(0, max(orbit)), xlab = 'time', ylab='', col='steelblue', lwd=2)
grid()
lines(orbit$y, col = 'tomato', lwd=2)
abline(h = parameters[4]/parameters[3], lty=2, col = 'steelblue', lwd = 2) # fixed point for prey
abline(h = parameters[1]/parameters[2], lty=2, col = 'tomato', lwd = 2) # fixed point for predator
legend('topleft',c('employment rate','labour share of income'), col = c('steelblue','tomato'), pch=16, bty='n')
alpha <- 0.25 # labor productivity growth rate
sigma <- 0.75 # capital/output ratio
beta <- 0.25 # population growth rate
gamma <- 0.25 # intercept of the wage curve close to full employment
parameters <- c(1/sigma - alpha - beta, # alpha*
1/sigma, # beta*
sigma, # delta*
alpha + gamma) # gamma*
orbit <- lotka.volterra(pars = parameters, time.length = 30, interval = 0.0001, init = c(0.5, 0.75))
plot(orbit$x, type='l' , xlim = c(0, nrow(orbit)), ylim=c(0, max(orbit)), xlab = 'time', ylab='', col='steelblue', lwd=2)
grid()
lines(orbit$y, col = 'tomato', lwd=2)
abline(h = parameters[4]/parameters[3], lty=2, col = 'steelblue', lwd = 2) # fixed point for prey
abline(h = parameters[1]/parameters[2], lty=2, col = 'tomato', lwd = 2) # fixed point for predator
legend('topleft',c('employment rate','labour share of income'), col = c('steelblue','tomato'), pch=16, bty='n')
The Flaschel-Semmler Model
Classical Gravitation
The Classical political economists did not believe that this competitive process would lead to an actual equalization of realized or prospective profit rates at any moment in time. The movement of capital from one line of production to another would upset the conditions of other lines of production, which, together with disturbances from outside the national economy, would always prevent the realization of a state of equalization of profit rates. They expected to see a ceaseless fluctuation of prices and profit rates as the outcome of the competitive process, rather than the achievement of a state of “equilibrium” in which prices settled down to levels (“natural prices”) at which profit rates would be equalized. Foley (2003)
Nonetheless, the concept of this equilibrium state (which has come to be referred to as “long-period” equilibrium) plays a natural and important part in the analysis of the real economy. The competitive dynamic, even if it is not stable in the mathematical sense of pushing the system to an equilibrium of equal profit rates, will prevent prices and profit rates from wandering indefinitely far from their equilibrium values. This idea is expressed by arguing that observed market prices tend to gravitate around the natural prices at which profit rates would be equalized. The abstract concept of long-period equilibrium natural prices plays a crucial analytical role in understanding the concrete fluctuations of observable market prices. Foley (2003)
The Flaschel-Semmler model of multi-sector growth (Flaschel and Semmler 1987) with circulating capital and constant technology describes price (as prey) and quantity (as predator) oscillations over time around the equilibrium values of N prices p (a row vector) and N quantities x (a column vector) of the Sraffa-von Neumann system (Von Neumann 1945; Sraffa 1960).
The neoclassical adjustment process is generally regarded to be insufficient to ensure convergence to equilibrium on its own (Hahn 1970). As Flaschel and Semmler note, this led “neoclassical equilibrium theorists to avoid such problems completely by restricting their attention exclusively to equilibrium situations and assuming a ‘law of demand and supply’ in operation which works with infinite speed in an asymptotically stable manner” (Flaschel and Semmler 1987). As a consequence, adjustment to short-run equilibrium is critically taken as a matter of fact without characterizing the dynamical process by which it happens, thus treating the time dimension in purely hypothetical, yet not real, terms.
Equilibrium
Equilibrium values in quantities are obtained from the equality of supply and demand, \[\underbrace{B x^*}_{\text{outputs (supply)}} = R \underbrace{A x^*}_{\text{inputs (demand)}}\] while equilibrium values in prices from the uniformity of profit rates, \[\underbrace{p^* B}_{\text{unit revenues}} = R \underbrace{p^* A}_{\text{unit costs}}\] Balanced Growth : \(p^*\) and \(x^*\) are the positive eigenvectors of \(A\), the matrix of input-output coefficients, and \(R = 1+r^*\), the expansion factor, corresponds to the dominant eigenvalue.
In scalar notation, \[\sum_j b_{ij} x^*_j = R \sum_j a_{ij} x^*_j \quad i=1,...,N\] \[ \sum_i p^*_i b_{ij} = R \sum_i p^*_i a_{ij} \quad j=1,...,N\]
Constant-Technology Dynamics
The model is based on the dynamic cross-dual linear adjustment between prices and quantities.
- Law of Excess Demand
- Demand-supply imbalances trigger changes in prices \[ \left(\frac{\dot{p}}{p}\right)^T = - \delta_{p} [B-RA] x = \delta_p [\underbrace{R A x}_{\text{demand}} - \underbrace{Bx}_{\text{supply}}] \]
- Law of Excess Profitability
- Profit imbalances trigger changes in quantities
\[ \frac{\dot{x}}{x} = + \delta_x [B-RA]^T p^T = \delta_x [\underbrace{B^Tp^T}_{\text{revenue}} - R \underbrace{A^Tp^T}_{\text{cost}}] \] \(\delta_p\) and \(\delta_x\) are vectors of adjustment coefficients to estimate econometrically in order to calibrate empirically the model.
Econometrics
The empirical calibration involves the econometric estimation of vectors of country-, sector-specific adjustment coefficients \(\delta^i_p\) and \(\delta^i_x\) that relate growth rates in prices and quantities to supply-demand and revenue-cost imbalances, respectively (i.e. excess demand and excess profitability), in a linear form: \[ y_{a,it} = \delta^i_a x_{a,it} \qquad a = p,x; \quad i=1,...,N; \quad t=2000,...,2014 \] This problem can be addressed applying a hierarchical mixed-effects linear model with varying slopes with respect to the country sector.
Computation
First, call the required packages (only for visualization) and define the technological matrix \(A\). In this code, there are three different technological matrices of different size (number of sectors: 2, 4, and 7) from which we can compute its dynamical evolution in prices and quantities. While the first two are toy matrices, the latter corresponds to a simple version of the US economy with 7 sectors (Miller and Blair 2009). We compute the equilibrium triplet where \(r0\) is the equilibrium profit rate, \(p0\) is the equilibrium price vector, and \(x0\) is the equilibrium output vector. In order to do so, we have to be careful with the dimensions of our variables (whether we have column or row vectors). \(gam\) corresponds to the \(\gamma\) adjustment parameter.
# FLASCHEL-SEMMLER MODEL
# WITH ADJUSTMENT TO STABILIZE THE OSCILLATIONS
library(colorspace)
library(scales)
library(RColorBrewer)
timelength <- 1000
# A <- matrix(c(0.625,0.2,0.8,0.05), nrow = 2)
A <- matrix(c(0.1,0.4,0,0.3,0,0.1,0.2,0.5,0.3,0.2,0.3,0,0.1,0.1,0.5,0.1), nrow = 4)
A <- matrix(c(.2008,.001,.0034,.1247,.0855,.0897,.0093,0,.0658,.0002,.0684,.0529,.1668,.0129,.0011,.0035,.0012,.1801,.0914,.1332,.0095,.0338,.0219,.0021,.2319,.0952,.1255,.0197,.0001,.0151,.0035,.0339,.0645,.1647,.019,.0018,.0001,.0071,.0414,.0315,.2712,.0184,.0009,.0026,.0214,.0726,.0528,.1873,.0228), nrow = 7)
K <- nrow(A)
r_eq <- r0 <- Re(1/eigen(A)$values[1] - 1)
q0 <- Re(eigen(A)$vectors[,1])
leontief <- solve(diag(K) - A) # Leontief Requirements Matrix
y0 <- 100*t(t(q0)) # Final Demand: we multiply the vector by 50
dx <- rep(0.1, K) # Adjustment parameters
dp <- rep(0.01, K) # Adjustment parameters
gam <- 0.98 # Gamma
x0 <- leontief %*% y0 # Input Vector # (A %*% p0) capital/output # actually p0 %*% A
p0 <- Re(t(eigen(t(A))$vectors[,1]))
# x0<-c(0.45,0.4)
# p0<-c(0.6,0.6)
# (diag(K)-R*A) %*% x0 should be zero (i.e. there is no extra surplus)
R <- 1 + r_eqIn the following code, we initialize our simulations, where we create a 3-dimensional array called capital defined by number of sectors \(\times\) number of attributes \(\times\) number of iterations where we will save the synthetic data we compute. Initial conditions are the equilibrium values of prices and quantities plus an error. Attributes are the capital stock \(M\), the sector index \(k\), the technological inputs for sector \(k\), price \(p\), profit rate \(r\), quantity \(X\), sales \(S\), demand \(D\), profits \(pi\), capital/output ratio \(K\), and growth rate \(g\).
# INITIALIZATION
labels <- c('M','k',paste0('a',1:K),'p','r','X','S','D','pi','K','g')
K <- nrow(A)
capital <- array(0,dim=c(K,length(labels),timelength),dimnames=list(paste0('K',c(1:K)),labels,paste0('t',c(1:timelength))))
capital[,paste0('a',1:K),1]<-A
capital[,c('p','X'), 1] <- c(p0, x0)
capital[,'X', 1]<- x0 + 5*rep(c(.1,-.1),K)[1:K]
capital[,'p', 1]<- p0 - 0.01*rep(c(.1,-.1),K)[1:K]
capital[,'M',1] <- capital[,'p',1] %*% A * capital[,'X',1]
capital[,c('pi','g'),1] <- NA
capital[,'r',1] <- capital[,'p',1]*capital[,'X',1]/(capital[,'p',1] %*% A * capital[,'X',1])-1Once we have set up the initial conditions, we can start the simulation, following the dynamical system:
\[ x_{t+1} = x_t + \langle d_x \rangle \langle x \rangle \left( [B-RA]^T p^T + \gamma s \right) \] \[ p_{t+1} = p_t - \langle d_p \rangle \langle p_t \rangle (B - RA) x_t \]
where there is an additional stabilization term \(S\):
\[s = \frac{d}{dt} [B - RA]^T p^T \]
for (t in 1:(timelength-1))
{
# Compute price and quantity
capital[,'p',t+1] <- capital[,'p',t] - dp * t((diag(K) - R * A) %*% t(t(capital[,'X', t])))*capital[,'p', t]
s<-t(diag(K)-R*A) %*% t(t(capital[,'p',t+1]-capital[,'p',t]))
capital[,'X',t+1]<-capital[,'X',t]+dx*(t(diag(K)-R*A) %*% t(t(capital[,'p',t]))+gam*s)*capital[,'X',t]
# Auxiliary Variables
capital[,'S',t+1]<-capital[,'X',t+1]
capital[,'D',t+1]<-capital[,'X',t+1]-(A %*% t(t(capital[,'X',t+1])) + y0)
capital[,'M',t+1]<-capital[,'p',t] %*% A * capital[,'X',t+1]
capital[,'g',t]<-capital[,'M',t+1]/capital[,'M',t]-1
capital[,'pi',t+1]<-capital[,'p',t+1]*capital[,'X',t+1]
capital[,'K',t+1]<-capital[,'M',t+1]/capital[,'X',t+1]
capital[,'r',t+1]<-capital[,'p',t]*capital[,'X',t]/(capital[,'p',t] %*% A * capital[,'X',t])-1
}Once we have already computed the capital array, we can visualize it.
colores<-rainbow_hcl(K)
lab_<-c('capital stock (M)',"output (M')",'price (p)','output (X)','profit rate (r)', 'growth (g)')
names(lab_)<-c('M','pi','p','X','r','g')
tval<-1:timelength
par(mfrow=c(2,3),mar=c(3, 2.7, 1.6, 1.4),oma=c(0,0,0,0))
for(lab in names(lab_))
{
yli<-c(min(Re(capital[,lab,tval]),na.rm=T),max(Re(capital[,lab,tval]),na.rm=T))
# yli<-c(min(avg[,lab,tval],na.rm=T),max(avg[,lab,tval],na.rm=T))
if(lab=='x'){yli[1]<-0}
#if(lab=='p'){yli[1]<-4}
plot(NA,xlim=c(min(tval),max(tval)),ylim=yli,xlab='',ylab='',main=lab_[lab],font.lab=2,cex.lab=1)
grid()
#if(lab=='r' | lab=='p'){for(i in 1:N){points(tval,capital[i,lab,tval],col=alpha(colori[N-i+1],0.25),cex=0.8,pch=16)}} else{for(i in 1:N){lines(tval,capital[i,lab,tval],col=alpha(colori[N-i+1],0.25),lwd=2)}}
for(k in 1:K){lines(tval,capital[k,lab,tval],col=colores[k],lwd=2.25)}
# if(lab=='x'){for(k in 1:K){lines(tval,yli[2]*N_h[tval,k]/N,col=rgb(t(col2rgb(colores[k])/1.4),maxColorValue=255),lwd=1.9,lty=3)}}
# if(lab=='x'){axis(4,at=seq(0,yli[2],length.out=5),labels=seq(0,yli[2],length.out=5)*N/yli[2])}
#if(lab=='p'){for(k in 1:K){lines(tval,p0[k,tval],col='black',lwd=0.8,lty=2)}}
if(lab=='p'){abline(h=p0,col='black',lwd=0.8,lty=2)}
if(lab=='p'){mtext(paste0('p',1:K),side=4,at=p0,las=1,cex=0.8,font=2)}
if(lab=='X'){abline(h=x0,col='black',lwd=0.8,lty=2)}
if(lab=='X'){mtext(paste0('X',1:K),side=4,at=x0,las=1,cex=0.8,font=2)}
if(lab=='r'){abline(h=r_eq,lty=2)}
legend('topleft',paste0('K',1:K),col=colores,pch=16,pt.cex=2,cex=0.9,bty='n')
mtext('time',side=1,cex=0.8,line=2,font=2)
}
Early Agent-Based Models
Cellular Automata
A cellular automaton is a discrete model studied in computer science, mathematics, physics, complexity science, theoretical biology and microstructure modeling.
A cellular automaton consists of a regular grid of cells, each in one of a finite number of states, such as on and off (in contrast to a coupled map lattice). The grid can be in any finite number of dimensions. For each cell, a set of cells called its neighborhood is defined relative to the specified cell. An initial state (time t = 0) is selected by assigning a state for each cell. A new generation is created (advancing t by 1), according to some fixed rule (generally, a mathematical function) that determines the new state of each cell in terms of the current state of the cell and the states of the cells in its neighborhood. Typically, the rule for updating the state of cells is the same for each cell and does not change over time, and is applied to the whole grid simultaneously, though exceptions are known, such as the stochastic cellular automaton and asynchronous cellular automaton.
Hence, we can consider CA as dynamical systems with (often) deterministic rules \(R\), operating on discrete time \(T\), over discrete space \(S\).
The theory of cellular automata, the “canonical complex system” (Bossomaier 2008, 59), has its origin in the later work of Von Neumann, who was dissatisfied with the theory of games he had created with Oskar Morgenstern in order to understand strategic interaction (Mirowski 1992). Cellular automata were framed as a possible idealization of biological systems [von Neumann 1966], with the particular purpose of modeling biological self-reproduction [Wolfram 1983]. An automaton is a mathematically well-defined set of rules through which an entity interacts with its environment [Foley 2005]. A cellular automaton is a particular simple class of automata that lives on a graph or network. For instance, BTW’s sandpile model is an example of a cellular automaton on a two-dimensional regular lattice. Each node of the graph can take on a finite range of states.
The key features of CAs necessary for practical applications are the locality, an evolution rule, the set of states available, the configurations feasible, and boundary conditions (Bossomaier 2008). In almost all cases, cellular automaton evolution is irreversible.
In Stephen Wolfram’s classic classification, cellular automata fall qualitatively into four basic behavior classes (Wolfram 1984; Bossomaier 2008):
- Class 1 leading to stable, static configurations, akin to very short transients to mainly a point attractor: the evolution is dominated by a unique state of its alphabet for any random initial condition.
- Class 2 settling in oscillations, akin to a limit cycle: the evolution is dominated by blocks of cells which are periodically repeated for any random initial condition.
- Class 3 induce chaotic attractors, with aperiodic oscillations; if for a long time and for any random initial condition, the evolution is dominated by sets of cells without any defined pattern.
- Class 4 further display complex patterns, with effectively very long transients; structured, but endlessly changing. are the really interesting ones. They display complex patterns, structured, but endlessly changing. The evolution is dominated by non-trivial structures emerging and traveling along the evolution space where uniform, periodic, or chaotic regions can coexist with these structures.
In 1984, Wolfram conjectured Class 4 automata may be capable of universal computation, but a few have been proved to do so, notably by Langton Langton (1990).
Examples include: von Neumann was the first to use automata theory to study the logical intricacies of biological reproduction, detailing an infinite two-dimensional array or cellular space of identical finite automata supporting an activity interpreted as self-reproduction (Neumann 1966). Von Neumann later showed the computational power of the highly complex cellular space, with 29 states per cell, and the self-reproducing, computing construction including 40,000 cells or so. Such a cellular space is then computationally equivalent to any given Turing machine.
Wolfram 1-dimensional’s CA (Wolfram 2002) The CA are one-dimensional with the states of the cells (0 or 1) at a given time forming each row of the figure and with increasing row number going downwards showing the evolution in time according to 2N different rules where N is the number of cells.
In late 1940, John von Neumann defined life as a creation (as a being or organism) which can reproduce itself and simulate a Turing machine. Von Neumann’s cellular automata are the original expression of cellular automata, the development of which was prompted by suggestions made to John von Neumann by his close friend and fellow mathematician Stanislaw Ulam. Their original purpose was to provide insight into the logical requirements for machine self-replication, and they were used in von Neumann’s universal constructor.
Conway’s Game of Life
The Game of Life, also known simply as Life, is a cellular automaton devised by the British mathematician John Horton Conway in 1970. Conway died on April 11th 2020 due to the coronavirus.
The game is a zero-player game, meaning that its evolution is determined by its initial state, requiring no further input. One interacts with the Game of Life by creating an initial configuration and observing how it evolves. It is Turing complete and can simulate a universal constructor or any other Turing machine.
The universe of the Game of Life is an infinite, two-dimensional orthogonal grid of square cells, each of which is in one of two possible states, alive or dead, (or populated and unpopulated, respectively). Every cell interacts with its eight neighbours, which are the cells that are horizontally, vertically, or diagonally adjacent.
At each step in time, the following transitions occur:
- Any live cell with fewer than two live neighbours dies, as if by underpopulation.
- Any live cell with two or three live neighbours lives on to the next generation.
- Any live cell with more than three live neighbours dies, as if by overpopulation.
- Any dead cell with exactly three live neighbours becomes a live cell, as if by reproduction. An example of Conway’s Game of Life can be found here.
Self-Organized Criticality: the Sandpile Model
Following Anderson’s words, Bak and co-authors coined the notion of “self-organized criticality” from statistical mechanics (Bak, Tang, and Wiesenfeld 1987): they aimed at describing a class of dynamical systems which naturally drive them- selves to a state where interesting dynamics occurs on all scales, that is, at developing a possible explanation of the omnipresent multi-scale structures throughout the natural and social world (Creutz 1997).
In this context, complex behavior in nature reflects the tendency of large systems with many components to evolve into a poised, “critical” state, way out of balance, where minor disturbances may lead to events, called avalanches, of all sizes.
The sandpile model generates very interesting complex fractal patterns, conceived as “snapshots of systems operating at the self-organized critical state”. Crucially, Bak et al showed the distribution of avalanches followed a power law. Power laws are relevant because:
- they are scale-free, that is, scaling the function with a constant a returns the same function multiplied by a−k). In other words, there is no preferred length scale, so no dimension units are well-defined.
- they lack well-defined statistics (the average is well-defined only if k > 2 and the variance only if k > 3; most phenomena show 2 < k < 3).
- they are universal in the sense that they belong to a class of universality determined by their critical exponents, revealing fundamental dynamics underlying diverse systems approaching criticality.
Power laws are of great importance since they feature heavy or fat tails in contrast to the Gaussian distribution, so that there is an over-representation in the data of large-scale events, that is, “black swans”.
Power laws are ubiquitous in nature: for instance, the Stefan-Boltzmann law for the electromagnetic radiation of a black body, the Zipf law for the rank-frequency distribution for words in natural languages; the distribution of city sizes, of avalanches, wars, and acts of terrorism. In economics, the Pareto distribution of income [Dragulescu and Yakovenko 2001], the distribution of price changes in the stock market [Mandelbrot 1963; Farmer et al. 2005] or the Zipf-law distribution of firm sizes [Simon and Bonini 1958; Axtell 2001] are all power laws, while the stationary distribution of profit rates [Scharfenaker and Semieniuk 2016] is a double-exponential that displays a heavy tail.
Schelling’s Urban Segregation Model
The Schelling model of segregation is an agent-based model that illustrates how individual tendencies regarding neighbors can lead to segregation. The model is especially useful for the study of residential segregation of ethnic groups where agents represent households who relocate in the city. In the model, each agent belongs to one of two groups and aims to reside within a neighborhood where the fraction of ‘friends’ is sufficiently high: above a predefined tolerance threshold value F. It is known that depending on F, for groups of equal size, Schelling’s residential pattern converges to either complete integration (a random-like pattern) or segregation. An example here.
Schelling’s model roughly corresponds to an Ising model at temperature T = 0.
Ising Model
The Ising model is a mathematical model of ferromagnetism in statistical mechanics. The model consists of discrete variables that represent magnetic dipole moments of atomic “spins” that can be in one of two states (+1 or 1). The spins are arranged in a graph, usually a lattice (where the local structure repeats periodically in all directions), allowing each spin to interact with its neighbors. Neighboring spins that agree have a lower energy than those that disagree; the system tends to the lowest energy but heat disturbs this tendency, thus creating the possibility of different structural phases. The model allows the identification of phase transitions, as a simplified model of reality.
The two-dimensional square-lattice Ising model is one of the simplest statistical models to show a phase transition.
There is also a long tradition of using the Ising model and its extensions to represent social interactions and organization (Wiedlich, 1971; 1991; 2000; Callen and Shapero, 1974; Montroll and Badger, 1974; Galam et al., 1982; Orlean, 1995).
A large set of economic models can be mapped onto various versions of the Ising model to account for social influence in individual decisions (Phan et al.,2004).
The Ising model is indeed one of the simplest models describing the competition between the ordering force of imitation or contagion and the disordering impact of private information or idiosyncratic noise, which leads already to the crucial concept of spontaneous symmetry breaking and phase transitions (McCoy and Wu, 1973). It is therefore not surprising to see it appearing in one guise or another in models of social imitation (Galam and Moscovici, 1991) and of opinion polarization (Galam, 2004; Sousa et al., 2005; Stauffer, 2005; Weidlich and Huebner, 2008). An illuminating way to justify the use in social systems of the Ising model (and of its many generalizations) together with a statistical physics approach (in terms of the Boltzmann factor) derives from discrete choice models.
# INITIALIZATION
library(caTools)
N <- 100
timelength <- 1000
spin <- array(NA, dim = c(timelength+1, N, N))
spinU <- array(NA, dim = c(timelength+1, N, N))
spin[1, , ] <- matrix(sample(c(-1, 1), N^2, replace = TRUE, prob = c(0.5, 0.5)), nrow = N, ncol = N)
#spin[1, , ] <- matrix(rep(1, N^2), nrow = N, ncol = N)
number <- c() # frequencies
temperature <- 0.1
beta <- 1/(temperature*1.38*10^{-23})
# SIMULATION
for(t in 1:timelength)
{
# if((t-1)%%7 == 0)
# {
# beta = betavector[count]
# count <- count+1
# }
for(i in 1:N)
{
for(j in 1:N)
{
# Include Boundary Conditions
i0 <- i %% N + 1
j0 <- j %% N + 1
iF <- (i + N - 2) %% N + 1
jF <- (j + N -2) %% N + 1
old <- spin[t, i, j]
new <- -old
val <- spin[t, i, j0] + spin[t, i0, j] + spin[t, iF, j] + spin[t, i, jF]
oldU <- old*val
newU <- new*val
dU <- beta*(newU - oldU)
spinU[t, i, j] <- newU-oldU
if(!dU <= 0 | (exp(dU) > runif(1))) spin[t+1, i, j] <- new else spin[t+1, i, j] <- old
}
}
number[t] <- table(spin[t, , ])['1']/N^2
}
spin1 <- aperm(spin, c(2, 3, 1)) # we need to change dimensions for the gif
write.gif(spin1, 'ising.gif', col = "jet", delay = 0)
# banana <- image_read("ising.gif")
# banana <- image_scale(banana, "400%") # resize the gif
# image_write(banana, "ising2.gif")
More Oscillators
Van der Pol Oscillator
The Van der Pol oscillator is a non-conservative oscillator with non-linear damping. It evolves in time according to a second-order differential equation with μ as a scalar parameter indicating the nonlinearity and the strength of the damping. Van der Pol found stable oscillations which he subsequently called relaxation-oscillations in the form of a limit cycle. When electric circuits of vacuum tubes were driven near the limit cycle, they become entrained, i.e. the driving signal pulls the current along with it. Van der Pol and his colleague, van der Mark, reported in the September 1927 issue of Nature that at certain drive frequencies an irregular noise was heard, which was later found to be the result of deterministic chaos. The Van der Pol equation has a long history of being used in both the physical and biological sciences, for instance the FitzHugh-Nagumo model for action potentials of spiking neurons.
Kaldor’s nonlinear model of business cycles
Inspired by Keynes’ income theory and Kalecki’s model of investment (Kalecki 1937), Kaldor (1940) formulated the first nonlinear model of endogenous business cycles by considering the interactions between the investment I(Y) and the savings S(Y) functions (both dependent on income Y) and the existence conditions for self-sustaining limit cycles. By noting that the linear forms of I(Y ) and S(Y ) fail to produce cyclical motions, Kaldor proposed a S-shaped (sigmoid) nonlinear form for I(Y ) and a mirror-imaged S-shaped nonlinear form for S(Y ) (Gabisch and Lorenz 1989), which yields oscillatory motion of business cycles. Chang and Smyth (1971) reformulated Kaldors model of business cycles, given by the following coupled dynamical equations:
FitzHugh-Nagumo
The FitzHughNagumo model (FHN), named after Richard FitzHugh who suggested the system in 1961 and J. Nagumo et al. who created the equivalent circuit the following year, describes a prototype of an excitable system (e.g. a neuron). The FHN Model is an example of a relaxation oscillator because, if the external stimulus Iext exceeds a certain threshold value, the system will exhibit a characteristic excursion in phase space, before the variables v and w relax back to their rest values. Further, Iext can be considered as a Gaussian white noise that may occasionally induce the neuron to spike, as well as a coupling term with adjacent oscillators, generating interesting interactions between coupled elements.