% Copyright 2005 William S. Krasker

% W gives the probability (of scoring or winning) at the moment the whistle blows, if clock is running.
% S gives the probability (of scoring or winning) at the moment the ball is snapped for a 
%   real play (not a spike). This is also the win probability when the whistle blows if the clock stops.

% scorefactor=0  if offense trails by 4 points or more. The outputs of the model are then
%    the probability of scoring a touchdown.
% scorefactor=1  if offense trails by 2 points or less. The outputs of the model are then
%    the probability of winning.
% scorefactor=.5 if offense trails by 3 points. The outputs of the model are then
%    the probability of winning.
scorefactor=1;

options = optimset('Display','off');

% probability of making a field goal from each field position.
fgprob=[.99 .98 .93 .85 .8 .74 .63 .5 .05 .01 zeros(1,11)];

spiketime=10;        % Seconds required, after whistle blows, to spike the ball.
timeouttime=2;       % Seconds required, after whistle blows, to call timeout.
letclockruntime=16;  % Seconds required, after whistle blows, to run another real play.

probob=.2;           % Probability that runner gets out of bounds after a completed pass.

maxgain=10;          % This is in "units." Gain of 1 unit is 5 yards.

psacklose0= .20;     %Probability that a sack loses 0 yards.
psacklose5= .39;     %Probability that a sack loses 5 yards.
psacklose10=.41;     %Probability that a sack loses 10 yards.

pstar=zeros(1,maxgain);
pstar(1)=.88;
pstar(5)=.35;
theta=1.6;
gamma=(pstar(1)/pstar(5)-1)/4^theta;
for gain=1:maxgain
   pstar(gain)=pstar(1)/(1+gamma*(gain-1)^theta);
end

alpha=zeros(1,maxgain);
alpha(1)=.16;
alpha(maxgain)=.42;
for gain=2:maxgain-1
   alpha(gain)=alpha(1)+((alpha(maxgain)-alpha(1))/(maxgain-1))*(gain-1);;
end
beta=.11;   % Probability of an interception divided by the probability of an incompletion.
totalcoverage=7.5;  % The constraint variable X.

psack=.0683;                       % Probability of a sack.
playtime=[7 7 7 7 7 7 7 7 7 7];    % Seconds used by a play, up until the whistle blows.

% Boundary conditions. If offense reaches goal line, they score.
S=zeros(21,121,5,4,4);
W=zeros(21,121,5,4,4);

for timeouts=1:4
   for down=1:5
      for togo=1:4
         for t=1:121
            S(1,t,down,togo,timeouts)=1;
            W(1,t,down,togo,timeouts)=1;
         end
      end
   end
end

% Backward induction

for t=2:121
   for yards=2:20
      for down=1:4
         for togo=1:4
            for timeouts=1:4
               
numgains=min(maxgain,yards-1);
               
sack=psacklose0*W(yards, max(t-playtime(1),1), down+1, togo, timeouts) ...
       +psacklose5*W(min(yards+1,21), max(t-playtime(1),1), down+1, min(togo+1,4), timeouts) ...
           +psacklose10*W(min(yards+2,21), max(t-playtime(1),1), down+1, min(togo+2,4), timeouts);
               
for gain=1:numgains
   I(gain)=S(yards,max(t-playtime(gain),1),down+1,togo,timeouts);
   Ibeta=I(gain)/(1+beta);
   if gain >= togo
      C(gain)=probob*S(yards-gain,max(t-playtime(gain),1),1,2,timeouts) ...
         +(1-probob)*W(yards-gain,max(t-playtime(gain),1),1,2,timeouts);      
   else
      C(gain)=probob*S(yards-gain,max(t-playtime(gain),1),down+1,togo-gain,timeouts) ...
         +(1-probob)*W(yards-gain,max(t-playtime(gain),1),down+1,togo-gain,timeouts);          
   end
   A1(gain)=pstar(gain)*(C(gain)-Ibeta)+Ibeta;
   A2(gain)=alpha(gain)*Ibeta;
   pbar(gain)=pstar(gain)/(1+alpha(gain)*(totalcoverage+0.1));
   A3(gain)=pbar(gain)*(C(gain)-Ibeta)+Ibeta;
end
v(1)=max( C(1:numgains) );
v(2)=max( A3(1:numgains) );

if v(1)==0
   S(yards,t,down,togo,timeouts)=fgprob(yards)*scorefactor;
else
   [val,fvalfromfzero,exitflag]=fzero('sumfunc',v,options,numgains,A1,A2,alpha,totalcoverage);
   if exitflag<0
      exitflag
   end
   S(yards,t,down,togo,timeouts)=max( psack*sack+(1-psack)*val , fgprob(yards)*scorefactor );
end

spike=S(yards,max(t-spiketime,1),down+1,togo,timeouts);

if timeouts>1
   calltime=S(yards,max(t-timeouttime,1),down,togo,timeouts-1);
else
   calltime=0;
end

letclockrun=S(yards,max(t-letclockruntime,1),down,togo,timeouts);

W(yards,t,down,togo,timeouts)=max([spike calltime letclockrun]);

            end
         end
      end
   end
   t
end