plus1sp

[1]Jason T. Maddox

# Bayesian estimation of in-game home team win probability for Division-I FBS
college football

Ryan Sides  Texas Woman's University, Mathematics and Computer Science,
Denton, TX, U.S., e-mail<EMAIL_ADDRESS>Jane L. Harvill  Baylor University,
Statistical Science, Waco, TX, U.S., e-mail<EMAIL_ADDRESS>
###### Abstract

Maddox et al., 2022a [9], Maddox et al., 2022b [10] establish Bayesian methods
for estimating home-team in-game win probability for college and NBA
basketball. This paper introduces a Bayesian approach for estimating in-game
home-team win probability for Division-I FBS college (American) football that
uses expected number of remaining possessions and expected score as two
predictors. Models for estimating these are presented and compared. These,
along with other predictors are introduced into two Bayesian approaches for
the final estimate of in-game home-team win probability. To illustrate
utility, methods are applied to the 2021 Big XII Conference Football
Championship game between Baylor and Oklahoma State.

##  1 Introduction

Sports analytics has become a well-established area of research. Work spans
more than 30 years and a large range of difficulty. Since the early 2000s,
research in statistical methods for sports analytics has risen dramatically.
The review articles of Kubatko et al., [6], Santos-Fernandez et al., [16], and
Terner and Franks, [17] provide a fairly comprehensive review for sports
analytics for a wide variety of sports, including football. One problem of
interest is predicting the probability that the home team wins during the
course of the game, or predicting ``in-game win probability.''

Speaking broadly, models for predicting the outcome of a sporting event can be
classified based on two objectives: (1) pregame prediction or (2) in-game, or
in-play, prediction. Pregame prediction involves determining the outcome of a
game before play begins. Once play begins, the process of predicting the
outcome ends. In contrast, in-game prediction attempts to use the progress
during a game to determine win probabilities that vary as a function of in-
game variables, for example, elapsed game time or score difference. The focus
of this paper is in-game prediction for college football.

One early paper on predicting in-game win probability for baseball is
attributed to Lindsey, [7], who used a conditional maximum likelihood
estimator to determine the probability the home team wins given the current
inning the game and the home team's lead. Lindsey, used historical information
to determine percentages the home team had won in a variety of scenarios, and
applied those to the current game to determine current team's estimated win
probability as the game progressed. Up until this work, baseball decisions
were based on what would maximize the scoring output for an inning. In
contrast, Lindsey,'s groundbreaking research focused on determining how each
decision would affect the probability that a team wins instead of only the
change in expected score. A more recent development in estimating in-game win
probability was approached by Benz, [1] who modeled college basketball in-game
win probability by using a series of logistic regressions at different times
throughout games on score differential and pregame win probability. He then
smooths the multiple logistic regression models into a single smooth function.
Maddox et al., 2022b [10] propose a Bayesian model for the National Basketball
Association (NBA) based on time and score differential as the two predictors.
Much of their methodology can be extended to college football, but different
predictors must be considered. In football, score and time are not as
informative for estimating win probability. A simple consideration of the two
sports makes clear a different model is required for each. For example, in
basketball, teams score often and quickly – often within seconds of each other
– and one possession can result in zero up to five points. However, it is
unusual for a team to score in seconds in a football game, and the possible
values for scores are not sequential integers. Additionally, in football, many
other variables contribute to the likelihood of a score, such as field
position or down-and-distance.

Only a few attempts have been proposed for in-game prediction of win
probability for (American) football. Lock and Nettleton, [8] use a random
forest of regression trees with up to ten predictor variables to combine pre-
play variables to estimate in-game win probability before any play of an
National Football League (NFL) game. Pro Football Reference, [12] create a
quasi ``black box'' model, where they go into some detail about creating their
win probability model using expected points along with pregame win probability
and the known standard deviation of the end of the game score differential.
However, many of the details used for their model are not provided. Later,
Ruscio and Brady, [14] compare the performance of the random forest model by
Lock and Nettleton, and the model put forth by Pro Football Reference, when
applied to the NFL. Their findings were that there was no discernable
difference between the two models. Ruscio and Brady, were able to obtain the
Pro Football Reference, model for their paper to reproduce the results. For
Australian rules football, Ryall, [15] used play-by-play data with pregame Elo
rankings to develop a model for Australian Rules football. In what follows, a
new approach that uses Bayesian methods is proposed. Explanatory variables for
predicting in-game win probability include expected number of possessions
remaining and expected score differential and with pregame power rankings.

The remainder of the paper is organized as follows. Section 2 details the
process of gathering and cleaning the data. Section 3 presents a method for
modeling the expected number of possession remaining in a college football
game. In Section 4, multiple models are proposed for best estimation of the
expected score and compared to each other. The resulting estimates of expected
remaining possessions and expected score are used as predictors for the
overall win probability model described in Section 5. The models for win
probability are evaluated and applied to a specific game in Section 6. Section
7 provides closing remarks.

##  2 Data Collection

The primary goal of the proposed models is effective practical prediction of
in-game home team win probability for a single college football game or a
collection of college football games during the regular season. The data
collected for investigating the proposed models performances were taken from
ESPN. Specifically, play-by-play data from ESPN was scraped using _R_ [13] and
the package rvest [19]. For football games, ESPN does not display all plays on
each game's webpage by default. Instead, they store the data in a back end
server that can be accessed through the ESPN Developer Center. 111The
description of the ESPN Developer Center may be found at
www.espn.com/apis/devcenter/overview.html. There is a general form of the URL
for college football back-end data providing access to the play-by-plays. For
each game, the game id changes in the URL. The game ids are scraped for all
games in each season, then input into this URL to scrape the play-by-plays for
each game. The data was collected starting with the beginning of the 2004
college football season through the 2021 season, excluding the 2020 season due
to the uniqueness of that season resulting from the COVID-19 pandemic.

There were some issues that occurred working with the raw data from ESPN's
back-end server. As when working with any raw data, there can be typographical
errors or mislabels that must be evaluated and corrected, if correction is
possible. Clear identifier were not always reliable; for example, errors were
found on which team had possession of the ball. The data was combed through
many times to ensure that specific identifiers for the game state were
correctly labeled. For some games, play-by-play data was not available in any
capacity. This was true especially for games during or close to the 2004
season. There was no clear pattern to which games had no play-by-play data.
However, enough games were collected on all season that the missing games
should not have a negative impact on the analysis.

##  3 Possessions Remaining Model

Within the game of college football, pace or tempo of play has been a key
discussion point for the last 20 years. Mike Leach is considered a modern day
football pioneer. While he was the head football coach at Texas Tech
University, Leach ushered in a new offensive play style which he carried with
him to Washington State University and then Mississippi State University. The
style of play is commonly called the ``air-raid'' or ``up-tempo'' offense. A
primary feature of the play style is that as little time as possible is used
between successive offensive plays. The philosophy is that minimizing time
between snaps prevents the opponent from successfully setting up their
defense.

The most common measure of a team's pace is the team's average number of plays
per game. However, Troch, [18] introduces time between consecutive offensive
plays as a viable alternative. He argues that average number of plays per game
does not take into account the tempo of the opponent or the number of run
versus pass plays a team calls. Because of incomplete passes, pass plays will
stop the clock more frequently than run plays. Therefore, the more pass plays
a team attempts, the more plays that team will run during a game, without it
necessarily being due to that team's tempo.

For the purposes of the win probability model, interest lies more in the
number of possessions that remain in the game than the time between
consecutive offensive plays. The more possessions that remain, the more snaps,
and the more potential points are left for teams to score in the game.
Troch,'s critique of plays per game not accounting for the opponent is valid.
To account for this, a new method is proposed that is similar to the method
put forth by Pomeroy, [11] for college basketball. He calculates pace based on
the number of possessions a team would expect to have in a game against a team
that plays with average tempo.

In what follows, the term ``pace'' is defined as the expected number of
possessions against an opponent that plays with average tempo. Pace is
calculated by a recursive algorithm as described with the following steps for
each season. FAT iteration $m$, for team $k=1,2,\ldots,n$, let $\xi_{k,m}$
represent pace. For a given season, $n$ is the number of Football Bowl
Subdivision (FBS) teams.

  1. 1.

For all teams $k=1,\ldots,n$, initialize $\xi_{k,0}\equiv 0$.

  2. 2.

Choose $\delta>0$ to be small.222The value $\delta=0.0001$ was used. While
$\max_{k}\\{|\xi_{k,m}-\xi_{k,m-1}|\\}>\delta$,

    1. (a)

Calculate the average pace for all teams,

$\displaystyle\mu_{m}=\frac{1}{n}\sum_{k=1}^{n}\xi_{k,m-1}.$

    2. (b)

For each team $k=1,\ldots,n$, let $\kappa_{k}$ represent the collection of
indices corresponding to opponent teams played during the regular season. For
each team $k$, sum the paces of the opponent teams,

$\displaystyle\psi_{k,m}=\sum_{j\in\kappa_{k}}\xi_{j,m-1},\quad k=1,\ldots,n.$

    3. (c)

Find the average difference between the total number of possessions played
$x_{k,m}$ and the expected possessions played $\psi_{k,m}$,

$\displaystyle\epsilon_{k,m}=\frac{x_{k,m}-\psi_{k,m}}{w_{k}},\quad
k=1,\ldots,n,$

where $w_{k}$ is the number of games the team $k$ played in the given season.

    4. (d)

Assign the updated pace for teams based on their number of possessions above
expected,

$\displaystyle\xi_{k,m}=\mu_{m}+\epsilon_{k,m},\quad k=1,\ldots,n.$

  3. 3.

When $\max_{k}\\{|\xi_{k,m}-\xi_{k,m-1}|\\}\leq\delta$, the convergence
criteria is satisfied and the pace $\xi_{k}$ for each team is taken to be the
final iterative value of $\xi_{k.m}$.

The fastest and slowest team paces for Power 5 FBS teams for the 2021 season
are shown in Table 1. The Power 5 teams with the highest pace are teams that
are known for having up-tempo, pass heavy offenses. For example, Pitt and
Tennessee had two of the highest scoring offenses in the nation due to their
fast pace. On the other end, teams known for their slow-it-down, run-based
offensive style have the lower pace. Kansas State and Boston College are known
for their grind-it-out style that does not translate into many possessions or
points.

Tab. 1: Fastest and slowest paces for Power 5 teams for the 2021 regular season. Rank | Team | Pace ($\xi_{k}$)
---|---|---
[]11 | Oklahoma State Cowboys | 30.11
16 | Duke Blue Devils | 29.69
21 | Pitt Panthers | 29.13
22 | Colorado Buffaloes | 29.04
23 | Tennessee Volunteers | 28.96
[]124 | Boston College Eagles | 23.96
125 | Kentucky Wildcats | 23.76
127 | Oregon Ducks | 23.19
128 | BYU Cougars | 22.90
129 | Kansas State Wildcats | 22.42

To compute the value of expected possessions remaining after a play, denoted
$\tau$, the two teams' paces are averaged and the average is weighted with a
decreasing linear function of the time left in the game. If $\xi_{1}$ and
$\xi_{2}$ represent the pace for the two teams in the game, and $t$ the number
of seconds elapsed during the game, then $\tau$ is given by

$\tau=\left(\frac{3600-t}{3600}\right)\left(\frac{\xi_{1}+\xi_{2}}{2}\right).$
(1)

##  4 Point Value Model

There have been many attempts to accurately predict the point value of a drive
in football. Many models are explained in generalities, but the details are
vague. For example, one of the more notable expected point models was posted
by a user named ``Mike'' and is featured on Sports-Reference.com [12]. A
linear model is used to predict the average point value for a possession given
down, distance, and yards to end zone. However the website does not disclose
the specifics of the model. The author of this paper contacted
SportsReference.com. They did not wish to provide any details regarding the
model.

Many of the models predict the point value of the current drive. However, they
fail to account for the dependency of the success of the next drive on the
outcome of the current drive. Consider, for example, the situation when a team
has a fourth down 99 yards away from the end zone. That team will almost
certainly punt, which will typically result in the opposing team receiving the
ball with good field position, making that opponent more likely to score. On
the other hand, if a team has the ball just a few yards from the end zone,
they are likely to score. They will then kick the ball off and the opponent
would likely start their following possession about 75 yards from the end
zone. These two contrasting scenarios illustrate the effect of a prior
possession on the likelihood of the subsequent possession resulting in a score
and thus on the in-game win probability. Therefore, instead of modeling the
point value of only the current possession, the proposed model is built on the
combined score for both the current possession and the following opponent
possession.

Four competing models are used to predict the point value of the current and
ensuing possessions using the predictor variables in Table 2. The first model
is a linear regression model with no interactions. The second is a linear
regression that includes all $2^{6}=64$ interactions between the nine
predictors, with the exception of any interactions involving offensive pace,
defensive pace, or number of possessions. The third and fourth models are a
random forest model and an extreme gradient boosting (XGBoost) model,
respectively. Literature on linear regression is prevalent, and the technique
is commonly applied. However, random forests and the XGBoost models are more
current, lesser known techniques. Random forest and XGBoost models have become
more important with advances in computation. Sections 4.1 and 4.2 provide
overviews of these techniques. For more detail on random forest see James et
al., [5], and for more detail on XGBoost, see Chen and Guestrin, [2] and Jain,
[4].

Tab. 2: Predictor variables in competing models for predicting the point values of the current and ensuing possessions. Time remaining ($3600-t$) | Offensive score | Defensive score
---|---|---
Down | Distance | Yards to end zone
Number of possessions played by time $t$ | Offensive pace | Defensive pace

###  4.1 Random Forest Model

The random forest model is built on many independent decision trees. A
decision tree is a machine learning algorithm that takes the data as a whole,
finds the ``best'' predictor/value of predictor pair to split the data into
two groups, or child nodes. What is ``best'' is determined by the predictive
power, that is, the one that would minimize root mean square error (RMSE) or
mean absolute error (MAE) of the model, using the average response in each
node as the prediction for the observations that fall in that node. Each of
these child nodes are then split again by the same process, with splitting
continuing to occur until some stopping criteria is reached. This stopping
criteria could be the maximum number of splits in a tree, the maximum depth of
the tree, or the maximum number of terminal nodes.

To build a random forest, a sample with replacement is taken from the data,
typically the same size as the data. This sample is used to build a decision
tree, then another sample, independent from the previous, is taken and another
tree is built. This is repeated many times to gather many independent trees.
To make a prediction, an observation is inserted into each tree, with the
prediction being the average response from its terminal node. These
predictions from each of the trees are then averaged to obtain the overall
prediction for the random forest.

Machine learning models, such as random forests, do not require assumptions
such as statistical models do. Instead, they have parameters that require
``tuning.'' For example, the maximum number of splits in a tree is a parameter
that could be tuned for each random forest model. The process would be to run
the model on several different values for the maximum number of splits, then
the number that minimizes the MAE of the test data would be chosen for this
model. The fewer number of splits in the tree, the less fit each tree becomes
to the training data, causing a greater risk of underfitting the data.
Likewise, the larger number of splits in each tree, the greater the risk of
overfitting the training data. Parameters for random forests can be split into
two categories, tree-based parameters and model-based parameters. Tree-based
parameters are applied to each tree individually and include, but are not
limited to maximum number of splits, minimum number of observations in a
terminal node, and maximum parameters considered at each split. Model-based
parameters are values pertaining to the random forest as a whole and include,
but are not limited to number of trees in the model and probabilities attached
to each observation of being included in each sample for building the tree.

###  4.2 XGBoost Model

An extreme gradient boosting model is similar to a random forest in that it is
a machine learning algorithm that is built on numerous decision trees. The
difference lies in that instead of the trees being independent from another,
each tree is built one at a time on the errors of the previous tree in
XGBoost. Each tree is also weighted down by a shrinkage factor,
$\eta,\leavevmode\nobreak\ 0<\eta<1$, to help prevent the tree from learning
too quickly and over-fitting the training data. In the training set, denote
${\mathbf{x}}$ as the collection of all predictor variables in the training
set and $x_{i}$ as the vector of predictor values of observation
$i,\leavevmode\nobreak\ y_{i}$ be the response variable value, and
$\hat{f}(x_{i})$ be predicted value from the model $\hat{f}({\mathbf{x}})$ at
$x_{i}$. The XGBoost algorithm can be set up as follows:

  1. 1.

Initialize $\hat{f}^{0}(x_{i})\equiv 0$ and $r_{i}^{0}=y_{i}$ for all $i$ in
the training set, where $r_{i}^{0}$ is the residuals for each observation
before any iterations in the algorithm.

  2. 2.

For $b=1,2,\ldots,B$, where $B$ is the total number of trees in the algorithm,
repeat the following:

    1. (a)

Fit a tree $\hat{f}^{b}(\mathbf{x})$ to $r^{b-1}$, the residuals from the most
recently updated model.

    2. (b)

Update $\hat{f}^{b}(\mathbf{x})$ by adding a shrunken version of the new tree,

$\displaystyle\hat{f}(\mathbf{x})\leftarrow\hat{f}(\mathbf{x})+\eta\hat{f}^{b}(\mathbf{x}).$

    3. (c)

For each $i$, update the residuals,

$\displaystyle r^{b}_{i}\leftarrow r^{b-1}_{i}-\eta\hat{f}^{b}(x_{i}).$

  3. 3.

Output the boosted model,

$\displaystyle\hat{f}(\mathbf{x})=\sum_{b=1}^{B}\eta\hat{f}^{b}(\mathbf{x}).$

Choice of $B$ is critical, since for larger values of $B$, that is, a model
with more trees, $\eta$ must be closer to zero to prevent over-fitting the
training data. If $\eta$ is large, say close to 1, then the first tree built
will start to fit the training data well, and trees later on in the algorithm
will be fitting residuals that are more likely due to noise than to an actual
signal.

Tuning of the XGBoost model can often be time consuming due to the number of
parameters to tune around and the amount of time it can take to run each
iteration of the model for large data sets. XGBoost has the same tree-based
parameters that are present in random forest models, most notably the depth of
the tree, because the depth of each tree is the largest possible interaction
depth the model can detect. The additional model-based parameters that need
tuning are $B$ and $\eta$ which are tuned simultaneously. There are many
methods for tuning parameters of the XGBoost model, one of which is performed
using the following steps.

  1. 1.

Choose a relatively high $\eta$, somewhere between 0.05 and 0.2.

  2. 2.

Optimize $B$ for this shrinkage value, keeping $B$ to a value where a machine
can run the model relatively quickly.

  3. 3.

Tune tree-based parameters using the values of $\eta$ and $B$ obtained in the
first two steps.

  4. 4.

Decrease $\eta$ and increase $B$ proportionally until the model's performance
improves minimally. Measures of the model's performance include statistics
like MAE or RMSE.

###  4.3 Test and Training Data for Models

To fit the point value model, and later the win-probability model, the data
was separated into a test data set and two training data sets. In particular,
the last five of the 17 seasons were designated as the test data for the win
probability model described in Section 5 and evaluated in Section 6. For the
remaining 12 seasons, half of the data were randomly selected to act as the
overall data for the point value model. This data will be referred to as the
``point value data.'' The other half were used to build the win probability
model.

###  4.4 Point Value Model Selection

To determine which of the four point value models is best, the point value
data was further separated into test and training data sets. Each of the four
point value models were built on the training set that was randomly pulled
from the point value data and then applied to the point value test data. The
models were compared using mean absolute error (MAE). The results are shown in
Table 3. For both linear regression models, the average difference between the
predicted point difference for the current and following possession and the
actual point difference was slightly more than three points per play. The
random forest model slightly outperforms the regression models, lowering the
MAE to slightly below three points. The XGBoost model performs the best, with
the lowest MAE by a margin at least 0.2949 points. Therefore, the XGBoost
model is adopted to predict the point differential for the current and
subsequent possessions. Conclusions about model performance are similar based
on root mean square error.

Tab. 3: Performance of point value models using test MAE. Model | MAE
---|---
[]Linear regression | 3.0805
Regression w/ inter | 3.0614
Random forest | 2.9751
XGBoost | 2.6802

##  5 Win Probability Model

For a specific game, consider the random process that is the expected lead for
the home team following the current and succeeding possession, which will now
be referred to as expected score, $\omega$. The expected score will shift as
time in the game moves forward, or equivalently as the expected possessions
remaining in the game, $\tau$, decreases, where $\tau$ is calculated as in
(1). For specific values of $\tau$ and $\omega$, let $p_{\tau,\omega}$
represent the in-game probability that the home team will win the game. When
considering multiple games $i=1,2,\ldots,M$, let $Y_{i}=1$ if the home team
wins game $i$ and 0 otherwise. Consider $p_{\tau,\omega}$ as a continuous
function of $\tau$ and $\omega$. Maddox et al., 2022a [9] introduce a Bayesian
estimator of in-game win probability, $p_{t,\ell}$, based on current time $t$
and score differential $\ell$, in college basketball. Their methods are
extended and adapted to the NBA in Maddox et al., 2022b [10]. In both papers,
Maddox et al., 2022a argue that the nature of a basketball game makes time and
score differential excellent predictors for in-game win probability. For
college football, instead of using time and score differential as predictors,
$\tau$ and $\omega$ are preferred for two reasons. First, and most obviously,
the nature of the games of basketball and football are inherently different.
As opposed to basketball, the predictors $\tau$ and $\omega$ contain more
information about current the state of the game, and account for which team
has the ball, how physically close (on the field) that team is to scoring,
etc. We now extend the methods of Maddox et al., 2022b [10] to estimate
$p_{\tau,\omega}$ for college football.

###  5.1 Naive Estimator of In-game Win Probability

For each combination of $\tau$ and $\omega$ rounded to the nearest whole
number, consider the $(\tau,\omega)$ ``cell.'' On a specific cell, the number
of wins by the home team, $n_{\tau,\omega}$, follows a
binomial$(N_{\tau,\omega},p_{\tau,\omega})$ distribution, where
$N_{\tau,\omega}$ is the total number of games that are observed in the
$(\tau,\omega)$ cell. Due to the total number of possible cells, one for each
combination of $\tau$ and $\omega$, some cells may have a small value of
$N_{\tau,\omega}$, resulting in a large standard error for any estimator for
$p_{\tau,\omega}$. Deshpande and Jensen, [3] and Maddox et al., 2022a [9],
Maddox et al., 2022b [10] suggest a binning approach to address the small
sample sizes. Windows centered around $(\tau,\omega)$ can be created in such a
way that the in-game win probability remains relatively constant across the
window. For creating the interval around expected score, since scores were
considered for each two possessions, a shift of two possessions will rarely
have a major affect on win probability. Moreover within college football for
similar score differentials, two or fewer possessions should not have a large
effect on the win probability, especially early in the game. Therefore a
reasonable interval for expected possessions remaining is defined as
$[\tau-2,\tau+2]$. To determine an interval for the expected score, note that
when an offense scores, the fewest points attainable is three points from a
field goal. Therefore, the interval on expected score is taken to be
$[\omega-3,\omega+3]$. The same notation will be adopted for any
$[\tau-2,\tau+2]\times[\omega-3,\omega+3]$ window; that is, $N_{\tau,\omega}$
is the number of games in the window in which the home team has an expected
lead by any value in $[\omega-2,\omega+2]$ points with any expected
possessions remaining in $[\tau-2,\tau+2]$, and given the specific value of
$N_{\tau,\omega},\leavevmode\nobreak\
n_{\tau,\omega}=\sum_{i=1}^{N_{\tau,\omega}}Y_{i}$, distributed as a
binomial($N_{\tau,\omega},p_{\tau,\omega})$ random variable. Based on the
binomial distribution, a simple estimator for in-game home team win
probability for for each $(\tau,\omega)$ window is the maximum likelihood
estimator

$\bar{p}_{\tau,\omega}=\frac{n_{\tau,\omega}}{N_{\tau,\omega}}.$ (2)

As a given game approaches the end of regulation, each individual point and
possession will have a larger impact on the in-game win probability. Therefore
the windows should be modified at the end of the game to reflect this.
Starting from fifteen expected possessions remaining, the proposed method
shortens window lengths and widths; that is, the length of interval around
each of expected score and expected number of possessions remaining will
gradually decrease until there is an expected two possesses remaining, when
the intervals' widths become zero.

###  5.2 Dynamic Bayesian Estimator

To elicit a prior distribution for in-game win probability, Maddox et al.,
2022b [10] suggest polling a sample of industry experts. The same can be done
for college football. The authors contacted a panel of 14 college football
experts, including coaches from a major Division I college football team and
respected media personnel. Each provided their estimate of the probability a
team wins for each combination of score differential, $\ell$, and time
elapsed, $t$, in Table 4, regardless of which team is the home team. Note that
$t$ and $\ell$ are the actual time and score differential in the hypothetical
game and are not the same as the previously defined $\tau$ and $\omega$. To
reduce the level of complexity and increase the likelihood of response, the
authors choose to ask the question of win probability in terms of $t$ and
$\ell$ as opposed to $\tau$ and $\omega$.

The Bayesian prior parameters have an interesting interpretation, first noted
by Deshpande and Jensen, [3]. The parameter $\alpha_{t,\ell}$ can be
interpreted as the number of ``pseudo-wins'' in the $(t,\ell)$ cell; likewise
$\beta_{t,\ell}$ as the number of ``pseudo-losses.'' Through this
interpretation, the two parameters can be seen as a way of increasing the
number of games in a specific $(t,\ell)$ cell. If the home team is ahead, then
first scale parameter being large effectively acts to increase the number of
wins in that cell. On the other hand, if the home team is behind, the second
scale parameter is large, and acts to increase the number of losses.

The sample mean $\tilde{p}_{t,\ell}$ and sample variance $s^{2}_{t,\ell}$, of
the probabilities were computed. The two scale parameters were estimated via a
method-of-moments type approach. The system of equations

$\displaystyle\tilde{p}_{t,\ell}$
$\displaystyle=\frac{\alpha_{t,\ell}}{\alpha_{t,\ell}+\beta_{t,\ell}},$
$\displaystyle s^{2}_{t,\ell}$
$\displaystyle=\frac{\alpha_{t,\ell}\beta_{t,\ell}}{\left(\alpha_{t,\ell}+\beta_{t,\ell}\right)^{2}\left(\alpha_{t,\ell}+\beta_{t,\ell}+1\right)},$

is solved for $\alpha_{t,\ell}$ and $\beta_{t,\ell}$, yielding

$\displaystyle\alpha_{t,\ell}$
$\displaystyle=-\frac{\tilde{p}_{t,\ell}\left(\tilde{p}_{t,\ell}^{2}-\tilde{p}_{t,\ell}+s_{t,\ell}^{2}\right)}{s^{2}_{t,\ell}}$
$\displaystyle\beta_{t,\ell}$
$\displaystyle=\frac{\left(\tilde{p}_{t,\ell}-1\right)\left(\tilde{p}_{t,\ell}^{2}-\tilde{p}_{t,\ell}+s^{2}_{t,\ell}\right)}{s_{t,\ell}^{2}},$

as long as
$\left(\tilde{p}_{t,\ell}-1\right)\tilde{p}_{t,\ell}\left(\tilde{p}_{t,\ell}^{2}-\tilde{p}_{t,\ell}+s^{2}_{t,\ell}\right)\neq
0$.

In Table 4 the score differential is presented when the home team has the
lead. The larger the lead, the more the prior is left-skewed. On the other
hand, if the visiting team has the lead, then the roles of $\alpha_{t,\ell}$
and $\beta_{t,\ell}$ are reversed, causing the prior to be right-skewed. At
any time, if the game is tied ($\ell=0$), or is sufficiently close in score
for the amount of time remaining, the prior distribution is a diffuse
beta$(1,1)$ prior, allowing the likelihood to drive the results in the
posterior.

Tab. 4: Imputed parameters for beta prior. Elapsed Time (t) | Home Team |  |
---|---|---|---
(sec.) | Lead ($\ell$) | $\alpha_{t,\ell}$ | $\beta_{t,\ell}$
[]$[0,900]$ | $[0,7]$ | 1 | 1
$[0,900]$ | $(7,14]$ | 21 | 10
$[0,900]$ | $(14,28]$ | 16 | 3
$[0,900]$ | $(28,\infty)$ | 59 | 1
$(900,1800]$ | $[0,7]$ | 1 | 1
$(900,1800]$ | $(7,14]$ | 15 | 7
$(900,1800]$ | $(14,28]$ | 11 | 2
$(900,1800]$ | $(28,\infty)$ | 44 | 1
$(1800,2700]$ | $[0,7]$ | 1 | 1
$(1800,2700]$ | $(7,14]$ | 14 | 4
$(1800,2700]$ | $(14,28]$ | 13 | 1
$(1800,2700]$ | $(28,\infty)$ | 126 | 1
$(2700,3300]$ | $[0,3]$ | 1 | 1
$(2700,3300]$ | $(3,7]$ | 28 | 16
$(2700,3300]$ | $(7,14]$ | 27 | 8
$(2700,3300]$ | $(14,21]$ | 11 | 1
$(2700,3300]$ | $(21,\infty)$ | 98 | 1
$(3300,3480]$ | 0 | 1 | 1
$(3300,3480]$ | $(0,3]$ | 29 | 18
$(3300,3480]$ | $(3,7]$ | 17 | 6
$(3300,3480]$ | $(7,10]$ | 13 | 2
$(3300,3480]$ | $(10,14]$ | 16 | 1
$(3300,3480]$ | $(14,\infty)$ | 98 | 1
$(3480,3600)$ | 0 | 1 | 1
$(3480,3600)$ | $(0,3]$ | 21 | 10
$(3480,3600)$ | $(3,7]$ | 18 | 5
$(3480,3600)$ | $(7,10]$ | 16 | 1
$(3480,3600)$ | $(10,\infty)$ | 91 | 1

The maximum likelihood estimator in (2) is based on $\tau$ and $\omega$.
However, the elicited prior is based on $t$ and $\ell$. Consequently, the in-
game win probability will be written $p_{t,\ell,\tau,\omega}$. Since the beta
family of distributions is a conjugate prior for the binomial distribution,
the beta-binomial connection is used to estimate $p_{t,\ell,\tau,\omega}$ with
the mean of the posterior beta distribution, specifically

$\hat{p}_{t,\ell,\tau,\omega}=\frac{n_{\tau,\omega}+\alpha_{t,\ell}}{N_{\tau,\omega}+\alpha_{t,\ell}+\beta_{t,\ell}}.$
(3)

###  5.3 Adjusted Dynamic Bayesian Estimator

During a game, the in-game win probability is clearly a function of expected
possessions remaining and expected score. However, it is also affected by the
overall skill of the teams playing the game. The skill level can be
incorporated using pregame win probabilities for each game. The normal
distribution quantile function was used to convert the pregame point spread to
a pregame home team win probability $\hat{p}_{p}$ for each game using the
method outlined by Maddox et al., 2022a [9].

Maddox et al., 2022b [10] introduce three different functions to incorporate
pregame probabilities. The first is a linear function of only time remaining.
The second is a linear function of time remaining and score differential. The
third is linear in time, but quadratic in score differential. Other variations
to the three weight functions were considered. For example, including a
quadratic term for time was attempted. However, the more complicated models
would not converge appropriately for _R_ to optimize the performance
accurately. The three weight functions specifically considered here are

$\displaystyle D_{1}$ $\displaystyle=bt,$ $\displaystyle D_{2}$
$\displaystyle=c_{0}+c_{1}t+c_{2}|\ell|,$ $\displaystyle D_{3}$
$\displaystyle=d_{0}+d_{1}t+d_{2}|\ell|+d_{3}\ell^{2}.$

Each of the weight functions yields a competing model for including pregame
win probabilities, which can be represented as

$p^{*}_{t,\ell,\tau,\omega,j}=\begin{cases}\hat{p}_{p},&D_{j}\leq 0\\\
\left(1-D_{j}\right)\hat{p}_{p}+D_{j}\hat{p}_{t,\ell,\tau,\omega},&0<D_{j}<1\\\
\hat{p}_{t,\ell,\tau,\omega},&D_{j}\geq 1\end{cases},\qquad j=1,2,3,$

where $p^{*}_{t,\ell,\tau,\omega,j}$ is the final predicted win probability
associated with weight function $D_{j}$.

Values for $b,c_{0},\ldots,d_{3}$ are estimated by minimizing Brier score for
each $p^{*}_{t,\ell,\tau,\omega,j}$ computed from the holdout test data. The
Brier scores for each $p^{*}_{t,\ell,\tau,\omega,j}$ are shown in Table 5. The
model $p^{*}_{t,\ell,\tau,\omega,2}$ is the most accurate predictor for the
test data. The fitted model for $D_{2}$ that provided the minimal Brier score
is

$\displaystyle D_{2}=0.09589+0.00018t+0.02523\ell.$

Using this expression for $B_{2}$, the model

$\displaystyle
p^{*}_{t,\ell,\tau,\omega,2}=\begin{cases}\hat{p}_{p},&D_{2}\leq 0\\\
\left(1-D_{2}\right)\hat{p}_{p}+D_{2}\hat{p}_{t,\ell,\tau,\omega},&0<D_{2}<1\\\
\hat{p}_{t,\ell,\tau,\omega},&D_{2}\geq 1\end{cases}$

is the ``adjusted dynamic Bayesian estimator.'' The percent of in-game win
probability accounted for in the final predicted win probability is displayed
in Figure 1. As desired, early on in any game the majority of the final
predicted win probability comes from the pregame win probability. However, as
the game draws to an end or one team takes a large lead, the final predicted
win probability gets closer and closer to the dynamic Bayesian estimator.

Tab. 5: Brier scores for models determining in-game probability proportion. Proportion Model | Brier Score
---|---
[]Linear Time ($D_{1}$) | 0.1272
Linear Time & Score ($D_{2}$) | 0.1250
Quadratic ($D_{3}$) | 0.1265
Fig. 1: Graphical function of $D_{2}$.

##  6 Model Evaluation and Application

For each play from each game from the 2017 through 2021 seasons, excluding the
2020 season, Brier score was computed. The Brier scores of the three models,
dynamic Bayesian estimator, adjusted dynamic Bayesian estimator, and the
random forest model proposed by Lock and Nettleton, [8] are shown in Table 6.
Brier's score is a statistic used to compare the performance of different
methods for estimating probabilities. Brier's score is the average of the
square of the difference between the estimated probability and the observed
binary outcome. In the context of in-game home team win probabilities, this
observed binary is $Y_{i}$ as defined in Section 5. To interpret Brier's
score, if $Y_{i}=1$ for all $i$, and the predicted probability is also one for
every $i$, then Brier's score will be zero, indicating perfect prediction. On
the other hand, if for all $i,Y_{i}=0$, and the estimated probability is one,
then Brier's score will be one, the worst possible Brier's score. Both
Bayesian models outperform the random forest model. The dynamic Bayesian model
with the pregame adjustment is the model that performs best overall.

Tab. 6: Brier scores for predictive performances for 2017 through 2021 seasons. Model | Brier Score
---|---
[]Dynamic Bayes | 0.1453
Adjusted dynamic Bayes | 0.1250
Random forest | 0.1705

To observe the performance of each model further, each model can be applied to
a specific game, observing the features of the models as the game progresses.
The models are run for the 2021 Big 12 Championship game between Baylor
University and Oklahoma State University. The results can be seen in Figure 2.

Fig. 2: In game win probability for 2021 Big 12 Championship.

On December 4, 2021, the Baylor Bears and Oklahoma State Cowboys met at AT&T
Stadium in Dallas, Texas for the Big 12 Championship game. Oklahoma State was
the home team. During the regular season, Oklahoma State had beaten Baylor,
and were considered the favorite to win this game. The probability traces for
the game for each of the three models is seen in Figure 2. The adjusted
dynamic Bayes model (red curve) shows that Oklahoma State is favored by
starting at well over 50% chance that Oklahoma State would win. Early in the
game, one of the features that leads to the Bayesian models outperforming the
random forest (blue line) is can be seen. The random forest model is too quick
to jump to a win probability close to 0 or 1, especially early in the game.
The game got off to a relatively slow start, but by halftime Baylor had
surprised many by jumping out to a 21-6 lead. At that time, all three models
predicted Baylor was most likely to win. In the second half, Oklahoma State
started to make a come back. Halfway through the third quarter, Oklahoma State
was able to score a touchdown and make the score 21-13, significantly raising
their probability of winning. The Cowboys then had another significant uptick
in their win probability with ten minutes to go when they earned a first-and-
goal from the Baylor one-yard line, appearing to be on the verge of scoring a
touchdown with the potential to tie the game. However, Baylor was able to hold
Oklahoma State to just a field goal, maintaining their lead and their edge in
win probability. The climax of the game came with three minutes left when
Oklahoma State started a drive from their own 10 yard line. They methodically
drove the ball down the field and once again ended up with a first-and-goal,
this time from the two-yard line with 80 seconds left. Being down by 5 points,
they were not going to kick a field goal. Instead, they had four attempts to
score a touchdown to win the game. With so many attempts from so close, it
appeared likely that Oklahoma State would be able to score, driving their win
probability in the two Bayesian models over 50% despite trailing. Once again
the Baylor defense made an impressive goal line stand, keeping Oklahoma State
out of the end zone, ultimately less than half of a yard short, sealing the
victory for Baylor.

##  7 Conclusion

Two new methods are proposed for estimating in-game win probability for
college football games. Both are an extension and enhancement of the methods
in Maddox et al., 2022a [9], Maddox et al., 2022b [10], which provide a number
of models for estimating in-game home-team win probabilities for NCAA
basketball and NBA basketball respectively. The first proposed ``dynamic
Bayesian estimator,'' uses expected score differential and expected
possessions remaining as predictors, which are modeled from the data
themselves, and a prior that has been calculated based on the distribution of
predicted win probabilities from 14 college football field experts. The second
method, referred to as the ``adjusted dynamic Bayesian estimator,'' adjusts
the dynamic Bayesian estimator based on pregame win probabilities obtained
from TeamRankings.com. The adjustment is optimized over a function of both
time and score so that as the game moves on or the score differential
increases, the adjusted dynamic Bayesian estimator will begin to approach the
dynamic Bayesian estimator rather than the pregame win probability. These two
methods are then compared to the random forest win probability model proposed
by Lock and Nettleton, [8]. Both new models outperform the standard random
forest model, with the adjusted dynamic Bayesian estimator performing the best
out of all of the models.

## Bibliography

  * Benz, [2019] Benz, L. (2019).  A new ncaahoopr win probability model.  [Online; posted December 26, 2019].
  * Chen and Guestrin, [2016] Chen, T. and Guestrin, C. (2016).  XGBoost: A scalable tree boosting system.
  * Deshpande and Jensen, [2016] Deshpande, S. K. and Jensen, S. T. (2016).  Estimating an NBA player’s impact on his team’s chances of winning.  Journal of Quantitative Analysis in Sports, 12(2):51–72.
  * Jain, [2022] Jain, A. (2022).  Complete guide to parameter tuning in XGBoost.
  * James et al., [2013] James, G., Witten, D., Hastie, T., and Tibshirani, R. (2013).  An Introduction to Statistical Learning: with Applications in R.  Springer.
  * Kubatko et al., [2007] Kubatko, J., Oliver, D., Pelton, K., and Rosenbaum, D. T. (2007).  A starting point for analyzing basketball statistics.  Journal of Quantitative Analysis in Sports, 3(3).
  * Lindsey, [1963] Lindsey, G. R. (1963).  An investigation of strategies in baseball.  Operations Research, 11(4):477–501.
  * Lock and Nettleton, [2014] Lock, D. and Nettleton, D. (2014).  Using random forests to estimate win probability before each play of an NFL game.  Journal of Quantitative Analysis in Sports, 10(2):197–205.
  * [9] Maddox, J. T., Sides, R., and Harvill, J. L. (2022a).  Bayesian estimation of in-game home team win probability for college basketball.  Journal of Quantitative Analysis in Sports.
  * [10] Maddox, J. T., Sides, R., and Harvill, J. L. (2022b).  Bayesian estimation of in-game home team win probability for national basketball association games.  Journal of Sport Analytics.
  * Pomeroy, [2012] Pomeroy, K. (2012).  Ratings glossary: The kenpom.com blog.
  * Pro Football Reference, [2012] Pro Football Reference (2012).  The p-f-r win probability model.  https://www.sports-reference.com/blog/2012/03/features-expected-points/.
  * R Core Team, [2021] R Core Team (2021).  R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing, Vienna, Austria.
  * Ruscio and Brady, [2021] Ruscio, J. and Brady, K. (2021).  Estimating win probability for nfl games.
  * Ryall, [2011] Ryall, R. (2011).  Predicting Outcomes in Australian Rules Football.  PhD thesis, Royal Melbourne Institute of Technology University.
  * Santos-Fernandez et al., [2019] Santos-Fernandez, E., Wu, P., and Mengersen, K. L. (2019).  Bayesian statistics meets sports: a comprehensive review.  Journal of Quantitative Analysis in Sports, 15(4):289–312.
  * Terner and Franks, [2021] Terner, Z. and Franks, A. (2021).  Modeling player and team performance in basketball.  Annual Review of Statistics and Its Applications, 8(1):1–23.
  * Troch, [2016] Troch, J. (2016).  A new look at college football tempo.
  * Wickham, [2021] Wickham, H. (2021).  rvest: Easily Harvest (Scrape) Web Pages.  R package version 1.0.2.


# Measuring periods in aperiodic light curves – Applying the GPS method to
infer rotation periods of solar-like stars

Timo Reinhold Max-Planck-Institut für Sonnensystemforschung, Justus-von-
Liebig-Weg 3, 37077 Göttingen, Germany  Alexander I. Shapiro Max-Planck-
Institut für Sonnensystemforschung, Justus-von-Liebig-Weg 3, 37077 Göttingen,
Germany  Sami K. Solanki Max-Planck-Institut für Sonnensystemforschung,
Justus-von-Liebig-Weg 3, 37077 Göttingen, Germany  Gibor Basri Department of
Astronomy, University of California, Berkeley, CA 94720, USA

###### Abstract

Light curves of solar-like stars are known to show highly irregular
variability. As a consequence, standard frequency analysis methods often fail
to detect the correct rotation period. Recently, Shapiro et al. (2020) showed
that the periods of such stars could still be measured by considering the
Gradient of the Power Spectrum (GPS) instead of the power spectrum itself. In
this study, the GPS method is applied to model light curves of solar-like
stars covering all possible inclination angles and a large range of
metallicities and observational noise levels. The model parameters are chosen
such that they resemble those of many stars in the Kepler field. We show that
the GPS method is able to detect the correct rotation period in $\approx 40\%$
of all considered cases, which is more than ten times higher than the
detection rate of standard techniques. Thus, we conclude that the GPS method
is ideally suited to measure periods of those Kepler stars lacking such a
measurement so far. We also show that the GPS method is significantly superior
to auto-correlation methods when starspot lifetimes are shorter than a few
rotation periods. GPS begins to yield rotation periods that are too short when
dominant spot lifetimes are shorter than one rotation period. We conclude that
new methods are generally needed to reliably detect rotation periods from
sufficiently aperiodic time series — these periods will otherwise remain
undetected.

stars: rotation — stars: solar-type — Kepler

††journal: ApJL

##  1 Introduction

Stellar rotation periods are usually measured by searching for periodic
patterns in long-term (photometric) time series caused by magnetic features
transiting the stellar disk several times. The high-precision, long-term
photometry of the Kepler space telescope has proven to be an ideal data set
for large rotation period surveys (e.g., Reinhold et al. 2013; McQuillan et
al. 2014; Santos et al. 2019, etc.). In particular, these surveys have shown
that rotation periods of stars with regular light curve patterns can reliably
be detected by standard frequency analysis methods such as the Lomb-Scargle
periodogram (Reinhold et al., 2013; Reinhold & Gizon, 2015), the auto-
correlation function (McQuillan et al., 2014), or wavelet power spectra
(Santos et al., 2019, 2021).

These methods, however, reach their limits when applied to stars around solar
spectral type that show comparably small and irregular variability. The Sun
also shows an irregular light curve pattern, which is mostly caused by the
relatively short lifetimes of the sunspots (days to weeks, see Petrovay & van
Driel-Gesztelyi 1997) compared to the solar rotation period of roughly 27
days. A similar behavior is seen for stars around solar spectral type, which
often show highly irregular light curves. As a consequence, rotation periods
could only be detected for 16% of the G-type dwarfs in the Kepler field (see
Table 1 and 2 in McQuillan et al. 2014). Basri et al. (2022) showed that the
stars they were unable to find reliable periods for are the ones whose
starspot lifetimes are shorter than a couple of rotation periods.

Recently, Reinhold et al. (2021) tested the performance of the auto-
correlation function (ACF) on simulated light curves with solar variability
patterns. These authors showed that the ACF is not well-suited for the period
detection in such aperiodic stars because the method (by design) searches for
a repeating variability pattern. Reinhold et al. (2021) further showed that
the period detection rate is as low as $\approx 3\%$ when setting similar
constraints on the auto-correlation peak heights as used in McQuillan et al.
(2014).

Thus, measuring periodic signals in a non- or hardly periodic time series is
challenging and requires novel techniques. Recently, Shapiro et al. (2020)
showed that the periods of such stars could still be measured by considering
the Gradient of the (global wavelet) Power Spectrum (GPS) instead of the power
spectrum itself. This so-called GPS method is able to detect the correct
rotation period in an aperiodic time series. In particular, it was shown that
the GPS method is able to determine the correct rotation period of the Sun
(Amazo-Gómez et al., 2020b). Furthermore, the GPS method has been successfully
applied to Kepler stars with previously determined rotation periods (Amazo-
Gómez et al., 2020a).

In contrast to classical period analysis methods (such as the Lomb-Scargle
periodogram, auto-correlation functions, or power spectra), the GPS method
does not require a repeatable pattern in the time series. Instead, the GPS
method is mostly sensitive to the shape and in particular the width of the
light curve profile of individual starspot transits averaged over time (see
Sect. 3.4). Consequently, the GPS method even works in cases when the magnetic
features live less long than the stellar rotation period such that no
reoccurring transits of the same magnetic features are needed. These
characteristics makes it a powerful tool for measuring periods in aperiodic
time series.

Here, we test the GPS method on simulated light curves with irregular
variability, and show that it reaches a detection rate of $\approx 40\%$,
which is more than ten times higher than the detection rates obtained by
classical methods (such as auto-correlation functions). In a forthcoming
publication, we will apply this method to all Kepler stars with near-solar
effective temperatures with the goal of measuring periods for the large sample
of stars lacking this information so far.

##  2 Data & Methods

###  2.1 The model light curves

In this study, we primarily apply the GPS method to simulated light curves
with solar-like variability. The light curves rely on models of solar
brightness variability, SATIRE (Spectral And Total Irradiance REconstruction,
Fligge et al. 2000; Krivova et al. 2003; Solanki et al. 2013). These models
build on the solar paradigm to simulate light curves with solar effective
temperature, rotation period, and activity level, and have recently been
extended to modeling stars observed at different inclination angles (Nèmec et
al., 2020b, a) and metallicities (Witzke et al., 2020; Reinhold et al., 2021).

The models considered here span 61 years of solar data between March 1948 and
May 2009, i.e., the data almost cover 6 full solar activity cycles with epochs
of high and low activity. The models were calculated for ten inclination
angles between pole-on ($i=0^{\circ}$) and edge-on ($i=90^{\circ}$) views, and
nine different metallicities $\rm-0.4\leq[Fe/H]\leq 0.4$ dex (in steps of 0.1
dex). In particular, the rotation period of all model light curves is equal to
the solar one, which roughly equals 27 days as seen from Earth.

To be comparable to Kepler observations later on, we chose a Monte Carlo
approach to account for a wide range of possible stellar parameters such as
inclination angles, metallicities, and noise levels. Following the strategy
described in Reinhold et al. (2021), the input model parameters inclination
and metallicity are chosen such that they resemble the distribution in the
Kepler field. From these distributions, we chose a random combination of
inclination and metallicity, and pick the model closest to the chosen values.
From this model, we randomly pick a 4-yr segment from the whole 61-years time
series, and remove the long-term variability of the activity cycle (see
Reinhold et al. 2021 for details). The GPS method is then applied to these
Keplerized light curves as described in the following.

Since the light curves are noise-free, we add white noise with zero mean and
standard deviation $\sigma$ to the data. The noise level $\sigma$ is
determined by the apparent magnitude of the observed star. We adopt the
distribution of apparent magnitudes (Kp) of solar-like stars in the Kepler
field, and compute different noise realizations following the procedure
described in Reinhold et al. (2020). In total, we have considered 5000 samples
each for the noise-free and the noisy light curves. For more details, we refer
the reader to Reinhold et al. (2021), who used the exact same light curves as
we do.

In the following, the GPS method will primarily be applied to the physics-
based SATIRE model light curves. In order to directly test the effect of
starspot lifetimes on the efficacy of the GPS and auto-correlation methods, we
also utilized some of the models described in Basri & Shah 2020 (see Sect.
3.4). These do not necessarily mimic stars similar to the Sun, but instead
were computed to include the general behavior of light curves of Kepler stars
whose periods have been determined. The models do not simulate physical flux
emergence processes, instead sampling light curves produced by many randomized
manifestations of starspots whose general parameters are specified and varied.
In particular, a set of models was employed with an average of six spots (all
with the same maximum size) present at a time, distributed randomly over the
stellar sphere, and viewed at an inclination of 60 degrees. Starspot lifetimes
for these noise-free model sets were set at values of 0.25, 0.5, 1, 2, 3, and
4 rotation periods. Each model set contains 1000 different instances of random
spot placements and birth dates and each run has 50 rotations. To better
illustrate the effect of different spot lifetimes, we refrain from adding
noise to these models.

###  2.2 The GPS method

Fig. 1 shows an example from our set of (noisy) model light curves (top
panel). Individual spot transits are clearly visible but without any obvious
periodicity. This observation is confirmed by the global wavelet111Here, we
use a 6th order Paul wavelet (e.g., Torrence & Compo, 1998). power spectrum of
the light curve (middle panel). The power spectrum increases toward longer
periods but does not show a distinct peak. Instead, it shows a plateau shape
at the model rotation period of 27 days. Thus, the rotation period would not
have been detected by this method.

Instead of searching for peaks of the power spectrum, we compute the Gradient
of the Power Spectrum (GPS), which is shown in the lower panel of Fig. 1. The
maximum of the GPS corresponds to the position of the inflection point (IP),
i.e., the point where the curvature of the high-frequency tail of the power
spectrum changes its sign. In this example, the highest peak is found at the
inflection point period $P_{\rm IP}\approx 5.41$ days. Shapiro et al. (2020)
showed that this period can be used to infer the correct rotation period by
the simple equation

$P_{\rm{rot}}=P_{\rm IP}/\alpha,$ (1)

where $\alpha$ is a calibration factor. We will see below how $\alpha$ depends
on the model parameters, and what the best choice of $\alpha$ might be for
real observations where physical parameters cannot be controlled. For the
model in Fig. 1, we find that $\alpha=5.41/27=0.20$ is the best choice. The
main idea behind the GPS method is that the high-frequency tail of the power
spectrum is much less affected by the evolution of magnetic features than the
power spectrum peak associated with the rotation period (see Fig. 3 in Shapiro
et al. 2020). We note that the inflection point itself does not have a
physical meaning, as stated in Shapiro et al. (2020). However, we ascribe a
sort of physical meaning to the inflection point, which can be associated with
the typical dip duration in the light curves caused by active regions crossing
the visible disc (see Sect. 3.4).

For some models, the GPS shows more than one peak, i.e., more than one
inflection point. From visual inspection, we found that this mostly occurs for
cases where the light curve shows some quiet and some active segment. Such
models usually return two peaks in the gradient, with one period much smaller
and the other peak much larger than the ”correct” inflection period. Thus, we
discard all cases where more than one GPS peak was found. This criterion
removes $\approx 25\%$ of the noise-free cases but only $\approx 10\%$ of the
noisy models. The reason why the noisy models less often show two peaks is
that the noise partly washes out the impact of small spot transits that have
very short lifetimes (often less than one day). Consequently, more power of
the GPS peak can be associated with the profiles of the remaining larger
spots.

Figure 1: Upper panel: Keplerized model light curve (including noise) with an
inclination of $i=70^{\circ}$ and solar metallicity $[Fe/H]=0.0$ dex.
Individual spot transits are visible but there is no obvious periodicity.
Middle panel: global wavelet power spectrum. The vertical blue line indicates
the model rotation period of 27 days. The power spectrum does not show any
peak at this period. Lower panel: Gradient of the Power Spectrum (GPS). The
highest peak indicates the period at the inflection point $P_{\rm{IP}}=5.41$
days. This period can be used to infer the rotation period.

##  3 Results

###  3.1 The calibration factor $\alpha$

To eventually determine the rotation period, we need to know the calibration
factor $\alpha$. According to eq. 1, this factor equals the measured
inflection period divided by the model rotation period of 27 days. Following
our Monte Carlo approach, we draw random samples of model light curves, and
measure the inflection point periods. Fig. 2 shows the distribution of
$\alpha$ for the noise-free (blue) and the noisy (red) models. The noise-free
distribution has a roughly Gaussian shape with a mean of
$\langle\alpha\rangle=0.165$ and standard deviation $\sigma_{\alpha}=0.032$.
The distribution for the noisy models is only approximately Gaussian with a
mean of $\langle\alpha\rangle=0.217$ and standard deviation
$\sigma_{\alpha}=0.047$. Generally, one can see that the noise shifts the
distribution to higher $\alpha$ values. We will see below (see Fig. 4) that
the noise-dominated models are mostly responsible for the long tail toward
larger $\alpha$ values (i.e. longer inflection periods). We like to emphasize
that the specific value of $\alpha$ is independent of the model rotation
period. In the following, we show how $\alpha$ depends on different model
parameters, and what value of $\alpha$ should be used for real observations
(see Sect. 3.3).

Figure 2: Distribution of the measured inflection point period divided by the
model rotation period of 27 days (i.e. $\alpha$) for the noise-free (blue) and
the noisy (red) models for 5000 random samples. Adding noise to the light
curves shifts the inflection point toward longer periods (i.e. larger values
of $\alpha$).

###  3.2 Dependence on model parameters

In Fig. 3, we can study the dependence of $\alpha$ on the model parameters in
more detail. The upper and middle panels show the noise-free models to study
the effects of inclination and metallicity independently of the noise, whereas
the lower panel shows the $\alpha$ dependence on apparent magnitude (i.e., the
effect of noise in the light curves). In all panels, the dashed blue line
shows $\langle\alpha\rangle$ and the red dots with error bars indicate the
mean and standard deviation in the chosen bins.

The upper panel shows the dependence of $\alpha$ on inclination for all
different metallicities drawn from the input distribution. One can see that
$\alpha$ does not show much dependence on inclination down to $i=40^{\circ}$.
Below that angle, the contribution of faculae to the overall variability
increases while the spot contribution shrinks (e.g., Knaack et al. 2001). The
faculae are well visible near the limb and leave a more sinusoidal variability
pattern in the light curves than the spots, leading to greater $\alpha$ values
eventually. The same is true for high latitude spots, which are seen much
longer than spots at lower latitudes. This interplay eventually leads to a
bimodal distribution of $\alpha$ for inclinations below $i=30^{\circ}$.

The middle panel shows that $\alpha$ does not depend on metallicity. As shown
in Witzke et al. (2020); Reinhold et al. (2021), a change in metallicity
mostly affects the facular contrasts, whereas the spots contrasts are almost
unchanged. This result is consistent with the observation that the GPS method
is most sensitive to the spot transit profiles.

The strongest dependence of $\alpha$ is found for observational noise such
that $\alpha$ strongly increases toward fainter stars. Since the majority
(more than 70%) of the (solar-like) stars in the Kepler field (and thus in our
models) are fainter than 14th magnitude, the noise is the dominant uncertainty
for $\alpha$.

Figure 3: Dependence of $\alpha$ on inclination (first panel), metallicity
(second panel), and Kepler magnitude Kp (third panel) for 5000 samples. The
dashed blue line indicates the median $\alpha$ of all models, and the red dots
and error bars indicate the mean and standard deviation of $\alpha$ in the
considered bins. We note that the upper and middle panels show the noise-free
models, whereas the lower panel shows the results for the noisy model light
curves.

Additionally to the model parameters, $\alpha$ shows some dependence on the
phase of the activity cycle. Fig. 4 shows the 61 years of the undetrended
model time series for the solar case, i.e., $i=90^{\circ}$ and $[Fe/H]=0$
(gray). Since we draw 4-year time series from the whole observing period, the
inflection periods are computed for different epochs of the solar cycle.
During activity maximum, the (rotational) variability is dominated by spots,
whereas faculae are the predominant type of activity during its minima. Thus,
the measured $\alpha$ values should depend on the phase of the cycle. This is
shown by the blue (noise-free) and red (noisy) lines, which depict the yearly
averages of $\alpha$ for all model realizations. To show both the light curve
and the $\alpha$ values on top of each other, the mean of the $\alpha$
distribution has been subtracted.

The red curve clearly peaks during activity minima. During these epochs the
noise dominates the light curve, which causes the very large $\alpha$ values.
During activity maxima, $\alpha$ lies mostly slightly below the average
(except for the last cycle). The $\alpha$ distribution of the noise-free
models also changes during the modeled period. However, this dependence is not
completely clear.

Figure 4: Dependence of $\alpha$ on the activity cycle. The gray data show the
underlying model light curve (with $90^{\circ}$ inclination and solar
metallicity) from which we picked 4-yr segments, computed the power spectrum,
its gradient and finally the inflection point period (i.e. $\alpha$). The blue
(noise-free) and red (noisy) curves shows the yearly averages of
$\alpha-\langle\alpha\rangle$ for 5000 randomly drawn samples. We subtracted
the mean value $\langle\alpha\rangle$ here to show the overlay with the
activity cycle.

###  3.3 Detection rate

We define the detection rate of a correct period detection as the number of
periods in the interval between 24–30 days divided by the number of samples
with a single inflection point. The question remains what value of $\alpha$
would be best for measuring the rotation periods, especially later in real
data. One choice would be to simply use the mean or median value for all
stars. For the noisy sample, we chose the median $\langle\alpha\rangle=0.209$
because the $\alpha$ distribution is not fully Gaussian. Another option would
be to use a magnitude-corrected value of $\alpha$. However, the functional
dependence of $\alpha$ on Kepler magnitude is not entirely clear (see bottom
panel in Fig. 3). First, we tried the simplest dependence, namely linear. We
also tested the robustness of the resulting period distribution and considered
an exponential fit, as the true dependence is not known. Eventually, we found
that the resulting rotation period distributions both are very similar. For
the fixed value $\alpha=0.209$, we find a detection rate of $39.7\pm 0.6\%$.
For the linear and exponential corrected $\alpha$ values instead, we find
slightly higher detection rates of $41.0\pm 0.6\%$ and $41.2\pm 0.6\%$,
respectively.

Given the fact that the same set of light curves has been analyzed in Reinhold
et al. (2021), this result can directly be compared to the detection rates of
auto-correlation functions. We find that the GPS method clearly outperforms
the auto-correlation functions, which yield very low detection rates from
$17\%$ for a local peak height (LPH) threshold of $LPH>0.1$ down to $1\%$ for
stricter criteria as $LPH>0.4$.

###  3.4 Effects of Spot Evolution

We now directly examine the effect of starspot lifetime on the efficacy of the
two methods for finding rotation periods. In our primary physical models spots
have a distribution of sizes and lifetimes (and smaller spots have shorter
lives, in agreement with sunspot observations, Solanki 2003). In contrast to
this, the models of Basri & Shah (2020) have spots with a single well-defined
lifetime. Spot areas grow linearly to a fixed maximum size over half the
lifetime then decay symmetrically in these latter models. We remind the reader
that these models have lifetimes in units of rotation period, so the rotation
period itself can be set to any arbitrary number and the inferred period will
scale accordingly. Here the period was arbitrarily set to 10 days simply for
display purposes. For each of the lifetimes tested we applied both period-
finding methods to each of the 1000 light curves produced. The histograms of
the derived periods for all cases are displayed in Fig. 5. The top panel shows
that for the GPS method, models where the spot lifetime was longer than one
rotation long did very well in finding the correct period. With $\alpha=0.183$
the mean period found for these models is 9.99 days with a 1$\sigma$
dispersion of less than 10% (see Table 1). The model with spot lifetime of one
rotation produced a mean period of 8.7 days at the same $\alpha$. This value
of $\alpha$ has been post-selected to match the model period of 10 days for
those models with lifetimes greater than 2 rotations. We emphasize that this
value is well within the $1\sigma$ range determined from the noise-free SATIRE
models (see Sect. 3.1).

Interestingly, the GPS-derived period (using the same $\alpha$) for a spot
lifetime of 0.5 rotations is a little under two-thirds of the actual rotation
period, and for a lifetime of 0.25 it is about one-third of the true period.
This demonstrates explicitly that the GPS result is sensitive to the typical
duration of a dip in the light curve. As a consequence, $\alpha$ does not
depend on the rotation period of the model light curve. When rotation is the
cause of the typical dip duration (spot evolution occurs on the rotational
timescale or longer) GPS delivers the actual rotation period. But when spot
evolution is the primary driver of the duration of the typical dip, GPS
responds to that more than to rotationally-driven features. Also of interest
is the fact that the GPS dispersion is twice as small (percentage-wise) for
the shortest lifetime than for the other cases. This may be because the
interaction between rotation and spot evolution is weakest in this case; spot
evolution completely dominates. This is most relevant for slowly rotating
stars, which are more likely to have spot lifetimes that are smaller fractions
of their rotation period. The implication is that some of the GPS-derived
periods for such stars could be shorter than the real periods. This analysis
suggests that for spot lifetimes of 2 periods and less previous analyses
(mostly using Kepler data) may not have determined the correct periods.

The results for the auto-correlation method are quite different. Although the
mean derived period is nearly correct for the shortest spot lifetime, that is
entirely coincidental given the nearly 40% 1$\sigma$ dispersion of the
individual periods (which range between 1 and 20 days). One can place no
confidence in any individual result under such circumstances. The mean derived
periods for lifetimes of a half and unit rotation periods are around twice
what they should be, with similarly large errors making them equally unusable.
At a spot lifetime of twice the rotation period the results for just under
half of the light curves begin to converge back to the correct period but the
majority of derived periods remain much too high and dispersed. It appears
that some of the light curves for lifetime two are periodic enough for auto-
correlation to detect it while more are not quite there, just due to random
variations in the spot patterns. We were unable to identify specific metrics
that distinguish the two groups (the light curves look very similar). The
auto-correlation method finally converges on just under the correct result
when the spot lifetimes are three or more rotation periods, and the dispersion
in the derived periods becomes very tight (significantly tighter than for
GPS). It can also be noted that the Basri & Shah (2020) models that have
lifetimes longer than 2 rotations produce light curves that are much like
those that McQuillan et al. (2014) were successful with. That this is true is
clear from the strong success of the auto-correlation method on them, and the
GPS method indeed found the same period.

Figure 5: Dependence of derived rotation periods on starspot lifetimes. Upper panel: histograms of derived rotation periods for a model set with different lifetimes utilizing the GPS method. Lower panel: the same but with periods derived using the auto-correlation method. The vertical dashed line in both panels indicates the model rotation period. Table 1: Results for Lifetime Models. Lifetime | P-GPS |  $\sigma$-GPS | %$\sigma$-GPS | P-ACH |  $\sigma$-ACH | %$\sigma$-ACH
---|---|---|---|---|---|---
0.25 | 3.35 | 0.16 | 4.8 | 9.84 | 3.89 | 39.5
0.5 | 6.12 | 0.55 | 9.0 | 18.3 | 6.08 | 33.2
1.0 | 8.73 | 0.84 | 9.6 | 25.5 | 7.25 | 28.4
2.0 | 9.91 | 0.85 | 8.6 | 19.9 | 10.42 | 52.4
3.0 | 10.08 | 0.81 | 8.0 | 9.72 | 1.51 | 15.5
4.0 | 9.99 | 0.78 | 7.8 | 9.81 | 0.19 | 1.9

##  4 Summary & Conclusions

Solar variability is known to be highly irregular on rotational timescales. In
this study, we employed physics-based models of stellar brightness variations
as these would be seen by the Kepler telescope for stars with different
inclination angles, metallicities, and noise levels. Previously, Reinhold et
al. (2021) used the same models to show that measuring accurate rotation
periods in such highly irregular light curves is a difficult task by itself,
and standard frequency analysis techniques such as auto-correlation functions
often fail to detect the correct period.

In this study, we took advantage of the GPS method (Shapiro et al., 2020) and
demonstrated that it is a powerful tool for detecting rotation periods in such
highly-irregular light curves. In particular, the detection rate of $\approx
40\%$ is much higher than the $\approx 3\%$ of the ACF method using common
peak height thresholds as $LPH>0.3$ (see Reinhold et al. 2021). We further
showed that the GPS method is largely insensitive to different inclinations
and metallicities but visibly reacts to noise, pushing the calibration factor
$\alpha$ toward larger values.

In summary, the distribution of $\alpha$ is broad which is only partly caused
by noise but also stems from the presence of the various manifestations of
magnetic activity, i.e. the interplay of dark spots and bright faculae (see
Fig. 4). Owing to the relatively large spread of the calibration factor
$\alpha$, the derived rotation periods are accompanied by errors on the order
of $\approx 20\%$.

Another direct comparison of the GPS and ACF methods showed that the ACF only
works in cases where the spot lifetimes are significantly longer than the
rotation periods (2 or more times). By contrast, the GPS method yields proper
rotation periods for spot lifetimes longer than one rotation period with
moderate relative uncertainties. For more periodic light curves, i.e. likely
those with long spot lifetimes, both methods converge and yield consistent
results, with ACF periods becoming increasingly accurate with increasing spot
lifetime (relative to the rotation period). Therefore, whereas for periodic or
quasi-periodic light curves, ACF may be the method of choice (giving somewhat
lower errors for the most periodic ones), while the GPS method is clearly
superior to the ACF method for finding periods in aperiodic light curves.

All these attributes make the GPS method a promising tool for determining
rotation periods of those Kepler stars lacking such a measurement so far. This
task will be addressed in a forthcoming publication. In particular, we will
apply the GPS method to all Kepler stars with near-solar effective
temperatures. This effort will generate the largest sample of solar-like stars
available so far. It will serve as test bed for various solar-stellar
comparison studies.

## References

  * Amazo-Gómez et al. (2020a) Amazo-Gómez, E. M., Shapiro, A. I., Solanki, S. K., et al. 2020a, A&A, 642, A225, doi: 10.1051/0004-6361/202038926
  * Amazo-Gómez et al. (2020b) —. 2020b, A&A, 636, A69, doi: 10.1051/0004-6361/201936925
  * Basri & Shah (2020) Basri, G., & Shah, R. 2020, ApJ, 901, 14, doi: 10.3847/1538-4357/abae5d
  * Basri et al. (2022) Basri, G., Streichenberger, T., McWard, C., et al. 2022, ApJ, 924, 31, doi: 10.3847/1538-4357/ac3420
  * Fligge et al. (2000) Fligge, M., Solanki, S. K., & Unruh, Y. C. 2000, A&A, 353, 380
  * Knaack et al. (2001) Knaack, R., Fligge, M., Solanki, S. K., & Unruh, Y. C. 2001, A&A, 376, 1080, doi: 10.1051/0004-6361:20011007
  * Krivova et al. (2003) Krivova, N. A., Solanki, S. K., Fligge, M., & Unruh, Y. C. 2003, A&A, 399, L1, doi: 10.1051/0004-6361:20030029
  * McQuillan et al. (2014) McQuillan, A., Mazeh, T., & Aigrain, S. 2014, ApJS, 211, 24, doi: 10.1088/0067-0049/211/2/24
  * Nèmec et al. (2020a) Nèmec, N. E., Işık, E., Shapiro, A. I., et al. 2020a, A&A, 638, A56, doi: 10.1051/0004-6361/202038054
  * Nèmec et al. (2020b) Nèmec, N. E., Shapiro, A. I., Krivova, N. A., et al. 2020b, A&A, 636, A43, doi: 10.1051/0004-6361/202037588
  * Petrovay & van Driel-Gesztelyi (1997) Petrovay, K., & van Driel-Gesztelyi, L. 1997, Sol. Phys., 176, 249, doi: 10.1023/A:1004988123265
  * Reinhold & Gizon (2015) Reinhold, T., & Gizon, L. 2015, A&A, 583, A65, doi: 10.1051/0004-6361/201526216
  * Reinhold et al. (2013) Reinhold, T., Reiners, A., & Basri, G. 2013, A&A, 560, A4, doi: 10.1051/0004-6361/201321970
  * Reinhold et al. (2020) Reinhold, T., Shapiro, A. I., Solanki, S. K., et al. 2020, Science, 368, 518, doi: 10.1126/science.aay3821
  * Reinhold et al. (2021) Reinhold, T., Shapiro, A. I., Witzke, V., et al. 2021, ApJ, 908, L21, doi: 10.3847/2041-8213/abde46
  * Santos et al. (2021) Santos, A. R. G., Breton, S. N., Mathur, S., & García, R. A. 2021, ApJS, 255, 17, doi: 10.3847/1538-4365/ac033f
  * Santos et al. (2019) Santos, A. R. G., García, R. A., Mathur, S., et al. 2019, ApJS, 244, 21, doi: 10.3847/1538-4365/ab3b56
  * Shapiro et al. (2020) Shapiro, A. I., Amazo-Gómez, E. M., Krivova, N. A., & Solanki, S. K. 2020, A&A, 633, A32, doi: 10.1051/0004-6361/201936018
  * Solanki (2003) Solanki, S. K. 2003, A&A Rev., 11, 153, doi: 10.1007/s00159-003-0018-4
  * Solanki et al. (2013) Solanki, S. K., Krivova, N. A., & Haigh, J. D. 2013, ARA&A, 51, 311, doi: 10.1146/annurev-astro-082812-141007
  * Torrence & Compo (1998) Torrence, C., & Compo, G. P. 1998, Bulletin of the American Meteorological Society, 79, 61, doi: 10.1175/1520-0477(1998)079<0061:APGTWA>2.0.CO;2
  * Witzke et al. (2020) Witzke, V., Reinhold, T., Shapiro, A. I., Krivova, N. A., & Solanki, S. K. 2020, A&A, 634, L9, doi: 10.1051/0004-6361/201936608