[1]Dimitris Karlis

[1]Department of Statistics, Athens University of Economics and Business,
Athens, Greece

2]Department of Business Administration and Economics, Bielefeld University,
Bielefeld, Germany

# Modelling handball outcomes using univariate and bivariate approaches

<EMAIL_ADDRESS>Rouven Michels<EMAIL_ADDRESS>Marius Ötting
<EMAIL_ADDRESS>*  [

###### Abstract

Handball has received growing interest during the last years, including
academic research for many different aspects of the sport. On the other hand
modelling the outcome of the game has attracted less interest mainly because
of the additional challenges that occur. Data analysis has revealed that the
number of goals scored by each team are under-dispersed relative to a Poisson
distribution and hence new models are needed for this purpose. Here we propose
to circumvent the problem by modelling the score difference. This removes the
need for special models since typical models for integer data like the Skellam
distribution can provide sufficient fit and thus reveal some of the
characteristics of the game. In the present paper we propose some models
starting from a Skellam regression model and also considering zero inflated
versions as well as other discrete distributions in $\mathbb{Z}$. Furthermore,
we develop some bivariate models using copulas to model the two halves of the
game and thus providing insights on the game. Data from German Bundesliga are
used to show the potential of the new models.

###### keywords:

Copulas, handball, inflation modelling, Skellam regression, sports analytics

##  1 Introduction

Handball is a dynamic sport, which is particularly popular in Europe and South
America. In January 2024, the European Men’s Handball Championship was held in
Germany, where more than 50,000 spectators attended the opening match.
Handball combines elements of football (soccer) and basketball, played on a
rectangular court with two teams of seven players each. A handball match
consists of a sequence of possessions, where each possession terminates with
either a score or a loss of possession.

Handball is an invasion game, played for fixed time. Each match lasts for 60
minutes (two halves of 30 minutes each), and is being played by two teams. The
ball possession ends with a field shot attempt, or a seven meter throw, that
is not rebounded by the offence (possession continues after an offensive
rebound), or in case of a turnover. The dynamics of the match process are
revealed due to the difference of the players’ attack and defence efficacy
during the match, which influence the scoring rate of each team.

While there is an increasing steam of research on handball, models for final
outcomes are less developed. Handball, contrary to football, is not a low-
score sport and hence modelling its outcome can be challenging. For football,
there are several proposals for modelling the score of a game ([17, 5, 19]),
mostly based on an assumption of a Poisson distribution for the number of
goals. In the case of handball, a Poisson distribution for the number of goals
for each team can be a starting assumption. Though, the number of goals often
has a variance smaller than the mean making the Poisson assumption
questionable. An alternative way is to assume a Gaussian distribution as an
approximation of the Poisson distribution, thereby ignoring the discrete
nature of the data.

There are only a few existing models on handball scores. Groll et al. [9] used
a regularised Poisson regression model together with over- and underdispersed
models. They also argue that since the Poisson mean is large enough, the model
can be approximated by a standard Gaussian model (with sparsity also). Smiatek
and Heuer [27] analyse team handball results of the German Bundesliga. Their
data set consists of ten seasons starting from 2001-02. While they do not have
a particular model, the descriptive statistics indicate a small but consistent
underdispersion. They also claim that the goal distribution can be described
as a Gaussian distribution. Recently, Felice and Ley [8] and Felice [7]
proposed a Conway-Maxwell-Poisson distribution. This is a generalisation of
the Poisson distribution that also allows for under-dispersion at the cost of
added computational cost and rather approximating probability calculations.
They also described some machine learning approaches to model the outcome and
not the score. Lago-Peñas et al. [16] tried to model the goal difference using
a simple regression model to examine how some technical characteristics affect
the final goal difference but also to examine the home advantage. A similar
approach has been used in Prieto et al. [24]. In contrast to modelling match
results, there is also some work on in-game modelling. Dumangane et al. [6]
mention that the dynamics of handball matches violate both the assumption of
independence and identical distribution, in some cases having a non-stationary
behaviour. Singh et al. [25] extended this model. They describe what is called
the Markov-match, where the possession of the ball alternates between
competitors and the probability of scoring, while different for each team, is
quite high — in fact, it can be up to $0.7$ for some teams. Thus, the number
of goals cannot be considered as a binomial random variable (further
approximated by a Poisson, if the goals are rare), since the individual
Bernoulli trials (the attacks) are not independent.

To model the winner of a handball match, we aim to overcome the problem of
underdispersion by proposing alternative approaches. In particular, we
consider the difference between the scores of the two competing teams. This
difference, which is a number in $\mathbb{Z}=\\{\ldots,-2,-1,0,1,2,\ldots\\}$,
can be modelled as a realisation of the Skellam distribution. The Skellam
distribution is a probability distribution that is frequently used in
analysing the difference between two independent Poisson-distributed random
variables. In this paper, we use the Skellam distribution to model the final
score difference of handball matches. However, as the final outcomes of
matches exhibit special patterns, we aim to explicitly address these in two
different ways. First, we consider a zero-inflated (ZI) Skellam distribution
which shifts probabilities to the zero from the other values, as well as other
discretised distributions defined in $\mathbb{Z}$ . Second, we investigate the
score difference of the first and second half separately as those differences
do not indicate systematic deviances. We model these differences jointly by
constructing bivariate distributions using copulas. Modelling the scores of
the first and second halves in such a way enables us to calculate conditional
probabilities of the final score based on the score at half time. Through this
approach, we aim to contribute to the evolving landscape of sports analytics,
providing a useful framework for understanding and predicting the outcomes in
handball.

The rest of the paper is structured as follows. Section 2 introduces the
univariate and bivariate distributions we use. These are applied in Section 3,
which presents models fitted to handball scores from the Bundesliga in
Germany. Concluding remarks can be found in Section 4.

##  2 Methods

This section first introduces the Skellam distribution, which is then extended
to obtain the zero-inflated Skellam distribution. Some other candidate models
like the discretized versions of normal and Laplace distributions are also
described. At the end of the section, we consider copula-based models with
Skellam marginals.

###  2.1 Skellam distribution

In handball, the score differences are represented as negative and positive
integers, that is values in $\mathbb{Z}$. Perhaps the most well-known
distribution in $\mathbb{Z}$ is the Skellam distribution. Irwin [11] derived
the distribution of the difference of two independent variables following the
Poisson distribution with the same mean, and Skellam [26] generalised this to
different means.

If the random variables $X$ and $Y$ follow independent Poisson distributions
with parameters $\theta_{1}>0$ and $\theta_{2}>0$, then the random variable
$Z=X-Y$ has probability mass function (pmf) given by

$P(Z=z|\theta_{1},\theta_{2})=e^{-(\theta_{1}+\theta_{2})}\left(\frac{\theta_{1}}{\theta_{2}}\right)^{z/2}I_{|z|}\left(2\sqrt{\theta_{1}\theta_{2}}\right),~{}~{}z\in\mathbb{Z},~{}~{}\theta_{1},\theta_{2}>0,$
(1)

where $I_{r}(x)$ is the modified Bessel function of order $r$ [1, pp. 375]

$I_{r}(x)=\left(\frac{x}{2}\right)^{r}\sum\limits_{m=0}^{\infty}\
\frac{\left(\frac{x^{2}}{4}\right)^{m}}{m!\Gamma(r+m+1)}.$

We will denote this distribution as the Skellam$(\theta_{1},\theta_{2})$
distribution. The Skellam distribution has also been called the Poisson
difference distribution. Its mean and variance are
${\mathrm{E}}(Z)=\theta_{1}-\theta_{2}$ and
$\mbox{Var}(Z)=\theta_{1}+\theta_{2}$. This can be used to reparameterize the
distribution; clearly, $\mbox{Var}(Z)\geq|{\mathrm{E}}(Z)|$. The skewness is
given as $(\theta_{1}-\theta_{2})/(\theta_{1}+\theta_{2})^{3/2}$. and it is
determined by the sign of $\theta_{1}-\theta_{2}$; in the case of
$\theta_{1}=\theta_{2}$, we come up with the symmetric Skellam distribution.
For large values of $\theta_{1}+\theta_{2}$, the distribution can be well
approximated by the normal distribution. If $\theta_{2}$ is very close to
zero, then the distribution tends to a Poisson distribution. If the parameter
$\theta_{1}$ approaches zero, then the distribution is the negative of a
Poisson distribution. The Skellam distribution is unimodal.

We can derive the Skellam distribution as the difference of other
distributions as well [14], which motivates its use in various applications.
In particular, it can be derived as the difference of two underdispersed
random variables as follows: Let $X_{1}$ and $X_{2}$ be two Poisson random
variables and $W$ be an underdispersed random variable, i.e. with variance
less that the mean, independently of $X_{1}$ and $X_{2}$. Consider the random
variables $Y_{1}=X_{1}+W$ and $Y_{2}=X_{2}+W$. Then, $Y_{1}$ and $Y_{2}$ will
be underdispersed. If we take their difference, we obtain
$Y_{1}-Y_{2}=X_{1}+W-X_{2}-W=X_{1}-X_{2}$ which is a Skellam random variable
being the difference of two Poisson variables. This makes the assumed Skellam
model as a reasonable approach for underdispersed data, such as the data on
handball considered in this work. Note that we do not need to assume some
parametric assumption for the two random variables, just for their difference,
as in the above definition $W$ can follow any underdispersed distribution.

To allow for better interpretation of the parameters but also more advanced
modelling approaches such as regression, Koopman et al. [15] use a
reparametrised version of the distribution with mean
$\mu=\theta_{1}-\theta_{2}$ and variance $\sigma^{2}=\theta_{1}+\theta_{2}$.
We will denote this by Skellam2$(\mu,\sigma^{2})$ in the following. Here, one
can add covariates in the mean as a simple linear model since the mean is
defined in $\mathbb{R}$, see Pelechrinis and Winston [23], Holmes and McHale
[10]. To cope with data for which the regular Skellam distribution lacks a
proper fit, the baseline model has been extended in various ways. Karlis and
Ntzoufras [14] consider a zero-inflated version, while a zero-deflated version
is given in Koopman et al. [15]. On a similar note, Ntzoufras et al. [21]
propose a truncated version of the Skellam distribution. A mixture of Skellam
distributions has been used by Jiang et al. [12] to cluster differentially
expressed genes. A modified version where probability mass is exchanged
between certain values is given in Koopman et al. [15]. For a recent review on
the Skellam distribution see Tomy and Veena [28], Karlis and Mamode Khan [13].

###  2.2 Zero-inflated Skellam distribution

From the extensions summarised above, we consider the zero-inflated version in
the following, since draws are overrepresented the handball data
(corresponding to a score difference of zero). The pmf of the zero-inflated
Skellam distribution is given by

$P_{Z}(Z=z|\theta_{1},\theta_{2},p)=\left\\{\begin{array}[]{cc}p+(1-p)P(Z=z|\theta_{1},\theta_{2})&\mbox{if}~{}~{}z=0,\\\
(1-p)P(Z=z|\theta_{1},\theta_{2})&\mbox{if}~{}~{}z\neq 0.\end{array}\right.$

with $p\in[0,1)$ and $P(\cdot)$ being the pmf of the Skellam distribution
given in (1). Thus, the zero-inflated Skellam distribution shifts some
probability from the other values towards zero.

The model can be fitted by the EM algorithm, which is a standard choice for
fitting zero inflated models. At first, we are not using covariates for
modelling the parameter $p$. We thus make a rather unnatural assumption that
each match, despite the difference between the teams’ abilities, has a
constant probability to end in a draw. However, we believe that this
assumption makes sense to avoid overparametrising the model.

A model that allows for covariates on the inflation parameter is also possible
by assuming that for match $i$

$p_{i}=\frac{\exp({\bf z}_{i}^{\prime}\boldsymbol{\gamma})}{1+\exp({\bf
z}_{i}^{\prime}\boldsymbol{\gamma})}$

where ${\bf z}_{i}$ is a vector with explanatory variables related to match
$i$ and $\boldsymbol{\gamma}$ a coefficient vector to be estimated.

###  2.3 Other univariate distributions in $\mathbb{Z}$

Alternative models for modelling data in $\mathbb{Z}$ can be also considered
either as the difference of discrete random variables or by discretizing
continuous ones (for an overview of such models see [13]). For example, a
discrete normal distribution can be considered based on a
$\mathcal{N}(\mu,\sigma^{2})$ distribution as

$\displaystyle P_{N}(Z=z|\mu,\sigma^{2})$ $\displaystyle=$
$\displaystyle\int_{z-0.5}^{z+0.5}\phi(t;\mu,\sigma^{2})dt$ $\displaystyle=$
$\displaystyle\Phi(z+0.5;\mu,\sigma^{2})-\Phi(z-0.5;\mu,\sigma^{2})$

where $\phi(z;\mu,\sigma^{2})$ and $\Phi(z;\mu,\sigma^{2})$ are the density
and the cumulative distribution function of a $\mathcal{N}(\mu,\sigma^{2})$,
respectively. An important remark, which is also a motivation for our
approach, is that the Skellam allows for skewness which is not the case for
the discrete normal.

Here, the same parameterisation as for the Skellam distribution is used for
$\mu$. In particular, $\mu_{i}$ is related to some covariates that reflect the
ability of the teams. Note that we here consider the discrete normal which is
a different approach than fitting a standard regression and then to simply
discretise the results of the regression.

Another option is to use the discretised version of the Laplace distribution
[2]. There are different ways to create a discretised version of the Laplace,
here we use the simple idea, as for the normal above, to take the cumulative
function of the continuous case at specific points. Thus, we have a pmf
defined as

$P_{L}(Z=z|\mu,\sigma^{2})=F(z+0.5;\mu,\sigma^{2})-F(z-0.5;\mu,\sigma^{2})$
(2)

where $F(x;\mu,\sigma^{2})$ is the cumulative distribution of a Laplace
distribution with location $\mu$ and scale $\sigma^{2}$. Such a model provides
more mass at zero, recall the shape of the Laplace distribution.

###  2.4 Bivariate distributions with copulas

Figure 1:  Marginal distributions are Skellam $\mu_{1}=0.5$ $\mu_{2}=0.5$,
$\sigma_{1}^{2}=15$ , $\sigma_{2}^{2}=14$, which are reasonable values for the
case study of handball scores, and varying copula parameter for the Frank
copula. The dependence parameter values are 10, 3, -3 and -10.

To provide a more thorough analysis of the results of handball matches, the
score difference can also be modelled for each half separately. In this case,
denote as $Y_{1}$ and $Y_{2}$ the score difference for the first and second
half. To jointly model such score differences, we need a bivariate
distribution defined in ${\mathbb{Z}}^{2}$. We can construct such a
distribution accounting for dependencies between both halves using copulas. In
particular, we assume that marginally both of the $Y_{1}$ and $Y_{2}$ follow a
Skellam distribution with parameters $\mu_{j},\sigma_{j}^{2},\ j=1,2$, coupled
with some copula $C(\cdot,\cdot;\theta)$, where $\theta$ is the dependence
parameter.

For the discrete case, we need to take differences to derive the joint pmf.
For example, the bivariate case with marginals $F(y_{1})$ and $G(y_{2})$ one
can derive the joint pmf as

$\displaystyle P(Y_{1}=y_{1},Y_{2}=y_{2})$ $\displaystyle=$ $\displaystyle
C\big{(}F(y_{1}),G(y_{2});\theta\big{)}-C\big{(}F(y_{1}-1),G(y_{2});\theta\big{)}-$
$\displaystyle
C\big{(}F(y_{1}),G(y_{2}-1);\theta\big{)}+C\big{(}F(y_{1}-1),G(y_{2}-1);\theta\big{)}$

We are going to use Skellam marginal distributions linked through the
bivariate Frank copula, which is given by

$C(u,v;\theta)=-\frac{1}{\theta}\log{\left[1+\frac{(\exp^{-\theta
u}-1)(\exp^{-\theta v}-1)}{(\exp^{-\theta}-1)}\right]}$

and a Gumbel copula with cdf

$C(u,v;\theta)=\exp\left(-((-\log(u))^{\theta}+(-\log(v))^{\theta})^{1/\theta}\right).$

We select the two copulas mainly to be as flexible as necessary since it
allows for both negative and positive correlation (Frank copula) and for tail
dependence (Gumbel copula). In our case, tail dependence means that for larger
differences we expect to occur both in the first and second half together,
implying also that they can occur more often than a model without (or small)
tail dependence.

However, of course any other copula could also be used. The approach further
applies for any marginal distribution in $\mathbb{Z}$. Marginal properties of
the model are based on the univariate Skellam. Figures 1 and 2 depict the
joint probability function with Skellam marginals and a Frank and Gumbel
copula respectively with various parameters. We have used marginal
distribution close to what we observe in handball and varying dependence
structure including both positive and negative dependence.

Figure 2:  Marginal distributions are Skellam $\mu_{1}=0.5$ $\mu_{2}=0.5$,
$\sigma_{1}^{2}=15$ , $\sigma_{2}^{2}=14$ and varying copula parameter for the
Gumbel copula. The values are 1 (independence),2,3 and 5.

For a broad review on bivariate Skellam and related distributions, see Karlis
and Mamode Khan [13].

##  3 Application to German Handball Bundesliga

In this section, we first introduce the data considered to fit the models
presented in the previous section. Afterwards, we report the estimation
results and subsequent analyses for the univariate and bivariate
distributions.

###  3.1 Data

Our data consist of the scores in the German Handball Bundesliga (HBL) for the
seasons 2017-18 up to 2022-23 which comprise, in total, 1844 matches. The HBL
is the top German professional handball league with 18 teams competing in a
round-robin format. The bottom two teams are relegated to the
2${}^{\text{nd}}$ Bundesliga, while the top two teams of the 2${}^{\text{nd}}$
Bundesliga will be directly promoted to the Bundesliga.

Figure 3:  Observed frequencies of the score differences in the season
2017-18-2022-23.

The observed means and variances for the home goals are $28.07$ and $19.67$,
and $26.92$ and $18.32$ for the away goals, thus implying a small home
advantage, but also underdispersion. For an attempt to model such
underdispersion see Felice [7], who uses a COM Poisson distribution fitted to
the home and away goals. In addition, the observed correlation between the
home and away goals in our data is relatively small (0.14). This is somewhat
counter-intuitive, since in handball the two teams are in offense
sequentially, hence one could expect to observe a higher correlation between
the number of goals scored. To investigate whether such patterns change over
time, and also to test the stability and robustness of the models, we
investigate the seasons considered in our analysis separately. Table 1
displays the season-by-season average goal difference per game, their standard
deviations as well as the skewness of the goal difference. The latter can be a
motivation for considering the Skellam model as we observe some skewness. Note
that season 2019-20 was cancelled after 240 matches due to COVID-19 and,
consequently, 20 teams participated in the 2020-2021 season instead of 18 (in
total 380 matches). For both seasons, we can observe a drop in the mean which
corresponds to a smaller margin of home wins on average. This coincides with
the well-known finding that the home advantage was mitigated during games
behind closed doors ([3]).

Season | Mean | SD | Skewness
---|---|---|---
2017-18 | 1.28 | 5.83 | -0.06
2018-19 | 1.21 | 5.75 | 0.04
2019-20 | 1.33 | 5.04 | 0.03
2020-21 | 1.07 | 5.65 | 0.19
2021-22 | 0.71 | 5.90 | 0.04
2022-23 | 1.35 | 5.96 | -0.10
Table 1:  Descriptive statistics for the score differences in the six years of
data.

###  3.2 Univariate analysis

In this section, we consider the mean/variance parametrization of the Skellam2
model such that we are able to add covariates. We model the mean goal
difference $\mu_{jk}$ between the two teams, i.e.

$\displaystyle\mu_{jk}$ $\displaystyle=$
$\displaystyle\alpha+\beta_{j}+\gamma_{k},$

where team $j$ plays at home against team $k$. In this model formulation,
$\beta_{j}$ is the ability of the $j$-th team at home while $\gamma_{j}$ is
the ability away. The team “Bergischer HC” has been used as the baseline team.
So the interpretation of $\beta_{j}$ is that when team $j$ plays at home with
Bergischer the expected score is $\alpha+\beta_{j}$ while if they play away
(so at home of Bergischer) it is $\alpha+\gamma_{j}$. Their difference
$\beta_{j}-\gamma_{j}$ measures the home advantage for the $j$-th team. The
model explicitly assumes a different home advantage for each team. Home
advantage has been discussed for handball in many instances, see for example
the work in [30]. We employ the same parametrisation of $\mu_{jk}$ also for
the zero-inflated version as well as for the discrete normal and discrete
Laplace distributions.

We have considered a regular Skellam model, followed by the zero-inflated
version. The latter has an extra additional parameter $p$ that reflects the
extra zeros we have observed. In the case of handball, the idea behind zero
inflation is that teams that are close to the score at the end of the game do
not want to risk and prefer the safe path of a draw. Recall that draws are not
common in handball as they are in football but it is still a possible outcome.
Moreover, we try other discrete distributions as well, such as the discrete
normal or the Laplace distribution. Finally we have also considered models
with covariates in $p$, i.e. the zero inflation parameter. One reason for that
is that using a common $p$ for all matches is somewhat unrealistic. An
interpretation of zero inflated models is some additional zeros, i.e. draws in
our case, that cannot be predicted by the model. It seems irrelevant to assume
this when we have teams of very different strength to be the same as for the
case of equal strength teams. For this reason some covariate may differentiate
this inflation probability. We have used as a covariate the betting odds at
the beginning of the game as an indication of the strength difference of the
two teams. The results did not show any statistical significance and the log-
likelihood improved very little. For this reason we skip them.

####  3.2.1 Model fit

| 2017-18 | 2018-19 | 2019-20 | 2020-21 | 2021-22 | 2022-23
---|---|---|---|---|---|---
Skellam | 1755.78 | 1743.45 | 1387.50 | 2245.23 | 1854.39 | 1813.42
ZI Skellam | 1751.07 | 1745.45 | 1389.50 | 2244.96 | 1850.77 | 1815.42
Discrete normal | 1756.25 | 1743.92 | 1386.47 | 2245.31 | 1855.83 | 1814.38
Discrete Laplace | 1765.62 | 1758.51 | 1376.45 | 2253.07 | 1854.42 | 1820.88
Table 2:  AIC results from fitting several models to the data from different
periods as well as the resulting ZI parameter in the ZI Skellam model.

We have fitted the Skellam model, the zero-inflated Skellam, the discrete
normal and the discrete Laplace. When fitting these models to the data, we
consider a season-by-season approach to check whether the discrepancy seen in
Figure 3 is not primarily driven by a single season. Table 2 displays AIC
values for the models in the different seasons. The zero-inflated model
provides slightly better log-likelihood at the cost of an additional
parameter, so judging from AIC it is not a preferable model in every season.
The discrete normal model provides slightly worst fit due to the small
skewness present in our data. Note that in 2019-20 season we observe a
restricted number of observations due to the cancelled seasons, which is the
reason why the structure of that season is least reliable. In all the other
seasons, either the Skellam or the zero-inflated Skellam are chosen, while the
latter one automatically includes the general Skellam distribution which is
why from a log-likelihood point of view, the zero-inflated Skellam is chosen
for every full played season.

####  3.2.2 Goodness-of-fit

In this subsection, we evaluate the goodness-of-fit. For this, we consider the
first season that was held without any COVID-19 restrictions, i.e. the 2021-22
season.

| home wins | draw | away wins
---|---|---|---
Observed | 153.00 | 29.00 | 124.00
Bets implied | 160.63 | 28.22 | 117.14
Skellam | 157.45 | 20.66 | 127.89
ZI Skellam | 152.02 | 30.39 | 123.59
Discrete Normal | 157.55 | 20.40 | 128.05
Discrete Laplace | 159.95 | 21.35 | 124.70
Table 3:  Expected outcomes over the 306 matches of the 2021-22 season based
on different models and bookmaker’s odds.

To evaluate the goodness-of-fit of the models, we proceed as follows: Based on
the fitted model, we predict for each match the probabilities of a home win, a
draw and an away win. We sum these predicted probabilities over all matches
and compare them to the corresponding actual frequencies found in the data. We
further consider betting odds111Betting odds were taken from
https://www.betexplorer.com/handball/germany/bundesliga-2017-2018/results/,
which were transformed to probabilities by assuming an equal vig for the three
outcomes. Table 3 shows the results. The baseline Skellam model predicts ten
fewer draws than we observe in the data, which is not negligible. Thus, we can
summarise that the fit of the Skellam is sufficient but it has some failure
for the draws, which can be corrected with the zero-inflated Skellam model.
The discrete normal is somewhat inferior as we have seen. This is confirmed
when comparing the models to the predictive performance of bookmakers: here,
the zero-inflated Skellam can not only keep up with the bookmaker’s accuracy
for draws (which is kind of the targeting outcome for the ZI Skellam) but is
also able to predict home- and away wins more precisely than bookmakers. This
is particularly astonishing considering that bookmakers are able to include
short-term news such as injuries of important players into their odds which we
did not include in the ZI Skellam.

Figure 4:  The observed numbers of points (blank points), the expected number
of points through simulations (filled points) and 95% confidence intervals
around.

As we have now seen that the zero-inflated Skellam model is superior in terms
of the AIC in the most seasons and the prediction of draws in the 2021-22
season, we perform the remaining goodness-of-fit checks using this model.

Figure 4 shows the observed number of points for each team, the expected
number of points based on 10,000 simulations of the league and the 95%
confidence interval as obtained via parametric bootstrap and use of the
percentile method. According to our results, SC Magdeburg does not seem to be
a clear champion, since Kiel also had also reasonable chances to win the
championship — however, SC Magdeburg won the championship by a margin of 6
points.

Figure 5:  Empirical distribution function of the observed score differences
and a 95% confidence interval under the ZI Skellam model.

As another goodness-of-fit check, we again consider the 10,000 simulated
championship results and investigate the corresponding empirical distribution
function. Based on all the replications, we estimated a 95% confidence
interval for $P(Y\leq y)$, which is indicated by the gray area in Figure 5.
The actual empirical cumulative distribution function (cdf) is indicated by
the black dots, demonstrating that there is only a minor deviation compared to
what we expect under the fitted model.

####  3.2.3 Out-of-sample prediction

Figure 6:  The expected points (filled) together with 95% confidence intervals
in the COVID-19 season if the season would have been finished under the
prediction using the zero-inflated Skellam. The blank points indicate the
number of points when the season was interrupted.

In this section, we present an analysis showcasing the feasibility of out-of-
sample prediction, using the context of the disrupted COVID-19 season as an
illustrative case study. Following the approach outlined by Csató [4] and Van
Eetvelde et al. [29], we employ a Monte Carlo simulation comprising 10,000
iterations to compute the probabilistic end-of-season standings under the
assumption of a completed regular season. The model was fitted to the data
available until the stop and then the estmated parameters were used for the
rest of the season. Figure 6 depicts the hypothetical end-of-season table,
while Figure 4 illustrates the expected remaining points to be accrued.
Notably, despite a significant number of fixtures remaining unplayed, the
expected season winner aligns with the outcome derived from terminating the
season prematurely. This observation is congruent with the findings presented
in Table 4, which displays the likelihood of winning the league and being
relegated following the season’s interruption, as determined by the Monte
Carlo simulation.

Furthermore, our analysis extends to examining the hypothetical relegation
scenarios (notwithstanding the absence of relegation in the respective year).
Of particular interest is the favourable outcome for Nordhorn-Lingen, who
consistently faced relegation across all simulation iterations, but were not
relegated at the end.

Moreover, the league administration faced the task of determining eligibility
for participation in international competitions in the subsequent year.
Initially, all teams up to Fuchse Berlin received international tickets,
presenting a unique circumstance for MT Melsungen, whose projected point came
close to the ones of Fuchse Berlin. However, an opportune withdrawal by
Hannover-Burgdorf facilitated MT Melsungen’s ascension to an international
slot and, thus prevented discussions about the fairness of the league’s
decision.

team | Champion | Relegation
---|---|---
Kiel | 0.972 | 0.000
Flensburg-H. | 0.028 | 0.000
Stuttgart | 0.000 | 0.001
Erlangen | 0.000 | 0.002
Minden | 0.000 | 0.006
HBW Balingen-Weilstetten | 0.000 | 0.273
Ludwigshafen | 0.000 | 0.719
Nordhorn-Lingen | 0.000 | 1.000
Table 4:  Probabilities for winning the championship/being relegated if the
2019-20 season would have been continued. For all other teams both
probabilities were zero.

###  3.3 Bivariate analysis

In sports analytics, copulas have been used to model jointly the number of
goals scored in soccer [see, e.g. 18]. Here we apply them for modelling the
difference in the number of goals between the two halves.

An interesting result about the Skellam distribution is an additive property:
the sum of Skellam random variables is also a Skellam. In a similar manner,
for the Skellam distribution we have that if a random variable $Y$ follows a
Skellam in the unit of time, then it follows a Skellam for any other unit of
time. In our bivariate analysis this simply implies that we may assume that
each half of the match the score difference also follows a Skellam
distribution. In the following, we investigate whether the scoring abilities
remain the same across the two halves of the game.

Figure 7:  Correlation between the score differences in the two separated
halves for all seasons. We do observe a slightly positive correlation (0.13),
rendering the Frank copula as suited. Figure 8:  Histogram of the score
differences in the two separated halves for all seasons. We do not observe any
structural errors, thus, the Skellam seems to be suited as the marginal
distribution for the Copula.

We consider the data from the German Bundesliga 2021-22 to model the score
differences in the first and second half. For the copulas, we consider a Frank
and a Gumbel copula.

As for the univariate case, we extend the basic bivariate model to include
covariates. Again, we work with two random variables $Y_{j1}$ and $Y_{j2}$,
where each follows marginally a Skellam2 distribution with mean $\mu_{j}$ and
variance $\sigma_{j}^{2}$, $j=1,2.$ Mean parameters are related to some
covariate vector ${\bf Z}_{ji}$ with $p$ elements. Clearly each variable can
be related to a different set of covariates. We also assume a linear
relationship,

$\displaystyle\mu_{1i}$ $\displaystyle=$ $\displaystyle\beta_{1}^{T}{\bf
Z}_{1i}$ $\displaystyle\mu_{2i}$ $\displaystyle=$
$\displaystyle\beta_{2}^{T}{\bf Z}_{2i},$

where $\beta_{j}$ are the regression coefficient vectors $j=1,2$. We further
model any dependence between the two variables using a copula (Frank and
Gumbel in our case) with parameter $\theta$, which may also depend on
covariates. Based on this parametrization, we need to estimate
$\Theta=(\beta_{1},\beta_{2},\sigma_{1}^{2},\sigma_{2}^{2},\theta)$.

Such bivariate models offer a great flexibility for modelling handball
results: we can, for example, estimate the probability that a team will win,
the score difference at the end of the match or even conditional probabilities
that a team wins the match given a score difference at half time.

We have fitted different models, thereby gradually increasing the model’s
complexity. Model A assumes independence, i.e. the two variables are fitted
separately without a copula. For each one we assume different effects, so we
assume that the abilities of the teams are different. Model B assumes
dependence through a Frank and Gumbel copula but the abilities are the same
for the two halves. Model C is the most flexible one, as it allows for
differences between the two halves resulting to a much larger number of
parameters. Table 5 displays the results. Model B, i.e. the model with a Frank
copula and the same parameters in the two halves is the preferred one. The
Frank copula may be preferred over the Gumbel as we do not observe any tail
dependence in the data.

Model |  | log-likelihood | #par | AIC
---|---|---|---|---
Model A |  | -1590.047 | 72 | 3324.09
Model B |  |  |  |
| Gumbel | -1603.235 | 38 | 3282.47
| Frank | -1601.736 | 38 | 3279.47
Model C |  |  |  |
| Gumbel | -1587.511 | 73 | 3321.02
| Frank | -1586.922 | 73 | 3319.84
Table 5:  Log-likelihood of fitted models of different complexity.

The results support the parsimonious model that the parameters representing
the abilities of the teams in home and away matches are the same. To further
investigate this finding, Figure 9 shows the estimated parameters for the
first and the second half. One can see that for most of the teams the
parameters for the first and second half are very close supporting that we can
be more parsimonious by considering common parameters.

Figure 9:  Estimated parameters for the 1${}^{\text{st}}$ and
2${}^{\text{nd}}$ half for home and away parameters. For most of the teams,
the estimated parameters are quite close. First half estimates are indicated
by the solid points, 2${}^{\text{nd}}$ half estimates are given by the open
points.

###  3.4 Conditional distributions

In this section, we consider the case that we know the half time result and we
want to find the probabilities for the final score difference. Based on the
fitted bivariate model, we can calculate the conditional probabilities

$P(Y=y|X=x)=\frac{P(Y=y,X=x)}{P(X=x)},$

where the nominator is a bivariate Skellam through a copula and the
denominator is the marginal Skellam distribution. From this we calculate the
probability for winning the game at the end. So if team, whose goals are
denoted by $Y$, is trailing by say $x$ goals at half time, then the
probability of winning is

$P(\text{Win})=P(Y>-x|X=x)=\sum\limits_{y=-x+1}^{\infty}P(Y>x|X=x)$

Figure 10 presents the conditional probabilities that the home team will win
at the end of the game conditional on the score of the first half, for the
first matchday of the 2021-2022 season. For the prediction of these
probabilities, we have used the model that was selected in Table 5. So, for
example, looking at the top left plot, we see that at the half time the score
between Erlangen and Leipzig was simply 8-8, and based on the team abilities
estimated we calculated that the probability that Erlangen will win is 0.49,
while the entire conditional distribution can be also seen.

Figure 10:  Conditional distribution for $Y|X=x$, for all the matches of the
first day. The plots show the conditional distribution given the half time
score for the difference in the final score. We also report the probability
that the home team to win. The half time score can be read in the caption of
each plot. The vertical line refers to the probability of a draw.

##  4 Discussion

In this study, the outcomes of handball results were modelled in several ways.
Initially, we decomposed the bivariate pair of a result by transforming it
into a univariate statistic, i.e. the score difference. We used the Skellam
distribution as the baseline model for modelling this difference. Considering
the patterns found in the data — an excess of draws — we considered a zero-
inflated version of the baseline model. Additionally, we compared both models
with a discretised normal and Laplace distributions. Modelling the score
difference overcomes the problem to specify the distributions for the number
of goals, which is not easy since we need underdispersed discrete
distributions.

In a second step, we addressed the problem of modelling the score differences
of the two halves in a handball match separately. Here, we employed copulas to
model the dependence between both halves’ score differences. Such bivariate
model formulations can be used to predict the final outcome as a conditional
probability, i.e. based on the score difference at half time.

The models provided in the present paper can serve as the basis for potential
further model extensions. For example, they can be the basis for modelling
certain characteristics of the game, such as the “7 vs 6” rule introduced in
2016, where a team can substitute the goalkeeper to have more players in
offence. This is a risky situation that has not been examined in detail in the
literature. To investigate the effect of replacing the goalkeeper for an
additional outfield player, one could consider the models developed here — in
particular, to estimate the effect on the score difference associated with
such strategy.

Finally, note that betting, and in particular in-game betting, is also an
important component of sports [20, 22]), also in handball. Since in handball
the score changes too often we need a model that can adjust on the score
differential and the remaining time in a flexible way. The Skellam
distribution shares such a property in the sense that changing the time domain
the distribution is still in the same family with altered parameters and hence
it can be very helpful for predicting in-game winning probabilities.

##  5 Competing Interests

The authors have no competing interests to declare that are relevant to the
content of this article.

## References

  *   * Abramowitz and Stegun [1974] Abramowitz, M. and I.A. Stegun. 1974.  Handbook of Mathematical Functions.  New York: Dover.
  * Barbiero [2014] Barbiero, A. 2014.  An alternative discrete skew Laplace distribution.  Statistical Methodology 16: 47–67 .
  * Bryson et al. [2021] Bryson, A., P. Dolton, J.J. Reade, D. Schreyer, and C. Singleton. 2021.  Causal effects of an absent crowd on performances and refereeing decisions during Covid-19.  Economics Letters 198: 109664 .
  * Csató [2021] Csató, L. 2021.  Coronavirus and sports leagues: obtaining a fair ranking when the season cannot resume.  IMA Journal of Management Mathematics 32(4): 547–560 .
  * Dixon and Coles [1997] Dixon, M.J. and S.G. Coles. 1997.  Modelling association football scores and inefficiencies in the football betting market.  Journal of the Royal Statistical Society: Series C (Applied Statistics) 46(2): 265–280 .
  * Dumangane et al. [2009] Dumangane, M., N. Rosati, and A. Volossovitch. 2009.  Departure from independence and stationarity in a handball match.  Journal of Applied Statistics 36(7): 723–741 .
  * Felice [2023] Felice, F. 2023.  Ranking handball teams from statistical strength estimation.  arXiv preprint arXiv:2307.06754 .
  * Felice and Ley [2023] Felice, F. and C. Ley. 2023.  Prediction of handball matches with statistically enhanced learning via estimated team strengths.  arXiv preprint arXiv:2307.11777 .
  * Groll et al. [2020] Groll, A., J. Heiner, G. Schauberger, and J. Uhrmeister. 2020.  Prediction of the 2019 IHF World Men’s Handball Championship–a sparse Gaussian approximation model.  Journal of Sports Analytics 6(3): 187–197 .
  * Holmes and McHale [2024] Holmes, B. and I.G. McHale. 2024.  Forecasting football match results using a player rating based model.  International Journal of Forecasting 40(1): 302–312 .
  * Irwin [1937] Irwin, J.O. 1937.  The frequency distribution of the difference between two independent variates following the same Poisson distribution.  Journal of the Royal Statistical Society 100(3): 415–416 .
  * Jiang et al. [2014] Jiang, L., K. Mao, and R. Wu. 2014.  A Skellam model to identify differential patterns of gene expression induced by environmental signals.  BMC Genomics 15(1): 1–9 .
  * Karlis and Mamode Khan [2023] Karlis, D. and N. Mamode Khan. 2023.  Models for integer data.  Annual Review of Statistics and its Application 10: 297–323 .
  * Karlis and Ntzoufras [2006] Karlis, D. and I. Ntzoufras. 2006.  Bayesian analysis of the differences of count data.  Statistics in Medicine 25(11): 1885–1905 .
  * Koopman et al. [2017] Koopman, S.J., R. Lit, and A. Lucas. 2017.  Intraday stochastic volatility in discrete price changes: the dynamic Skellam model.  Journal of the American Statistical Association 112(520): 1490–1503 .
  * Lago-Peñas et al. [2013] Lago-Peñas, C., M.A. Gomez, J. Viaño, and I. González-García. 2013\.  Home advantage in elite handball: the impact of the quality of opposition on team performance.  International Journal of Performance Analysis in Sport 13(3): 724–733 .
  * Maher [1982] Maher, M.J. 1982.  Modelling association football scores.  Statistica Neerlandica 36(3): 109–118 .
  * McHale and Scarf [2007] McHale, I. and P. Scarf. 2007.  Modelling soccer matches using bivariate discrete distributions with general dependence structure.  Statistica Neerlandica 61(4): 432–445 .
  * Michels et al. [2023] Michels, R., M. Ötting, and D. Karlis. 2023.  Extending the Dixon and Coles model: an application to women’s football data.  arXiv preprint arXiv:2307.02139 .
  * Michels et al. [2023] Michels, R., M. Ötting, and R. Langrock. 2023.  Bettors’ reaction to match dynamics: Evidence from in-game betting.  European Journal of Operational Research 310(3): 1118–1127 .
  * Ntzoufras et al. [2021] Ntzoufras, I., V. Palaskas, and S. Drikos. 2021.  Bayesian models for prediction of the set-difference in volleyball.  IMA Journal of Management Mathematics 32(4): 491–518 .
  * Ötting et al. [2024] Ötting, M., R. Michels, R. Langrock, and C. Deutscher. 2024.  Demand for live betting: An analysis using state-space models.  Applied Stochastic Models in Business and Industry .
  * Pelechrinis and Winston [2021] Pelechrinis, K. and W. Winston. 2021.  A Skellam regression model for quantifying positional value in soccer.  Journal of Quantitative Analysis in Sports 17(3): 187–201 .
  * Prieto et al. [2016] Prieto, J., M.Á. Gómez, A. Volossovitch, and J. Sampaio. 2016.  Effects of team timeouts on the teams’ scoring performance in elite handball close games.  Kinesiology 48(1.): 115–123 .
  * Singh et al. [2023] Singh, A., P. Scarf, and R. Baker. 2023.  A unified theory for bivariate scores in possessive ball-sports: The case of handball.  European Journal of Operational Research 304(3): 1099–1112 .
  * Skellam [1946] Skellam, J. 1946.  The frequency distribution of the difference between two Poisson variates belonging to different populations.  Journal of the Royal Statistical Society. Series A 109(3): 296–296 .
  * Smiatek and Heuer [2012] Smiatek, J. and A. Heuer. 2012.  A statistical view on team handball results: home advantage, team fitness and prediction of match outcomes.  arXiv preprint arXiv:1207.0700 .
  * Tomy and Veena [2022] Tomy, L. and G. Veena. 2022.  A retrospective study on Skellam and related distributions.  Austrian Journal of Statistics 51(1): 102–111 .
  * Van Eetvelde et al. [2023] Van Eetvelde, H., L.M. Hvattum, and C. Ley. 2023.  The probabilistic final standing calculator: a fair stochastic tool to handle abruptly stopped football seasons.  AStA Advances in Statistical Analysis 107(1-2): 251–269 .
  * Volossovitch and Debanne [2021] Volossovitch, A. and T. Debanne. 2021.  Home advantage in handball.  Home Advantage in Sport: Causes and the Effect on Performance: 220–227 .


# Low-Complexity Detection for Faster-than-Nyquist Signaling based on
Probabilistic Data Association

Michel Kulhandjian, Ebrahim Bedeer, Hovannes Kulhandjian, Claude D’Amours, and
Halim Yanikomeroglu M. Kulhandjian and C. D’Amours are with the School of
Electrical Engineering, & Computer Science, University of Ottawa, Ottawa,
Canada, e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>Bedeer is with the
Department of Electrical & Computer Engineering, University of Saskatchewan,
Saskatoon, Canada, e-mail<EMAIL_ADDRESS>Kulhandjian is with the
Department of Electrical & Computer Engineering, California State University,
Fresno, U.S.A., e-mail<EMAIL_ADDRESS>Yanikomeroglu is with
the Department of Systems & Computer Engineering, Carleton University, Ottawa,
Canada, e-mail<EMAIL_ADDRESS>
###### Abstract

In this paper, we investigate the sequence estimation problem of faster-than-
Nyquist (FTN) signaling as a promising approach for increasing spectral
efficiency (SE) in future communication systems. In doing so, we exploit the
concept of Gaussian separability and propose two probabilistic data
association (PDA) algorithms with polynomial time complexity to detect binary
phase-shift keying (BPSK) FTN signaling. Simulation results show that the
proposed PDA algorithm outperforms the recently proposed SSSSE and SSSgb$K$SE
algorithms for all SE values with a modest increase in complexity. The PDA
algorithm approaches the performance of the semidefinite relaxation (SDRSE)
algorithm for SE values of $0.96$ bits/sec/Hz, and it is within the $0.5$ dB
signal-to-noise ratio (SNR) penalty at SE values of $1.10$ bits/sec/Hz for the
fixed values of $\beta=0.3$.

###### Index Terms:

Faster-than-Nyquist (FTN) signaling, intersymbol interference (ISI), Mazo
limit, sequence estimation.

##  I Introduction

Improving spectral efficiency (SE) is one of the key requirements of future
communication systems. Faster-than-Nyquist (FTN) signaling is a promising
physical layer transmission technique for increasing SE of future
communication systems compared to conventional Nyquist signaling [1]. The
basic idea of FTN signaling is to transmit the time domain pulses with a rate
that exceeds the Nyquist limit, which intentionally creates inter-symbol
interference (ISI) at the receiver side. Interestingly, J. Mazo in 1975 showed
that transmitting time domain pulses every $0.802\>T$, where $T$ is the symbol
duration, results in the same minimum Euclidean distance (and, hence, the same
error probability) as conventional Nyquist signaling [2]. This is equivalent
to a $25\%$ increase in SE for the same energy per symbol and transmission
bandwidth. The concept of FTN signaling has been extended to different
domains, e.g., frequency-domain [3], non-binary modulation [4], other pulse
shapes [5], to name a few. Recent research results extended the concept of FTN
signaling to different domains and suggested that it can be a promising
candidates for high capacity peer-to-peer links.

Given that optimal detection of FTN signaling is shown to be NP-hard [6],
several linear/nonlinear detection techniques have been developed to remove
ISI with different degrees of computational complexity [1]. For high spectral
efficiency systems, i.e., severe ISI, M-Bahl, Cocke, Jelinek, and Raviv
(M-BCJR) algorithm [7] and the truncated Viterbi algorithm [8] are examples of
nonlinear detectors that have shown SE gains; however, this has included
exponential complexity in ISI length. It is worth mentioning that iterative
detection of FTN signaling was shown to be effective in approaching the zero-
ISI performance at high signal-to-noise ratio (SNR) values in [1, 9]. To
achieve moderate SE gains with reasonable detection complexity, nonlinear FTN
signaling detection algorithms based on semi-definite relaxation were proposed
in [10] for high-order modulations with polynomial time complexity. A low-
complexity nonlinear frequency domain equalizer was proposed in [11]; however,
it is found that the insertion of the cyclic prefix could result in
performance loss for short block transmissions. The authors in [12] proposed a
linear algorithm that exploits the structure of ISI at the receiver and
detects FTN signaling on a symbol-by-symbol basis. It was shown in [12] that
the proposed algorithm is suitable for low ISI scenarios.

In this paper, we exploit the concept of Gaussian separability to propose
relatively low-complexity probabilistic data association (PDA) algorithms to
detect binary phase-shift keying (BPSK)111The proposed algorithms can be
straightforwardly extended to quadrature-phase-shift keying (QPSK). FTN
signaling at polynomial time complexity. The proposed algorithms are
iterative, i.e., at each iteration a given FTN signaling symbol is detected to
improve the overall error performance. Simulation results show the
effectiveness of the proposed algorithms and their merit over other competing
schemes in the literature. The proposed algorithms outperform the successive
symbol-by-symbol sequence estimation (SSSSE) and successive symbol-by-symbol
with go-back-$K$ sequence estimation (SSSgb$K$SE) algorithms in [12] for all
SE values. Additionally, the proposed algorithms approach the performance of
the semidefinite relaxation (SDRSE) algorithm proposed in [10] for SE values
of $0.96$ bits/sec/Hz, and they are within $0.5$ dB SNR penalty at SE values
of $1.10$ bits/sec/Hz.

The remainder of this paper is organized as follows. Section II presents the
FTN signaling system model and the linear separability concept. The concept of
Gaussian separability and the proposed PDA algorithms are presented in Section
III, while the simulation results are discussed in Section IV. Finally, the
paper is concluded in Section V.

The following notations are used in this paper. All boldface lower case
letters indicate column vectors, and upper case letters indicate matrices.
Additionally, $(\cdot)^{H}$ denotes Hermitian operation, $\mathsf{sgn}(\cdot)$
denotes the sign function, $|\cdot|$ is the scalar magnitude,
$\mathsf{var}[\cdot]$ is variance, $\mathsf{cov}[\cdot]$ denotes covariance,
$\mathsf{diag}[\cdot]$ denotes a matrix with diagonal elements, and
$\mathbb{E}\\{\cdot\\}$ denotes the expected value.

##  II FTN Signaling System Model

At the FTN signaling transmitter, the BPSK data symbols are shaped by a unit-
energy pulse $p(t)$. We transmit a total of $N$ symbols at every $\tau T$,
where $\tau\in(0,1]$ is the time packing parameter. The received signal, after
being affected by an additive white Gaussian noise (AWGN) channel, passes
through a filter matching the transmit pulse and is given as

$\displaystyle{y}(t)=\sqrt{\tau\>E_{s}}\>\sum\nolimits_{n=1}^{N}{a}_{n}g(t-n\tau
T)+{q}(t),$ (1)

where ${a}_{n},\>\>n=1,\ldots,N,$ is the BPSK data symbol, $E_{s}$ is the data
symbol energy, $g(t)=\int\nolimits p(x)p(x-t)dx$, ${q}(t)=\int\nolimits
n(x)p(x-t)dx$, where $n(t)$ is the AWGN with zero mean and variance
$\sigma^{2}$, and $1/(\tau T)$ is the signaling rate. The output signal of the
matched-filter is sampled at every $\tau T$ and can be expressed as

$\displaystyle{y}_{k}$ $\displaystyle=$
$\displaystyle\sqrt{\tau\>E_{s}}\sum\nolimits_{n=1}^{N}{a}_{n}g(k\tau T-n\tau
T)+{q}(k\tau T),$ $\displaystyle=$
$\displaystyle\underbrace{\sqrt{\tau\>E_{s}}\>{a}_{k}\>g(0)}_{\textup{desired
symbol}}$
$\displaystyle+\underbrace{\sqrt{\tau\>E_{s}}\>\sum\nolimits_{n=1,\>n\neq
k}^{N}{a}_{n}\>g((k-n)\tau T)}_{\textup{ISI}}+\>{q}(k\tau T).$

We rewrite (II) in a vector form as

$\displaystyle{\mathbf{y}}_{\textup{c}}$ $\displaystyle=$
$\displaystyle\sqrt{\tau\>E_{s}}\>{\mathbf{a}}\ast\mathbf{g}+{\mathbf{q}}_{\textup{c}},$
(3)

where $\mathbf{a}$, $\mathbf{g}$, $\mathbf{q}_{c}$, and $\ast$ are the
transmit data symbol vector, ISI vector, colored noise vector, and convolution
operator, respectively. Since the noise samples are correlated, we design a
whitening matched filter by applying spectral factorization to the z-transform
of $\mathbf{g}$ to obtain the approximate whitening filter coefficients. For
more details on the design of such filters for FTN receivers, we refer the
reader to [8, 10]. Hence, after passing (3) through the whitening filter, we
have

$\displaystyle{\mathbf{y}}$ $\displaystyle=$
$\displaystyle\sqrt{\tau\>E_{s}}\>{\mathbf{a}}\ast\mathbf{v}+{\mathbf{n}_{\rm{w}}},$
(4)

where ${\mathbf{n}}_{\rm{w}}$ is AWGN $\sim\mathcal{N}(0,\sigma^{2})$ and
$\mathbf{v}$ is the causal ISI vector and is constructed from $\mathbf{g}$ as
$\mathbf{v}[n]\ast\mathbf{v}[-n]=\mathbf{g}$. Equation (4) can be rewritten as

${\mathbf{y}}=\mathbf{G}\mathbf{a}+{\mathbf{n}_{\rm{w}}},$ (5)

where $\mathbf{G}$ is an $N\times N$ Toeplitz Gram matrix that represents the
ISI. The maximum likelihood sequence estimation (MLSE) problem for detecting
the FTN signaling in (5) is formally expressed as

$\displaystyle\widehat{\bm{a}}$ $\displaystyle=$
$\displaystyle\operatornamewithlimits{argmin}_{\mathbf{a}\in\\{\pm
1\\}^{N\times
1}}\mathbf{a}^{H}\mathbf{G}\mathbf{a}-2{\mathbf{y}}^{H}\mathbf{a}.$ (6)

It is known that obtaining the MLSE solution is generally NP-hard [6], which
may be prohibitively complex for use in practical detectors. That being said,
finding low-complexity optimal detectors of FTN signaling is important for
improving SE in future communication systems.

In this work, we aim to identify what properties the ISI matrix $\mathbf{G}$
needs to have such that the sequence estimation can achieve asymptotically
zero probability of error detection. In [12], the authors exploit the idea of
linear separability to identify an operating region where perfect
reconstruction of FTN signaling is guaranteed for noise-free transmission. It
was found that such a linear separability operating region depends on the
raised cosine pulse shape, its roll-off factor $\beta$, and the FTN signaling
time acceleration parameter $\tau$. The perfect reconstruction conditions of
the BPSK FTN signaling in [12] can be summarized in the following lemma.

###### Lemma 1

[12] We assume perfect estimation conditions for BPSK FTN signaling over the
noise-free transmission. Regardless of the value of the current data symbol
$a_{k}$, the upcoming $L-1$ data symbols $a_{k+1},...,a_{k+L-1}$, and the
value of $L$, the following inequality holds for a certain range of $\tau$ and
$\beta$:

$|G_{1,1}a_{k}|>|G_{1,2}a_{k+1}+\ldots+G_{1,L}a_{k+L-1}|,$ (7)

where $G_{i,j}$ is the $i$-th row and $j$-th column value of the ISI matrix
$\mathbf{G}$. Following (7), in this work, we define the conditionally linear
margin as

$\delta_{k}=|G_{1,1}a_{k}|-|G_{1,2}a_{k+1}+\ldots+G_{1,L}a_{k+L-1}|,$ (8)

which can be seen that $\delta_{k}>0$ is equivalent to Lemma 1 and it is the
sufficient and necessary condition for linear separability. Therefore, the
linear margin, $\delta_{k}$, represents a distance measure for the
separability condition of Lemma 1 to hold. Obviously, the greater the linear
margin $\delta_{k}$, the better the error performance.

##  III Proposed Detection Algorithms

In this section, low-complexity iterative detectors are proposed to detect FTN
signaling on the basis of Gaussian separability [13]. We show that Gaussian
separability has a greater classification margin compared to the linear
separability in [12], which thereby improves the detection error performance
of FTN signaling. To illustrate the idea, let us define a subset
$\mathcal{F}\subset\mathcal{U}=\\{1,2,\ldots,N\\}$ of the symbols that are
already detected, i.e., $\\{a_{l}=a_{l}^{*},\>l\in\mathcal{F}\\}$. To detect
the $k$-th symbol $k\in\mathcal{U}-\mathcal{F}$, we should be able to separate
the set of undetected symbols in $\mathcal{U}-\mathcal{F}$ in two groups
corresponding to $a_{k}=1$ and $a_{k}=-1$, respectively. Following this
notation, the observation vector ${\mathbf{y}}$ in (5) is rewritten as

$\displaystyle{\mathbf{y}}$ $\displaystyle=$
$\displaystyle\mathbf{g}_{k}a_{k}+\bm{\mu}_{\mathcal{F}}^{*}+\bm{\mu}_{\mathcal{F}_{k}}+\mathbf{n}_{\rm{w}},$
(9)

where $\mathbf{g}_{k}$ is the $k$-th column of ISI matrix
$\mathbf{G}=[\mathbf{g}_{1},\mathbf{g}_{2},\dots,\mathbf{g}_{N}]$,
$\bm{\mu}_{\mathcal{F}}^{*}=\sum_{l\in\mathcal{F}}\mathbf{g}_{l}a_{l}^{*}$ and
$\bm{\mu}_{\mathcal{F}_{k}}=\sum_{j\in\mathcal{U}-\mathcal{F}-k}\mathbf{g}_{j}a_{j}$.
Linear detectors can be viewed as hyperplanes on the subspace spanned by
$\mathbf{g}_{k}$ for $1\leq k\leq N$ columns of matrix $\mathbf{G}$. A linear
detector that achieves asymptotic negligible error probability as
$\sigma^{2}\rightarrow 0$ is said to be asymptotically efficient.
Mathematically, it can be expressed as, $\forall a_{j}\in\\{\pm 1\\}$,
$j\in\mathcal{U}-\mathcal{F}-k$,

$\displaystyle\mathbf{c}^{H}\mathbf{y}$ $\displaystyle\geq
0\>\>\text{if}\>\>a_{k}=1$ $\displaystyle\mathbf{c}^{H}\mathbf{y}$
$\displaystyle<0\>\>\text{if}\>\>a_{k}=-1.$

From this definition of linear separability, the following theorem can be
proved.

###### Theorem III.1

[14] The $k$-th symbol is asymptotically linearly separable conditionally on
$\mathcal{F}$ if for $\sigma^{2}\rightarrow 0$, an $N$-dimensional linear
filter $\mathbf{c}$,

$|\mathbf{c}^{H}\mathbf{g}_{k}|>\sum_{j\in\mathcal{U}-\mathcal{F}-k}|\mathbf{c}^{H}\mathbf{g}_{j}|.$
(10)

It can be shown that Lemma 1 is a special case of Theorem III.1 when the
linear filter $\mathbf{c}=\mathbf{e}_{k}$ in (10), where $\mathbf{e}_{k}$ is a
vector with value $1$ at position $k$ and $0$ in all other positions under the
constraint that data symbols $a_{k}\in\\{\pm 1\\}$, $\forall k$.

TABLE I: Symbol-by-symbol Separability ($\beta=0.3$)  $\tau$ | 0.6 | 0.7 | 0.8 | 0.9
---|---|---|---|---
$\delta_{max},\delta_{ave}$ | $0.97,0.11$ | $0.97,0.20$ | $0.97,0.38$ | $0.97,0.53$
TABLE II: Gaussian Separability ($\beta=0.3$)  SNR (dB) | $\tau$ = 0.6 | 0.7 | 0.8 | 0.9
---|---|---|---|---
$0$ | $0.00,0.00$ | $1.16,0.49$ | $1.08,0.58$ | $1.02,0.73$
2 | $2.01,0.68$ | $1.84,0.82$ | $1.71,0.96$ | $1.62,1.21$
4 | $3.19,1.12$ | $2.92,1.35$ | $2.71,1.57$ | $2.56,2.00$
6 | $5.06,1.73$ | $4.62,2.16$ | $4.29,2.55$ | $4.06,3.30$
8 | $8.02,2.71$ | $7.32,3.45$ | $6.80,4.13$ | $6.43,5.42$

Recall that if the condition densities for all the symbols
$\mathbf{f}(\mathbf{y}|a_{j})$ are Gaussian then the optimal filter for
user-$k$ is a scaled version of e.g.,
$\mathbf{c}_{k}\triangleq\mathbf{R}_{\mathcal{F}_{k}}^{-1}\mathbf{g}_{k}$. The
solution of the optimal separating hyperplane in a linear detector
(classifier), despite the simple formulation, may be very difficult when the
classes do not follow a Gaussian distribution. If we assume that the ISI is
Gaussian, then the linear classifier becomes the optimal maximum likelihood
(ML) decision with hyperplanes dividing the $N$-dimensional space.

###### Theorem III.2

[14] The $k$-th symbol is a Gaussian separable conditionally on $\mathcal{F}$
if for $\sigma^{2}\rightarrow 0$,

$\mathbf{g}_{k}^{H}\mathbf{R}_{\mathcal{F}_{k}}^{-1}\mathbf{g}_{k}>\sum_{j\in\mathcal{U}-\mathcal{F}-k}|\mathbf{g}_{k}^{H}\mathbf{R}_{\mathcal{F}_{k}}^{-1}\mathbf{g}_{j}|,$
(11)

where

$\mathbf{R}_{\mathcal{F}_{k}}=\sum_{j\in\mathcal{U}-\mathcal{F}-k}\mathbf{g}_{j}\mathbf{g}_{j}^{H}+\sigma^{2}\mathbf{I}_{N},$
(12)

and the parameter

$\Delta_{\mathcal{F}_{k}}=\mathbf{g}_{k}^{H}\mathbf{R}_{\mathcal{F}_{k}}^{-1}\mathbf{g}_{k}-\sum_{j\in\mathcal{U}-\mathcal{F}-k}|\mathbf{g}_{k}^{H}\mathbf{R}_{\mathcal{F}_{k}}^{-1}\mathbf{g}_{j}|$
(13)

is called conditionally Gaussian margin. Unlike the linear margin
$\delta_{k}$, the Gaussian margin $\Delta_{\mathcal{F}_{k}}$ value depends on
the noise variance $\sigma^{2}$. We measured conditionally margins,
$\delta_{k}$ and $\Delta_{\mathcal{F}_{k}}$, in the operating region, where
ISI matrix $\mathbf{G}$ is formed with roll-off factor $\beta=0.3$ and
$\tau\in\\{0.6,0.7,0.8,0.9\\}$. At each pair of $\beta,\tau$, we generate
matrix $\mathbf{G}$ and then compute $\delta_{k}$ and
$\Delta_{\mathcal{F}_{k}}$ for $1\leq k\leq N$ using equations (8) and
(13)222We note that $\Delta_{\mathcal{F}_{k}}$ is computed using the following
set $\mathcal{F}=\\{i_{1},i_{2},\dots,i_{k-1}\\}$ for $1\leq k\leq N$. The
order of $\\{i_{1},i_{2},\dots,i_{k-1}\\}$ depends on the ISI matrix
$\mathbf{G}$ and $i_{t}$’s are precisely computed in our PDA algorithm at line
$6$., respectively. We evaluate their maximum and average values as
$\delta_{max}=\operatornamewithlimits{argmax}_{k}\delta_{k}$,
$\delta_{ave}=\frac{1}{N}\sum_{k=1}^{N}\delta_{k}$,
$\Delta_{max}=\operatornamewithlimits{argmax}_{k}\Delta_{\mathcal{F}_{k}}$,
$\Delta_{ave}=\frac{1}{N}\sum_{k=1}^{N}\Delta_{\mathcal{F}_{k}}$.

Since the linear margin $\delta_{k}$ does not depend on the noise variance,
$\sigma^{2}$, there is only one row with values ($\delta_{max}$,
$\delta_{ave}$) as shown in Table I. On the other hand, Table II shows the
($\Delta_{max}$, $\Delta_{ave}$) for different SNR values at each row, where
we vary the noise variance, $\sigma^{2}$, and keep the matrix $\mathbf{G}$
unchanged. We observe from Table I and Table II that the Gaussian
classification margins have greater maximum and average values than their
linear counterparts. Therefore, we can theorize that using a Gaussian
classification margin can potentially improve the performance in terms of bit
error rate (BER).

The main idea of the PDA detection algorithm is based on conditionally
Gaussian separability criterion and the fact that the Gaussian margin
$\Delta_{\mathcal{F}_{k}}$ is greater than the linear margin $\delta_{k}$, as
shown in Tables I and II. One can re-express (9) as

$\displaystyle{\mathbf{y}}$ $\displaystyle=$
$\displaystyle{\mathbf{g}_{k}a_{k}+\mathbf{w}_{k}},$ (14)

where $\mathbf{w}_{k}=\mathbf{G}_{k}\mathbf{a}_{k}+\mathbf{n}_{\rm{w}}$,
$\mathbf{G}_{k}$ is the $N\times N-1$ resulting matrix after excluding the
$k$-th column of $\mathbf{G}$, and $\mathbf{a}_{k}$ is $N-1\times 1$ vector
resulting after excluding $k$-th element of $\mathbf{a}$.

We associate the probabilities $P_{a}(k)$ and $1-P_{a}(k)$ with the estimate
of the symbols $a_{k}=1$ and $a_{k}=-1$, respectively. Now, for an arbitrary
symbol $a_{k}$, we treat the other symbols $a_{j}$, $\forall j\neq k$ as
binary random variables and $\mathbf{w}_{k}$ as a Gaussian noise vector. Then
the likelihood ratio for the $a_{k}$ symbol can be written as

$\displaystyle\Lambda(a_{k})$ $\displaystyle=$
$\displaystyle\frac{P\\{a_{k}=1|\mathbf{y},\\{P_{a}(j)\\}_{j\neq
k}\\}}{P\\{a_{k}=-1|\mathbf{y},\\{P_{a}(j)\\}_{j\neq k}\\}},$ $\displaystyle=$
$\displaystyle\frac{f_{\mathbf{y}}(\mathbf{y}|a_{k}=1|,\\{P_{a}(j)\\}_{j\neq
k})}{f_{\mathbf{y}}(\mathbf{y}|a_{k}=-1,\\{P_{a}(j)\\}_{j\neq k})}$
$\displaystyle=$
$\displaystyle\frac{\mathcal{N}(\mathbf{y}|\mathbf{g}_{k}+\bm{\mu}_{k},\mathbf{C}_{k})}{\mathcal{N}(\mathbf{y}|\mathbf{g}_{k}-\bm{\mu}_{k},\mathbf{C}_{k})}$
$\displaystyle=$
$\displaystyle\exp\left[2(\bm{y}-\bm{\mu}_{k})^{H}\mathbf{C}_{k}^{-1}\mathbf{g}_{k}\right],$
(16)

where $P_{a}(k)=P\\{a_{k}=1|\mathbf{y},\\{P_{a}(j)\\}_{j\neq k}\\}$ and
$f_{\mathbf{y}}$ is the multivariate Gaussian distribution conditioned on
$\\{\mathbf{y},\\{P_{a}(j)\\}_{j\neq k}\\}$. The covariance can be calculated
as

$\displaystyle\mathbf{C}_{k}$ $\displaystyle=$
$\displaystyle\mathsf{cov}[\mathbf{w}_{k}]$ (17) $\displaystyle=$
$\displaystyle\mathbf{G}_{k}\mathbb{E}[\mathbf{a}_{k}\mathbf{a}_{k}^{H}]\mathbf{G}_{k}^{H}-\mathbf{G}_{k}\mathbb{E}[\mathbf{a}_{k}]\mathbb{E}[\mathbf{a}_{k}^{H}]\mathbf{G}_{k}^{H}+\mathbf{R}_{n}$
$\displaystyle=$
$\displaystyle\mathbf{G}_{k}\mathbf{C}_{\mathbf{a}_{k}}\mathbf{G}_{k}^{H}+\mathbf{R}_{n},$

where $\mathsf{var}[a_{k}]=4P_{a}(k)(1-P_{a}(k))$,
$\mathbf{C}_{\mathbf{a}_{k}}=\mathsf{cov}[\mathbf{a}_{k}]=\mathsf{diag}(\mathsf{var}[a_{1}],\cdots,\mathsf{var}[a_{k-1}],\mathsf{var}[a_{k+1}],\cdots,\mathsf{var}[a_{n}])$
and $\mathbf{R}_{n}=\mathbb{E}[\mathbf{n}_{\rm{w}}\mathbf{n}_{\rm{w}}^{H}]$,
are the covariance matrix of $\mathbf{a}_{k}$ and the autocorrelation of noise
vector, respectively. The knowledge about the $a_{k}$ symbol is carried over
in each iteration through a dynamic update of means and covariances as follows

$\displaystyle\bm{\mu}_{k}$ $\displaystyle=$
$\displaystyle\mathbb{E}[\mathbf{w}_{k}]=\mathbb{E}[\mathbf{G}_{k}\mathbf{a}_{k}+\mathbf{n}_{\rm{w}}]=\mathbf{G}_{k}\mathbb{E}[\mathbf{a}_{k}]$
(18) $\displaystyle=$
$\displaystyle\mathbf{G}_{k}\begin{bmatrix}2P_{a}(1)-1\\\ \vdots\\\
2P_{a}(k-1)-1\\\ 2P_{a}(k+1)-1\\\ \vdots\\\ 2P_{a}(N)-1\end{bmatrix}$
$\displaystyle=$ $\displaystyle\mathbf{G}_{k}\mathbf{p}_{a,k},$

where $\mathbf{p}_{a,k}$ is the estimated $\mathbf{a}$ signal vector without
the $k$-th signal. As we can see, the likelihood ratio in (III) is in the form
of $P_{a}(k)/(1-P_{a}(k))=\exp(x)$. Hence, the posteriori probability
$P_{a}(k)$ can be calculated as

$\displaystyle
P_{a}(k)=1/\left(1+\exp\left[-2(\mathbf{y}-\bm{\mu}_{k})^{H}\mathbf{C}_{k}^{-1}\mathbf{g}_{k}\right]\right).$
(19)

That being said, the proposed PDA algorithm for detecting FTN signaling can be
summarized in the next page, where $M$ is the number of iterations and
$D_{k}=\mathbf{g}_{k}^{H}\mathbf{C}_{k}^{-1}\mathbf{g}_{k}$.

Proposed PDA Algorithm
---
Input: $\mathbf{y}$, $\mathcal{F}\leftarrow\\{1,2,\dots,N\\}$, $M$
1: Initialize : $P_{a}(k)=1/2$, $k=1,\dots,N$
2: Compute $D_{k}$, $k=1,\cdots,N$
3: for $i=1$ to $M$
4: $\mathcal{U}\leftarrow\mathcal{F}$
5: while ($\mathcal{U}\neq\emptyset$)
6: Update $P_{a}(k_{m})$,
$k_{m}=\operatornamewithlimits{argmax}_{k\in\mathcal{U}}(D_{k})$
7: Subtract $k_{m}$ from $\mathcal{U}$
8: Recompute $D_{k}$, $\forall k\in\mathcal{U}$
Output: $\mathbf{\hat{a}}=\mathsf{sgn}\left(\mathbf{p}_{a}\right)$

###  III-A Modified PDA algorithm

The complexity of the PDA algorithm can be further reduced if the detected
data symbols whose probabilities lie within a confidence interval are removed
in each iteration. In other words, if the probability of the symbol $k$,
$P_{a}(k)$ lies in confidence region $\epsilon$, if $|P_{a}(k)-d|<\epsilon$,
where $d=\\{0,1\\}$. All the symbols that satisfy the confidence interval
define the set $\mathcal{C}$. We assign
$P_{a}(k)=\\{d:|P_{a}(k)-d|<\epsilon\\}$ for $\forall k\in\mathcal{C}$. Next,
we remove symbols $\mathcal{C}$ from $\mathcal{F}$,
$\mathcal{F}=\mathcal{F}-\mathcal{C}$. The modified PDA algorithm is presented
next.

Modified PDA Algorithm
---
Input: $\mathbf{y}=\mathbf{G}\mathbf{a}+\mathbf{n}$,
$\mathcal{F}\leftarrow\\{1,2,\dots,N\\}$, $M$
1: Initialize $P_{a}(k)=1/2$, $k=1,\cdots,N$
2: Compute $D_{k}$, $k=1,\dots,N$
3: for $i=1$ to $M$
4: Find $k$ that satisfy $|P_{a}(k)-d|<\epsilon$
5: if $\mathcal{C}\neq\emptyset$
6: $P_{a}(k)=d$, $\forall k\in\mathcal{C}$
7: $\mathcal{F}\leftarrow\mathcal{F}-\mathcal{C}$
8: $\mathcal{U}\leftarrow\mathcal{F}$
9: while ($\mathcal{U}\neq\emptyset$)
10: Update $P_{a}(k_{m})$,
$k_{m}=\operatornamewithlimits{argmax}_{k\in\mathcal{U}}(D_{k})$
11: Subtract $k_{m}$ from $\mathcal{U}$
12: Recompute $D_{k}$, $\forall k\in\mathcal{U}$
Output: $\mathbf{\hat{a}}=\mathsf{sgn}\left(\mathbf{p}_{a}\right)$

###  III-B Complexity analysis

The main complexity of the proposed PDA algorithm lies in updating
$\mathbf{C}_{k}^{-1}$ in (19). However, efficient numerical techniques such as
the Durbin-Levinson algorithm can exploit the Toeplitz structure of $C_{k}$
and compute its inverse, i.e., $\mathbf{C}_{k}^{-1}$, with a complexity of
$\mathcal{O}(N^{2})$ [15]. However, since $C_{k}$ is not always guaranteed to
be Toeplitz, the calculations of $C_{k}^{-1}$ can be in the order of
$\mathcal{O}(N^{3})$. Since evaluating (19) is expected to occur at most $N$
times and the algorithm has deterministic execution steps, the total average
and worst case computational complexity of the proposed PDA algorithm can be
in the order of $\mathcal{O}(N^{4})$. The PDA algorithm has the same
polynomial complexity order (i.e., both complexities scale in polynomial time
with the number of transmit symbols) as the SDRSE algorithm in [10] (please
note that complexity measure based on the number of operations depends mainly
on hardware implementations and is outside the scope of this work), while it
has higher complexity compared to the work in [12].

##  IV Simulation Results

In this section, we assess the performance of our proposed algorithms for
detecting BPSK FTN signaling with a root-raised cosine (rRC) pulse and roll-
off factor of $\beta=0.3$. Simulations were performed over the AWGN channel,
both with and without channel encoding. The iteration number for proposed PDA
is $8$ for all cases and the number of iterations we averaged the simulations
over are $50,000$ and $5,000$ for uncoded and coded cases, respectively.

Figure 1: BER performance of BPSK FTN signaling as a function of
$\frac{E_{b}}{N_{o}}$ using the proposed PDA algorithm, SSSSE [12], SSSgb$K$SE
[12], and the SDRSE [10], at $\beta=0.3$ and $\tau=0.8$, SE $=0.96$
bits/sec/Hz.

Figure 1 depicts the BER of BPSK FTN signaling as a function of
$\frac{E_{b}}{N_{o}}$ for the proposed PDA algorithm, comparing it with the
SSSE and SSSgb$K$SE ($K=3$) algorithms presented in [12], the SDRSE algorithm
in [10], and the frequency domain equalizer (FDE) in [11] for $\beta=0.3$ and
$\tau=0.8$, i.e., a SE of $0.96$ bits/sec/Hz. The SE is calculated as
$\frac{\log_{2}{M}}{(1+\beta)\tau}$, where $M$ is the constellation size and
the SE loss due to adding the cyclic prefix in the FDE in [11] is not
considered. As we can see, the proposed PDA algorithm approaches the
performance of the SDRSE algorithm in [10] and outperforms the works in [12,
11] for SE value of $0.96$ bits/sec/Hz.

Figure 2: BER performance of BPSK FTN signaling as a function of
$\frac{E_{b}}{N_{o}}$ using the proposed PDA algorithm, SSSSE [12], SSSgb$K$SE
[12], and the SDRSE [10], at $\beta=0.3$ and $\tau=0.7$, SE $=1.10$
bits/sec/Hz. Figure 3: BER performance for confidence region.

Figure 2 plots the BER of the BPSK FTN signaling as a function of
$\frac{E_{b}}{N_{o}}$ for the proposed PDA algorithm, the SSSSE and SSSgb$K$SE
($K=3$) algorithms in [12], and the SDRSE algorithm in [10] for $\beta=0.3$
and $\tau=0.7$, i.e., a SE of $1.10$ bits/sec/Hz. As we can see in Fig. 2, the
proposed PDA algorithm’s performance is within almost $0.5$ dB of the SDRSE
performance in [10], and it extensively outperforms the works in [12] that
fail at such value of SE. We observe in Fig. 3 of our simulation for the
confidence region $\epsilon=0.4$, the average number of iteration steps
reduces by around $30\%$. At higher SNR values, we observe higher percentage
of reduction in the iteration steps.

We also performed the simulation using a low-density parity-check (LDPC) and
turbo codes. For implementation, we used a custom parity based matrix
generator for LDPC and the long-term evolution (LTE) turbo channel codes
described in [16]. For both cases, we used $1/3$ code rate channel coding with
an input message length of $320$ and an output code length of $972$ for both
LDPC and LTE turbo, respectively. The construction of the LTE interleaver was
based on the quadratic permutation polynomial (QPP) scheme [16]. Figure 4
shows the BER performance for LDPC and turbo encoded with $\tau=0.7$ and
$\tau=0.8$. The BER performance of turbo coded is better compared to LDPC in
both cases. The SSSE and SSSgb$K$SE are not included in the simulation results
since they do not produce soft outputs required by the LDPC and turbo soft
decoders.

Figure 4: BER performance for LDPC and turbo encoding.

##  V Conclusion

Faster-than-Nyquist (FTN) signaling is a promising candidate for improving
spectral efficiency (SE) in future communication systems. In this paper, we
showed that Gaussian separability has a greater margin compared to linear
separability, which has been used recently in the detection of FTN. In
exploiting Gaussian separability, we proposed two probabilistic data
association (PDA) algorithms of polynomial time computational complexity to
detect binary phase-shift keying (BPSK) FTN signaling. Simulation results
showed that the proposed algorithms outperform the successive symbol-by-symbol
sequence estimation (SSSSE) and successive symbol-by-symbol with go-back-$K$
sequence estimation (SSSgb$K$SE) algorithms in [12] for all SE values.
Additionally, they approach the performance the semidefinite relaxation
(SDRSE) algorithm in [10] for SE values of $0.96$ bits/sec/Hz and they are
within the $0.5$ dB signal-to-noise ratio (SNR) penalty at SE values of $1.10$
bits/sec/Hz.

## References

  * [1] J. B. Anderson, F. Rusek, and V. Öwall, “Faster-than-Nyquist signaling,” _Proc. IEEE_ , vol. 101, no. 8, pp. 1817–1830, Aug. 2013\.
  * [2] J. Mazo, “Faster-than-Nyquist signaling,” _Bell Syst. Tech. J._ , vol. 54, no. 8, pp. 1451–1462, Oct. 1975.
  * [3] F. Rusek and J. B. Anderson, “The two dimensional Mazo limit,” in _Proc. IEEE Int. Symp. Inf. Theory (ISIT)_ , Adelaide, Australia, Sep. 2005, pp. 970–974.
  * [4] ——, “Non binary and precoded faster than Nyquist signaling,” _IEEE Trans. Commun._ , vol. 56, no. 5, pp. 808–817, May 2008.
  * [5] A. D. Liveris and C. N. Georghiades, “Exploiting faster-than-Nyquist signaling,” _IEEE Trans. Commun._ , vol. 51, no. 9, pp. 1502–1511, Sep. 2003.
  * [6] R. Lupas and S. Verdu, “Linear multiuser detectors for synchronous code-division multiple-access channels,” _IEEE Trans. Inf. Theory_ , vol. 35, no. 1, pp. 123–136, Jan. 1989.
  * [7] J. B. Anderson, A. Prlja, and F. Rusek, “New reduced state space BCJR algorithms for the ISI channel,” in _Proc. IEEE Int. Symp. Inf. Theory (ISIT)_ , Seoul, Korea, Jun. 2009, pp. 889–893.
  * [8] A. Prlja, J. B. Anderson, and F. Rusek, “Receivers for faster-than-Nyquist signaling with and without turbo equalization,” in _Proc. IEEE Int. Symp. Inf. Theory (ISIT)_ , Toronto, Canada, Jul. 2008, pp. 464–468.
  * [9] A. Barbieri, D. Fertonani, and G. Colavolpe, “Time-frequency packing for linear modulations: spectral efficiency and practical detection schemes,” _IEEE Trans. Commun._ , vol. 57, no. 10, pp. 2951–2959, Oct. 2009\.
  * [10] E. Bedeer, M. H. Ahmed, and H. Yanikomeroglu, “Low-complexity detection of high-order QAM faster-than-Nyquist signaling,” _IEEE Access_ , vol. 5, pp. 14 579–14 588, Jul. 2017.
  * [11] S. Sugiura, “Frequency-domain equalization of faster-than-Nyquist signaling,” _IEEE Wireless Commun. Lett._ , vol. 2, no. 5, pp. 555–558, Oct. 2013.
  * [12] E. Bedeer, M. H. Ahmed, and H. Yanikomeroglu, “A very low complexity successive symbol-by-symbol sequence estimator for faster-than-Nyquist signaling,” _IEEE Access_ , vol. 5, no. 1, pp. 2169–3536, Mar. 2017.
  * [13] J. Luo, K. R. Pattipati, and P. K. Willett, “Near-optimal multiuser detection in synchronous CDMA using probabilistic data association,” _IEEE Commun. Lett._ , vol. 5, no. 9, pp. 361–366, Sept. 2001.
  * [14] G. Romano, F. Palmieri, P. K. Willett, and D. Mattera, “Separability and gain control for overloaded CDMA,” in _Proc. of the Conf. on Inf. Sci. and Sys. (CISS)_ , Princeton, NJ, USA, Mar. 2005, pp. 1–6.
  * [15] T. Wysocki, H. Razavi, and B. Honary, _Digital signal processing for communication systems_. Springer Science & Business Media, 2013.
  * [16] H. Zarrinkoub, _Understanding LTE with MATLAB: From Mathematical Modeling to Simulation and Prototyping_ , 1st ed. Wiley Publishing, 2014.