% Likelihood analysis of all possible outcomes of election based on
% State-by-state polls from http://www.realclearpolitics.com/bush_vs_kerry_sbys.html
% Calculations by Sam Wang, Princeton University.
%
% The point of this calculation is to take into account the variability of polls.
% It uses all data for the last N polls. Typically N=3.

bushsafe=142; kerrysafe=149; % states other than HI, NJ, CO, NC, TN, VA, and the RealClearPolitics Battleground States list
polls.state=['AR AZ CO FL IA ME MI MN MO NC NV NH NM OH OR PA TN VA WA WV WI NJ HI '];
polls.EV=[ 6 10 9 27 7 4 17 10 11 15 5 4 5 20 7 21 11 13 11 5 10 15 4];
numpolls=size(polls.EV,2);
uncertainvotes=sum(polls.EV);

% checksum to make sure no double assignment or missed assignment
if (kerrysafe+bushsafe+uncertainvotes~=538)
    warning('Electoral votes do not sum to 538!')
    [kerrysafe bushsafe uncertainvotes]
end

% Data as of Friday, October 29, 2004 4:00PM EDT.
% 999 indicates no more polls for that state
% The rule is: last three polls, or all polls in last seven days, whichever is greater.
% Second rule: only one poll from Rasmussen at a time (because of overlapping samples).
 polls.margin=[-7, 0, -1.3, 999, 0, 0, 0, 0, 0, 0, 0, 0, 0 ; ... %AR 1    
    -5, -5, -7, -11, 999, 0, 0, 0, 0, 0, 0, 0, 0; ... %AZ 2  
    4.1, 1, 0, 999, 0, 0, 0, 0, 0, 0, 0, 0, 0; ... %CO 3       
    -3, -2.7, -3, -8, 3, -2, 2, -8, 0, 999, 0, 0, 0; ... %FL 4     
    -1, -0.4, -4, 999, 0, 0, 0, 0, 0, 0, 0, 0, 0; ... %IA 5       
    11, 6, 2, 999, 0, 0, 0, 0, 0, 0, 0, 0, 0; ... %ME 6
    2, 4, -0.6, 0, 5, 10, 999, 0, 0, 0, 0, 0, 0; ... %MI 7  
    0, -2, 2.6, 7, -3, 1.2, 999, 0, 0, 0, 0, 0, 0; ... %MN 8  
    -3, -4, -7, 999, 0, 0, 0, 0, 0, 0, 0, 0, 0; ... %MO 9 
    -9, -10, -8, -6, 999, 0, 0, 0, 0, 0, 0, 0, 0; ... %NC 10   
    -7.3, -2, -2, 999, 0, 0, 0, 0, 0, 0, 0, 0, 0; ... %NV 11   
    4, 8.8, 2, -1, 5.1, 999, 0, 0, 0, 0, 0, 0, 0; ... %NH 12    
    -3.2, -4, -5.3, 999, 0, 0, 0, 0, 0, 0, 0, 0, 0; ... %NM 13     
    -3, 1.6, 6, 3, 2, -1, -1, -5, 999, 0, 0, 0, 0; ... %OH 14      
    6, 3, 4, -2, 999, 0, 0, 0, 0, 0, 0, 0, 0; ... %OR 15        
    3, 2.8, -2, 3, 0, 8, 3, 3, 2, 2, 2, 999, 0; ... %PA 16      
    -12, -17, -22, 999, 0, 0, 0, 0, 0, 0, 0, 0, 0; ... %TN 17  
    -6, -5, -4, 999, 0, 0, 0, 0, 0, 0, 0, 0, 0; ... %VA 18   
    6, 6, 6, 999, 0, 0, 0, 0, 0, 0, 0, 0, 0; ... %WA 19  
    -2.8, -5, -2, 999, 0, 0, 0, 0, 0, 0, 0, 0, 0; ... %WV 20
    1, -4, -5, -3.2, 999, 0, 0, 0, 0, 0, 0, 0, 0; ... %WI 21  
    0, 0, 7, 6, 8, 999, 0, 0, 0, 0, 0, 0, 0; ... %NJ 22  
    4, -0.6, -0.7, 999, 0, 0, 0, 0, 0, 0, 0, 0, 0]' %HI 23  

% Find out how many polls are on each list
lastpollindex=zeros(numpolls,1)+12;
[foo foo2]=find(polls.margin==(999));
lastpollindex(foo2)=foo-1;
%lastpollindex=zeros(numpolls,1)+3; % force consideration of last three polls only

% Calculate mean and SEM of margin, all polls weighted equally.
% THIS INCLUDES THE UNDECIDED ALLOCATION OF 2.5 +/- 2.0 (mean +/- SD)
% OCTOBER 28 REVISED TO 1.0 +/- 0.8 (mean +/- SD)
for i=1:numpolls
    polls.averagemargin(i)=mean(polls.margin(1:lastpollindex(i),i))+1.0;
    polls.SEM(i)=max(std(polls.margin(1:lastpollindex(i),i))/sqrt(lastpollindex(i)),2);
% switch to median and estimated SD.
%    polls.medianmargin(i)=median(polls.margin(1:lastpollindex(i),i))+1.0;
%    polls.mad(i)=median(abs(polls.margin(1:lastpollindex(i),i)-polls.medianmargin(i)));
%    polls.SEM(i)=max(polls.mad(i)/0.6745/sqrt(lastpollindex(i)),2);
    polls.SEM(i)=sqrt(polls.SEM(i)*polls.SEM(i)+0.8*0.8);
end

% Calculate z-score and convert to probability, assuming normal distribution.
polls.z=polls.averagemargin./polls.SEM;
polls.probkerrywin=(erf(polls.z/sqrt(2))+1)/2;
polls.probbushwin=1-polls.probkerrywin;

hold on;
xlabel('Electoral votes');
ylabel('Probability of Kerry getting fewer than # votes');
plot([270 270],[0 1],'r-')

%%%%%%%%%%%%%%%%%%%%%%%%%%% The calculation %%%%%%%%%%%%%%%%%%%%%%%%
%%%% Can assign bias and then cut-and-paste this to get result %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate z-score and convert to probability, assuming normal distribution.
polls.z=(polls.averagemargin+bias)./polls.SEM;
polls.probkerrywin=(erf(polls.z/sqrt(2))+1)/2;
polls.probbushwin=1-polls.probkerrywin;

polls.probkerrywin(4)=0;
polls.probkerrywin(14)=0;

% Show probabilities to user.
polls.state
probs=round(polls.probkerrywin*100);
[probs(1:4) probs(23) probs(5:10) probs(22) probs(11:21)]

EVpoly=[polls.probkerrywin(1) zeros(1, polls.EV(1)-1) 1-polls.probkerrywin(1)];
for i=2:numpolls
    nextEV=[polls.probkerrywin(i) zeros(1, polls.EV(i)-1) 1-polls.probkerrywin(i)];
    EVpoly=conv(EVpoly,nextEV);
end
histogram=fliplr(EVpoly);

% Cumulative histogram of all possibilities
cumulative_prob=cumsum(histogram);
electoralvotes=kerrysafe:kerrysafe+uncertainvotes;
plot(electoralvotes,cumulative_prob,'k-')

% Calculate statistics from cumulative histogram
confidence_low=electoralvotes(max(find(cumulative_prob<0.025))); % 95-pct lower limit
KerrymedianEV=electoralvotes(min(find(cumulative_prob>0.5)));  % 50-pct outcome
confidence_high=electoralvotes(min(find(cumulative_prob>0.975))); % 95-pct upper limit
probability_Kerry_win=1-cumulative_prob(min(find(electoralvotes>=269)));
BushmedianEV=538-KerrymedianEV;

% Read out results
bias
[KerrymedianEV BushmedianEV probability_Kerry_win]
[confidence_low confidence_high]

%%%%%%%%%%%%%%% Display probabilities in rank order %%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[sortedpolls,ix]=sort(1-polls.probkerrywin); %sort states by decreasing probability

% Kerry winning: greater than 95% probability
wins95=find(sortedpolls<=0.05);
if wins95
    foo=[];
    for i=1:max(wins95)
        foo=[foo ix(i)*3-2 ix(i)*3-1 ix(i)*3];
    end
    Kerry_wins95 = polls.state(foo);
    Kerry_wins95 = strrep(Kerry_wins95,' ','/');
    Kerry_wins95 = Kerry_wins95(1:size(Kerry_wins95,2)-1)
else
    wins95=0;    
end

% Kerry likely winning: 80-94.99% probability
wins80_95=find(sortedpolls<=0.20);
if (max(wins80_95)>max(wins95))
    foo=[];
    for i=(max(wins95)+1):max(wins80_95)
        foo=[foo ix(i)*3-2 ix(i)*3-1 ix(i)*3];
    end
    Kerry_wins80_95 = polls.state(foo);
    Kerry_wins80_95 = strrep(Kerry_wins80_95,' ',' / ');
    Kerry_wins80_95 = Kerry_wins80_95(1:size(Kerry_wins80_95,2)-2)
end

% Kerry win probability 20-79.99%
wins20_80=find(sortedpolls<=0.80);
if (max(wins20_80)>max(wins80_95))
    foo=[];
    for i=(max(wins80_95)+1):max(wins20_80)
        foo=[foo ix(i)*3-2 ix(i)*3-1 ix(i)*3];
    end
    Kerry_wins20_80 = polls.state(foo);
    Kerry_wins20_80 = strrep(Kerry_wins20_80,' ',' / ');
    Kerry_wins20_80 = Kerry_wins20_80(1:size(Kerry_wins20_80,2)-2)
end

% Bush likely winning: Kerry to win at 5-19.99% probability
wins5_20=find(sortedpolls<=0.95);
if (max(wins5_20)>max(wins20_80))
    foo=[];
    for i=(max(wins20_80)+1):max(wins5_20)
        foo=[foo ix(i)*3-2 ix(i)*3-1 ix(i)*3];
    end
    Kerry_wins5_20 = polls.state(foo);
    Kerry_wins5_20 = strrep(Kerry_wins5_20,' ',' / ');
    Kerry_wins5_20 = Kerry_wins5_20(1:size(Kerry_wins5_20,2)-2)
end

% Bush certain wins: Kerry win probability less than 5%
if (size(ix,2)>max(wins5_20))
    foo=[];
    for i=(max(wins5_20)+1):size(ix,2)
        foo=[foo ix(i)*3-2 ix(i)*3-1 ix(i)*3];
    end
    Bush_wins95 = polls.state(foo);
    Bush_wins95 = strrep(Bush_wins95,' ','/');
    Bush_wins95 = Bush_wins95(1:size(Bush_wins95,2)-1)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
