%Francisco M. Torralba
%June 2003
%MATLAB CODE FOR THE LUCAS-PRESCOTT MODEL OF SEARCH
%This program calculates the equilibrium theta and the equilibrium unemployment
%rate for different values of p, which measures the persistence of the shock z, 
%and for different values of N, which is the size of the labor force
clear all;
%BENCHMARK PARAMETER VALUES
alpha=0.64;            %labor share
A=1;                   %productivity
beta=0.995;            %discount factor
I=2000;                %number of grid points of labor 
J=7;                   %number of grid points of the shock
sigma=1;               %stdev of log(z)
mu=-(sigma^2)/2;       %mean of log(z)
N=100;                 %size of the labor force
epsilon=0.0001;        %precision level in the search of v
epsilon2=0.5;          %precision level in the search of theta
T=200000;              %Number of simulations
M=(J-1)/2;
P=zeros(99);
for y=1:99;
    P(y)=y/100;
end
for l=1:99;
    p=P(l)
%PARAMETERIZATION OF THE PROBLEM
for i=1:J;             %These loops generate the elements of z and Q, where
for j=1:J;             %z is the vector of shocks and Q is the transition matrix
    z(j)=exp(mu+((4*sigma)/(J-1))*(j-M-1));
    if i==j ;
        Q(i,j)=p;
    else 
        Q(i,j)=(1-p)/(2*M);
    end
end
end
theta0=z(J)*A*alpha*(N^(alpha-1))/(1-beta);     %Theta0 is the initial guess for theta, which 
                                                %we will use to define the range of x      
thetahigh=theta0;                               %The loop for the search of theta starts here  
thetalow=0;

D(l)=0;
while abs(D(l)-N)>epsilon2;
    thetahigh;
    thetalow;
theta=(thetahigh+thetalow)/2;
x(1)=(1.01*(1-beta)*theta/(z(1)*A*alpha))^(1/(alpha-1));
x(I)=(0.99*(1-beta)*theta/(z(J)*A*alpha))^(1/(alpha-1));
step=(x(I)-x(1))/(I-1);                        %step is the increment of x 
for  q=1:I;                                    %This loop generates the x grid    
    x(q)=step*(q-1)+x(1);
end
%COMPUTATIONAL STEP 1: Find the fixed point for v(theta)
%Creating the initial matrix of values v0
for r=1:I;
    for c=1:J;
    f(r,c)=z(c)*A*alpha*(x(r)^(alpha-1));    
    v0(r,c)=max([theta f(r,c)+beta*theta]);        
    end
end
%Iterating on v until convergence
vnew=v0;                                       %vnew is v(n+1)
diff=1;                                        %diff will collect the differences between the      
                                               %elements of two successive iterations on the matrix v
while diff>epsilon ;               
    v=vnew;                                    
for r=1:I;
    for c=1:J;  
        vnew(r,c)=max([theta f(r,c)+min([theta beta*v(r,:)*(Q(c,:)')])]);
    end
end
diff=max(max(abs(vnew-v))');   %max(max(diff)) ensures that the biggest difference across 
end                            %elements in the matrix is smaller than epsilon when we converge  
vconv=vnew;
%COMPUTATIONAL STEP 2: Find gmax(theta),gmin(theta) and g(theta)for r=1:I;
for r=1:I;
      for c=1:J;           
          k2(r,c)=abs(beta*vconv(r,:)*(Q(c,:)')-theta);          %This comes from equation (1)
          k3(r,c)=abs(f(r,c)+beta*vconv(r,:)*(Q(c,:)')-theta);   %This comes from equation (2) 
      end 
  end
           
for c=1:J;
    [k2min,a(c)]=min(k2(:,c));    %This recovers gmin and gmax from the x vector           
    [k3min,b(c)]=min(k3(:,c));
    gmin(c)=x(a(c));
    gmax(c)=x(b(c));                                            
end
for r=1:I;
    for c=1:J;
        g(r,c)=min([gmax(c) max([x(r) gmin(c)])]);
    end
end
%COMPUTATIONAL STEP 3: Find the labor size consistent with theta
%Simulation of zs and xs (simulated time series)
a=zeros(T,1);
zs=zeros(T,1);
xs=zeros(T,1);
xs(1)=(x(1)+x(I))/2;                         
zs(1)=z(M+1);
a(1)=M+1;
for t=2:T;
    w(t)=rand(1);
    wnew(t)=rand(1);
    for h=1:a(t-1)-1;
    if w(t)<=p;
        zs(t)=zs(t-1);
        a(t)=a(t-1);
    elseif (wnew(t)>=((h-1)/(2*M))) & (wnew(t)<=(h/(2*M)));
            zs(t)=z(h);
            a(t)=h;
        end
    end        
    for h=a(t-1)+1:J;
    if w(t)<=p;
       zs(t)=zs(t-1);
       a(t)=a(t-1);
   elseif (wnew(t)>=(h-2)/(2*M)) & (wnew(t)<=((h-1)/(2*M)));
           zs(t)=z(h);
           a(t)=h;
       end
   end
xs(t)=min([gmax(a(t-1)) max([xs(t-1) gmin(a(t-1))])]);
us(t-1)=max([(xs(t-1)-xs(t)) 0]);
end
D(l)=mean(xs);
E(l)=mean(us);
urate(l)=E(l)/N;
%The updating of theta in each state is done according to the method of bisection
if D(l)-N>0;
 thetalow=theta;
else
    thetahigh=theta;
end                       
end                    %The loop for the search of theta ends here    
thetastar(l)=theta
end     
thetastar
E
urate
for j=1:N;                          %Using the fact that theta(N*)=theta(N)*((N*/N)^(alpha-1)),
    d(j)=(1/(N+1-j))^(alpha-1);     %I compute thetastar for N=1 to N=100. 
end
d;
B=d'*thetastar;                     %Here I use the fact that the unemployment rate does not  
C=ones(100,1)*urate;                %change with N 
subplot(1,2,1);
surf(C);
axis([1 100 1 100 0 .3])
hold on;
subplot(1,2,2);
surf(B);
axis([1 100 1 100 0 300])
hold off;