clear

%%%%%%%%%%%%%%
% USER INPUT %
%%%%%%%%%%%%%%
I = 1;   % number of private goods markets
J = 1;   % number of public goods markets

% private goods market parameters
a = [ 150 ]';
b = [ 3 ]';
g = [ 30 ]';
d = [ 1 ]';

% public goods market parameters
A = [ 160 ]';
B = [ 8 ]';
G = [ 48 ]';
D = [ 0 ]';

Lmax = 2;   density = 3001;
Lgrid = linspace(0,Lmax,density)';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SEARCH FOR EQUILIBRIUM LAMBDA %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

surplus = zeros(density,1);   revenue = zeros(density,1);   expenditure = zeros(density,1);
for Lindex = 1:density
    L = Lgrid(Lindex);
    shadow = 1 + L;
    
    % private goods markets
    R = zeros(I,1);   t = zeros(I,1);   x = zeros(I,1);   
    for i = 1:I
        t(i) = L*( a(i) - g(i) ) / ( 2*L + 1 );
        x(i) = ( a(i) - g(i) - t(i) ) / ( b(i) + d(i) );
        R(i) = t(i) * x(i);
    end
    sumR = sum(R);
    
    % public goods markets
    y = zeros(J,1);   E = zeros(J,1);
    for j = 1:J
        y(j) = ( A(j) - shadow*G(j) ) / ( B(j) + shadow*D(j) );
        E(j) = G(j) * y(j) + .5 * D(j) * y(j)^2;
    end
    sumE = sum(E);
    
    revenue(Lindex) = sumR;
    expenditure(Lindex) = sumE;
    surplus(Lindex) = sumR - sumE;
end
[num,Leqmindex] = min( abs(surplus) );
L = Lgrid(Leqmindex);
shadow = L + 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CALCULATIONS USING EQUILIBRIUM LAMBDA %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% private goods markets
R = zeros(I,1);   t = zeros(I,1);   x = zeros(I,1);   
for i = 1:I
    t(i) = L*( a(i) - g(i) ) / ( 2*L + 1 );
    x(i) = ( a(i) - g(i) - t(i) ) / ( b(i) + d(i) );
    R(i) = t(i) * x(i);
end
sumR = sum(R);

% public goods markets
y = zeros(J,1);   E = zeros(J,1);
for j = 1:J
    y(j) = ( A(j) - shadow*G(j) ) / ( B(j) + shadow*D(j) );
    E(j) = G(j) * y(j) + .5 * D(j) * y(j)^2;
end
sumE = sum(E);

if I == J
    report = [ a b g d A B G D x y t R E ]';
end