% This is a matlab file for Theorem 8 in Sugaya(2008)
clc
clear all

global a b r
a=2/3;
b=2/3;
r=2;

beta = .8;
lambda = .2;
c = 3;

gridv = .015;
V = gridv:gridv:1;
[nv kv] = size(V);
VV = zeros(1,kv);
Theta = zeros(1,kv);
U = zeros(1,kv);

maxx = 3;
maxitertheta = 100000;
maxitervalue = 1000;
maxiterh     = 1000;
grid = .005;
X = 0:grid:maxx;
[n k] = size(X);

F  = zeros(1,k);
for iF = 1:k
        F(iF) = k^2-iF^2;
end
F = F/sum(F);

G  = zeros(1,k);
for iG = 1:k
        G(iG) = k^2-iG^2;
end
G = G/sum(G);

for vi=1:kv
    v = V(vi);
    h0u = .25;
    h1u = .25;
    H1e = .25*G;
    H2e = .25*G;
    theta0 = theta(.5);
    Hnotbuystay = H1e;
    Hnotbuyquit = zeros(1,k);
    for i=1:maxitertheta
        V0u = zeros(1,k);
        V1u = zeros(1,k);
        V1e = zeros(1,k);
        V2e = zeros(1,k);
        for j=1:maxitervalue
            V0u1 = beta*theta0*(max(V1e-X,V0u)*F')*ones(1,k)+beta*(1-theta0)*V0u;
            V1u1 = v*(c*ones(1,k)+beta*theta0*(max(V1e-X,V0u)*F')*ones(1,k)+beta*(1-theta0)*V0u)+...
                (1-v)*(beta*theta0*(max(V2e-X,V1u)*F')*ones(1,k)+beta*(1-theta0)*V1u);        
            V1e1 = v*(c*ones(1,k)+beta*lambda*(max(V1e-X,V0u)*G')*ones(1,k)+...
                beta*(1-lambda)*(max(V1e-X,V0u)))+...
                (1-v)*(beta*lambda*(max(V2e-X,V1u)*G')*ones(1,k)+...
                beta*(1-lambda)*(max(V2e-X,V1u)));
            V2e1 = c*ones(1,k)+beta*lambda*(max(V2e-X,V1u)*G')*ones(1,k)+...
                beta*(1-lambda)*(max(V2e-X,V1u));
            distance = max(max(max(abs(V0u1-V0u),abs(V1u1-V1u)),...
                max(abs(V1e1-V1e),abs(V2e1-V2e))));
            if distance<.001
                break
            else
                V0u=V0u1;
                V1u=V1u1;
                V1e=V1e1;
                V2e=V2e1;
            end
            if j == maxitervalue
                disp('not converge')
            end            
        end

        [R0 r0] = min(abs(V1e-X-V0u));
        [R1 r1] = min(abs(V2e-X-V1u));
        F0      = [F(1:r0) zeros(1,k-r0)];
        F1      = [F(1:r1) zeros(1,k-r1)];      
        G0      = [G(1:r0) zeros(1,k-r0)];
        G1      = [G(1:r1) zeros(1,k-r1)];
        h0u1 = v*lambda*(1-sum(G0))*sum(H1e)+... % h1e buy and quit
            (1-theta0)*h0u+theta0*(1-sum(F0))*h0u+... % h0u cannot find
            v*(1-theta0)*h1u+v*theta0*(1-sum(F0))*h1u; % h1u buy and cannot find
        h1u1 = (1-v)*((1-lambda)*sum(Hnotbuyquit)+lambda*(1-sum(G1))*sum(H1e))+... %h1e not buy and quit
            lambda*(1-sum(G1))*sum(H2e)+... %h2e buy and quit
            (1-v)*(1-theta0)*h1u+(1-v)*theta0*(1-sum(F1))*h1u; %h1u not buy and find
        H1e1 = v*(1-lambda)*H1e + v*lambda*G0*sum(H1e) +... % h1e buy and stay
            v*theta0*F0*h1u + theta0*F0*h0u; % h1u buy and find and h0u find 
        H2e1 = (1-v)*(1-lambda)*Hnotbuystay + (1-v)*lambda*G1*sum(H1e)+... % h1e not buy and stay
            (1-lambda)*H2e + lambda*G1*sum(H2e) +... % h2e buy and stay
            (1-v)*theta0*F1*h1u; % h1u not buy and find
        U1 = h0u1+h1u1;
        theta1 = theta(U1);
        if max(max(abs(h1u1-h1u),abs(h0u1-h0u)),max(max(abs(sum(H2e1)-sum(H2e))),max(abs(sum(H1e1)-sum(H1e)))))<0.0001
            VV(vi)=H1e*V1e'+H2e*V2e'+h0u*V0u(1,1)+h1u*V1u(1,1);
            Theta(vi)=theta0;
            U(vi)=U1;
            break
        else
            h0u=h0u1;
            h1u=h1u1;
            H1e = H1e1;
            H2e = H2e1;
            theta0=theta1;
            Hnotbuystay = [H1e(1:r1) zeros(1,k-r1)];
            Hnotbuyquit = [zeros(1,r1) H1e(r1+1:k)];
        end
        if i == maxitertheta
            disp('distribution does not converge')
            vi
            VV(vi)=H1e*V1e'+H2e*V2e'+h0u*V0u(1,1)+h1u*V1u(1,1);
            Theta(vi)=theta0;
            U(vi)=U1;
        end 
    end
end
figure (1)
plot(V,VV)
title('the value')
xlabel('excess demand')
ylabel('value')

figure (2)
plot(V,Theta)
title('\theta')
xlabel('excess demand')
ylabel('\theta')

figure (3)
plot(V,U)
title('Unemployment')
xlabel('excess demand')
ylabel('U')