clc;
clear all;
%Chapitre 1
%Processus de Wiener
%%
%1.1.1. Simulation du mouvement Brownien standard à 4 dimensions
clc;
clear all;

T=365; % Instants du processus
experiences=10000; % Nombre d'expériences indépendantes

position=zeros(4,T+1,experiences); % Matrice pour stocker les expériences : (dimension, temps, expérience)
for i=1:experiences
    for t=1:T
        position(:,t+1,i)=position(:,t,i)+(1/sqrt(T))*randn(4,1);
    end
end

means=zeros(4,3);
stds=zeros(4,3);
range_of_k=[10 100 300];
correlations=zeros(4,4,3);
for i=1:3;
    k=range_of_k(i);
    means(:,i)=[mean(position(1,k,:)); mean(position(2,k,:)); mean(position(3,k,:)); mean(position(4,k,:))];
    stds(:,i)=[std(position(1,k,:)); std(position(2,k,:)); std(position(3,k,:)); std(position(4,k,:))];
    for m=1:4
        for n=m:4
    correlation_mn=corrcoef(squeeze(position(m,k,:)),squeeze(position(n,k,:)));
    correlations(m,n,i)=correlation_mn(1,2);
    correlations(n,m,i)=correlation_mn(1,2);
        end
    end
end


%1.1.2. Simulation du mouvement Brownien standard à 4 dimensions corrélées

Gamma = [1 0.2 0.8 0.5; 0.2 1 0.3 0.2; 0.8 0.3 1 0.9; 0.5 0.2 0.9 1];
choleski_G=chol(Gamma);
position_corrigee=zeros(4,T+1,experiences);
for i=1:experiences
    for t=1:T+1
        position_corrigee(:,t,i)=choleski_G*position(:,t,i);
    end
end
means_corrigee=zeros(4,3);
stds_corrigee=zeros(4,3);
range_of_k=[10 100 300];
correlations_corrigee=zeros(4,4,3);
for i=1:3;
    k=range_of_k(i);
    means_corrigee(:,i)=[mean(position_corrigee(1,k,:)); mean(position_corrigee(2,k,:)); mean(position_corrigee(3,k,:)); mean(position_corrigee(4,k,:))];
    stds_corrigee(:,i)=[std(position_corrigee(1,k,:)); std(position_corrigee(2,k,:)); std(position_corrigee(3,k,:)); std(position_corrigee(4,k,:))];
    for m=1:4
        for n=m:4
    correlation_mn=corrcoef(squeeze(position_corrigee(m,k,:)),squeeze(position_corrigee(n,k,:)));
    correlations_corrigee(m,n,i)=correlation_mn(1,2);
    correlations_corrigee(n,m,i)=correlation_mn(1,2);
        end
    end
end

%%
%1.2.1 Processus de Poisson
lambda=150; %intervalle moyen entre deux appels
N=10000; %nombre de clients
[t_attente, t_fin]=attente(gen_exp(N,lambda),normrnd(120,20,1,N));

hist(t_attente)
close all;

indices = find(t_fin>36000);
t_fin(indices) = [];
n=length(t_fin); %nombre de clients en 10 heures
t_attente=t_attente(1:n); %temps d'attente de ces clients

%Nous ne considérons que les gens qui ont été effecitvement traités dans la
%journée.

hist(t_attente)
proportion_dattente_10minutes=length(find(t_attente>600))/n;

%%
%1.3.1 Chaîne de Markov
close all;
clc;
X0=randi([-42 42]);
N=100000;
Z=2*randi([0 1],1,N)-ones(1,N);

X=zeros(1,N);
X(1)=X0;
for k=1:N-1
    X(k+1)=X(k)+Z(k);
end
plot(X)

%%
%1.4.1 Processus de diffusion
close all;
clc;
dt=0.01;
x0=1;
N=1000;

parameters=[1 1 1 0.01 1;-1 1 1 0.01 1;1 1 0.1 0.01 1;1 2 0.3 0.01 1];
%De la forme a b sigma dt x0
X=zeros(1,N);
X(1)=x0;
i=4;
for k=1:N-1
    X(k+1)=X(k)+(parameters(i,1)*(parameters(i,2)-X(k))*dt+parameters(i,3)*normrnd(0,dt));
end
%plot(X)

%1.4.2 deuxième equation

%Lorsque x est proche de 0.5, le facteur ALEATOIRE sera dominant. Lorsque x
%est proche de 0 ou 1, le facteur DE GAUCHE le rapproche de 0.5 => On
%reste entre 0 et 1
%L'amplitude du bruit augmente lorsque x s'approche de 0.5

dt=0.01;
N=1000;
X=zeros(1,N);
X(1)=rand;
for k=1:N-1
    X(k+1)=X(k)+(0.5-X(k))*dt+sqrt(X(k)*(1-X(k)))*normrnd(0,dt);
end
plot(X)


function [time] = gen_exp(n,lambda)
time=exprnd(lambda,1,n);