%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%			Mortensen-Pissarides 
%
%
%   George J Hall
%   Brandeis University
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
! rm mp.out
diary mp.out; 
disp('Mortensen-Pissarides');
disp('');
t1 = clock;
%
%  set parameter values
%
c      = 0.05;            % cost to posting vacancy
beta   = 0.98;            % discount factor 
delta  = 0.50;            % sharing parameter
A      = 0.20;            % match parameter
alpha  = 0.50;            % match parameter
Ustart = 27.33;            % initial value for U
g      = 0.95;            % relaxation parameter
%
%   form productivity grid
%   
maxY = 0.8;                         % maximum value of Y grid   
incY = .8/40;                       % size of Y grid increments
minY = .8/40;                       % minimum value of Y grid
NY   = round((maxY-minY)/incY+1);   % number of grid points
ygrid = [minY:incY:maxY];
%
%   form transition probability matrix
%
prob = zeros(NY,NY);
for i = 1:NY;
  for j = 1:NY;
    if j <= i; prob(i,j) = (j/i)^10; end;
    if j >  i; prob(i,j) = (1 -.5*((j-i)/(NY-i)))^10; end
  end
end
psum = sum(prob');
for i = 1:NY;
  prob(i,:) = prob(i,:)/psum(i);
end
%
%   loop to find fixed point for U
%
liter   = 1;
maxiter = 100;
toler   = 1e-8;
metric  = 10;
U = Ustart;
disp('ITERATING ON U');
disp('');
disp('    liter     metric     Uprime      Uold');
while  (metric > toler) & (liter <= maxiter);
  %
  %  initialize some variables
  %
  S       = zeros(NY,1);
  TS      = zeros(NY,1);
  decis   = zeros(NY,1);
  test    = 10;
  %
  %  iterate on Bellman's equation and get the decision 
  %  rules and the value function at the optimum         
  %
  while test > 1e-8;
    for i=1:NY;
      y = minY + (i-1)*incY;
      [TS(i), decis(i)]=max([y-(1-beta)*U+beta*prob(i,:)*S, 0 ]);
    end;
    test=max((TS-S)./(S+eps));
    S=TS;
  end;
  %  
  %  compute new value of U
  %
  J = (1-delta)*S;
  W = delta*S + U;
  q = c/J(NY);
  voveru = (q/A)^(inv(-alpha));
  Uprime = voveru*q*W(NY) + beta*(1-voveru*q)*U;
  %
  %  check for convergence
  %
  Unew = g*Uprime + (1-g)*U;
  metric = abs((Uprime-U)/U);
  disp([ liter metric Uprime U ]);
  U = Unew;
  liter = liter+1;
end;
%
%  compute distributions
%
I = (decis == 1);
u = 0.05;
psi = inv(sum(I))*I;
while (abs(sum(psi) + u - 1.0) > 1e-8)
  psi = (1-u)*inv(sum(I))*I;
  test=1;
  while test > 1e-8;
    psi1       = (prob'*psi).*I;
    gamma      = (prob'*psi).*(ones(NY,1)-I);
    rho        = sum(gamma);
    u          = (u+rho)*(1-voveru*q);
    v          = voveru*u;
    psi1(NY,1) = prob(:,NY)'*(psi.*I) + v*q + rho*q*(voveru);
    test = max(abs(psi1-psi));
    psi  = psi1;
  end;
end;
%
%   print out results
%
disp('PARAMETER VALUES');
disp('');
disp('       c       beta      delta       A      alpha'); 
disp([ c beta delta A alpha ]);
disp(''); 
disp('EQUILIBRIUM RESULTS ');
disp('');
disp('      U         u        v      q');
disp([ U u  v q ]);
%
%   plot some graphs
%
figure(1)
plot(psi)
print psi.ps