% R(t,d,y,g) = prob that we win if in regular (scrimmage) mode % K(t,d,g) = prob that we win if in kickoff mode % E(t,d,g) = prob that we win if in extra-point mode % t is time period, with t=1 the final period % d indexes point differential; maxd is largest differential, and d runs 1:2*maxd+1 % y is yardage marker; there are numy of them % g indicates possession; g=1 if we have the ball, g=2 if they have the ball % if onside kick, success with probability ponside % Copyright 2004 by William S. Krasker % Takes 30 minutes to run when T=181, numy=6, maxd=100, and tdrivemax=50. % We require that the probability of winning in OT be 1/2, to justify symmetry. pot=0.5; p1pt=0.985; p2pt=0.4; ponside=0.17; timeperiodspermin=6; numy=6; maxd=100; T=181; numg=2; % For the hurry-up offense, the scoring probabilities are decreased by the factor % hurryimpairment. There is additional reduction in the probabilities if the time % remaining in the game is less than hurryclockusemin miutes. hurryimpairment=0.8; hurryclockusemin=2.0; hurryclockuse=round(hurryclockusemin*timeperiodspermin); dtop=2*maxd+1; % The following three matrices quantify the state transitions. The transition % probabilities are given by ptrans. The conditional means and standard deviations % for elapsed time (in minutes) are given by tmeanminutes and tstddevminutes. % These are converted to mean and standard deviation for elapsed time in time % periods, tmean and tstddev. % 1 2 3 4 5 6 TD FG TDOPP ptrans=[0.03 0.2 0.27 0.2 0.08 0.02 0.11 0.07 0.02; ... 0.05 0.26 0.26 0.09 0.05 0.01 0.16 0.1 0.02; ... 0.08 0.32 0.16 0.05 0.02 0.01 0.22 0.13 0.01; ... 0.07 0.22 0.13 0.01 0.01 0.01 0.27 0.27 0.01; ... 0.04 0.11 0.04 0 0 0 0.4 0.4 0.01; ... 0.02 0.02 0 0 0 0 0.77 0.18 0.01 ]; % 1 2 3 4 5 6 TD FG TDOPP tmeanminutes=[4.5 3.6 2.5 1.8 1.9 1.8 4.7 5.5 2.0; ... 3.8 2.9 2.1 1.8 1.2 1.2 3.8 4.9 2.0; ... 2.6 2.3 2 1.3 0.9 1 3.6 4 2.0; ... 2 1.8 1.6 1.6 1.6 1.6 2.5 2.9 1.6; ... 1.8 1.4 1.4 1.4 1.4 1.4 1.4 1.9 1.5; ... 1.5 1.3 1.2 1.1 1.0 0.9 0.7 1 0.8 ]; % 1 2 3 4 5 6 TD FG TDOPP tstddevminutes=[1.4 1.6 1.3 1 1 1.0 2.7 1.8 1.5; ... 1.7 1.5 1.3 1.1 1 0.8 2.1 2.1 1.3; ... 1.3 1.2 1.3 0.9 0.7 0.7 1.8 1.6 1.3; ... 1.2 1 1 0.9 0.9 0.9 1.4 1.2 1.2; ... 1 1.1 0.9 0.9 0.9 0.9 1.1 0.9 1.0; ... 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.6 0.6 ]; tmean=tmeanminutes*timeperiodspermin; tstddev=tstddevminutes*timeperiodspermin; % The next array gives transition probabilities when we kick off. % They get ball We get ball % ----------------------------------- --------------------------------------- % 1 2 3 4 5 6 TD 1 2 3 4 5 6 TD ptranskickoff=[.01 .64 .30 .02 .01 .005 .005 .001 .001 .001 .001 .002 .002 .002 ]; % We now construct the conditional probabilities for elapsed time on a drive, % assuming the time is lognormal. tdrivemax=50; tprobs=zeros(numy,numy+3,tdrivemax); for i=1:numy for j=1:numy+3 [mu, sigma]=lognorparams(tmean(i,j),tstddev(i,j)); tprobtemp=lognordensity(mu,sigma,tdrivemax); tprobs(i,j,:)=tprobtemp; end end % Initialize arrays R=zeros(T,dtop,numy,numg); K=zeros(T,dtop,numg); E=zeros(T,dtop,numg); Oneortwo=zeros(T,dtop); Doonside=zeros(T,dtop); Dohurry=zeros(T,dtop,numy); % Boundary conditions for g=1 for y=1:numy for d=1:maxd R(1,d,y,1)=0; end for d=maxd+2:dtop R(1,d,y,1)=1; end R(1,maxd+1,y,1)=pot; end for d=1:maxd K(1,d,1)=0; end for d=maxd+2:dtop K(1,d,1)=1; end K(1,maxd+1,1)=pot; for d=1:maxd-2 E(1,d,1)=0; end E(1,maxd-1,1)=p2pt*pot; E(1,maxd,1)=max(p1pt*pot,p2pt); E(1,maxd+1,1)=p1pt+(1-p1pt)*pot; for d=maxd+2:dtop E(1,d,1)=1; end % fill in boundary conditions for g=2 by symmetry for d=1:dtop K(1,d,2)=1-K(1,dtop+1-d,1); E(1,d,2)=1-E(1,dtop+1-d,1); for y=1:numy R(1,d,y,2)=1-R(1,dtop+1-d,y,1); end end % backward induction for t=2:T for d=1:dtop kickaway=ptranskickoff(numy+1)*E(t-1,max(d-6,1),2)+ptranskickoff(2*numy+2)*E(t-1,min(d+6,dtop),1); for y2=1:numy kickaway=kickaway+ptranskickoff(y2)*R(t-1,d,y2,2)+ptranskickoff(y2+numy+1)*R(t-1,d,y2,1); end kickonside=ponside*R(t-1,d,3,1)+(1-ponside)*R(t-1,d,4,2); K(t,d,1)=max(kickaway,kickonside); if kickonside>kickaway Doonside(t,d)=1; else Doonside(t,d)=0; end end for d=1:dtop gofor1=p1pt*K(t,min(d+1,dtop),1)+(1-p1pt)*K(t,d,1); gofor2=p2pt*K(t,min(d+2,dtop),1)+(1-p2pt)*K(t,d,1); E(t,d,1)=max(gofor1,gofor2); if gofor2>gofor1 Oneortwo(t,d)=2; else Oneortwo(t,d)=1; end end for d=1:dtop for y=1:numy rvalue1=0; for k=1:min(t-1,tdrivemax) rvalue1=rvalue1+ptrans(y,numy+1)*tprobs(y,numy+1,k)*E(t-k,min(d+6,dtop),1) ... +ptrans(y,numy+2)*tprobs(y,numy+2,k)*K(t-k,min(d+3,dtop),1) ... +ptrans(y,numy+3)*tprobs(y,numy+3,k)*E(t-k,max(d-6,1),2); for y2=1:numy rvalue1=rvalue1+ptrans(y,y2)*tprobs(y,y2,k)*R(t-k,d,y2,2); end end for k=min(t-1,tdrivemax)+1:tdrivemax rvalue1=rvalue1+ptrans(y,numy+1)*tprobs(y,numy+1,k)*R(1,d,y,2) ... +ptrans(y,numy+2)*tprobs(y,numy+2,k)*R(1,d,y,2) ... +ptrans(y,numy+3)*tprobs(y,numy+3,k)*R(1,d,y,2); for y2=1:numy rvalue1=rvalue1+ptrans(y,y2)*tprobs(y,y2,k)*R(1,d,y2,2); end end elapsedtime=min(hurryclockuse,round(t/2)); thurry=t-elapsedtime; totalimpairment=hurryimpairment*min((t-1)/hurryclockuse,1); rvalue2=totalimpairment*ptrans(y,numy+1)*E(thurry,min(d+6,dtop),1) ... +totalimpairment*ptrans(y,numy+2)*K(thurry,min(d+3,dtop),1) ... +(1-totalimpairment*ptrans(y,numy+1)-totalimpairment*ptrans(y,numy+2)) ... *R(thurry,d,numy+1-y,2); R(t,d,y,1)=max(rvalue1,rvalue2); if rvalue2>rvalue1 Dohurry(t,d,y)=1; else Dohurry(t,d,y)=0; end end end for d=1:dtop K(t,d,2)=1-K(t,dtop+1-d,1); E(t,d,2)=1-E(t,dtop+1-d,1); for y=1:numy R(t,d,y,2)=1-R(t,dtop+1-d,y,1); end end end % Fill in the array needed for boundary condition if we want to do the first half. B=zeros(dtop,numg); for g=1:2 for d=1:dtop B(d,g)=K(T,d,g); end end