%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Mortensen (1998) avec capital humain 
%
%  François Langot -- 08/2002
%
%   simplification r = 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all

% paramétrisation Mortensen 2002
s     = 0.287;  % sorties (e->u) + (e->e) mortensen p.54  E(w)=497
s     = 0.3189; % les femmes E(w)=4773
lamb  = 0.142;  % mortensen p.54
psi   = 0.5;

b     = 0.3;
alpha = .6;

%%%%%%%%%%%%%%%%%%%
%
%  équilibre
%
%%%%%%%%%%%%%%%%%%%

v = lamb^(1/psi);

kmin = alpha/(s+lamb);
p    = kmin^(-alpha);

kmax = (p*alpha/s)^(1/(1-alpha));

ymin = 1;
ymax = p*(kmax^alpha);

wmax = ymax - s*kmax - ((s/(s+lamb))^2)*(ymin - b - kmin*(s+lamb));

gamma = (lamb/v)*(s/(s+lamb))*((ymin - b - kmin*(s+lamb))/(s+lamb));

incw   = (wmax-b)/50;                 
gridw  = b:incw:wmax;	    
nw     = length(gridw); 

Fw = zeros(1,nw);
fw = zeros(1,nw);
Gw = zeros(1,nw);
gw = zeros(1,nw);

for i=1:nw
    w      = gridw(i);
    % amorce pour alpha =.1 b=.7
    Fw0    = 1;    
    k0     = .55;
    % amorce pour alpha =.3 b=.7
    Fw0    = 1;    
    k0     = 10;
    % amorce pour alpha =.3 b=.3
    Fw0    = 1;    
    k0     = 10;    
    P      = [s,lamb,alpha,ymin,b,p,w,kmin];
    % amorce pour alpha =.1 b=.3
    X      = fsolve('distwkh',[Fw0 k0],optimset('fsolve'),P);  
    Fw(i)  = X(1);
    k(i)   = X(2);
end

for i=1:nw
    Gw(i) = (s*Fw(i))/(s + lamb*(1-Fw(i)));
end

fw(1) = Fw(1);
gw(1) = Gw(1);
for j=2:nw
    fw(j)  = Fw(j)-Fw(j-1);
    gw(j)  = Gw(j)-Gw(j-1);    
end

w = gridw;
debg = 1;

hw = (lamb/v).*(s./(s+lamb*(1-Fw)));
dw = s+lamb*(1-Fw);

%
% Graphiques
%
subplot(221),plot(w(1:end),hw,'-',w(1:end),dw,'--')
legend('h(w)','d(w)')
xlabel('w')
ylabel('h(w),d(w)')

subplot(222),plot(w(1:end),(p.*(k.^alpha)-w(1:end))./w(1:end))
legend('(pf(k)-w)/w')
xlabel('w')
ylabel('(pf(k)-w)/w')

subplot(223),plot(w(1:end),k(1:end),'-')
legend('k(w)')
xlabel('w')
ylabel('k(w)')

subplot(224),plot(w(2:end),gw(2:end),'-')
legend('g(w)')
xlabel('w')
ylabel('g(w)')


pause
print mortks0
close

plot(w(debg+1:end),fw(debg+1:end),'-o',w(debg+1:end),gw(debg+1:end),'-*')
grid on;
axis([b 1.05*max(w) 0 1.05*max([gw(debg:end) fw(debg:end)])]);
title('DISTRIBUTIONS : DENSITES');
ylabel('offres, salaires');
xlabel('w');
legend('f(w)','g(w)',2)
pause
print mortks1
close

plot(gridw(debg:end),Fw(debg:end),'-o',gridw(debg:end),Gw(debg:end),'-*')
grid on;
axis([b 1.05*max(w) 0 1.05*max([Gw(debg:end) Fw(debg:end)])]);
title('DISTRIBUTIONS : REPARTITIONS');
ylabel('offres, salaires');
xlabel('w');
legend('F(w)','G(w)',2)
pause
print mortks2
close

% sauvegarde
Fw0  = Fw;
Gw0  = Gw;
w0   = w;
fw0  = fw;
gw0  = gw;
hw0  = hw;
dw0  = dw;
k00  = k;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% variante
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

b     = 0.28;

alpha = 0.61;

%%%%%%%%%%%%%%%%%%%
%
%  équilibre
%
%%%%%%%%%%%%%%%%%%%

kmin = alpha/(s+lamb);
p    = kmin^(-alpha);

kmax = (p*alpha/s)^(1/(1-alpha));

ymin = 1;
ymax = p*(kmax^alpha);

v    = ((1/gamma)*(s/(s+lamb))*((ymin - b - kmin*(s+lamb))/(s+lamb)))^(1/(1-psi));
lamb = v^psi;

wmax = ymax - s*kmax - ((s/(s+lamb))^2)*(ymin - b - kmin*(s+lamb));

incw   = (wmax-b)/50;                 
gridw  = b:incw:wmax;	    
nw     = length(gridw); 

Fw = zeros(1,nw);
fw = zeros(1,nw);
Gw = zeros(1,nw);
gw = zeros(1,nw);

for i=1:nw
    w      = gridw(i);
    % amorce pour alpha =.1 b=.7
    Fw0    = 1;    
    k0     = .55;
    % amorce pour alpha =.3 b=.7
    Fw0    = 1;    
    k0     = 10;
    % amorce pour alpha =.3 b=.3
    Fw0    = 1;    
    k0     = 10;    
    P      = [s,lamb,alpha,ymin,b,p,w,kmin];
    % amorce pour alpha =.1 b=.3
    X      = fsolve('distwkh',[Fw0 k0],optimset('fsolve'),P);  
    Fw(i)  = X(1);
    k(i)   = X(2);
end

for i=1:nw
    Gw(i) = (s*Fw(i))/(s + lamb*(1-Fw(i)));
end

fw(1) = Fw(1);
gw(1) = Gw(1);
for j=2:nw
    fw(j)  = Fw(j)-Fw(j-1);
    gw(j)  = Gw(j)-Gw(j-1);    
end

w = gridw;
debg = 1;

hw = (lamb/v).*(s./(s+lamb*(1-Fw)));
dw = s+lamb*(1-Fw);

%
% Graphiques
%
subplot(221),plot(w(1:end),hw,'-o',w(1:end),dw,'-*',w0(1:end),hw0,'-',w0(1:end),dw0,'--')
legend('h(w)','d(w)','h0(w)','d0(w)')
xlabel('w')
ylabel('h(w),d(w)')

subplot(222),plot(w(1:end),(p.*(k.^alpha)-w(1:end))./w(1:end),'-o',w0(1:end),(p.*(k00.^alpha)-w0(1:end))./w0(1:end),'--')
legend('(pf(k)-w)/w','(pf(k0)-w0)/w0')
xlabel('w')
ylabel('(pf(k)-w)/w')

subplot(223),plot(w(1:end),k(1:end),'-o',w0(1:end),k00(1:end),'--')
legend('k(w)','k0(w)')
xlabel('w')
ylabel('k(w)')

subplot(224),plot(w(2:end),gw(2:end),'-o',w0(2:end),gw0(2:end),'--')
legend('g(w)','g0(w)')
xlabel('w')
ylabel('g(w)')


pause
print mortksvar
close

plot(w(debg+1:end),fw(debg+1:end),'-o',w(debg+1:end),gw(debg+1:end),'-*')
grid on;
axis([b 1.05*max(w) 0 1.05*max([gw(debg:end) fw(debg:end)])]);
title('DISTRIBUTIONS : DENSITES');
ylabel('offres, salaires');
xlabel('w');
legend('f(w)','g(w)',2)
pause
close

plot(gridw(debg:end),Fw(debg:end),'-o',gridw(debg:end),Gw(debg:end),'-*')
grid on;
axis([b 1.05*max(w) 0 1.05*max([Gw(debg:end) Fw(debg:end)])]);
title('DISTRIBUTIONS : REPARTITIONS');
ylabel('offres, salaires');
xlabel('w');
legend('F(w)','G(w)',2)
pause
close

