% 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=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 Monday, November 1, 2004 11:00PM EST.
% 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 and other rolling polls at a time 
%    (because of overlapping sample problem).
 polls.margin=[-3.2, 0, -8, -1.5, 999, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0; ... %AR 1    
    -15, -5, -5, 999, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0; ... %AZ 2  
    -1.7, -6, -5.2, -7, -9, 1, 999, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0; ... %CO 3       
    0, 5, -4, -3, -3, 2, -4, 0, -4, -2, 1.6, -3, -8, 3, 999, 0, 0; ... %FL 4     
    4.9, -4, 0, -2, -1, 3, -5, 0, 1, -3, -1, -1, -4, 999, 0, 0, 0; ... %IA 5       
    8, 11, 6, 999, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0; ... %ME 6
    5.4, 3, 4, 2, -2, 4, 0, 5, 999, 0, 0, 0, 0, 0, 0, 0, 0; ... %MI 7  
    6.5, 1, 8, 6, 8, -1, -2, -7, -3, 999, 0, 0, 0, 0, 0, 0, 0; ... %MN 8  
    -4.2, -5, -5, -5, -3, -4, -7, 999, 0, 0, 0, 0, 0, 0, 0, 0, 0; ... %MO 9 
    -8, -11, -6, -9, 999, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0; ... %NC 10   
    -4.7, -6, 0, -1, -2, 999, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0; ... %NV 11   
    4.5, 0, 1, 2, 3, 1, 8.8, 999, 0, 0, 0, 0, 0, 0, 0, 0, 0; ... %NH 12    
    2.6, 1, -3, -4, -9, -4, -5.3, 999, 0, 0, 0, 0, 0, 0, 0, 0, 0; ... %NM 13     
    -5.6, -4, -3, -2, -0.9, -4, 0, -2, -1, -1, 1.6, 6, 3, 2, -1, -1, -5; ... %OH 14      
    7, 6, 10.1, 6, 3, 4, -2, 999, 0, 0, 0, 0, 0, 0, 0, 0, 0; ... %OR 15        
    3.7, 2, 4, 1, -4, 3, 2, 5.3, 1, 0, -2, 0, 3, 8, 3, 3, 999; ... %PA 16      
    -3.3, -12, -17, 999, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0; ... %TN 17  
    -5, -4, -9, -6, 999, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0; ... %VA 18   
    11.2, 9, 6, 6, 6, 999, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0; ... %WA 19  
    -8, -8, -2.8, 999, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0; ... %WV 20
    5.5, -3,  -8, 7.6, 2, 5, -3, -3, 1, -5, -3.2, 999, 0, 0, 0, 0, 0; ... %WI 21  
    0, 4, 12, 7, 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, 0, 0, 0, 0]' %HI 23  

% 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)+17;
[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
%    original calculation of aerage
    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);
% switch to median and estimated SD.
%    polls.medianmargin(i)=median(polls.margin(1:lastpollindex(i),i));
%    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);
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)=1;

% 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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%