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
- System parameters
- Feedback control will be p=pref+gain(R_+ + R_-)
- nruns runs of nsteps time stepper steps
- scan through all values of reference parameter
- perform nsteps time steps starting from initial condition x=[0,...,0] or from previous run
- 1d bifurcation diagram with histogram coloring and time series of input and output
- Densities of agents' internal states (in [-1,1]) along branch
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');