% Macro2 HW2 P4 Anton Cheremukhin, Paulina Restrepo Echavarria, Hisayuki Yoshimoto
clear;
%define values
c=0.1:0.1:2.0;      %i
a=[0.2 0.4];        %j
q=[1 2];            %k
P=[0.8 0.2 0.2 0.8];
w0=0:(1.1/19):1.1;

%define indicator matrices and payoffs
MAX=80;
U=zeros(1,MAX);
Ij=zeros(2,MAX);
Ik=zeros(2,MAX);
Il=zeros(4,MAX);
f=zeros(MAX,1);
for i=1:20
for j=1:2
for k=1:2
    l=2*(j-1)+k;
    m=4*(i-1)+l;
    U(m)=sqrt(c(i))*(1-a(j));
    f(m)=q(k)-c(i);
    if j==1
       Ij(1,m)=1;
    elseif j==2
       Ij(2,m)=1;
    end
    
    if k==1
       Ik(1,m)=1;
    elseif k==2
       Ik(2,m)=1;
    end
    
    if l==1
       Il(1,m)=1;
    elseif l==2
       Il(2,m)=1;
    elseif l==3
       Il(3,m)=1;
   elseif l==4
       Il(4,m)=1;
    end
end
end
end

%define matrices for constraints
Aeq=[zeros(4,MAX);ones(1,MAX)];
A=[zeros(2,MAX);-U];
for i=1:20
for j=1:2
for k=1:2
    l=2*(j-1)+k;
    m=4*(i-1)+l;
    Aeq(1,m)=(P(1)-Ik(1,m))*Ij(1,m);
    Aeq(2,m)=(P(2)-Ik(2,m))*Ij(1,m);
    Aeq(3,m)=(P(3)-Ik(1,m))*Ij(2,m);
    Aeq(4,m)=(P(4)-Ik(2,m))*Ij(2,m);

    A(1,m)=((P((3-j-1)*2+k)*U(m+6-4*j)/P((j-1)*2+k))-U(m))*Ij(1,m);
    A(2,m)=((P((3-j-1)*2+k)*U(m+6-4*j)/P((j-1)*2+k))-U(m))*Ij(2,m);
end
end
end

%have IC constraints
lb=zeros(MAX,1);
ub=ones(MAX,1);
beq=[0;0;0;0;1];
x=zeros(20,MAX);
for h=1:20
    b=[0;0;-w0(h)];
    x0=linprog(-f,A,b,Aeq,beq,lb,ub);
    x(h,:)=x0';
end

% no IC constraints
y=zeros(20,MAX);
for h=1:20
    b=[-w0(h)];
    y0=linprog(-f,A(3,:),b,Aeq,beq,lb,ub);
    y(h,:)=y0';
end
%plot
for i=1:10
    hold on;
    figure(i)
    for j=1:2;
        m=2*(i-1)+j;
        subplot(2,2,2*j-1);
        plot(c,x(m,1:4:77),'-',c , x(m,2:4:78),'--',c,x(m,3:4:79),':',c,x(m,4:4:80),'-.',...
                                                                                    'LineWidth',2)
        legend('0.2-1','0.2-2','0.4-1','0.4-2',2);
        title(['have ICC, w0=' num2str(w0(2*(i-1)+j))]);
        subplot(2,2,2*j);
        plot(c,y(m,1:4:77),'-',c,y(m,2:4:78),'--',c,y(m,3:4:79),':',c,y(m,4:4:80),'-.',...
                    'LineWidth',2);
        legend('0.2-1','0.2-2','0.4-1','0.4-2',2);
        title(['no ICC, w0=' num2str(w0(2*(i-1)+j))]);
    end
    %print
    if i==1
        hold off; print -dbmp graphc1.bmp
    elseif i==2
        hold off; print -dbmp graphc2.bmp
    elseif i==3
        hold off; print -dbmp graphc3.bmp
    elseif i==4
        hold off; print -dbmp graphc4.bmp
    elseif i==5
        hold off; print -dbmp graphc5.bmp
    elseif i==6
        hold off; print -dbmp graphc6.bmp
    elseif i==7
        hold off; print -dbmp graphc7.bmp
    elseif i==8
        hold off; print -dbmp graphc8.bmp
    elseif i==9
        hold off; print -dbmp graphc9.bmp
    elseif i==10
        hold off; print -dbmp graphc10.bmp
    end
end
%effort probabilities
prob0x=zeros(2,20);
prob0y=zeros(2,20);
for m=1:20
    prob0x(:,m)=[x(m,1:4:77)*ones(20,1)+x(m,3:4:79)*ones(20,1); x(m,2:4:78)*ones(20,1)+x(m,4:4:80)*ones(20,1)];
    prob0y(:,m)=[y(m,1:4:77)*ones(20,1)+y(m,3:4:79)*ones(20,1); y(m,2:4:78)*ones(20,1)+y(m,4:4:80)*ones(20,1)];
end
hold on;
figure(11)
subplot(2,2,1)
plot(w0,prob0x(1,:),'-',w0,prob0x(2,:),':',...
                    'LineWidth',2)
legend('low effort','high effort',2);
title('effort probabilities private');
subplot(2,2,2)
plot(w0,prob0y(1,:),'-',w0,prob0y(2,:),':',...
                    'LineWidth',2)
legend('low effort','high effort',2);
title('effort probabilities public');

%output probabilities
prob1x=zeros(2,20);
prob1y=zeros(2,20);
for m=1:20
    prob1x(:,m)=[x(m,1:4:77)*ones(20,1)+x(m,2:4:78)*ones(20,1); x(m,3:4:79)*ones(20,1)+x(m,4:4:80)*ones(20,1)];
    prob1y(:,m)=[y(m,1:4:77)*ones(20,1)+y(m,2:4:78)*ones(20,1); y(m,3:4:79)*ones(20,1)+y(m,4:4:80)*ones(20,1)];
end
subplot(2,2,3)
plot(w0,prob1x(1,:),'-',w0,prob1x(2,:),':',...
                    'LineWidth',2);
legend('low q','high q',2);
title('output probabilities private');
subplot(2,2,4)
plot(w0,prob1y(1,:),'-',w0,prob1y(2,:),':',...
                    'LineWidth',2);
legend('low q','high q',2);
title('output probabilities public');

hold off; print -dbmp graphc11.bmp

%expected utilities by states
probx=zeros(4,20);
proby=zeros(4,20);
for m=1:20
    probx(:,m)=[x(m,1:4:77)*U(1:4:77)'; x(m,2:4:78)*U(2:4:78)'; x(m,3:4:79)*U(3:4:79)'; x(m,4:4:80)*U(4:4:80)'];
    proby(:,m)=[y(m,1:4:77)*U(1:4:77)'; y(m,2:4:78)*U(2:4:78)'; y(m,3:4:79)*U(3:4:79)'; y(m,4:4:80)*U(4:4:80)'];
end
hold on;
figure(12)
subplot(2,2,1)
plot(w0,probx(1,:),'-',w0,probx(2,:),'--',w0,probx(3,:),':',w0,probx(4,:),'-.',...
                    'LineWidth',2)
legend('0.2-1','0.2-2','0.4-1','0.4-2',2)
title('expected utilities by state private');
subplot(2,2,2)
plot(w0,proby(1,:),'-',w0,proby(2,:),'--',w0,proby(3,:),':',w0,proby(4,:),'-.',...
                    'LineWidth',2)
legend('0.2-1','0.2-2','0.4-1','0.4-2',2)
title('expected utilities by state public');

%expected utilities by effort
probx=zeros(2,20);
proby=zeros(2,20);
for m=1:20
    probx(:,m)=[x(m,1:4:77)*U(1:4:77)'*0.8+x(m,2:4:79)*U(2:4:79)'*0.2; x(m,3:4:78)*U(3:4:78)'*0.2+x(m,4:4:80)*U(4:4:80)'*0.8];
    proby(:,m)=[y(m,1:4:77)*U(1:4:77)'*0.8+y(m,2:4:79)*U(2:4:79)'*0.2; y(m,3:4:78)*U(3:4:78)'*0.2+y(m,4:4:80)*U(4:4:80)'*0.8];
end
subplot(2,2,3)
plot(w0,probx(1,:),'-',w0,probx(2,:),':',...
                    'LineWidth',2)
legend('low effort','high effort',2)
title('expected utilities by effort private');
subplot(2,2,4)
plot(w0,proby(1,:),'-',w0,proby(2,:),':',...
                    'LineWidth',2)
legend('low effort','high effort',2)
title('expected utilities by effort public');
hold off; print -dbmp graphc12.bmp

%firm's profit and people's utility
hold on;
figure(13)
plot(w0,f'*x','-',w0,f'*y','--',w0,U*x',':',w0,U*y','-.',...
                    'LineWidth',2)
legend('fx private','fx public','Ux private','Ux public',1)
title('firms profit and consumers utility');
hold off; print -dbmp graphc13.bmp