%Définition de "gamma"

function g=gamma(x,m,s,c,a)

g1=s/sqrt(2*pi)*exp(-(x-m-c/a)^2/(2*s^2))+(m+c/a-x)/2*erfc((x-(m+c/a))/(sqrt(2)*s));
g=a*x-g1;

%Méthode par dichotomie :
function  g=zgamma(n,m,s,c,a)

g0=0;
g1=m+2*c/a;

for i=1:n
    g2=(g1+g0)/2;
    k0=gamma(g0,m,s,c,a);
    k2=gamma(g2,m,s,c,a);
    
    if k0*k2<0 g1=g2;
    else g0=g2;
    end
end
g=g2;

%Premier calcul :
clear
m=200000;s=0.05*m;c=200;a=0.001;
P=zgamma(100,m,s,c,a)-c/a
g=zgamma(100,m,s,c,a);
p=1-normcdf(g-c/a,m,s);
T=(1-p)/p
e=sqrt(1-p)/p
Pam=g*(1+a)-c/a
c/a



%Simulation du flux actualisé à m,s,c,a fixés (le temps d'attente est-il
%réaliste?)
%Simulation des 3 situations d'offres de prix
clear
m=200000;s=0.1*m;c=200;a=0.001;
P=zgamma(100,m,s,c,a)-c/a;
g=zgamma(100,m,s,c,a);
p=1-normcdf(g-c/a,m,s);
T=(1-p)/p;
t=floor(6*T);

Offres=zeros(3,t);
for i=1:3
    Offres(i,:)=m*ones(1,t)+randn(1,t)*s;
end

%Moment d'atteinte de la valeur P
tp=ones(1,3);
for i=1:3
    j=1;
    while Offres(i,j)-P<0
        tp(i)=tp(i)+1;
        j=j+1;
    end
end

%Moment d'atteinte de la valeur m
tm=ones(1,3);
for i=1:3
    j=1;
    while Offres(i,j)-m<0
        tm(i)=tm(i)+1;
        j=j+1;
    end
end


Pa=zeros(3,t);
A=(1+a)^(-1)*ones(1,t);
for i=1:3
    for j=1:t
        A1=cumsum(cumprod(A(1,1:j)));
        Pa(i,j)=Offres(i,j)*(1+a)^(-j)-c*A1(1,j);
    end
end


subplot(3,1,1)
hold on
title('Evolution temporelle du prix actualisé')
xlabel('t')
ylabel('Pa')
plot(Pa(1,:));

Tp1=Pa(1,tp(1))*ones(1,t);
plot(Tp1,'r');
X=[tp(1),tp(1)];Y=[0,Pa(1,tp(1))];
plot(X,Y,'r');

plot(Pa(1,tm(1))*ones(1,t),'g');
X=[tm(1),tm(1)];Y=[0,Pa(1,tm(1))];
plot(X,Y,'g');

X=[T];Y=[0];
plot(X,Y,'ro');
hold off



subplot(3,1,2)
hold on
xlabel('t')
ylabel('Pa')
plot(Pa(2,:));

Tp2=Pa(2,tp(2))*ones(1,t);
plot(Tp2,'r');
X=[tp(2),tp(2)];Y=[0,Pa(2,tp(2))];
plot(X,Y,'r');

plot(Pa(2,tm(2))*ones(1,t),'g');
X=[tm(2),tm(2)];Y=[0,Pa(2,tm(2))];
plot(X,Y,'g');

X=[T];Y=[0];
plot(X,Y,'ro');
hold off



subplot(3,1,3)
hold on
xlabel('t')
ylabel('Pa')
plot(Pa(3,:));

Tp3=Pa(3,tp(3))*ones(1,t);
plot(Tp3,'r');
X=[tp(3),tp(3)];Y=[0,Pa(3,tp(3))];
plot(X,Y,'r');

plot(Pa(3,tm(3))*ones(1,t),'g');
X=[tm(3),tm(3)];Y=[0,Pa(3,tm(3))];
plot(X,Y,'g');

X=[T];Y=[0];
plot(X,Y,'ro');
hold off
print -dps flux.ps




%Elaboration des histogrammes
clear
m=200000;s=0.3*m;c=250;a=0.001;
P=zgamma(20,m,s,c,a)-c/a;
g=zgamma(20,m,s,c,a);
p=1-normcdf(g-c/a,m,s);
T=(1-p)/p;

N=200;
t=floor(10*T);

%Simulation de N situations
clear
m=200000;s=0.1*m;c=200;a=0.001;
P=zgamma(20,m,s,c,a)-c/a;
g=zgamma(20,m,s,c,a);
p=1-normcdf(g-c/a,m,s);
T=(1-p)/p;

N=300;
t=floor(10*T);


Offres=zeros(N,t);
for i=1:N
    Offres(i,:)=m*ones(1,t)+randn(1,t)*s;
end

%Moment de dépassement de m
tm=ones(1,N);
for i=1:N
    j=1;
    while Offres(i,j)-m<0
        tm(i)=tm(i)+1;
        j=j+1;
    end
end

%Moment de dépassement du seuil P
tp=ones(1,N);
for i=1:N
    j=1;
    while Offres(i,j)-P<0 
        tp(i)=tp(i)+1;
        j=j+1;
    end
end


%Calcul du prix actualisé
Pa=zeros(N,t);
A=(1+a)^(-1)*ones(1,t);
for i=1:N
    for j=1:t
        A1=cumsum(cumprod(A(1,1:j)));
        Pa(i,j)=Offres(i,j)*(1+a)^(-j)-c*A1(1,j);
    end
end

%Prix actualisé obtenu avec la stratégie P
Sp=zeros(1,N);
for i=1:N
    Sp(1,i)=Pa(i,tp(i));
end

%Prix actualisé obtenu avec la stratégie P
Sm=zeros(1,N);
for i=1:N
    Sm(1,i)=Pa(i,tm(i));
end

%Prix actualisé maximal
M=zeros(1,N);
for i=1:N
    for j=1:t
        if Pa(i,j)>M(1,i) M(1,i)=Pa(i,j);
        end
    end
end

subplot(2,2,1)
hold on
title('Répartition pour la stratégie P');
xlabel('Pa');
hist(Sp,30);
hold off

subplot(2,2,2)
hold on
title('Répartition pour la stratégie m');
xlabel('Pa')
hist(Sm,20);
hold off

subplot(2,2,4)
hold on
title('Comparaison des stratégies P et m');
xlabel('Pa(P)-Pa(m)')
hist(Sp-Sm,30);
hold off

subplot(2,2,3)
hold on
title('Comparaison stratégie P et flux maximal');
xlabel('M-Pa(P)');
hist(M-Sp,30);
hold off

print -dps hist.ps






%Impact de la valeur de s :
clear
for i=1:50
    s(i)=i*0.01;
end
ls=length(s);
m=200000;c=200;a=0.001;

Pam=zeros(1,ls);
for i=1:ls
    Pam(i)=zgamma(20,m,s(i)*m,c,a)*(1+a)-c/a;
end


T=zeros(1,ls);
for i=1:ls
    g(i)=zgamma(20,m,s(i)*m,c,a);
    p(i)=1/2*erfc((g(i)-m-c/a)/sqrt(2*(s(i)*m)^2));
    T(i)=(1-p(i))/p(i);
end

hold on
subplot(2,1,1)
hold on
title('influence de s sur le prix actualisé moyen')
xlabel('s/m')
ylabel('Pam')
plot(s,Pam)
hold off

subplot(2,1,2)
hold on
title('influence de s sur le temps d attente moyen')
xlabel('s/m')
ylabel('T')
plot(s,T)
hold off
hold off
print -dps infs.ps



%Impact du prix des charges
clear
for i=1:2000
    c(i)=10*i;
end
lc=length(c);
m=200000;s=0.05*m;a=0.001;

Pam=zeros(1,lc);
for i=1:lc
    Pam(i)=zgamma(200,m,s,c(i),a)*(1+a)-c(i)/a;
end


T=zeros(1,lc);
for i=1:lc
    g(i)=zgamma(20,m,s,c(i),a);
    p(i)=1/2*erfc((g(i)-m-c(i)/a)/sqrt(2*s^2));
    T(i)=(1-p(i))/p(i);
end

subplot(2,1,1)
hold on
title('influence de c sur le prix actualisé moyen')
xlabel('c')
ylabel('Pam')
plot(c,Pam)
hold off

subplot(2,1,2)
hold on
title('influence de c sur le temps d attente moyen')
xlabel('c')
ylabel('T')
plot(c,T)
hold off
print -dps infc.ps


%Impact du taux de revenu
clear
for i=1:8
    a(i)=i*0.001;
end
la=length(a);
m=200000;s=0.05*m;c=200;

Pam=zeros(1,la);
for i=1:la
    Pam(i)=zgamma(20,m,s,c,a(i))*(1+a(i))-c/a(i);
end

T=zeros(1,la);
for i=1:la
    g(i)=zgamma(20,m,s,c,a(i));
    p(i)=1/2*erfc((g(i)-m-c/a(i))/sqrt(2*s^2));
    T(i)=(1-p(i))/p(i);
end

subplot(2,1,1)
hold on
title('influence de a sur le prix actualisé moyen')
xlabel('a')
ylabel('Pam')
plot(a,Pam)
hold off

subplot(2,1,2)
hold on
title('influence de a sur le temps d attente')
xlabel('a')
ylabel('T')
plot(a,T)
hold off
print -dps infa.ps