%plot different random walks associated to different noise realizations
N=10000;
nseed=10;
%draw random Gaussian variables
x=randn(N,nseed);
figure(1)
subplot(2,1,1)
plot(x)
title('driving noise');
subplot(2,1,2)
plot(cumsum(x))
title('random walk');
%%
%plot different random walks associated to different noise realizations but
%with different timestep in the discretization scheme
%simulation duration
T=100;
%time steps
dt=[0.001,0.005,0.025,0.125];
%number of distinct time steps
nd=length(dt);
figure(2)
hold on
for i=1:nd
t=0:dt(i):T;
N =length(t);
%deterministic Euler discretization scheme
x=dt(i)*cumsum(randn(N,1));
length(x)
plot(t,x)
end
title('dt scaling');
%%
%simulation duration
T=100;
%time steps
dt=[0.001,0.005,0.025,0.125];
%number of distinct time steps
nd=length(dt);
figure(3)
hold on
for i=1:nd
t=0:dt(i):T;
N =length(t);
%stochastic Euler discretization scheme
x=sqrt(dt(i))*cumsum(randn(N,1));
length(x)
plot(t,x)
end
title('sqrt(dt) scaling');
%%
%simulation of a stochastic ODE with a deterministic part (relaxation to zero)
%and a stochastic part (Gaussian white noise). This noisy relaxation to
%zero is called an Ornstein-Uhlenbeck process
%simulation duration
T=100;
%time steps
dt=[0.001,0.005,0.025,0.125];
%number of distinct time steps
nd=length(dt);
figure(3)
hold on
for i=1:nd
t=0:dt(i):T;
N=length(t);
x=zeros(1,N);
x(1)=10;
for j=1:N-1
%Euler scheme for the Ornstein Uhlenbeck process (deterministic part"dt(i)*...", stochastic part "sqrt(dt(i))*...")
x(j+1)=x(j)-dt(i)*x(j)+sqrt(dt(i))*randn();
end
length(x)
plot(t,x)
end
title('Ornstein-Uhlenbeck processes');
%%
%simulation of coupled stochastic ODEs: n is an Ornstein-Uhlenbeck process
%that should be seen as a noise source perturbing a bistable system
%number of time steps
ntimes=1000000;
nconds=1;
%data vector (voltages)
v=zeros(nconds,ntimes);
n=zeros(nconds,ntimes);
r=randn(nconds,ntimes);
%initial condition
for i=1:nconds
v(i)=(i-nconds/2)/nconds;
end
%time step
dt=0.01;
%amplitude of the noise (the smaller the value of amp the least likely are
%transitions between stable wells)
amp=0.2;
for i=1:ntimes-1
v(:,i+1)=v(:,i)+dt*(v(:,i)-v(:,i).^3)+amp*n(:,i);
n(:,i+1)=n(:,i)-dt*n(:,i)+sqrt(dt)*r(:,i);
end
figure(4)
plot(v')
%%
%script computing the statisitics of switch times from one well to another
%detect switch to the "up" state i.e 1
up=(diff(v>0)==1);
%detect switch to the "down" state i.e -1
down=(diff(v<0)==1);
%bring together upward switches and downward switches
ind=up+down;
%find time index of switch events
time=find(ind==1);
%plot histogram of "dominance durations"
figure(5)
hist(diff(time),100)