Demo for simplified Omurtag Sirovich model

see related papers by Omurtag & Sirovich and Siettos et al

Contents

Background

This is a discrete-time model. N traders trade the same asset. Each trader's state of mind is a number $x\in[-1,1]$, expressing his/her inclination to buy or sell. At every time step the following happens.

The trader's state relaxes with rate $\gamma$ back to zero (forgetting). The trader receives good news and bad news according to a Poisson distribution with mean rate $\mu_\pm$. For each piece of good news his/her state jumps $\epsilon_+$ up and for each piece of bad news it jumps down by $\epsilon_-$. 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

$$\mu_\pm=\nu_\pm(1+g R_\pm)$$

where $\nu_\pm$ is the base arrival rate, $g$ is a parameter controlling the feedback effect, and $R_\pm$ is the buying/selling rate averaged over the last nT steps:

$R_+=$ number of buys during last navg steps divided by time unit tstep (same for $R_-$ and sells).

When $\epsilon_+=\epsilon_-$ and $\nu_+\nu_-$ we get a subcritical pitchfork bifurcation at some critical value of $g$. In the two-parameter plane $(g,\epsilon_+-\epsilon_-)/2$ 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 ($\epsilon_+>\epsilon_-$).

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 $\epsilon_+=\epsilon_-$ 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: $R_++R_-$ (note that $R_-$ is negative).

pitchfork bifurcation diagram

Two-parameter bifurcation diagram in $g$ and $(\epsilon_+-\epsilon_-)/2$

(high-resolution and matlab figure in folder traders) When $\epsilon_+$ is different from $\epsilon_-$ then the symmetry is broken. In this picture the mean of $\epsilon_+$ and $\epsilon_-$ 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 $g=0$ in a cusp. The black dots are the intersection of the solution surface with $g=0$ (corresponding to the one-parameter figure).

two-parameter bifurcation diagram

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 $(\epsilon_+-\epsilon_-)/2,g)$ (left) and parameter-solution projection $(g,R_++R_-)$; 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.

two-parameter bifurcation diagram