%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 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
lamb  = 0.143;  % mortensen p.54
psi   = .5;

b     = 0.3;
alpha = 0.4;

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

v  = lamb^(1/psi);
hr = (lamb/v)*(s/(s+lamb));

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

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

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

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

wmax = ymax + gamma - s*((s+(lamb/v))/(lamb/v))*kmax - ...
       hr*((s+(lamb/v))/(lamb/v))*((ymin - b + gamma)/(s+lamb+hr))  + kmin*s*((s+(lamb/v))/(lamb/v));

if wmax<b
    disp('probleme wmax<b')
    pause
end
   
valr = hr*((ymin-b+gamma)/(s+lamb+hr)) - s*kmin;
   
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);
    Fw0    = .01;    
    k0     = 1;    
    P      = [valr,s,lamb,alpha,p,w,v,gamma];
    X      = fsolve('distwkg',[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)')
pause
print mortkg0
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 mortkg1
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 mortkg2
close


