%histogramme.m


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
Offres=zeros(N,t); for i=1:N
    Offres(i,:)=m*ones(1,t)+randn(1,t)*s;
end

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

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


%Calcul du flux
Fa=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)));
        Fa(i,j)=Offres(i,j)*(1+a)^(-j)-c*A1(1,j);
    end
end

%Flux actualisé pour la stratégie m
S1=zeros(1,N); for i=1:N
    S1(1,i)=Fa(i,v1(i));
end

%Flux actualisé pour la stratégie P
S2=zeros(1,N); for i=1:N
    S2(1,i)=Fa(i,v(i));
end

%Valeur maximale du flux actualisé
M=zeros(1,N); for i=1:N
    for j=1:t
        if Fa(i,j)>M(1,i) M(1,i)=Fa(i,j);
        end
    end
end
