Demo for simplified Omurtag Sirovich model
see related papers by Omurtag & Sirovich and Siettos et al
Contents
- Background
- Set basic system parameters
- Initialize loop
- loop over all parameter values
- Each run starts from the last state of the previous run
- Plot results
- Bifurcation diagrams
- One-parameter bifurcation diagram (high-resolution in folder traders)
- Two-parameter bifurcation diagram in
and
- Zoom-in of two-parameter diagram near the cusp
Background
This is a discrete-time model. N traders trade the same asset. Each trader's state of mind is a number , expressing his/her inclination to buy or sell. At every time step the following happens.
The trader's state relaxes with rate back to zero (forgetting). The trader receives good news and bad news according to a Poisson distribution with mean rate
. For each piece of good news his/her state jumps
up and for each piece of bad news it jumps down by
. When the state of the trader crosses +1 in a time step s/he "buys", and the state jumps back to zero. When the state of the trader crosses -1 in a time step s/he "sells", and the state jumps back to zero.
A positive feedback effect is introduced in the following way: the mean of the Poisson distribution depends on the recent number of buys and sells via
where is the base arrival rate,
is a parameter controlling the feedback effect, and
is the buying/selling rate averaged over the last nT steps:
number of buys during last navg steps divided by time unit tstep (same for
and sells).
When and
we get a subcritical pitchfork bifurcation at some critical value of
. In the two-parameter plane
we get a cusp.
The demo below performs a simple parameter sweep, gradually increasing g from 2 until run-away buying is observed. We set the parameters slightly asymmetric ().
Set basic system parameters
clear addpath('../utilities/'); N=10000; % number of traders x=zeros(N,1); % initial state of all traders navg=5; % number of time steps over which buying and selling % rates should be averaged par=struct(... 'tstep',5e-2,... % time steps size (only used for calculating buy/sell/forget rates 'nup',20,... % default Poisson rate of arrival of good news 'nun',20,... % default Poisson rate of arrival of bad news 'gamma',1,... % default rate of forgetting 'g',2,... % amplification: how much change buy/sell the news? 'epsplus', 0.071,...% size of jump up for good news 'epsminus',0.069); % size of jump down for bad news nsteps=2000; % number of time steps per run threshold=40; % stop when |buyingrate - selling rate|> threshold gvals=2:0.1:7; % values for g to try nruns=length(gvals); % expected number of runs
Initialize loop
Rp=zeros(nsteps,nruns); % storage array for buying rates Rn=Rp; % storage array for selling rates xm=Rp; % storage array for mean agent states evp=zeros(navg,1); % initial history for number of buying events evn=evp; % initial history for number of selling events esc=false; % flag signaling if runaway has occured
loop over all parameter values
for k=1:nruns
par.g=gvals(k);
Each run starts from the last state of the previous run
for i=1:nsteps [xnew,evp,evn]=OmurtagSirovich(x,par,evp,evn); x=xnew; Rpcur=mean(evp,1)/par.tstep; Rncur=mean(evn,1)/par.tstep; Rp(i,k)=Rpcur; Rn(i,k)=Rncur; xm(i,k)=mean(x,1); if Rp(i,k)>threshold || Rn(i,k)<-threshold esc=true; break end end fprintf('k=%d,g=%g,xm=%g,Rp+Rn=%g\n',... k,par.g,mean(xm(:,k),1),mean(Rp(:,k)+Rn(:,k),1)); if esc break end
k=1,g=2,xm=0.0519675,Rp+Rn=0.0063854
k=2,g=2.1,xm=0.0543985,Rp+Rn=0.0065386
k=3,g=2.2,xm=0.0559977,Rp+Rn=0.006753
k=4,g=2.3,xm=0.0576206,Rp+Rn=0.007209
k=5,g=2.4,xm=0.0582166,Rp+Rn=0.0072216
k=6,g=2.5,xm=0.0603645,Rp+Rn=0.0072238
k=7,g=2.6,xm=0.061528,Rp+Rn=0.0075164
k=8,g=2.7,xm=0.0630029,Rp+Rn=0.0076118
k=9,g=2.8,xm=0.0660094,Rp+Rn=0.0081458
k=10,g=2.9,xm=0.0683721,Rp+Rn=0.0086804
k=11,g=3,xm=0.0707366,Rp+Rn=0.008863
k=12,g=3.1,xm=0.0730904,Rp+Rn=0.009303
k=13,g=3.2,xm=0.0763318,Rp+Rn=0.0097474
k=14,g=3.3,xm=0.0808009,Rp+Rn=0.0105666
k=15,g=3.4,xm=0.0824016,Rp+Rn=0.0105674
k=16,g=3.5,xm=0.0873056,Rp+Rn=0.0113974
k=17,g=3.6,xm=0.0931323,Rp+Rn=0.0125292
k=18,g=3.7,xm=0.098363,Rp+Rn=0.0133318
k=19,g=3.8,xm=0.111415,Rp+Rn=0.015852
k=20,g=3.9,xm=0.112361,Rp+Rn=0.087366
end
Plot results
The figure below shows the result of the parameter scan, showing histograms of |buying and selling rates The color coding shows the distribution of values, giving an idea about the fluctuations around the mean (black dot).
if esc k=k-1; end Rp=Rp(:,1:k); Rn=Rn(:,1:k); xm=xm(:,1:k); plot_histograms(gvals(ones(nsteps,1),1:k),Rp+Rn,... 'ay',[min(Rp(:)+Rn(:)),max(Rp(:)+Rn(:))]); xlabel('g'); ylabel('R_++R_-');

Bifurcation diagrams
Included are bifurcation diagrams from some preliminary computations using N=2000 traders and, apart from g, epsplus and epsminus the same paramerters as in the demo.
One-parameter bifurcation diagram (high-resolution in folder traders)
For the system is symmetric (selling and buying are treated the same). When varying g we obtain a subcritical pitchfork bifurcation. The branches of approximate unstable equilibria are the non-zero solutions. x-axis: g, y-axis:
(note that
is negative).
Two-parameter bifurcation diagram in
and 
(high-resolution and matlab figure in folder traders) When is different from
then the symmetry is broken. In this picture the mean of
and
is kept constant (0.07), but the difference between them is varied, along with g.The dots are results of the computations. The blue surface is an interpolation/regression of these points. The red curve is the curve of saddle-node bifurcations, intersecting the plane
in a cusp. The black dots are the intersection of the solution surface with
(corresponding to the one-parameter figure).
Zoom-in of two-parameter diagram near the cusp
(high-resolution and matlab figure in folder traders) Grey dots: points at which equilibrium was evaluated; red curve: fold in two-parameter plane (left) and parameter-solution projection
; pink curves quantify uncertainty: assuming that the disturbances of the measurements of equilibria follows a Gaussian distribution (the standard deviation was estimated) then the linear corrections of the fold curve in the Newton iteration have a disturbance with a 2.58 times standard deviation (0.99 quantile) as indicated by the pink curve.