% R(t,d,g) = prob that team 1 wins if in regular (scrimmage) mode
% E(t,d,g) = prob that team 1 wins 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
% g indicates possession; g=1 if team 1 has the ball, g=2 if team 2 has the ball

% probability that favored team 1 wins in OT
pot=0.5;

spread=7;      % point spread for favored team 1
totpoints=42;   % expected total points
theta=0.6;      % expectd ratio of FGs to TDs

% probabilities of a TD or FG on a drive, for the two teams
ptd1=(spread+totpoints)/(24*(7+3*theta));
ptd2=(ptd1*12*(7+3*theta)-spread)/(12*(7+3*theta));
pfg1=ptd1*theta;
pfg2=ptd2*theta;

% low and high clock usage if using hurry-up offense
thurrylow=6;
thurryhigh=18;
% low and high clock usage if using slow-down offense
tslowlow=9;
tslowhigh=27;

% our probabilities for 1- and 2-point conversions
p1pt1=0.985;
p2pt1=0.4;
% opponent's probabilities for 1- and 2-point conversions
p1pt2=0.985;
p2pt2=0.4;

maxd=100;
T=361;
numg=2;

dtop=2*maxd+1;  

% Initialize arrays

R=zeros(T,dtop,numg);
E=zeros(T,dtop,numg);
Oneortwo=zeros(T,dtop,numg);
Dohurry=zeros(T,dtop,numg);

% Boundary conditions for g=1

for g=1:2
   for d=1:maxd
      R(1,d,g)=0;
   end
   for d=maxd+2:dtop
      R(1,d,g)=1;
   end
   R(1,maxd+1,g)=pot;
end

for d=1:maxd-2
   E(1,d,1)=0;
end
E(1,maxd-1,1)=p2pt1*pot;
E(1,maxd,1)=max(p1pt1*pot,p2pt1);
E(1,maxd+1,1)=p1pt1+(1-p1pt1)*pot;
for d=maxd+2:dtop
   E(1,d,1)=1;
end

for d=1:maxd
   E(1,d,2)=0;
end
E(1,maxd+1,2)=(1-p1pt2)*pot;
E(1,maxd+2,2)=min(1-p1pt2+p1pt2*pot,1-p2pt2);
E(1,maxd+3,2)=1-p2pt2+p2pt2*pot;
for d=maxd+4:dtop
   E(1,d,2)=1;
end

% backward induction

for t=2:T

   for d=1:dtop
      gofor1=p1pt1*R(t-1,min(d+1,dtop),2)+(1-p1pt1)*R(t-1,d,2);
      gofor2=p2pt1*R(t-1,min(d+2,dtop),2)+(1-p2pt1)*R(t-1,d,2);
      E(t,d,1)=max(gofor1,gofor2);
      if gofor2>gofor1
         Oneortwo(t,d,1)=2;
      else
         Oneortwo(t,d,1)=1;
      end
   end
   
   for d=1:dtop
      gofor1=p1pt2*R(t-1,max(d-1,1),1)+(1-p1pt2)*R(t-1,d,1);
      gofor2=p2pt2*R(t-1,max(d-2,1),1)+(1-p2pt2)*R(t-1,d,1);
      E(t,d,2)=min(gofor1,gofor2);
      if gofor2<gofor1
         Oneortwo(t,d,2)=2;
      else
         Oneortwo(t,d,2)=1;
      end
   end
  
   for d=1:dtop
      
%                                              Insufficient time even for hurry-up offense
      if thurrylow > t-1
         R(t,d,1)=R(1,d,1);
         R(t,d,2)=R(1,d,2);
%                                              Only time for hurry offense, if low time outcome
      elseif tslowlow > t-1
         rvalue1=0.5*(ptd1*E(t-thurrylow,min(d+6,dtop),1)+pfg1*R(t-thurrylow,min(d+3,dtop),2) ...
            +(1-ptd1-pfg1)*R(t-thurrylow,d,2))+0.5*R(1,d,2);
         rvalue2=R(1,d,2);
         R(t,d,1)=max(rvalue1,rvalue2);
         if rvalue1 > rvalue2
            Dohurry(t,d,1)=1;
         else
            Dohurry(t,d,1)=0;
         end
         rvalue1=0.5*(ptd2*E(t-thurrylow,max(d-6,1),2)+pfg2*R(t-thurrylow,max(d-3,1),1) ...
            +(1-ptd2-pfg2)*R(t-thurrylow,d,1))+0.5*R(1,d,1);
         rvalue2=R(1,d,1);
         R(t,d,2)=min(rvalue1,rvalue2);
         if rvalue1 < rvalue2
            Dohurry(t,d,2)=1;
         else
            Dohurry(t,d,2)=0;
         end
%                                              Time for either offense if low time outcome
      elseif thurryhigh > t-1
         rvalue1=0.5*(ptd1*E(t-thurrylow,min(d+6,dtop),1)+pfg1*R(t-thurrylow,min(d+3,dtop),2) ...
            +(1-ptd1-pfg1)*R(t-thurrylow,d,2))+0.5*R(1,d,2);
         rvalue2=0.5*(ptd1*E(t-tslowlow,min(d+6,dtop),1)+pfg1*R(t-tslowlow,min(d+3,dtop),2) ...
            +(1-ptd1-pfg1)*R(t-tslowlow,d,2))+0.5*R(1,d,2);
         R(t,d,1)=max(rvalue1,rvalue2);
         if rvalue1 > rvalue2
            Dohurry(t,d,1)=1;
         else
            Dohurry(t,d,1)=0;
         end
         rvalue1=0.5*(ptd2*E(t-thurrylow,max(d-6,1),2)+pfg2*R(t-thurrylow,max(d-3,1),1) ...
            +(1-ptd2-pfg2)*R(t-thurrylow,d,1))+0.5*R(1,d,1);
         rvalue2=0.5*(ptd2*E(t-tslowlow,max(d-6,1),2)+pfg2*R(t-tslowlow,max(d-3,1),1) ...
            +(1-ptd2-pfg2)*R(t-tslowlow,d,1))+0.5*R(1,d,1);
         R(t,d,2)=min(rvalue1,rvalue2);
         if rvalue1 < rvalue2
            Dohurry(t,d,2)=1;
         else
            Dohurry(t,d,2)=0;
         end
%                                              Time for hurry up, and for slowdown if low time outcome
      elseif tslowhigh > t-1
         rvalue1=0.5*(ptd1*E(t-thurrylow,min(d+6,dtop),1)+pfg1*R(t-thurrylow,min(d+3,dtop),2) ...
            +(1-ptd1-pfg1)*R(t-thurrylow,d,2))+0.5*(ptd1*E(t-thurryhigh,min(d+6,dtop),1) ...
            +pfg1*R(t-thurryhigh,min(d+3,dtop),2)+(1-ptd1-pfg1)*R(t-thurryhigh,d,2));
         rvalue2=0.5*(ptd1*E(t-tslowlow,min(d+6,dtop),1)+pfg1*R(t-tslowlow,min(d+3,dtop),2) ...
            +(1-ptd1-pfg1)*R(t-tslowlow,d,2))+0.5*R(1,d,2);
         R(t,d,1)=max(rvalue1,rvalue2);
         if rvalue1 > rvalue2
            Dohurry(t,d,1)=1;
         else
            Dohurry(t,d,1)=0;
         end
         rvalue1=0.5*(ptd2*E(t-thurrylow,max(d-6,1),2)+pfg2*R(t-thurrylow,max(d-3,1),1) ...
            +(1-ptd2-pfg2)*R(t-thurrylow,d,1))+0.5*(ptd2*E(t-thurryhigh,max(d-6,1),2) ...
            +pfg2*R(t-thurryhigh,max(d-3,1),1)+(1-ptd2-pfg2)*R(t-thurryhigh,d,1));
         rvalue2=0.5*(ptd2*E(t-tslowlow,max(d-6,1),2)+pfg2*R(t-tslowlow,max(d-3,1),1) ...
            +(1-ptd2-pfg2)*R(t-tslowlow,d,1))+0.5*R(1,d,1);
         R(t,d,2)=min(rvalue1,rvalue2);
         if rvalue1 < rvalue2
            Dohurry(t,d,2)=1;
         else
            Dohurry(t,d,2)=0;
         end
%                                                 Time for either offense
      else
         rvalue1=0.5*(ptd1*E(t-thurrylow,min(d+6,dtop),1)+pfg1*R(t-thurrylow,min(d+3,dtop),2) ...
            +(1-ptd1-pfg1)*R(t-thurrylow,d,2))+0.5*(ptd1*E(t-thurryhigh,min(d+6,dtop),1) ...
            +pfg1*R(t-thurryhigh,min(d+3,dtop),2)+(1-ptd1-pfg1)*R(t-thurryhigh,d,2));
         rvalue2=0.5*(ptd1*E(t-tslowlow,min(d+6,dtop),1)+pfg1*R(t-tslowlow,min(d+3,dtop),2) ...
            +(1-ptd1-pfg1)*R(t-tslowlow,d,2))+0.5*(ptd1*E(t-tslowhigh,min(d+6,dtop),1) ...
            +pfg1*R(t-tslowhigh,min(d+3,dtop),2)+(1-ptd1-pfg1)*R(t-tslowhigh,d,2));
         R(t,d,1)=max(rvalue1,rvalue2);
         if rvalue1 > rvalue2
            Dohurry(t,d,1)=1;
         else
            Dohurry(t,d,1)=0;
         end
         rvalue1=0.5*(ptd2*E(t-thurrylow,max(d-6,1),2)+pfg2*R(t-thurrylow,max(d-3,1),1) ...
            +(1-ptd2-pfg2)*R(t-thurrylow,d,1))+0.5*(ptd2*E(t-thurryhigh,max(d-6,1),2) ...
            +pfg2*R(t-thurryhigh,max(d-3,1),1)+(1-ptd2-pfg2)*R(t-thurryhigh,d,1));
         rvalue2=0.5*(ptd2*E(t-tslowlow,max(d-6,1),2)+pfg2*R(t-tslowlow,max(d-3,1),1) ...
            +(1-ptd2-pfg2)*R(t-tslowlow,d,1))+0.5*(ptd2*E(t-tslowhigh,max(d-6,1),2) ...
            +pfg2*R(t-tslowhigh,max(d-3,1),1)+(1-ptd2-pfg2)*R(t-tslowhigh,d,1));
         R(t,d,2)=min(rvalue1,rvalue2);
         if rvalue1 < rvalue2
            Dohurry(t,d,2)=1;
         else
            Dohurry(t,d,2)=0;
         end
      end
      
   end  
end