
%simulation.m

%Notations
%m=prix moyen de la maison, fixé par l'agence
%s=erreur commise sur ce prix

%a=taux de revenu de l'argent placé (par semaine)
%c=manque à gagner (par semaine)


%Calcul de "gamma"
%%On remarque qu'il suffit d'annuler la fonction suivante :

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;

%Pour cela, on va utiliser une 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;



%On peut alors calculer le prix auquel on va vendre P et le temps moyen de
%vente T :
P=zgamma(20,m,s,c,a)-c/a
%Calcul du paramètre de la loi géométrique suivie par T :
g=zgamma(20,m,s,c,a); p=1-normcdf(g-c/a,m,s); T=(1-p)/p



%Premier calcul :
clear m=200000;s=m;c=5000;a=0.006; 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
Fa=g*(1+a)-c/a




%Simulation du flux actualisé à m,s,c,a fixés (le temps d'attente est-il
%réaliste?)
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
t=floor(3*floor(T));

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

v=ones(1,3); for i=1:3
    j=1;
    while Offres(i,j)-P<0
        v(i)=v(i)+1;
        j=j+1;
    end
end

v1=ones(1,3); for i=1:3
    j=1;
    while Offres(i,j)-m<0
        v1(i)=v1(i)+1;
        j=j+1;
    end
end

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


subplot(3,1,1) hold on title('flux actualisé') xlabel('temps')
ylabel('flux actualisé') plot(F(1,:));

V1=F(1,v(1))*ones(1,t); plot(V1,'r');
X=[v(1),v(1)];Y=[0,F(1,v(1))]; plot(X,Y,'r');

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

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

subplot(3,1,2) hold on xlabel('temps') ylabel('flux actualisé')
plot(F(2,:));

V2=F(2,v(2))*ones(1,t); plot(V2,'r')
X=[v(2),v(2)];Y=[0,F(2,v(2))]; plot(X,Y,'r');

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


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

hold off

subplot(3,1,3) hold on xlabel('temps') ylabel('flux actualisé')
plot(F(3,:));

V3=F(3,v(3))*ones(1,t); plot(V3,'r')
X=[v(3),v(3)];Y=[0,F(3,v(3))]; plot(X,Y,'r');

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

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




%Elaboration des histogrammes
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

S1=zeros(1,N); for i=1:N
    S1(1,i)=Fa(i,v1(i));
end

S2=zeros(1,N); for i=1:N
    S2(1,i)=Fa(i,v(i));
end

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

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

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

subplot(2,2,4) hold on title('Comparaison des stratégies');
hist(S2-S1,30); hold off

subplot(2,2,3) hold on title('Comparaison stratégie P et flux
maximal'); hist(M-S2,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=250;a=0.001; P=zeros(1,ls); for
i=1:ls
    P(i)=zgamma(20,m,s(i)*m,c,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)=floor(1/p(i))+1;
end

hold on subplot(2,1,1) hold on title('influence de s sur le prix
de vente') xlabel('s/m') ylabel('P') plot(s,P) hold off

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



%Impact du prix des charges
clear for i=0:200
    c(i+1)=10*i;
end lc=length(c); m=200000;s=0.3*m;a=0.001; P=zeros(1,lc); for
i=1:lc
    P(i)=zgamma(20,m,s,c(i),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)=floor(1/p(i))+1;
end


subplot(2,1,1) hold on title('influence de c sur le prix de
vente') xlabel('c') ylabel('P') plot(c,P) hold off

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


%Impact du taux de revenu
clear for i=1:6
    a(i)=i*0.001;
end la=length(a); m=200000;s=0.3*m;c=250; P=zeros(1,la); for
i=1:la
    P(i)=zgamma(20,m,s,c,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)=floor(1/p(i))+1;
end

subplot(2,1,1) hold on title('influence de a sur le prix de
vente') xlabel('a') ylabel('P') plot(a,P) 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
