Sequence of 1-parameter bifurcation diagrams for OS model with feedback control

We create a 1-parameter bifurcation diagrams in the bias parameter p=(eps_+ - eps_-)/2. This is also a good parameter to use as feedback control input.

Contents

clear
clf
addpath('../utilities/');

System parameters

N=2000;                 % number of traders
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',4,...           % amplification: how much change buy/sell the news?
    'epsplus', 0.075,...% size of jump up for good news
    'epsminus',0.065);  % size of jump down for bad news
navg=5;                 % number of time steps over which buying and selling
                        % rates should be averaged
%
nsteps=2000;            % number of time steps OS model is run
t=(1:nsteps)*par.tstep; % time points
histbins=20;            % number of bins for histograms of agents distribution
fdeco={'fontweight','bold','fontsize',14};
ldeco={'linewidth',2};

Feedback control will be p=pref+gain(R_+ + R_-)

where p is the bias: (par.epsplus-par.epsminus)/2 The control appears to be fairly globally stabilizing, so we can start bifurcation diagram from trader state equal zero.

pref=1e-3*(-15:15);
pscale=1e-2;
gain=-3e-1;
Rscale=pscale/abs(gain);
feedback_control=@(par,output,gain,ref)...
    setfield(setfield(par,...
    'epsplus',(par.epsplus+par.epsminus)/2+ref+gain*output),...
    'epsminus',(par.epsplus+par.epsminus)/2-(ref+gain*output)); %#ok<SFLD>

nruns runs of nsteps time stepper steps

nruns=length(pref);
pvals=NaN(nsteps,nruns); % array storing all values of p during all runs
Rp=pvals;                   % array of R_+ outputs
Rn=pvals;                   % array of R_- outputs
xm=pvals;                   % array of mean agent states
tail=round(nsteps/4):nsteps;% discard first quarter of time series for plotting
xdist=NaN(histbins,length(tail),nruns);  % histograms of agent distributions
x=zeros(N,1);         % initial state of all traders
evp=zeros(navg,nruns);% initialization for number of past buying events
evn=evp;              % initialization for number of past selling events

scan through all values of reference parameter

for k=1:length(pref)
    if k>1
        evp(k)=evp(k-1);
        evn(k)=evn(k-1);
    end

perform nsteps time steps starting from initial condition x=[0,...,0] or from previous run

    for i=1:nsteps
        Rpcur=mean(evp(k),1)/par.tstep;  % # buying events avg over past navg steps
        Rncur=mean(evn(k),1)/par.tstep;  % # selling events avg over past navg steps
        par=feedback_control(par,Rpcur+Rncur,gain,pref(k)); % feedback control applied to p (bias)
        Rp(i,k)=Rpcur;              % store R_+
        Rn(i,k)=Rncur;              % store R_-
        xm(i,k)=mean(x,1);          % store mean agent state
        pvals(i,k)=(par.epsplus-par.epsminus)/2;   % store current parameter value
        if i>=tail(1)
            [xdist(:,i-tail(1)+1,k),xcenter]=agentdist(x,histbins);
        end
        [xnew,evp(k),evn(k)]=OmurtagSirovich(x,par,evp(k),evn(k)); % single time step
        x=xnew;
    end
    fprintf('k=%d of %d, mean(p)=%g, mean(R- + R+)=%g\n',...
        k,length(pref),mean(pvals(tail,k)),mean(Rp(tail,k)+Rn(tail,k)));
k=1 of 31, mean(p)=-4.59694e-05, mean(R- + R+)=-0.0498468
k=2 of 31, mean(p)=-0.000283145, mean(R- + R+)=-0.0457229
k=3 of 31, mean(p)=-0.000466356, mean(R- + R+)=-0.0417788
k=4 of 31, mean(p)=-0.00046569, mean(R- + R+)=-0.0384477
k=5 of 31, mean(p)=-0.000740839, mean(R- + R+)=-0.0341972
k=6 of 31, mean(p)=-0.00091006, mean(R- + R+)=-0.0302998
k=7 of 31, mean(p)=-0.00102931, mean(R- + R+)=-0.026569
k=8 of 31, mean(p)=-0.000960693, mean(R- + R+)=-0.0234644
k=9 of 31, mean(p)=-0.00102598, mean(R- + R+)=-0.0199134
k=10 of 31, mean(p)=-0.000949367, mean(R- + R+)=-0.0168354
k=11 of 31, mean(p)=-0.000908728, mean(R- + R+)=-0.0136376
k=12 of 31, mean(p)=-0.000852099, mean(R- + R+)=-0.010493
k=13 of 31, mean(p)=-0.00064557, mean(R- + R+)=-0.0078481
k=14 of 31, mean(p)=-0.000443038, mean(R- + R+)=-0.00518987
k=15 of 31, mean(p)=-0.000218521, mean(R- + R+)=-0.00260493
k=16 of 31, mean(p)=4.59694e-05, mean(R- + R+)=-0.000153231
k=17 of 31, mean(p)=0.00022052, mean(R- + R+)=0.00259827
k=18 of 31, mean(p)=0.000465023, mean(R- + R+)=0.00511659
k=19 of 31, mean(p)=0.000707528, mean(R- + R+)=0.00764157
k=20 of 31, mean(p)=0.000826116, mean(R- + R+)=0.0105796
k=21 of 31, mean(p)=0.00096469, mean(R- + R+)=0.013451
k=22 of 31, mean(p)=0.000991339, mean(R- + R+)=0.0166955
k=23 of 31, mean(p)=0.000956029, mean(R- + R+)=0.0201466
k=24 of 31, mean(p)=0.00103065, mean(R- + R+)=0.0232312
k=25 of 31, mean(p)=0.000965356, mean(R- + R+)=0.0267821
k=26 of 31, mean(p)=0.000830113, mean(R- + R+)=0.0305663
k=27 of 31, mean(p)=0.000740839, mean(R- + R+)=0.0341972
k=28 of 31, mean(p)=0.00055563, mean(R- + R+)=0.0381479
k=29 of 31, mean(p)=0.000436376, mean(R- + R+)=0.0418787
k=30 of 31, mean(p)=0.000205197, mean(R- + R+)=0.0459827
k=31 of 31, mean(p)=5.996e-05, mean(R- + R+)=0.0498001
end

1d bifurcation diagram with histogram coloring and time series of input and output

subplot(2,1,1);
plot_histograms(pvals(tail,1:k),Rp(tail,1:k)+Rn(tail,1:k),'nb',6,...
    'ax',prctile(pvals(:),[1,99]),'ay',prctile(Rp(:)+Rn(:),[1,99]));
xlabel('p=(\epsilon_+-\epsilon_-)/2')
title({'bifurcation diagram, histograms and mean values';'of p=(\epsilon_+-\epsilon_-)/2 and R_++R_-'});
subplot(2,1,2);
ind=1:k*nsteps;
plot(ind*par.tstep,pvals(ind)/pscale,'.',ind*par.tstep,(Rp(ind)+Rn(ind))/Rscale,'.');
xlabel('t');
ylabel(sprintf('p/%3.1g, (R_++R_-)/%3.1g',pscale,Rscale));
title(sprintf('time series: p/%3.1g and (R_++R_-)/%3.1g',pscale,Rscale));
set(gca,fdeco{:},ldeco{:})

Densities of agents' internal states (in [-1,1]) along branch

pmean=mean(pvals(tail,:),1);
Rmean=mean(Rp(tail,:)+Rn(tail,:));
xdens=squeeze(mean(xdist,2));
figure(2);
subplot(1,2,1)
plot_histograms(pvals(tail,:),Rp(tail,:)+Rn(tail,:),'nb',6,'lw',3,...
    'ax',[min(pmean),max(pmean)],'ay',[min(Rmean),max(Rmean)]);
xlabel('p=(\epsilon_+-\epsilon_-)/2')
ylabel('R_++R_-');
pos1=get(gca,'position');
subplot(1,2,2);
contourf(xcenter,Rmean,xdens');
pos2=get(gca,'position');
set(gca,'position',[pos2(1),pos1(2),pos2(3),pos1(4)],'ytick',[],fdeco{:},ldeco{:})
xlabel('internal state')
ylabel('R_++R_-');
cb=colorbar;
title(cb,'density',fdeco{:})
set(cb,fdeco{:});
colormap(jet)
save('bif1d_sweep.mat');