% Likelihood analysis of all possible outcomes of election based on
% state polls from race2004.net, realclearpolitics.com, and
% davidwissing.com
% 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=153; % states other than 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 '];
polls.EV=[ 6 10 9 27 7 4 17 10 11 15 5 4 5 20 7 21 11 13 11 5 10 15];
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 Monday, October 25, 2004 6: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=[0, -1.3, -5, 999, 0, 0, 0, 0, 0, 0 ; ... %AR 1    
    -7, -11, -5, 999, 0, 0, 0, 0, 0, 0; ... %AZ 2  
    3.6, -7, -6, -5, 999, 0, 0, 0, 0, 0; ... %CO 3       
    0, -2.1, 0, 1, 0, 0, -2, -1, -1.2, 999; ... %FL 4     
    -2.3, -2, 0, -6, 1, -6, 3.2, 999, 0, 0; ... %IA 5       
    11, 6, 2, 999, 0, 0, 0, 0, 0, 0; ... %ME 6
    10, 5, 6, 2, 7, -4.2, 1, 6.7, 999, 0; ... %MI 7  
    1.2, -2, 11.1, -2, 999, 0, 0, 0, 0, 0; ... %MN 8  
    -5, -6, -3.1, 999, 0, 0, 0, 0, 0, 0; ... %MO 9 
    -8, -6, -8, -3, 999, 0, 0, 0, 0, 0; ... %NC 10   
    -3.7, -2, -1, -4, -7, -3.9, 999, 0, 0, 0; ... %NV 11   
    8.8, 2, -1, 5.1, 999, 0, 0, 0, 0, 0; ... %NH 12    
    -5.3, -5, 2, 9.5, 999, 0, 0, 0, 0, 0; ... %NM 13     
    -5, -5, 4, 1, -2, 2, -5, -3, 999, 0; ... %OH 14      
    6, 7, 1, 8, 13.1, 999, 0, 0, 0, 0; ... %OR 15        
    2, 1.9, 5, 2, 5, 1, 1, 5.7, 999, 0; ... %PA 16      
    -12, -17, -22, -2.5, 999, 0, 0, 0, 0, 0; ... %TN 17  
    -5, -4, -2.2, 999, 0, 0, 0, 0, 0, 0; ... %VA 18   
    8, 5, 10.4, 999, 0, 0, 0, 0, 0, 0; ... %WA 19  
    -2.8, -5, -2, 999, 0, 0, 0, 0, 0, 0; ... %WV 20
    -3.2, -3, -1, -6, 0, 0, 3.8, 999, 0, 0; ... %WI 21  
    7, 6, 8, 1, 8, 999, 0, 0, 0, 0]'; %NJ 22  

% Note on dealing with Maine: one way would be to create two virtual states, ME 3 with EV and ME2 with 1 EV.
% This would treat them as nominally independent even though the outcomes are linked - movement in one poll
% affects the probability of the other event.

% Find out how many polls are on each list
lastpollindex=zeros(numpolls,1)+10;
[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.
for i=1:numpolls
    polls.averagemargin(i)=mean(polls.margin(1:lastpollindex(i),i));
    polls.SEM(i)=max(std(polls.margin(1:lastpollindex(i),i))/sqrt(lastpollindex(i)),2);
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;

% Show probabilities to user.
polls.state
probs=round(polls.probkerrywin*100);
[probs(1: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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%