% Chapter 3 Simulations A %
% Jamus Jerome Lim %

--------------------------------------------------------------------------------------------------------------------------------------------------------

function [fval] = effort(e)

global xhat xbar xi ebar J;
fval = diag(xi(1,1:J))*0.5*(1/ebar)*(((e(1:J,1)/sum(e(1:J,1)))'*eye(J)*xhat(1,1:J)' - xhat(1,1:J)').^2 + (xbar*ones(J,1) - xhat(1,1:J)').^2) - e(1:J,1);

--------------------------------------------------------------------------------------------------------------------------------------------------------

clear;

global xhat xbar xi ebar J;

% Initially homogeneous groups %

xhatJ = 2;

while xhatJ < 101; % Loop to vary preferred policy of treatment
      xhat = [2 xhatJ];
      xbar = 1;
      xi = [0.1 0.1];
      ebar = 10;
      J = 2;
      row = xhatJ-1;
      subrow = (row-1)*4+1;

      xiJ = 0.1;

      while xiJ < 1.01; % Loop to vary size of shock
            xi = [0.1 xiJ];
            e = broyden('effort',[0.2;0.2]);
            K = J - 1;
            J = K; % K used as J to simplify computation of etilde
            etildeJ = [broyden('effort',[0.2]); 0];
            xi = [xiJ];
            etilde1 = [0; broyden('effort',[0.2])];
            xi = [0.1 xiJ];
            etilde = [etilde1 etildeJ];
            J = J + 1;
            column = xiJ*100-9;
            eJ(round(row),round(column)) = e(J);
            e1(round(row),round(column)) = e(1);
            PC = -((sum(e(1:J,1))/ebar)*(((e(1:J,1)/sum(e))'*xhat(1,1:J)'*ones(J,1) - xhat(1,1:J)').^2 - (xbar*ones(J,1) - xhat(1,1:J)').^2) - diag(sum(etilde)./ebar)*(((etilde*diag(1./sum(etilde)))'*xhat(1,1:J)' - xhat(1,1:J)').^2 - (xbar*ones(J,1) - xhat(1,1:J)').^2) + e(1:J,1).^2./xi(1,1:J)'); % Participate iff > 0
            PC1(round(row),round(column)) = PC(1);
            PCJ(round(row),round(column)) = PC(J);
            H(round(subrow):round(subrow+1),round(column)) = sign(PC); % H < 0 => Do not participate
            h(round(row),round(column)) = 0;
            for i = 0:1
                if H((round(subrow)+i),round(column)) > 0
                   h(round(row),round(column)) = h(round(row),round(column)) + 1;
                end;
            end;
            xiJ = xiJ + 0.01;
      end;
      xhatJ = xhatJ + 1;
end;

% Control group with divergent preferred policy %

xhat = [2 20];
xbar = 19;
xi = [0.1 0.1];
ebar = 10;
J = 2;
row = 1;
subrow = 1;

xiJ = 0.1;

while xiJ < 1.01; % Loop to vary size of shock
    xi = [0.1 xiJ];
    e = broyden('effort',[0.2;0.2]);
    K = J - 1;
    J = K; % K used as J to simplify computation of etilde
    etildeJ = [broyden('effort',[0.2]); 0];
    xi = [xiJ];
    etilde1 = [0; broyden('effort',[0.2])];
    xi = [0.1 xiJ];
    etilde = [etilde1 etildeJ];
    J = J + 1;
    column = xiJ*100-9;
    eprimeJ(round(row),round(column)) = e(J);
    eprime1(round(row),round(column)) = e(1);
    PC = -((sum(e(1:J,1))/ebar)*(((e(1:J,1)/sum(e))'*xhat(1,1:J)'*ones(J,1) - xhat(1,1:J)').^2 - (xbar*ones(J,1) - xhat(1,1:J)').^2) - diag(sum(etilde)./ebar)*(((etilde*diag(1./sum(etilde)))'*xhat(1,1:J)' - xhat(1,1:J)').^2 - (xbar*ones(J,1) - xhat(1,1:J)').^2) + e(1:J,1).^2./xi(1,1:J)'); % Participate iff > 0
    PCprime1(round(row),round(column)) = PC(1);
    PCprimeJ(round(row),round(column)) = PC(J);
    Hprime(round(subrow):round(subrow+1),round(column)) = sign(PC); % H < 0 => Do not participate
    hprime(round(row),round(column)) = 0;
    for i = 0:1
        if Hprime((round(subrow)+i),round(column)) > 0
        hprime(round(row),round(column)) = hprime(round(row),round(column)) + 1;
        end;
    end;
    xiJ = xiJ + 0.01;
end;

% Control group with initially larger shocks and stronger preferred policy %

xhat = [2.5 2];
xbar = 1;
xi = [0.5 0.1];
ebar = 10;
J = 2;
row = 1;
subrow = 1;

xiJ = 0.1;

while xiJ < 1.01; % Loop to vary size of shock
    xi = [0.5 xiJ];
    e = broyden('effort',[0.2;0.2]);
    K = J - 1;
    J = K; % K used as J to simplify computation of etilde
    etildeJ = [broyden('effort',[0.2]); 0];
    xi = [xiJ];
    etilde1 = [0; broyden('effort',[0.2])];
    xi = [0.5 xiJ];
    etilde = [etilde1 etildeJ];
    J = J + 1;
    column = xiJ*100-9;
    eprimeprimeJ(round(row),round(column)) = e(J);
    eprimeprime1(round(row),round(column)) = e(1);
    PC = -((sum(e(1:J,1))/ebar)*(((e(1:J,1)/sum(e))'*xhat(1,1:J)'*ones(J,1) - xhat(1,1:J)').^2 - (xbar*ones(J,1) - xhat(1,1:J)').^2) - diag(sum(etilde)./ebar)*(((etilde*diag(1./sum(etilde)))'*xhat(1,1:J)' - xhat(1,1:J)').^2 - (xbar*ones(J,1) - xhat(1,1:J)').^2) + e(1:J,1).^2./xi(1,1:J)'); % Participate iff > 0
    PCprimeprime1(round(row),round(column)) = PC(1);
    PCprimeprimeJ(round(row),round(column)) = PC(J);
    Hprimeprime(round(subrow):round(subrow+1),round(column)) = sign(PC); % H < 0 => Do not participate
    hprimeprime(round(row),round(column)) = 0;
    for i = 0:1
        if Hprimeprime((round(subrow)+i),round(column)) > 0
        hprimeprime(round(row),round(column)) = hprimeprime(round(row),round(column)) + 1;
        end;
    end;
    xiJ = xiJ + 0.01;
end;

% Graph participation constraints with shock %

%{
title('Participation decisions by magnitude of shock');
xlabel('Magnitude of shock');
ylabel('Participation constraint');
hold on;
plot(shock,PC1(1,:),shock,PCJ(1,:),shock,0,'r-'); % xhatJ = 3
text('Interpreter','latex','String','PC_J (\hat{x}_J = 3)','Position',[0.35 0.0025]);
text('Interpreter','latex','String','PC_1 (\hat{x}_1 = 2)','Position',[0.5 0.002]);
plot(shock,PC1(J,:),shock,PCJ(J,:),shock,0,'r:'); % xhatJ = 6
text('Interpreter','latex','String','PC_J (\hat{x}_J = 6)','Position',[0.3 0.007]);
text('Interpreter','latex','String','PC_1 (\hat{x}_1 = 6)','Position',[0.5 0.0025]);
plot(shock,PC1(9,:),shock,PCJ(9,:),shock,0,'r:'); % xhatJ = 11
text('Interpreter','latex','String','PC_J (\hat{x}_J = 9)','Position',[0.7 0.01]);
text('Interpreter','latex','String','PC_1 (x_J = 9)','Position',[0.5 -0.0275]);
hold off;
%}

% Graph effort, participation constraints, and movement size with shock %

shock = [0.1:0.01:1];
figure(1);
set(gcf,'Position',[30 60 1025 850]);

% xhatJ = 2 %
subplot(4,3,1), plot(shock,e1(1,:),'--',shock,eJ(1,:),'linewidth',1.5);
title({'Optimal effort levels by magnitude of shock';''});
xlabel('Magnitude of shock','FontSize',8);
ylabel('Optimal effort','FontSize',8);
legend('e_1','e_2',2,'FontSize',8);
legend boxoff
subplot(4,3,2), plot(shock,PC1(1,:),'--',shock,PCJ(1,:),shock,0,'r:','linewidth',1.5);
title({'Participation decisions by magnitude of shock';''});
xlabel('Magnitude of shock','FontSize',8);
ylabel('Participation constraint','FontSize',8);
legend('PC_1','PC_2','Threshold',2,'FontSize',8);
legend boxoff
subplot(4,3,3), bar(shock,h(1,:));
axis([0,1,0,3]);
colormap(gray);
title({'Movement formation by magnitude of shock';''});
xlabel('Magnitude of shock','FontSize',8);
ylabel('Number of groups','FontSize',8);
%{
% xhatJ = 3 %
subplot(4,3,4), plot(shock,e1(2,:),'--',shock,eJ(2,:),'linewidth',1.5);
xlabel('Magnitude of shock','FontSize',8);
ylabel('Optimal effort','FontSize',8);
subplot(4,3,5), plot(shock,PC1(2,:),'--',shock,PCJ(2,:),shock,0,'r:','linewidth',1.5);
title({'Participation decisions by magnitude of shock';''});
xlabel('Magnitude of shock','FontSize',8);
ylabel('Participation constraint','FontSize',8);
subplot(4,3,6), bar(shock,h(2,:));
axis([0,1,0,3]);
colormap(gray);
xlabel('Magnitude of shock','FontSize',8);
ylabel('Number of groups','FontSize',8);
%}
% xhatJ = 4 %
subplot(4,3,4), plot(shock,e1(3,:),'--',shock,eJ(3,:),'linewidth',1.5);
xlabel('Magnitude of shock','FontSize',8);
ylabel('Optimal effort','FontSize',8);
subplot(4,3,5), plot(shock,PC1(3,:),'--',shock,PCJ(3,:),shock,0,'r:','linewidth',1.5);
xlabel('Magnitude of shock','FontSize',8);
ylabel('Participation constraint','FontSize',8);
subplot(4,3,6), bar(shock,h(3,:));
axis([0,1,0,3]);
colormap(gray);
xlabel('Magnitude of shock','FontSize',8);
ylabel('Number of groups','FontSize',8);
%{
% xhatJ = 20 %
subplot(4,3,4), plot(shock,e1(19,:),'--',shock,eJ(19,:),'linewidth',1.5);
xlabel('Magnitude of shock','FontSize',8);
ylabel('Optimal effort','FontSize',8);
subplot(4,3,5), plot(shock,PC1(19,:),'--',shock,PCJ(19,:),shock,0,'r:','linewidth',1.5);
xlabel('Magnitude of shock','FontSize',8);
ylabel('Participation constraint','FontSize',8);
subplot(4,3,6), bar(shock,h(19,:));
axis([0,1,0,3]);
colormap(gray);
xlabel('Magnitude of shock','FontSize',8);
ylabel('Number of groups','FontSize',8);
% xhatJ = 40 %
subplot(4,3,7), plot(shock,e1(39,:),'--',shock,eJ(39,:),shock,0,'r:','linewidth',1.5);
xlabel('Magnitude of shock','FontSize',8);
ylabel('Optimal effort','FontSize',8);
subplot(4,3,8), plot(shock,PC1(39,:),'--',shock,PCJ(39,:),shock,0,'r:','linewidth',1.5);
xlabel('Magnitude of shock','FontSize',8);
ylabel('Participation constraint','FontSize',8);
subplot(4,3,9), bar(shock,h(39,:));
axis([0,1,0,3]);
colormap(gray);
xlabel('Magnitude of shock','FontSize',8);
ylabel('Number of groups','FontSize',8);
% xhatJ = 100 %
subplot(4,3,7), plot(shock,e1(99,:),'--',shock,eJ(99,:),'linewidth',1.5);
xlabel('Magnitude of shock','FontSize',8);
ylabel('Optimal effort','FontSize',8);
subplot(4,3,8), plot(shock,PC1(99,:),'--',shock,PCJ(99,:),shock,0,'r:','linewidth',1.5);
xlabel('Magnitude of shock','FontSize',8);
ylabel('Participation constraint','FontSize',8);
subplot(4,3,9), bar(shock,h(99,:));
axis([0,1,0,3]);
colormap(gray);
xlabel('Magnitude of shock','FontSize',8);
ylabel('Number of groups','FontSize',8);
%}
% xi1 = 0.5 %
subplot(4,3,7), plot(shock,eprimeprime1(1,:),'--',shock,eprimeprimeJ(1,:),'linewidth',1.5);
xlabel('Magnitude of shock','FontSize',8);
ylabel('Optimal effort','FontSize',8);
subplot(4,3,8), plot(shock,PCprimeprime1(1,:),'--',shock,PCprimeprimeJ(1,:),shock,0,'r:','linewidth',1.5);
xlabel('Magnitude of shock','FontSize',8);
ylabel('Participation constraint','FontSize',8);
subplot(4,3,9), bar(shock,hprimeprime(1,:));
axis([0,1,0,3]);
colormap(gray);
xlabel('Magnitude of shock','FontSize',8);
ylabel('Number of groups','FontSize',8);
% xhatJ = 20, xbar = 19 %
subplot(4,3,10), plot(shock,eprime1(1,:),'--',shock,eprimeJ(1,:),'linewidth',1.5);
xlabel('Magnitude of shock','FontSize',8);
ylabel('Optimal effort','FontSize',8);
subplot(4,3,11), plot(shock,PCprime1(1,:),'--',shock,PCprimeJ(1,:),shock,0,'r:','linewidth',1.5);
xlabel('Magnitude of shock','FontSize',8);
ylabel('Participation constraint','FontSize',8);
subplot(4,3,12), bar(shock,hprime(1,:));
axis([0,1,0,3]);
colormap(gray);
xlabel('Magnitude of shock','FontSize',8);
ylabel('Number of groups','FontSize',8);

print shocks.eps

% Graph participation constraints with homogeneity and variation of shock and policy %

figure(2);
set(gcf,'Position',[30 502 1025 418]);
title('Participation decisions by magnitude of shock and intensity of policy preferences');
colormap default;
subplot(1,2,1), surf(PCJ);
shading interp;
xlabel('Magnitude of shock','FontSize',8);
ylabel('Intensity of policy preference','FontSize',8);
zlabel('Participation constraint','FontSize',8);
subplot(1,2,2), surf(PC1);
xlabel('Magnitude of shock','FontSize',8);
ylabel('Intensity of policy preference','FontSize',8);
zlabel('Participation constraint','FontSize',8);
shading interp;
view(-37.5,10);
print participate.eps