# Boole’s probability bounding problem, linear programming aggregations, and
nonnegative quadratic pseudo-Boolean functions

Endre Boros and Joonhee Lee
RUTCOR and MSIS Department
Rutgers Business School
Rutgers University

###### Abstract

George Boole first proposed the union bounding problem, which is a class of
probabilistic satisfiability problems. Based on Hailperin’s linear programming
formulation of the problem, a very large body of research was conducted to
provide good bounds by considering a relaxation that is either polynomial
dimension in $n$ or is at least polynomially computable. We study monotone
linear mappings $\mathcal{L}$ to a lower $N$ dimensional space such that
$\mathcal{L}(\mathbb{R}_{+}^{2^{n}-1})\subseteq\mathbb{R}_{+}^{N}$ and provide
a complete characterization of the cone
$L_{\mathcal{E}}(\mathbb{R}_{+}^{\Omega})$ for three aggregation models from
the literature. This characterization leads to strictly improved bounds in two
of these cases as a result. We also obtain an $O(n^{2})$ dimensional linear
program that is equivalent with Hailperin’s model. An interesting side result
is that there exists a connection between Hailperin’s model and the cut-
polytope.

##  1 Boole’s Problem and our results

In his seminal works, George Boole [1, 2] (see also [3]) introduced a
difficult family of problems that turned out to be very important in numerous
applications and became widely studied. In its simplest and most frequently
studied form, we are given reals $p_{i}$, $i=1,...,n$, and $p_{i,j}$ for
$1\leq i<j\leq n$, for some positive integer $n$, and the problems is to
determine tight lower and upper bounds for the probability of the union of $n$
events, $A_{i}$, $i=1,..,n$, assuming the equalities $\mathbb{P}(A_{i})=p_{i}$
for $i=1,...,n$ and $\mathbb{P}(A_{i}\cap A_{j})=p_{i,j}$ for all $1\leq
i<j\leq n$. This problem, known as the _union bounding problem_ , turns out to
be quite difficult.

Hailperin [17] proposed a linear programming formulation. Note that the events

$\omega_{S}~{}=~{}\bigcap_{i\in S}A_{i}\cap\bigcap_{i\not\in S}\bar{A}_{i}$

for $S\subseteq\\{1,...,n\\}$ partition the probability space, and the events
$A_{i}$ or $A_{i}\cap A_{j}$ or the union of all these events can be viewed as
unions of some of these _atomic_ events. Hailperin proposed to view

$x_{S}~{}=~{}\mathbb{P}(\omega_{S})~{}~{}~{}\text{ for all
}~{}~{}~{}S\subseteq\\{1,...,n\\}$

as unknown variables, and formulated the union bounding problem as the linear
optimization problem of the form

$\begin{array}[]{rrcll}\max/\min&\displaystyle\sum_{S\subseteq\\{1,...,n\\}\atop
S\neq\emptyset}x_{S}&&&\\\ s.t.&&&&\\\
&\displaystyle\sum_{S\subseteq\\{1,...,n\\}\atop S\supseteq
Q}x_{S}&=&p_{Q}&~{}~{}\forall~{}~{}\begin{array}[]{c}Q\subseteq\\{1,...,n\\},\\\
|Q|\leq 2,\end{array}\\\\[19.91692pt]
&x_{S}&\geq&0&~{}~{}\forall~{}~{}S\subseteq\\{1,...,n\\}.\end{array}$ (1)

Let us denote by $LB(H)$ and $UB(H)$ the optimum values of the minimization
and maximization, respectively, of above linear program.

Let us add here a few remarks, to enhance clarity of our following
discussions:

  * (R1)

Note that $LB(H)$ and $UB(H)$ are _sharp bounds_ , meaning that for an
arbitrary real $Z$ with $LB(H)\leq Z\leq UB(H)$ there exists a probability
space with events $A_{i}$, $i=1,...,n$ such that
$Z=\mathbb{P}\left(\displaystyle\bigcup_{i=1}^{n}A_{i}\right)$.

  * (R2)

From a computational point of view, the input size depends only on $n$. Thus,
for a reasonably bounded rational input the input size is polynomially bounded
by $n$, while the number of variables in (1) is exponential in $n$. This in
itself of course is not an indication of computational hardness. However, it
was shown that trying to use column generation leads, in every iteration, to a
binary unconstrained quadratic optimization, that in general is known to be
NP-hard [19]. It was also shown that the feasibility of (1) belongs to the
family of probabilistic satisfiablity problems, that are also known to be NP-
hard, in general [20]. These indications are not a definitive proof of
computational hardness. Nevertheless, the perceived hardness of (1) motivated
a large amount of research to provide bounds that are close to $LB(H)$ from
below and to $UB(H)$ from above and are computable in polynomial time in terms
of the input size.

  * (R3)

We can observe that variable $x_{\emptyset}$ appears only once, in the
equality corresponding to $p_{\emptyset}=1$. Thus, we can eliminate both
$x_{\emptyset}$ and the first equality from our formulation yielding a
somewhat simplified formulation:

$\begin{array}[]{rrcll}\max/\min&\displaystyle\sum_{S\subseteq\\{1,...,n\\}\atop
S\neq\emptyset}x_{S}&&&\\\ s.t.&&&&\\\
&\displaystyle\sum_{S\subseteq\\{1,...,n\\}\atop S\supseteq
Q}x_{S}&=&p_{Q}&~{}~{}\forall~{}~{}\begin{array}[]{c}Q\subseteq\\{1,...,n\\},\\\
1\leq|Q|\leq 2,\end{array}\\\\[19.91692pt]
&x_{S}&\geq&0&~{}~{}\forall~{}~{}\emptyset\neq
S\subseteq\\{1,...,n\\}.\end{array}$ (2)

This simplification may lead to a feasible formulation (2), even when problem
(1) is infeasible [8]. Still, denoting by $Z^{\max}$ and $Z^{\min}$ the
optimum values in (2), we have the relations that (1) is infeasible if and
only if $Z^{\min}>1$, and for feasible inputs we have $LB(H)=Z^{\min}$ and
$UB(H)=\min\\{1,Z^{\max}\\}$. Thus, problems (1) and (2) are equivalent.

There is a very large body of research on providing good bounds for (2),
mostly by considering a relaxation that is either polynomial dimension in $n$
or is at least polynomially computable, due to some structural property. We
refer the reader to [8] for a concise summary of these results. Many of these
results are based on relaxations of (2) obtained by aggregating some of the
variables of the original formulation (and sometimes, aggregating also some of
the equalities of (2)). In fact such aggregations are monotone linear mappings
$\mathcal{L}$ to a lower $N$ dimensional space such that
$\mathcal{L}(\mathbb{R}_{+}^{2^{n}-1})\subseteq\mathbb{R}_{+}^{N}$. A commonly
overlooked fact is that one may have strict containment in this relation. In
our paper we study such mappings, and provide a complete polyhedral
characterization of $\mathcal{L}(\mathbb{R}_{+}^{2^{n}-1})$ for a large family
of aggregations. Our results allow us to provide new polynomially computable
probability bounds that strictly improve on existing results of e.g., [27,
28]. We also obtain an $O(n^{2})$ dimensional linear program that is
equivalent with (2). An interesting side result is that there exists a
somewhat surprising connection between problem (2) and the cut-polytope (a
very different connection from the one pointed out by [15]).

##  2 Notation and terminology

For a positive integer $n$ let us define $V=[n]=\\{1,2,...,n\\}$, and define
$\Omega=2^{V}\setminus\\{\emptyset\\}$.

In this paper we focus on aggregations of the variables of (2) of the
following type: To a hypergraph $\mathcal{E}=\\{E_{1},\ldots,E_{N}\\}\subseteq
2^{\Omega}$ we associate

$y^{k}~{}=~{}\sum_{S\in E_{k}}x_{S}~{}~{}~{}\text{ for }k=1,...,N.$ (3)

This defines a monotone linear mapping
$L_{\mathcal{E}}:\mathbb{R}^{\Omega}\mapsto\mathbb{R}^{N}$. In each of the
following special cases, we choose $\mathcal{E}$ in such a way that the
objective function and the equality constraints in (2) are easily
representable (sometimes after applying additional aggregations of some of the
equalities).

From the above definition it is also immediate that
$L_{\mathcal{E}}(\mathbb{R}_{+}^{\Omega})\subseteq\mathbb{R}_{+}^{N}$. A
rarely recognized property in the literature is that many times we have
$L_{\mathcal{E}}(\mathbb{R}_{+}^{\Omega})\neq\mathbb{R}_{+}^{N}$. Thus, the
obtained relaxation could frequently be tightened with additional inequalities
corresponding to the facets of the cone
$L_{\mathcal{E}}(\mathbb{R}_{+}^{\Omega})$.

We consider three aggregation models from the literature and provide a
complete characterization of the cone
$L_{\mathcal{E}}(\mathbb{R}_{+}^{\Omega})$. This characterization leads to
strictly improved bounds in two of these cases.

In the subsequent sections variables will be associated to various subfamilies
of $\Omega$. One of the dividing parameters used many times is the cardinality
of the subsets. Accordingly, we introduce

$\Omega^{\ell}~{}=~{}\\{S\in\Omega\mid|S|=\ell\\}$

for $\ell=1,...,n$, and therefore we have the partition
$\Omega~{}=~{}\bigcup_{\ell=1}^{n}\Omega^{\ell}$.

Let us start first recalling some of the aggregations of this type from the
literature, and see if one can improve on them by describing the image of the
positive orthant.

##  3 The Binomial Moment Problem

Let us first consider the aggregation with the family
$\mathcal{E}^{0}=\\{E^{0}_{k}\mid k=1,...,n\\}$, where
$E^{0}_{k}=\\{S\subseteq V\mid|S|=k\\}$ for $k=1,...,n$ resulting in the so
called _binomial moment problem_ :

$\max\left\\{\sum_{k=1}^{n}y^{k}\left|\begin{array}[]{rcl}\displaystyle\sum_{k=1}^{n}ky^{k}&=&\displaystyle\sum_{i\in
V}p_{i}\\\
\displaystyle\sum_{k=1}^{n}\binom{k}{2}y^{k}&=&\displaystyle\sum_{i,j\in
V\atop i<j}p_{i,j}\\\ y^{k}&\geq&0~{}~{}\text{ for
}~{}~{}k=1,...,n.\end{array}\right.\right\\}$

Numerous papers dealt with this problem ([7], [9], [10], [11], [16], [18],
[21], [23], [24], [25], and [26]). As we noted above the nonnegativity of the
$x_{S}$, $S\in\Omega$ variables simply replaced by the nonnegativity of the
aggregated $y^{k}$, $k=1,...,n$ variables. It turns out that in this case this
is correct:

###### Theorem 1

For this aggregation we have

$L_{\mathcal{E}^{0}}\left(\mathbb{R}_{+}^{\Omega}\right)~{}=~{}\mathbb{R}_{+}^{n}.$

Proof: Assume $\alpha_{0},...,\alpha_{1}$ are real numbers such that

$\alpha_{0}+\sum_{k=1}^{n}\alpha_{k}y^{k}~{}\geq~{}0$

whenever $\mathbf{x}=(x_{S}\mid S\in\Omega)\in\mathbb{R}_{+}^{\Omega}$.

Since $\mathbf{x}=0$ is such a nonnegative vector, for which $\mathbf{y}=0$,
and thus $\alpha_{0}\geq 0$ follows immediately.

Furthermore we can consider the vectors
$\mathbf{x}\in\mathbb{R}_{+}^{\Omega}$, in which $x_{S}=0$ for all
$S\in\Omega$, $|S|\neq k$. For this case we have $y^{j}=0$ for all $j\neq k$,
and $y^{k}$ can take arbitrary positive values. Since we must have
$\alpha_{0}+\alpha_{k}\cdot y^{k}\geq 0$ for all these substitutions, we must
have $\alpha_{k}\geq 0$, too. Since all coefficients in the above inequality
must be nonnegative, it is a trivial consequence of the nonnegativity of the
$y^{k}$, $k=1,...,n$ variables. $\Box$

##  4 The Aggregation by Prékopa and Gao

Prékopa and Gao [27] introduced an aggregation of (2) that corresponds to the
case when

$\mathcal{E}^{1}=\\{E^{k}_{i}\mid i\in V,~{}k=1,...,n\\},$

where

$E^{k}_{i}=\\{S\in\Omega^{k}\mid i\in S\\}.$

Thus, $L_{\mathcal{E}^{1}}$ is a mapping from $\mathbb{R}_{+}^{\Omega}$ to the
space of the $n^{2}$ variables defined by

$y^{k}_{i}~{}=~{}\sum_{S\in\Omega^{k}\atop i\in S}x_{S}~{}~{}~{}~{}\text{ for
all }i\in V\text{ and }k=1,...,n.$ (4)

Using these definitions, they considered the aggregation

$\begin{array}[]{crcll}\max~{}(\min)&\displaystyle\sum_{k=1}^{n}\frac{1}{k}\sum_{i\in
V}y^{k}_{i}\\\\[14.22636pt] s.t.\\\
&\displaystyle\sum_{k=1}^{n}y^{k}_{i}&=&p_{\\{i\\}}&\forall~{}~{}i\in V,\\\
&\displaystyle\sum_{k=2}^{n}(k-1)\cdot
y^{k}_{i}&=&\displaystyle\sum_{R\in\Omega^{2}\atop i\in
R}p_{R}&\forall~{}~{}i\in V\\\\[19.91692pt]
&y^{k}_{i}&\geq&0&\begin{array}[]{c}\forall~{}~{}i\in V,\text{ and}\\\
k=1,...,n.\end{array}\end{array}$ (5)

They proved that this maximization problem has the same maximum value as the
upper bound by [24] and the minimization problem has the same value as the
lower bound by [22]. Let us denote by $PG$ the linear programming problem
formed by (5) and by $LB(PG)$ and $UB(PG)$ the corresponding minimum and
maximum values.

We can observe that in this case however
$L_{\mathcal{E}^{1}}(\mathbb{R}_{+}^{\Omega})\neq\mathbb{R}_{+}^{n^{2}}$, thus
formulation (5) could potentially be tightened up.

To see this, let us observe first that by (4) we have the equalities

$\sum_{{S\subseteq V\atop|S|=k}\atop i\not\in
S=\emptyset}x_{S}~{}~{}=~{}~{}\sum_{j\in
V\setminus\\{i\\}}y^{k}_{j}~{}~{}-~{}~{}(k-1)\cdot y^{k}_{i}$

for all $i\in V$ and $k=1,...,n$. Thus, the inequalities

$\sum_{j\in V\setminus\\{i\\}}y^{k}_{j}~{}~{}-~{}~{}(k-1)\cdot
y^{k}_{i}~{}~{}\geq~{}~{}0$ (6)

are valid inequalities for $L_{\mathcal{E}^{1}}(\mathbb{R}_{+}^{\Omega})$, for
all $i\in V$ and $1\leq k\leq n$, and they are not linear consequences of the
nonnegativities in (5).

Let us introduce the convex cone

$U^{1}~{}~{}~{}=~{}~{}~{}\left\\{\mathbf{y}\in\mathbb{R}_{+}^{n^{2}}\left|\begin{array}[]{rcl}y^{k}_{i}&\geq&0\\\\[14.22636pt]
\displaystyle\sum_{j\in V\setminus\\{i\\}}y^{k}_{j}-(k-1)\cdot
y^{k}_{i}&\geq&0\end{array}~{}~{}~{}\begin{array}[]{c}\forall~{}i\in V\text{
and}\\\ k=1,...,n.\\\ $~{}~{}~{}$\end{array}\right.\right\\}$

Then we have the relations

$L_{\mathcal{E}^{1}}(\mathbb{R}_{+}^{\Omega})~{}\subseteq~{}U^{1}~{}\subsetneq~{}\mathbb{R}_{+}^{n^{2}}.$
(7)

Let us pause for a moment, and note that by its definition
$L_{\mathcal{E}^{1}}(\mathbb{R}_{+}^{\Omega})$ is a convex cone, however, it
is far from obvious that it is a polyhedral convex cone. Our first result in
this section is to prove that in fact it is a polyhedral convex cone, and we
provide a full characterization of its facets.

For real numbers $\alpha_{0}$, $\alpha^{k}_{i}$, $i\in V$, $k=1,...,n$ we call
the inequality

$\alpha_{0}~{}+~{}\sum_{k=1}^{n}\sum_{i\in V}\alpha^{k}_{i}\cdot
y^{k}_{i}~{}~{}\geq~{}~{}0$ (8)

a _valid inequality for $L_{\mathcal{E}^{1}}(\mathbb{R}_{+}^{\Omega})$_ if it
holds for all $\mathbf{x}\in\mathbb{R}_{+}^{\Omega}$ under the substitutions
(4). Our first observation is that the tightest valid inequalities have a
special format.

###### Lemma 2

If (8) is a tightest valid inequality for
$L_{\mathcal{E}^{1}}(\mathbb{R}_{+}^{\Omega})$, then we must have
$\alpha_{0}=0$, and there exists an integer $1\leq k\leq n$ such that
$\alpha^{\ell}_{i}=0$ for all $i\in V$ and $\ell\neq k$.

Proof: Let us note first that for $\mathbf{x}=0\in\mathbb{R}_{+}^{\Omega}$ we
have $y^{k}_{i}=0$ for all $i\in V$ and $k=1,...,n$. Thus, $\alpha_{0}\geq 0$
follows.

We claim next, that for an arbitrary value of $k$ ($1\leq k\leq n$) the
inequality

$\sum_{i\in V}\alpha^{k}_{i}\cdot y^{k}_{i}~{}~{}\geq~{}~{}0$

is also a valid inequality for $L_{\mathcal{E}^{1}}(\mathbb{R}_{+}^{\Omega})$.
For this let us denote the left hand side of the above inequality as
$F^{k}(\mathbf{x})$ for $\mathbf{x}\in\mathbb{R}_{+}^{\Omega}$, $k=1,...,n$,
assuming the substitutions (4). Let us also note that this expression depends
only on variables $x_{S}$ with $|S|=k$. Consequently we can assume that
$x_{S}=0$ for all $S\in\Omega$ with $|S|\neq k$. Assume then indirectly that
the above inequality is not valid, that is that there exits a vector
$\mathbf{x}\in\mathbb{R}_{+}^{\Omega}$ with $x_{S}=0$ for all $S\in\Omega$
with $|S|\neq k$ such that $F^{k}(\mathbf{x})<0$. For such a vector, and for
any positive real $\lambda>0$ we have that the left hand side of (8) is equal
to

$\alpha_{0}+F^{k}(\lambda\cdot\mathbf{x})~{}=~{}\alpha_{0}+\lambda\cdot
F^{k}(\mathbf{x})$

which for large $\lambda$ values would turn negative, contradicting the
validity of (8). This contradiction, when applied for all values of $k$, and
the nonnegativity of $\alpha_{0}$ implies that (8) is a linear consequence of
the inequalities

$\sum_{i\in V}\alpha^{k}_{i}\cdot y^{k}_{i}~{}~{}\geq~{}~{}0~{}~{}~{}\forall
k=1,...,n.$

$\Box$

###### Lemma 3

Given reals $\alpha_{i}$, $i\in V$, and a integer $1\leq k\leq n$, the
inequality

$\sum_{i\in V}\alpha_{i}\cdot y^{k}_{i}~{}~{}\geq~{}~{}0$ (9)

is valid for $L_{\mathcal{E}^{1}}(\mathbb{R}_{+}^{\Omega})$ if and only if for
all subsets $S\in\Omega^{k}$ we have

$\sum_{i\in S}\alpha_{i}~{}~{}\geq~{}~{}0.$

Proof: By Lemma 2 all valid inequalities for
$L_{\mathcal{E}^{1}}(\mathbb{R}_{+}^{\Omega})$ have the form like in (9).
Viewing the left hand side as a linear combination of the $x_{S}$,
$S\in\Omega^{k}$ variables, it is easy to see that

$\sum_{i\in S}\alpha_{i}$

is the coefficient of $x_{S}$ in this linear combination, for all
$S\in\Omega^{k}$. Thus, if this quantity was negative for a subset
$S\in\Omega^{k}$, then changing the value of $x_{S}$ to a larger value, while
not changing the values of the other $x_{P}$, $P\in\Omega^{k}$, $P\neq S$
variables would strictly decrease (at a linear rate) the left hand side of
(9). Since $x_{S}$ could take an arbitrarily large value, the claim follows.
$\Box$

Our main result in this section is the following equality:

###### Theorem 4

$L_{\mathcal{E}^{1}}(\mathbb{R}_{+}^{\Omega})~{}=~{}U^{1}.$

Proof: According to (7), $U^{1}$ is not smaller than
$L_{\mathcal{E}^{1}}(\mathbb{R}_{+}^{\Omega})$. Assume indirectly that
$U^{1}\neq L_{\mathcal{E}^{1}}(\mathbb{R}_{+}^{\Omega})$, and choose a
feasible solution $\mathbf{y}\in U^{1}\setminus
L_{\mathcal{E}^{1}}(\mathbb{R}_{+}^{\Omega})$. This implies that $\mathbf{y}$
must violate one of the valid inequalities for
$L_{\mathcal{E}^{1}}(\mathbb{R}_{+}^{\Omega})$. According to Lemmas 2 and 3,
this (violated) inequality should look like

$\sum_{i\in V}\alpha_{i}\cdot y^{k}_{i}~{}~{}<~{}~{}0.$ (10)

According to Lemma 3, permuting the coefficients in this inequality also
yields a valid inequality for $L_{\mathcal{E}^{1}}(\mathbb{R}_{+}^{\Omega})$.
Thus we can assume a permutation that yields the lowest left hand side in
(10), i.e., we can assume that

$\begin{array}[]{c}y^{k}_{1}~{}\geq~{}y^{k}_{2}~{}\geq~{}\cdots~{}\geq~{}y^{k}_{n}\\\
\alpha_{1}~{}\leq~{}\alpha_{2}~{}\leq~{}\cdots~{}\leq~{}\alpha_{n}.\end{array}$
(11)

Note that by (10) we must have $\alpha_{1}<0$, and by Lemma 3
$\alpha_{1}+\cdots+\alpha_{k}\geq 0$, implying $\alpha_{k}>0$.

Note that after the sorting of the components as in (11), if we decrease
$\alpha_{i}$ by epsilon, and increase $\alpha_{j}$, $j>i$ by the same
$\epsilon$, then the left hand side of (10) can only decrease. Let us perform
this operations, decreasing $\alpha_{1}$ and increasing $\alpha_{j}$ for
$2\leq j<k$ until we arrive to the vector

$\begin{array}[]{rcl}\beta_{i}~{}=~{}\begin{cases}\alpha_{1}-(k-2)\alpha_{k}+\alpha_{2}+\cdots+\alpha_{k-1}&\text{
if }i=1,\\\ \alpha_{k}&\text{ if }2\leq i\leq k,\\\ \alpha_{i}&\text{ if
}i>k.\end{cases}\end{array}$

Note that we have $\beta_{1}\geq-(k-1)\alpha_{k}$, and
$\beta_{i}\geq\alpha_{k}$ for all $i=2,...,n$. Thus we can write the chain of
inequalities

$0>\sum_{i\in V}\alpha_{i}\cdot y^{k}_{i}\geq\sum_{i\in V}\beta_{i}\cdot
y^{k}_{i}\geq\alpha_{k}\left(-(k-1)y^{k}_{1}+\sum_{i=2}^{n}y^{k}_{i}\right)~{}\geq~{}0$

Here the last inequality follows since $\alpha_{k}>0$ and $\mathbf{y}\in
U^{1}$. The derived contradiction proves that
$U^{1}=L_{\mathcal{E}^{1}}(\mathbb{R}_{+}^{\Omega})$, as claimed. $\Box$

Now, we fully characterize facets of a polyhedral convex cone, and (5) is
tighten up with $U^{1}$.

$\begin{array}[]{crcll}\max~{}(\min)&\displaystyle\sum_{k=1}^{n}\frac{1}{k}\sum_{i\in
V}y^{k}_{i}\\\\[14.22636pt] s.t.\\\
&\displaystyle\sum_{k=1}^{n}y^{k}_{i}&=&p_{\\{i\\}}&~{}~{}\forall~{}~{}i\in
V,\\\ &\displaystyle\sum_{k=2}^{n}(k-1)\cdot
y^{k}_{i}&=&\displaystyle\sum_{R\in\Omega^{2}\atop i\in
R}p_{R}&~{}~{}\forall~{}~{}i\in V\\\\[19.91692pt] &\displaystyle\sum_{j\in
V\setminus\\{i\\}}y^{k}_{j}-(k-1)\cdot
y^{k}_{i}&\geq&0&\begin{array}[]{c}\forall~{}~{}i\in V,\text{ and}\\\
k=1,...,n.\end{array}\\\\[19.91692pt]
&y^{k}_{i}&\geq&0&\begin{array}[]{c}\forall~{}~{}i\in V,\text{ and}\\\
k=1,...,n.\end{array}\end{array}$ (12)

Let us denote by $IPG$ the linear programming problem formed by (12) and by
$LB(IPG)$ and $UB(IPG)$ the corresponding minimum and maximum values.

| $UB(PG)$ | $UB(IPG)$ | $LB(PG)$ | $LB(IPG)$
---|---|---|---|---
n | Mean | Std | Max | Mean | Std | Max | Mean | Std | Max | Mean | Std | Max
7 | 5.56 | 1.16 | 8.67 | 3.49 | 1.18 | 6.61 | 0.99 | 1.46 | 6.66 | 0.88 | 1.41 | 6.66
8 | 6.74 | 1.19 | 10.3 | 4.41 | 1.23 | 7.88 | 1.32 | 1.80 | 9.72 | 0.73 | 1.65 | 9.06
9 | 7.84 | 1.09 | 10.6 | 5.38 | 1.27 | 8.94 | 1.65 | 2.02 | 9.73 | 1.06 | 1.81 | 9.02
10 | 8.83 | 1.25 | 12.1 | 6.52 | 1.19 | 9.77 | 2.00 | 1.41 | 7.07 | 0.98 | 1.35 | 5.66
11 | 9.88 | 1.28 | 13.5 | 7.29 | 1.43 | 10.6 | 1.68 | 1.19 | 5.28 | 0.73 | 1.04 | 5.28
12 | 11.1 | 1.10 | 13.9 | 8.53 | 1.26 | 11.5 | 2.02 | 1.54 | 8.90 | 1.25 | 1.33 | 6.12
13 | 11.8 | 1.22 | 15.0 | 9.09 | 1.42 | 12.9 | 1.63 | 1.38 | 5.33 | 0.85 | 1.27 | 4.03
14 | 12.9 | 1.23 | 16.9 | 10.1 | 1.49 | 14.3 | 1.76 | 0.75 | 3.06 | 1.46 | 0.91 | 3.06
15 | 13.9 | 1.33 | 17.2 | 11.1 | 1.48 | 14.4 | 1.91 | 0.88 | 3.47 | 1.60 | 0.82 | 3.21
16 | 15.1 | 1.51 | 18.5 | 12.3 | 1.53 | 15.8 | 1.83 | 0.99 | 4.12 | 1.25 | 0.85 | 3.02
17 | 15.8 | 1.30 | 18.4 | 12.9 | 1.40 | 15.8 | 2.47 | 1.24 | 4.79 | 2.21 | 1.29 | 4.78
18 | 16.6 | 1.17 | 19.1 | 13.7 | 1.26 | 16.7 | 2.13 | 1.24 | 4.48 | 1.79 | 1.11 | 4.00
Table 1: Mean, standard deviation, and maximum value of the percentage
relative error.

In Table 1 we compare the bounds. For upper bounds, we generate 100 instances
for n from 7 to 15; for n from 16 to 18, we generate 50 instances. And, for
lower bounds, we generate 100 instances for n from 7 to 13; for n from 14 to
18, we generate 20 instances. Let _opt_ be the tightest possible bound, which
can be achieved by solving (1). Then we calculate the percentage relative
error by $100\times\lvert~{}\text{bound}~{}-~{}opt~{}\rvert~{}/~{}opt$. We can
see that in many examples the addition of $U^{1}$ improves strictly both the
lower and upper bounds.

##  5 The Aggregation by Yang et al.

In [28] a novel aggregation model was considered that in our terms and
notation corresponds to the case of

$\mathcal{E}^{2}~{}=~{}\\{E^{k}_{Q}\mid
Q\in\Omega^{1}\cup\Omega^{2},~{}~{}|Q|\leq k\leq n\\},$

where

$E^{k}_{Q}~{}=~{}\\{S\subseteq V\mid Q\subseteq
S,~{}~{}|S|=k\\}~{}~{}~{}\forall~{}Q\in\Omega^{1}\cup\Omega^{2},~{}~{}|Q|\leq
k\leq n.$

In other words, they consider the aggregation, defined by the following
mapping $L_{\mathcal{E}^{2}}(\mathbf{x})=\mathbf{y}=(y^{k}_{Q}\mid
Q\in\Omega^{1}\cup\Omega^{2},~{}~{}|Q|\leq k\leq n)$, where

$y^{k}_{Q}~{}=~{}\sum_{S\in\Omega^{k}\atop Q\subseteq
S}x_{S}~{}~{}~{}\forall~{}Q\in\Omega^{1}\cup\Omega^{2},~{}~{}|Q|\leq k\leq n.$
(13)

Note also that this notation is somewhat redundant, since we have the
equalities

$y^{k}_{\\{i\\}}~{}=~{}\frac{1}{k-1}\sum_{Q\in\Omega^{2}\atop i\in
Q}y^{k}_{Q}$ (14)

for all $i\in V$ and $2\leq k\leq n$. Note finally that in general there are
no other linear relations between these variables. Thus in fact, this
aggregation is a linear mapping $L_{\mathcal{E}^{2}}$ into the space of
dimension $N=n+(n-1)\cdot\binom{n}{2}$.

It is easy to verify that with this aggregation the objective function and the
left hand sides of the equality constraints in (2) can be represented as
linear expressions, Namely we have

$\displaystyle\displaystyle\sum_{S\in\Omega}x_{S}~{}=~{}$
$\displaystyle\displaystyle\sum_{Q\in\Omega^{1}}y^{1}_{Q}~{}+~{}\sum_{k=2}^{n}\frac{1}{\binom{k}{2}}\sum_{Q\in\Omega^{2}}y^{k}_{Q}$
$\displaystyle\rightarrow$ $\displaystyle~{}\max~{}$ ($E1$)
$\displaystyle\displaystyle\sum_{S\in\Omega\atop i\in S}x_{S}~{}=~{}$
$\displaystyle\displaystyle
y^{1}_{\\{i\\}}~{}+~{}\sum_{k=2}^{n}\frac{1}{k-1}\sum_{Q\in\Omega^{2}\atop
i\in Q}y^{k}_{Q}$ $\displaystyle~{}=~{}$ $\displaystyle p_{\\{i\\}}$
$\displaystyle\forall~{}i\in V$ ($E2$)
$\displaystyle\displaystyle\sum_{S\in\Omega\atop Q\subseteq S}x_{S}~{}=~{}$
$\displaystyle\displaystyle\sum_{k=2}^{n}y^{k}_{Q}$ $\displaystyle~{}=~{}$
$\displaystyle p_{Q}$ $\displaystyle\forall~{}Q\in\Omega^{2}$ ($E3$)
$\displaystyle y^{k}_{\\{i\\}}$ $\displaystyle~{}\geq~{}$ $\displaystyle 0$
$\displaystyle\forall~{}\begin{array}[]{c}i\in V\\\ 1\leq k\leq n\end{array}$
($E4$) $\displaystyle y^{k}_{Q}$ $\displaystyle~{}\geq~{}$ $\displaystyle 0$
$\displaystyle\forall~{}\begin{array}[]{c}Q\in\Omega^{2}\\\ 2\leq k\leq
n\end{array}$ ($E5$)

leading to a quite tight polynomial size aggregation of our exponential sized
original problem. Note that inequalities ($E4$) include redundant ones,
according to (14), but we keep this redundancy, since it will help us to
realize a more general structure to these inequalities.

The nonnegativity of the $y^{1}_{\\{i\\}}$, $i\in V$ and $y^{k}_{Q}$,
$Q\in\Omega^{2}$ variables are implied by their definition (13), since all
$x_{S}$ variables are nonnegative. The authors of this aggregation however
noticed that in this case
$L_{\mathcal{E}^{2}}(\mathbb{R}_{+}^{\Omega})\neq\mathbb{R}_{+}^{N}$. In fact,
several other nonnegative combination of the $x_{S}$ variables are linear
functions of $\mathbf{y}$. For instance, we can note that for $i\in V$ and
integers $2\leq k\leq n$ we have

$\sum_{S\in\Omega^{k}\atop i\not\in
S}x_{S}~{}~{}=~{}~{}\frac{1}{\binom{k}{2}}\sum_{Q\in\Omega^{2}}y^{k}_{Q}-\frac{1}{k-1}\sum_{Q\in\Omega^{2}\atop
i\in Q}y^{k}_{Q}~{}\geq~{}0$

Note that the nonnegativity of this expression is not a consequence of the
nonnegativity of $\mathbf{y}$. Thus the nonnegativity of this expression could
be added to the aggregation, tightening it up. In fact in [28] the authors
came up with the following additional nonnegative combinations of the
$\mathbf{x}$ variables that can be expressed as a linear function of the
$\mathbf{y}$ variables when we use substitutions (13).

$\sum_{S\in\Omega^{k}\atop R\cap
S=\emptyset}x_{S}~{}=~{}\frac{1}{\binom{k}{2}}\left(\sum_{Q\in\Omega^{2}}y^{k}_{Q}-\frac{k}{2}\sum_{Q\in\Omega^{2}\atop|R\cap
Q|=1}y^{k}_{Q}+\frac{k(k-3)}{2}y^{k}_{R}\right)~{}\geq~{}0$

for all $R\in\Omega^{2}$, and $2\leq k\leq n$.

$\sum_{{S\in\Omega^{k}}\atop R\cap
S=\\{i\\}}x_{S}~{}=~{}\frac{1}{k-1}\left(\sum_{Q\in\Omega^{2}\atop Q\cap
R=\\{i\\}}y^{k}_{Q}-(k-2)y^{k}_{R}\right)~{}\geq~{}0$

for all $i\in V$, $i\in R\in\Omega^{2}$, and $2\leq k\leq n$.

$\sum_{S\in\Omega^{k}\atop|R\cap
S|\in\\{0,3\\}}x_{S}~{}=~{}\frac{1}{\binom{k}{2}}\sum_{Q\in\Omega^{2}}y^{k}_{Q}-\frac{1}{k-1}\sum_{Q\in\Omega^{2}\atop|R\cap
Q|=1}y^{k}_{Q}+\frac{k-3}{k-1}\sum_{Q\in\Omega^{2}\atop Q\subseteq
R}y^{k}_{Q}~{}\geq~{}0$

for all $R\in\Omega^{3}$ and $2\leq k\leq n$.

$\sum_{{{S\in\Omega^{k}\atop R\cap S=\\{i\\}}\atop\textbf{or}}\atop R\setminus
S=\\{i\\}}x_{S}~{}=~{}\frac{1}{k-1}\sum_{Q\in\Omega^{2}\atop Q\cap
R=\\{i\\}}y^{k}_{Q}+y^{k}_{R\setminus\\{i\\}}-\frac{k-2}{k-1}\sum_{Q\in\Omega^{2}\atop
i\in Q\subseteq R}y^{k}_{Q}~{}\geq~{}0$

for all $i\in V$, $i\in R\in\Omega^{3}$, and $2\leq k\leq n$.

Let us remark that the polynomial sized linear programming aggregation ($E1$)
- (5) proposed by [28] provides a bound for our original problem that is not
weaker than the one proposed by [27] in the previous subsection. This is
because problem (5) and its strenghtening by $U^{1}$ are both aggregations of
problem ($E1$) - (5). To see this, note that $y^{1}_{i}$, $i\in V$ have
identical definitions in both aggregations (see (4) and (13)), and for
$y^{k}_{i}$, $i\in V$ and $k=2,...,n$ we have the equalities (14). Furthermore
the inequalities defining $U^{1}$ are consequences of ($E4$) - (5).

The question arise: do the inequalities ($E4$)-(5) describe the cone
$L_{\mathcal{E}^{2}}(\mathbb{R}^{\Omega}_{+})$?  The answer is no. In what
follows we provide an exact polyhedral description of
$L_{\mathcal{E}^{2}}(\mathbb{R}^{\Omega}_{+})$.

Let us start observing that the analogue of Lemma 2 can be shown for this
case, too.

For real numbers $\alpha_{0}$, $\alpha^{k}_{Q}$,
$Q\in\Omega^{1}\cup\Omega^{2}$, $k=1,...,n$ we call the inequality

$\alpha_{0}~{}+~{}\sum_{k=1}^{n}\sum_{Q\in\Omega^{1}\cup\Omega^{2}}\alpha^{k}_{Q}\cdot
y^{k}_{Q}~{}~{}\geq~{}~{}0$ (15)

a _valid inequality for $L_{\mathcal{E}^{2}}(\mathbb{R}_{+}^{\Omega})$_ if it
holds for all $\mathbf{x}\in\mathbb{R}_{+}^{\Omega}$ under the substitutions
(13).

###### Lemma 5

If (15) is a tightest valid inequality for
$L_{\mathcal{E}^{2}}(\mathbb{R}_{+}^{\Omega})$, then we must have
$\alpha_{0}=0$, and there exists an integer $1\leq k\leq n$ such that
$\alpha^{\ell}_{Q}=0$ for all $Q\in\Omega^{1}\cup\Omega^{2}$ and $\ell\neq k$.

Proof: Let us note first that for $\mathbf{x}=0\in\mathbb{R}_{+}^{\Omega}$ we
have $y^{k}_{Q}=0$ for all $Q\in\Omega^{1}\cup\Omega^{2}$ and $k=1,...,n$.
Thus, $\alpha_{0}\geq 0$ follows.

We claim next, that for an arbitrary value of $k$ ($1\leq k\leq n$) the
inequality

$\sum_{k=1}^{n}\sum_{Q\in\Omega^{1}\cup\Omega^{2}}\alpha^{k}_{Q}\cdot
y^{k}_{Q}~{}~{}\geq~{}~{}0$

is also a valid inequality for $L_{\mathcal{E}^{2}}(\mathbb{R}_{+}^{\Omega})$.
For this let us denote the left hand side of the above inequality as
$F^{k}(\mathbf{x})$ for $\mathbf{x}\in\mathbb{R}_{+}^{\Omega}$, $k=1,...,n$,
assuming the substitutions (13). Let us also note that this expression depends
only on variables $x_{S}$ with $|S|=k$. Consequently we can assume that
$x_{S}=0$ for all $S\in\Omega$ with $|S|\neq k$. Assume then indirectly that
the above inequality is not valid, that is that there exits a vector
$\mathbf{x}\in\mathbb{R}_{+}^{\Omega}$ with $x_{S}=0$ for all $S\in\Omega$
with $|S|\neq k$ such that $F^{k}(\mathbf{x})<0$. For such a vector, and for
any positive real $\lambda>0$ we have that the left hand side of (15) is equal
to

$\alpha_{0}+F^{k}(\lambda\cdot\mathbf{x})~{}=~{}\alpha_{0}+\lambda\cdot
F^{k}(\mathbf{x})$

which for large $\lambda$ values would turn negative, contradicting the
validity of (15). This contradiction, when applied for all values of $k$, and
the nonnegativity of $\alpha_{0}$ implies that (15) is a linear consequence of
the inequalities

$\sum_{k=1}^{n}\sum_{Q\in\Omega^{1}\cup\Omega^{2}}\alpha^{k}_{Q}\cdot
y^{k}_{Q}~{}~{}\geq~{}~{}0~{}~{}~{}\forall~{}k=1,...,n.$

$\Box$

Thus we can focus on subsets of a certain size. Let us also note that for sets
$S\in\Omega^{1}$ we have $x_{S}=y^{1}_{S}\geq 0$ by ($E4$), thus we can focus
on set sizes $2\leq k\leq n$. Note also that all inequalities ($E4$)-(5)
involve $y^{k}_{Q}$ variables for a fixed value of $k$.

For notational simplicity let us introduce $\Omega^{0}=\\{\emptyset\\}$, set
$\Sigma=\Omega^{0}\cup\Omega^{1}\cup\Omega^{2}$, and define for $2\leq k\leq
n$

$y^{k}_{\emptyset}~{}=~{}\sum_{S\in\Omega^{k}}x_{S}~{}=~{}\frac{1}{\binom{k}{2}}\sum_{Q\in\Omega^{2}}y^{k}_{Q},$
(16)

where the second equality follows by the definitions of the $\mathbf{y}^{k}$
variables, see (13) and (16).

Let us now fix a $2\leq k\leq n$ value, and introduce
$\mathbf{x}^{k}=(x_{S}\mid S\in\Omega^{k})$ and $\mathbf{y}^{k}=(y^{k}_{Q}\mid
Q\in\Sigma)$. Recall that the $y^{k}_{Q}$ variables are in fact (linear)
functions of $\mathbf{x}^{k}$ according to (13) and (16). Thus, for a real
vector $\alpha\in\mathbb{R}^{\Sigma}$ we say that the inequality

$\sum_{Q\in\Sigma}\alpha_{Q}\cdot y^{k}_{Q}~{}\geq~{}0$ (17)

is _valid_ if it holds for all nonnegative
$\mathbf{x}^{k}\in\mathbb{R}^{\Omega^{k}}_{+}$. We denote the left hand side
of inequality (17) by $L_{\alpha}(*)$, where $*$ indicates the set of
variables we want to use. Thus inequality (17), as it is written, is the
expression $L_{\alpha}(\mathbf{y}^{k})\geq 0$. We can also view the left hand
side as a function of $\mathbf{x}^{k}$, and then we write
$L_{\alpha}(\mathbf{x}^{k})$. Due to the relations (13) and (16), we can also
write $L_{\alpha}(y^{k}_{Q}\mid Q\in\Omega^{2})$. Note finally that the same
$\alpha$ vector can be used as coefficients for different $k$ values, e.g.,
$L_{\alpha}(\mathbf{y}^{\ell})~{}=~{}L_{\alpha}(\mathbf{x}^{\ell})$ for
$\ell\neq k$ are also well defined.

In what follows we characterize all valid inequalities (for all $2\leq k\leq
n$), and hence provide a complete description for
$L_{\mathcal{E}^{2}}(\mathbb{R}^{\Omega}_{+})$. To arrive to such a complete
description, we consider a (perhaps surprising) mapping into the space of
nonnegative quadratic pseudo-Boolean functions.

Introduce $\mathbf{Z}=(Z_{i}\mid i\in V)\in\\{0,1\\}^{V}$, and associate to
inequality (17) the quadratic pseudo-Boolean function (or QPBF in short)

$F_{\alpha}(\mathbf{Z})~{}=~{}\sum_{Q\in\Sigma}\alpha_{Q}\cdot\prod_{i\in
Q}Z_{i},$ (18)

where we have $\prod_{i\in\emptyset}Z_{i}~{}=~{}1$ by definition. For a subset
$S\subseteq V$ we denote by $\chi_{S}\in\\{0,1\\}^{V}$ its characteristic
vector.

###### Lemma 6

Given a real vector $\alpha\in\mathbb{R}^{\Sigma}$ and a subset $\emptyset\neq
S\subseteq V$, the real value $F_{\alpha}(\chi_{S})$ is the coefficient of
$x_{S}$ in $L_{\alpha}(\mathbf{x}^{|S|})$.

Proof: Observe that $x_{S}$ appears in $y^{k}_{Q}$ with coefficient $1$
exactly when $Q\subseteq S$, and we have

$F_{\alpha}(\chi_{S})~{}=~{}\sum_{Q\in\Sigma\atop Q\subseteq S}\alpha_{Q}$

by (18). $\Box$

###### Corollary 1

If $F_{\alpha}(\mathbf{Z})\geq 0$ for all $\mathbf{Z}\in\\{0,1\\}^{V}$, then
$L_{\alpha}(\mathbf{y}^{k})\geq 0$ is a valid inequality for all $k=2,...,n$.

Proof: By Lemma 6 the coefficients of the $x_{S}$, $S\in\Omega$ variables are
all nonnegative, since $F_{\alpha}$ is a nonnegative function. Thus
$L_{\alpha}(\mathbf{y}^{k})$ is a nonnegative real for all
$\mathbf{x}^{k}\in\mathbb{R}^{\Sigma}_{+}$. $\Box$

Since linear combinations of QPBF-s is again a QPBF, the set of nonnegative
QPBF-s form a convex cone (in the space of their coefficients, i.e., in
dimension $1+n+\binom{n}{2}$. It is well-known that this cone is polyhedral,
and its extremal rays are in a one-to-one correspondence with the facets of
the cut polytope [12]. Let us denote by $\mathcal{G}$ the finite set of
extremal nonnegative QPBF-s, and thus $cone(\mathcal{G})$ is the set of
nonnegative QPBF-s.

Let us now consider the polyhedral cone defined by the valid inequalities
corresponding to nonnegative QPBFs.

$U^{2}~{}=~{}\left\\{\mathbf{y}~{}\left|~{}\begin{array}[]{rlc}y^{1}_{\\{i\\}}&\geq~{}0&\forall~{}i\in
V,\\\\[14.22636pt] \displaystyle\sum_{Q\in\Sigma}\alpha_{Q}\cdot
y^{k}_{Q}&\geq~{}0&\begin{array}[]{c}\forall~{}\alpha\in\mathbb{R}^{\Sigma}\text{
and }2\leq k\leq n\\\ \text{such that
}F_{\alpha}(\mathbf{Z})\in\mathcal{G}\end{array}\end{array}\right.\right\\}$
(19)

Since by Corollary 1 all these inequalities are valid for
$L_{\mathcal{E}^{2}}(\mathbb{R}^{\Omega}_{+})$, we have the relation

$L_{\mathcal{E}^{2}}(\mathbb{R}^{\Omega}_{+})~{}\subseteq~{}U^{2}$

In what follows we show that in fact we have equality here. The main
difficulty in proving this claim stems from the fact that an inequality of the
form (17) may correspond to multiple QPBF-s $F_{\alpha}(\mathbf{Z})$. This is
because of the linear dependencies, like equations (14) and (16), we have in
the $\mathbf{y}^{k}$ space. In fact, to a given QPBF we have several
corresponding inequalities of the form (17), with different $k$ values. Some
may be valid for $L_{\mathcal{E}^{2}}(\mathbb{R}^{\Omega}_{+})$ and some may
not. More precisely, a real vector $\alpha$ as coefficients in such an
inequality may yield a valid one with variables $\mathbf{y}^{k}$, while if we
replace $\mathbf{y}^{k}$ with $\mathbf{y}^{\ell}$ for some $\ell\neq k$, then
the same expression may not be a valid inequality. Thus we need a mechanism
that can transform a QPBF into another one, without changing a corresponding
inequality (for a particular value of $k$).

A very helpful observation is the following claim:

###### Lemma 7

Given an integer $2\leq k\leq n$, let $\alpha\in\mathbb{R}^{\Sigma}$ be the
coefficient vector of the QPBF

$\left(k-\sum_{i\in V}Z_{i}\right)^{2}~{}=~{}k^{2}-(2k-1)\cdot\sum_{i\in
V}Z_{i}+2\cdot\sum_{1\leq i<j\leq n}Z_{i}Z_{j}.$

Then $L_{\alpha}(y^{k}_{Q}\mid Q\in\Omega^{2})$ is the identically zero
function.

Proof: Let us observe first that the equality in the above formula holds,
because for binary variables we have $Z_{i}^{2}=Z_{i}$ for all $i\in V$. Thus
we have

$L_{\alpha}(\mathbf{y}^{k})~{}=~{}k^{2}y^{k}_{\emptyset}-(2k-1)\sum_{Q\in\Omega^{1}}y^{k}_{Q}+2\sum_{Q\in\Omega^{2}}y^{k}_{Q}.$

The claim now follows by elementary algebra, using the equalities (14) and
(16).

$\displaystyle L_{\alpha}(y^{k}_{Q}\mid Q\in\Omega^{2})$
$\displaystyle=~{}k^{2}\left(\frac{1}{\binom{k}{2}}\sum_{Q\in\Omega^{2}}y^{k}_{Q}\right)-(2k-1)\sum_{i\in
V}\left(\frac{1}{k-1}\sum_{Q\in\Omega^{2}\atop i\in Q}y^{k}_{Q}\right)$
$\displaystyle~{}~{}~{}+2\sum_{Q\in\Omega^{2}}y^{k}_{Q}$

Thus, for an arbitrary $Q\in\Omega^{2}$ the coefficient of $y^{k}_{Q}$ in the
above expression is

$k^{2}\frac{1}{\binom{k}{2}}-(2k-1)\frac{2}{k-1}+2~{}~{}~{}=~{}~{}~{}0.$

$\Box$

Note that $L_{\alpha}(y^{\ell}_{Q}\mid Q\in\Omega^{2})$ may not be identically
zero if $\ell\neq k$.

Now, we are ready to claim our main result in this section.

###### Theorem 8

$L_{\mathcal{E}^{2}}(\mathbb{R}^{\Omega}_{+})~{}=~{}U^{2}$

Proof: For $k=1$ we have $y^{1}_{i}=x_{\\{i\\}}$ for all $i\in V$, and clearly
no other inequality of these variables can be minimally valid. For $k\geq 2$,
the containment“$\subseteq$” is implied by Corollary 1. To see the reverse
containment, let us assume that $\alpha^{*}\in\mathbb{R}^{\Sigma}$ is a real
vector for which the inequality

$\sum_{Q\in\Sigma}\alpha^{*}_{Q}y^{k}_{Q}~{}\geq~{}0$

is valid, and for which $F_{\alpha^{*}}(\mathbf{Z})$ is not a nonnegative
QPBF. Let us introduce

$\gamma~{}=~{}\min_{\mathbf{Z}\in\\{0,1\\}^{V}}F_{\alpha^{*}}(\mathbf{Z})~{}<~{}0.$

Then,

$F_{\alpha^{\prime}}(\mathbf{Z})~{}=~{}F_{\alpha^{*}}(\mathbf{Z})-\gamma\left(k-\sum_{i\in
V}Z_{i}\right)^{2}$

is a nonnegative QPBF, and $L_{\alpha^{\prime}}(y^{k}_{Q}\mid
Q\in\Omega^{2})\geq 0$ is identical to $L_{\alpha^{*}}(y^{k}_{Q}\mid
Q\in\Omega^{2})\geq 0$ by Lemma 7. $\Box$

Note that the inequalities introduced by [28] correspond to some known
extremal nonnegative QPBFs, see [5] or see Proposition 18 in [6]. Namely,
inequality ($E4$) corresponds to $Z_{i}\geq 0$. Inequality ($E5$) corresponds
to $Z_{i}Z_{j}\geq 0$. Inequality (5) corresponds to $Z_{i}(1-Z_{j})\geq 0$.
Inequality (5) corresponds to $(1-Z_{i})(1-Z_{j})\geq 0$. Inequality (5)
corresponds to

$Z_{i}Z_{j}Z_{\ell}+(1-Z_{i})(1-Z_{j})(1-Z_{\ell})~{}=~{}1-Z_{i}-Z_{j}-Z_{\ell}+Z_{i}Z_{j}+Z_{i}Z_{\ell}+Z_{j}Z_{\ell}~{}\geq~{}0.$

Finally, inequality (5) corresponds to

$Z_{i}(1-Z_{j})(1-Z_{\ell})+(1-Z_{i})Z_{j}Z_{\ell}~{}=~{}Z_{i}-Z_{i}Z_{j}-Z_{i}Z_{\ell}+Z_{j}Z_{\ell}~{}\geq~{}0.$

These nonnegative quadratic functions are known to belong to $\mathcal{G}$,
see [4]. There are however many more members in $\mathcal{G}$ that can be
added to tighten up the formulation. In fact the set $\mathcal{G}$ is
partitioned
$\mathcal{G}=\mathcal{G}^{2}\cup\mathcal{G}^{3}\cup\dots\cup\mathcal{G}^{n}$
into families of increasing complexity, and it is known that none of these is
empty, and together they contain exponentially many functions [4].

By the above analysis and results, the model proposed by [28] utilizes exactly
the members of $\mathcal{G}^{2}\cup\mathcal{G}^{3}$. Let us denote by $YAT$
their model and by $LB(YAT)$ and $UB(YAT)$ the corresponding minimum and
maximum values.

$\begin{array}[]{crcll}\max~{}(\min)&\displaystyle\sum_{Q\in\Omega^{1}}y^{1}_{Q}~{}+~{}\sum_{k=2}^{n}\frac{1}{\binom{k}{2}}\sum_{Q\in\Omega^{2}}y^{k}_{Q}\\\
s.t.\\\ &\displaystyle
y^{1}_{\\{i\\}}~{}+~{}\sum_{k=2}^{n}\frac{1}{k-1}\sum_{Q\in\Omega^{2}\atop
i\in Q}y^{k}_{Q}&~{}={}&p_{\\{i\\}}&\forall~{}i\in V\\\
&\displaystyle\sum_{k=2}^{n}y^{k}_{Q}&~{}={}&p_{Q}&\forall~{}Q\in\Omega^{2}\\\
&y^{k}_{\\{i\\}}&~{}\geq{}&0&\forall~{}\begin{array}[]{c}i\in V\\\ 1\leq k\leq
n\end{array}\\\\[14.22636pt] &\displaystyle\sum_{Q\in\Sigma}\alpha_{Q}\cdot
y^{k}_{Q}&\geq{}&0&\forall~{}\begin{array}[]{c}\alpha\in\mathbb{R}^{\Sigma}\text{
and }\\\ 2\leq k\leq n~{}\text{where}\\\
F_{\alpha}(\mathbf{Z})\in\mathcal{G}^{2}\cup\mathcal{G}^{3}\end{array}\end{array}$
(20)

The above theorem provides a surprising connection to the cone of nonnegative
quadratic pseudo-Boolean functions, which is known to be equivalent with the
cut-polytope (see [12]). It was already known that the dual of problem (2) has
strong connections to the cut-polytope [13, 14], but this connection to an
aggregation of the primal problem is different.

Let us denote by $QPB$ the linear programming problem formed by the objective
($E1$), the equalities ($E2$) and ($E3$), and the linear inequalites
(exponentially many) describing $L_{\mathcal{E}^{2}}(\mathbb{R}^{\Omega}_{+})$
as in Theorem 8, and denote by $LB(QPB)$ and $UB(QPB)$ the minimum and maximum
values of it.

$\begin{array}[]{crcll}\max~{}(\min)&\displaystyle\sum_{Q\in\Omega^{1}}y^{1}_{Q}~{}+~{}\sum_{k=2}^{n}\frac{1}{\binom{k}{2}}\sum_{Q\in\Omega^{2}}y^{k}_{Q}\\\
s.t.\\\ &\displaystyle
y^{1}_{\\{i\\}}~{}+~{}\sum_{k=2}^{n}\frac{1}{k-1}\sum_{Q\in\Omega^{2}\atop
i\in Q}y^{k}_{Q}&~{}={}&p_{\\{i\\}}&\forall~{}i\in V\\\
&\displaystyle\sum_{k=2}^{n}y^{k}_{Q}&~{}={}&p_{Q}&\forall~{}Q\in\Omega^{2}\\\
&y^{k}_{\\{i\\}}&~{}\geq{}&0&\forall~{}\begin{array}[]{c}i\in V\\\ 1\leq k\leq
n\end{array}\\\\[14.22636pt] &\displaystyle\sum_{Q\in\Sigma}\alpha_{Q}\cdot
y^{k}_{Q}&\geq{}&0&\forall~{}\begin{array}[]{c}\alpha\in\mathbb{R}^{\Sigma}\text{
and }\\\ 2\leq k\leq n~{}\text{where}\\\
F_{\alpha}(\mathbf{Z})\in\mathcal{G}\end{array}\end{array}$ (21)

###### Corollary 2

We have the equalities

$LB(H)~{}=~{}LB(QPB)~{}~{}~{}\text{ and }~{}~{}~{}UB(H)~{}=~{}UB(QPB).$

Proof: By Theorem 8 for every feasible solution $\mathbf{y}$ of problem $QPB$
there exists a corresponding feasible solution $\mathbf{x}$ of problem (2)
satisfying the equalities (13). $\Box$

In practice, we cannot use all inequalities described by Theorem 8. In this
paper we suggest to use some of the functions in
$\mathcal{G}^{4}\cup\dots\cup\mathcal{G}^{n}$ to tighten up the formulation
($E1$) - (5).

For instance, we recall the following result from [5]: given a subset of the
literals $W\subseteq\\{Z_{1},\dots,Z_{n}\\}\cup\\{1-Z_{1},\dots,1-Z_{n}\\}$ of
cardinality $|W|\geq 4$ and an integer $1\leq\gamma\leq|W|-2$, then the
function defined by following binomial expression

$G_{W,\gamma}(\mathbf{Z})~{}=~{}\binom{\displaystyle\sum_{w\in
W}w~{}-~{}\gamma}{2}$

is a member of $\mathcal{G}^{|W|}$. As an example, consider
$W=\\{Z_{1},Z_{2},1-Z_{3},Z_{4}\\}$ and $\gamma=2$. Then we have

$\displaystyle G_{W,\gamma}(\mathbf{Z})$
$\displaystyle=~{}\binom{Z_{1}+Z_{2}+(1-Z_{3})+Z_{4}-2}{2}$
$\displaystyle=~{}\frac{(Z_{1}+Z_{2}-Z_{3}+Z_{4}-1)(Z_{1}+Z_{2}-Z_{3}+Z_{4}-2)}{2}$
$\displaystyle=~{}\frac{(Z_{1}+Z_{2}-Z_{3}+Z_{4})^{2}-3(Z_{1}+Z_{2}-Z_{3}+Z_{4})+2}{2}$
$\displaystyle=~{}1-Z_{1}-Z_{2}+Z_{3}-Z_{4}+Z_{1}Z_{2}-Z_{1}Z_{3}+Z_{1}Z_{4}-Z_{2}Z_{3}+Z_{2}Z_{4}-Z_{3}Z_{4}$

Please recall that we have $Z_{i}^{2}=Z_{i}$ and $(1-Z_{i})^{2}=1-Z_{i}$ for
all $i\in V$, since these expressions take only binary values. Thus, from this
small example we get the inequality

$y^{k}_{\emptyset}-y^{k}_{\\{1\\}}-y^{k}_{\\{2\\}}+y^{k}_{\\{3\\}}-y^{k}_{\\{4\\}}+y^{k}_{\\{1,2\\}}-y^{k}_{\\{1,3\\}}+y^{k}_{\\{1,4\\}}-y^{k}_{\\{2,3\\}}+y^{k}_{\\{2,4\\}}-y^{k}_{\\{3,4\\}}~{}\geq~{}0.$

Let us introduce these inequalities with $|W|=4$ and $\gamma=1,2$. Since we
utilize some functions in $\mathcal{G}^{4}$, we call this weaker model as
problem $QPB^{-}$ (and denote by $LB(QPB^{-})$ and $UB(QPB^{-})$ the
corresponding minimum and maximum values). Note that $QPB^{-}$ is still a
polynomial sized formulation, and thus we can compute the corresponding lower
and upper bounds in polynomial time.

###### Corollary 3

We have the inequalities

$\begin{array}[]{rclll}LB(H)&=&LB(QPB)&\geq~{}LB(QPB^{-})&\geq~{}LB(YAT)~{}~{}~{}\text{
and }\\\\[8.53581pt]
UB(H)&=&UB(QPB)&\leq~{}UB(QPB^{-})&\leq~{}UB(YAT).\end{array}$

$\Box$

In fact, we do not generate all such inequalities, but add them one-by-one, as
needed. The following table with numerical results demonstrates that already
the addition of a few new inequalities provide strict improvements when
compared to (20).

| $UB(YAT)$ | $UB(QPB^{-})$ | $LB(YAT)$ | $LB(QPB^{-})$
---|---|---|---|---
n | Mean | Std | Max | Mean | Std | Max | Mean | Std | Max | Mean | Std | Max
8 | 0.15 | 0.15 | 0.71 | 0.10 | 0.13 | 0.65 | 1.31 | 1.63 | 6.83 | 0.62 | 1.04 | 3.28
9 | 0.32 | 0.23 | 1.12 | 0.31 | 0.23 | 1.08 | 1.51 | 1.31 | 3.77 | 1.26 | 1.29 | 3.73
10 | 0.57 | 0.29 | 1.24 | 0.55 | 0.28 | 1.22 | 1.57 | 1.59 | 5.44 | 1.57 | 1.59 | 5.44
11 | 0.91 | 0.29 | 1.47 | 0.58 | 0.32 | 1.21 | 1.14 | 1.00 | 3.88 | 1.09 | 1.01 | 3.88
12 | 1.29 | 0.5 | 2.45 | 1.17 | 0.48 | 2.35 | 1.02 | 0.86 | 3.28 | 0.99 | 0.89 | 3.28
13 | 1.58 | 0.33 | 2.33 | 1.12 | 0.50 | 1.85 | 1.83 | 0.99 | 3.31 | 1.44 | 1.16 | 3.20
14 | 1.76 | 0.42 | 2.62 | 1.54 | 0.39 | 2.28 | 1.35 | 1.52 | 5.05 | 1.28 | 1.55 | 5.05
15 | 2.47 | 0.50 | 3.21 | 1.81 | 0.58 | 3.03 | 1.58 | 1.51 | 4.23 | 1.47 | 1.43 | 3.93
Table 2: Mean, standard deviation, and maximum value of the percentage
relative error.

In Table 2 we compare the bounds provided by these models. For upper bounds,
we generate 100 instances from n from 8 to 10; for n from 11 to 15, we
generate 20 instances. For lower bounds, we generate 100 instances for n from
8 to 10; for n from 11 to 12, we generate 20 instances; for n from 13 to 15,
we generate 10 instances. Note that in many cases we get a strict improvement
by the addition of only a few members of $\mathcal{G}^{4}$ of the above type.

## References

  * [1] G. Boole.  Laws of thought.  American reprint of 1854 edition, Dover, 1854.
  * [2] G. Boole.  Of propositions numerically definite.  Transactions of Cambridge Philosophical Society, Part II, XI, 1868\.
  * [3] G. Boole.  Collected logical works, Vol I, Studies in Logic and Probability.  Open Court Publ. Co., 1952.
  * [4] E. Boros, Y. Crama, and P. L. Hammer.  Chvátal cuts and odd cycle inequalities in quadratic 0–1 optimization.  SIAM Journal on Discrete Mathematics, 5(2):163–177, 1992.
  * [5] E. Boros, Y. Crama, and P.L. Hammer.  Upper-bounds for quadratic 0–1 maximization.  Operations Research Letters, 9(2):73–79, 1990.
  * [6] Endre Boros and Peter L. Hammer.  Pseudo-boolean optimization.  Discrete Applied Mathematics, 123(1):155–225, 2002.
  * [7] Endre Boros and András Prékopa.  Closed form two-sided bounds for probabilities that at least r and exactly r out of n events occur.  Mathematics of Operations Research, 14(2):317–342, 1989.
  * [8] Endre Boros, Andrea Scozzari, Fabio Tardella, and Pierangela Veneziani.  Polynomially computable bounds for the probability of the union of events.  Mathematics of Operations Research, 39:1311–1329, 2014.
  * [9] K. Chung and P. Erdős.  On the application of the borel-cantelli lemma.  Transactions of the American Mathematical Society, 72:179–186, 1952\.
  * [10] D. A. Dawson and D. Sankoff.  An inequality for probabilities.  Proceedings of The American Mathematical Society, 18:504–504, 06 1967.
  * [11] D. DeCaen.  A lower bound on the probability of a union.  Discrete Math., 169:217–220, 1997.
  * [12] C. De Simone.  The cut polytope and the boolean quadric polytope.  Discrete Mathematics, 79:71–75, 1989-90.
  * [13] M. Deza and M. Laurent.  Facets for the cut cone i.  Math. Programming Ser. A, 56:121–160, 1992.
  * [14] M. Deza and M. Laurent.  Facets for the cut cone ii: clique-web inequalities.  Math. Programming Ser. A, 56:161–188, 1992.
  * [15] M. M. Deza and M. Laurent.  Geometry of Cuts and Metrics.  Springer-Verlag, Berlin, 1997.
  * [16] Janos Galambos.  Bonferroni Inequalities.  The Annals of Probability, 5(4):577 – 581, 1977.
  * [17] T. Hailperin.  Best possible inequalities for the probability of a logical function of events.  Amer. Math. Monthly, 72:343–359, 1965.
  * [18] F. Hoppe and E. Seneta.  Bonferroni-type inequalities and the methods of indicators and polynomials.  Advances in Applied Probability, 22:241–246, 1990.
  * [19] B. Jaumard, P. Hansen, and M. Poggi de Aragão.  Column generation methods for probabilistic logic.  J. Computing, 3:135–148, 1991.
  * [20] D. Kavvadias and C. H. Papadimitriou.  A linear programming approach to reasoning about probabilities.  Ann. Math. Artif. Intel., 1:189–205, 1990.
  * [21] H. Kuai, F. Alajaji, and G. Takahara.  A lower bound for the probability of a finite union of events.  Discrete Appl. Math., 215:147–158, 2000.
  * [22] H. Kuai, F. Alajaji, and G. Takahara.  A lower bound on the probability of a finite union of events.  Discrete Mathematics, 215(1):147–158, 2000.
  * [23] S. M. Kwerel.  Bounds on probability of a union and intersection of m events.  Adv. Appl. Probab., 7:431–448, 1975.
  * [24] S. M. Kwerel.  Most stingent bounds on aggregated probabilities of partially specified dependent probability systems.  J. Am. Statistical Association, 70:472–479, 1975.
  * [25] A. Prékopa.  Boole-bonferroni inequalities and linear programming.  Oper. Res., 36:145–162, 1988.
  * [26] András Prékopa.  Boole-bonferroni inequalities and linear programming.  Operations Research, 36(1):145–162, February 1988.
  * [27] András Prékopa and Linchun Gao.  Bounding the probability of the union of events by aggregation and disaggregation in linear programs.  Discrete Applied Mathematics, 145(3):444–454, 2005.
  * [28] J. Yang, F. Alajaji, and G. Takahara.  Linear programming bounds on the union probability.  Communications in Statistics - Simulation and Computation, 48:1–13, 10 2018.


# Learning Non-Stationary Time-Series with
Dynamic Pattern Extractions

Xipei Wang∗<EMAIL_ADDRESS>Haoyu Zhang
<EMAIL_ADDRESS>Yuanbo Zhang<EMAIL_ADDRESS>Meng
Wang<EMAIL_ADDRESS>Jiarui Song<EMAIL_ADDRESS>Tin
Lai<EMAIL_ADDRESS>Matloob Khushi∗ The authors are with The School of
Computer Science, The University of Sydney, Australia. ∗Corresponding author.
<EMAIL_ADDRESS>
###### Abstract

The era of information explosion had prompted the accumulation of a tremendous
amount of time-series data, including stationary and non-stationary time-
series data. State-of-the-art algorithms have achieved a decent performance in
dealing with stationary temporal data. However, traditional algorithms that
tackle stationary time-series do not apply to non-stationary series like Forex
trading. This paper investigates applicable models that can improve the
accuracy of forecasting future trends of non-stationary time-series sequences.
In particular, we focus on identifying potential models and investigate the
effects of recognizing patterns from historical data. We propose a combination
of the seq2seq model based on RNN, along with an attention mechanism and an
enriched set features extracted via dynamic time warping and zigzag peak
valley indicators. Customized loss functions and evaluating metrics have been
designed to focus more on the predicting sequence’s peaks and valley points.
Our results show that our model can predict 4-hour future trends with high
accuracy in the Forex dataset, which is crucial in realistic scenarios to
assist foreign exchange trading decision making. We further provide
evaluations of the effects of various loss functions, evaluation metrics,
model variants, and components on model performance.
Impact statement—We investigate pattern-matching techniques in non-stationary
time-series data, which closely resemble sensory and financial data in our
society. Neural networks learn to predict future trends based on patterns or
features that it discovers from historical data; therefore, pattern-matching
techniques for non-stationary series allow for a more expressiveness model. We
introduce three kinds of feature extraction techniques in this work, as a
generic framework that applies to any general time-series data. We focus our
attention on the Forex financial time series, and experimentally we illustrate
the best practical model architecture for financial data. We further
investigate the impact of loss functions that are tailored to financial series
to accelerate model convergence.

###### Index Terms:

Time-series, forecast, RNN, attention, seq2seq, GRU

##  I Introduction

Over the past decades, improvements in internet technology have led to the
explosion of information and massive data accumulation. Businesses,
technologies, transportation, media, studies, and many other sectors had
created a tremendous amount of data and became dependent on data-driven
approaches. However, in addition to extracting information from past data,
most businesses also want to forecast and infer future information. Through
the use of historical data, we can analyse and extract valuable knowledge,
which allows us to infer and predict future information. Extraction of
information will be a significant help for us to make further informed
decisions based on the gathered data.

Timeline based data, also known as time-series data, records the change of
attributes of certain events or objects over a period of time. Time series
data can be divided into stationary and non-stationary data. Heartbeat,
timestamps, train arrival, and departure can all be regarded as stationary
time-series data; their trends tend to be obvious and apparent. On the other
hand, non-stationary data refers to data that constantly change without
obvious regularity, such as stock, wind speed and rainfall. Non-stationary
data do not have any obvious trends, and they are often affected by numerous
internal or external factors beyond prediction. The major problem studied in
this paper is the question of _how to forecast future trends in non-stationary
time-series data_.

In the time-series literature, there are guidelines to classify whether a
time-series belongs to the stationary or non-stationary category. Previous
works [10, 1] point out that patterns extracted from historical time-series
data can often be used to improve the performance of future trend forecasting.
Bishop [2] further defines pattern recognition as the process of the automatic
regularity discovery in data by using algorithms.

In this paper, the proposing model can autonomously extract and recognises
existing patterns within the historical foreign exchange data. The extracted
patterns are then further utilised in a neural network model to forecast
future trends. Two main types of features are extracted from the input data:
(i) the similarity feature [18] that denotes the resemblance of the input data
against a pre-defined set of common patterns, and (ii) the zigzag indicators
[16] which extracts the peak and valley points from the input trend as an
expressive financial indicator. We focus on examining the predictability and
sensitivity of the proposing model with the additional autonomously discovered
features. This paper provides the basis for future research on investigating
the effectiveness of autonomous pattern extraction and pattern compositions in
the area of non-stationary time-series forecast.

Two major objectives are identified in our investigation. The first objective
is to structure methods that can recognise a series of typical patterns that
residue in the foreign exchange data and verify their corresponding
effectiveness with respect to the improvement towards forecasting future
trends. Two general categories of patterns, including similarity features and
zigzag peak valley indicators, as explained in detail in section IV-A1, will
be composed together or used separately. The second objective is to design a
forecasting model architecture, which can exploit the input non-stationary
time-series data stacked with the enriched features set. Our contributions are
as follows: (i) we propose a hybrid seq2seq model based on RNN with attention
mechanism and investigates its performance against other state-of-the-art
methods; (ii) we identify enrichment features in the non-stationary series and
propose automatic extraction of such features in historical data; and (iii) we
investigate the effect of custom loss functions that identify the peak and
valley points in the time-series trend. We report the best composition of
model components, such as loss functions, specific RNN variants, the best
combination of patterns and features that obtained the best performance to
forecast foreign exchange time series.

##  II Literature Review

A time-series consists of a sequence of observations of which are recorded at
some specific timestamps [4]. Time series can further be classified according
to their variability along time: if the mean and variance of a time series are
constants over time, it is identified as stationary; otherwise, the time
series is said to be non-stationary.

###  II-A Patterns in non-stationary time-series

Most time-series data consists of three main common patterns, including (i)
the trend which exists when there is a long-term increase or decrease, (ii)
the seasonality, which is a repetitive shape of the data and whose occurrence
is dependent on seasonal factors, such as a year, a week, an hour, etc., (iii)
and the cyclic pattern which refers to the phenomenon where the data rise and
fall without fixed frequencies [10]. For example, an important subcategory of
derived attributes in stock or foreign exchange data can be categorised by
chart patterns such as W, Flag, Wedge, ‘Head and Shoulder’, and ‘Cup and
Handle’, of which these patterns have long been utilised as indicators to
assist in making trading decisions.

There are various ways to extract patterns elaborated in related studies.
Spearman’s rank correlation, rule sets and sliding window are used in a real-
time hybrid pattern-matching algorithm [26]. The sliding window enables
pattern matching to be performed based on the latest received sub-sequence of
time-series data. Therefore, the proposed algorithm can be applied in real-
time application. Based on the sliding window, the rank correlation and rule
sets are used to distinguish patterns. In addition, the Gini function can be
used to determine each pattern’s beginning point and length [25]. Dynamic Time
Warping is another method that predicts the future trend of the financial
series by comparing the price sequence with the predefined patterns [12].
These results implies that we can often interpret the inherent patterns within
historical financial data, which would allow us extracts more information
about the pattern from the original data. If this information is input into
the model as additional features, it is possible to achieve better
performance.

###  II-B Time-series forecasting

Time-series forecasting is the procedure of predicting the future values of a
time-series according to the information extracted from its present and past
values.

####  II-B1 Auto-regressive Moving-Average (ARMA) models

ARMA models, first introduced by [23], have been proved to be effective in
short-term forecasting [9, 22]. An ARMA model, being the combination of an
Auto-Regressive (AR) model and a Moving-Average (MA) model, took advantage of
both of the two models by making use of the historical values and lagged
errors. However, since most models designed for stationary time-series
forecasting are based on the assumption of constant means and variance, they
cannot effectively capture the empirical features in time-series when time-
series become non-stationary, such as financial data [15]. Thus, more advanced
techniques for non-stationary time-series forecasting were introduced.

####  II-B2 ARMA Variants

An improved model had been developed to extend ARMA’s capability on processing
non-stationary data—AutoRegressive Integrated Moving Average (ARIMA) model,
which jointly utilises the ARMA model and differencing method. Differencing
can convert non-stationary values sequence to stationary residual sequence,
which is the strength of ARMA methods. Previous research has proved the
effectiveness of the ARIMA model in forecasting non-stationary time-series
such as farm price and oil price [14]. Another variant is called Generalized
AutoRegressive Conditional Heteroskedasticity (GARCH) model, which has shown
to be able to exploit the non-linear and non-stationary properties of time-
series data such as traffic load [27]. Although these methods had been applied
in specific domains, there exist critical limitations in applying to forecast
over-complicated financial data. For example, the ARIMA model is built upon
the assumption of constant variance, which does not conform to the high
volatility of most trading data [15].

####  II-B3 Recurrent Neural Network

Recurrent Neural Networks are suitable for handling temporal sequences due to
their nature of transferring data in-between cells among each layer, where the
order of the sequence can be memorized [17]. Studies have been done to compare
the performance of the ARIMA models and RNNs, and it is shown that RNNs are
able to provide higher satisfactory performances in comparison to the ARIMA
models [7]. Long Short-term Memory (LSTM), first proposed by [8], outperforms
vanilla RNN which has the problems of exploding gradient and vanishing
gradient [8, 19]. These authors prove that the exploding gradient problem can
cause the weights in the neurons being accumulated to an extremely large value
as to overflow, while the vanishing gradient problem results in the weights
being varied slowly; both of which will lead to the inability for the model to
continue to learn from the input data. LSTM, having a more complicated RNN
structure and more parameters than that of vanilla RNN, mitigates the above
problems by introducing three gates in each neuron that can decide what
information can be dropped and what information can be kept. Gated Recurrent
Units (GRU), introduced by [3], is relatively simpler than LSTM with less
complicated inner structure, fewer parameters and less computation. Both GRU
and LSTM work well in handling long sequence time-series data and there has
yet to be any consensus on which one is better among academics. Forecasting
future trends based on the past events of time-series data falls into the
problem of sequence to sequence prediction, which can be handled by a seq2seq
model [21]. The seq2seq model contains an encoder to obtain the input sequence
and a decoder to output the forecasting sequence. Both the encoder and decoder
can adopt the aforementioned recurrent neural networks to handle the temporal
sequences.

Although the seq2seq model based on RNNs has already achieved a decent
performance as evidential by numerous empirical studies, it remains to exhibit
the vanishing gradient problem in handling long sequences. The attention
mechanism is necessary to further mitigate this shortcoming. The seq2seq model
backed with an attention mechanism outweighs the vanilla seq2seq model which
follows the encoder-decoder paradigm in handling the long time-series sequence
[20, 11, 24]. The attention is a method that, during the decoding processes,
allows the decoder to directly connect to the encoder to focus on the most
relevant features from the input signals [13, 24].

##  III Forex Dataset

###  III-A Data Collection

The non-stationary time-series datasets used in this research are the ever-
changing price data of the Forex market over several years. These datasets
include attributes that describe the datetime, open, close, high, low price,
and volume of each time interval of a certain currency pair. All the datasets
are obtained with fxcmpy [6].

###  III-B Data Analysis

The data obtained from Forex are collected at different time intervals: 5, 15,
30 minutes and 1, 2, 4 hours intervals. For the practical reason of proposing
a realistic model that can forecast relatively long future trends, e.g.
several hours or one day, we had chosen 15-minute-interval data with 4 hours
as our forecast length (equivalent to 16 data points) of future trends. The
reason is to balance the trade-off between predicting an overly long
sequence—which will over complicates the problem for the model to handle, and
overly short sequence—which would under-represent the upcoming trend.

The collected data, as specified in Section III-A, consists of date, time,
various prices and volume. Within these attributes, the time information is
uninformative; the open price, which is the close price of the previous time
slot, is redundant information that can be disregarded. The close price, being
the final state of the current time point, is chosen as the primary input of
the model.

The data amount is sufficiently large since the data dates back to the 1990s,
which contains many representative features for the model to learn. There are
also data of various currency pairs, e.g. USD to JPY and GBP to AUD. In this
study, we focus on the currency pair of USD to JPY within the time frame of
2015 to 2019 is chosen for evaluation purposes, which contain over 120,000
records in total.

###  III-C Data Processing

The input features are first scaled to a common magnitude such that it enables
the model to utilise the learned features across different currency pairs with
varying scales. Such an operation enable us to utilise the same underlying
patterns across currency-pairs that have different relative scale. Then, the
data are split into a training set and a validation set with a rate of 9:1 for
cross-validation. In addition, the data from the training and validation set
will be cut by shifting windows into input and output sequence pairs for the
seq2seq model. The overall mechanism of data processing is illustrated in
Figure 1.

Figure 1: Data pre-processing

##  IV Methodology

###  IV-A Model Architecture

The overall model architecture is shown in Figure 2.

Figure 2: The overall model architecture

Firstly, the series of close price $D=[d_{0},d_{1},...,d_{T-1}]$ are extracted
from the original raw data where $T$ (over 120,000) is the total number of
data points. Each series is sliced as a sliding window with model input $X$,
where $|X|=672$, is the input feature in a 7 days period. The output sequence
$y$, where $|y|=16$, is the prediction size with a 4 hours period. Each
prediction takes a step size of 16. Then, we begin to extract the similarity
pattern and zigzag peak valley features based on each input sequence.
Subsequently, these features are stacked onto the input sequence. Finally, the
stacked data and features are input to the seq2seq model, where we predict a
sequence from the model and compare it against the ground truth $y$. We will
begin to detail the feature extraction and the structure of the seq2seq model
in the following sections.

####  IV-A1 Feature Extraction

This paper proposes two different approaches in pattern matching and feature
extraction for time series data. The first one is to directly match different
segments of the time series data against a predefined set of representative
patterns, then subsequently computes similarities score as an additional
feature. Hereafter we will refer to this as similarity feature. Since the
proposing model needs to predict the highest and lowest price point in a short
horizon to aid trading decisions, the second approach is to recognize the
oscillation in time series by marking out the peaks (local maxima) and valley
(local minima). Subsequently, we can use this peak-valley series as an
additional feature, which is we will refer to as zigzag feature.

Figure 3: Diagram of feature extraction process

Figure 3 illustrates the overview of the feature extraction process. For each
data pair that contains the input sequence and target sequence, its input
sequence will be divided into several sub-windows with equal length. Then,
each sub-window is matched against our pre-defined set of representative
patterns that are common in the financial market [12], which produces a
distance vector of length 13 that measures the similarities between each sub-
window and the patterns (Figure 4). The minimum value in the vector represents
the most similar representative pattern, and the vector is further converted
to a one-hot encoding form based on the minimum value. The same process is
repeated for each pattern to form a similarity feature matrix for the
corresponding input sequence. As shown in Figure 4, these representative
patterns can be used to represent the basic trend in the trading chart. Thus,
the various combination of these representative patterns can effectively
represent most of the common chart patterns in Forex or the stock market.

Figure 4: Structure of the 13 representative patterns [12]

As shown in Figure 3, several groups of sub-windows is used for similarity
feature generation, where each group consists of different window sizes. The
purpose of this strategy is to extract features with different granularity;
thus, features of various lengths can be recognized and extracted. To solve
the problem of distance calculation for sub-windows of different lengths,
Dynamic Time Warping (DTW) is introduced. DTW is a well-known algorithm for
matching two time series of different lengths. In this paper, the Euclidean
distance-based DTW algorithm is used to calculate the distance between each
sub-window and each representative pattern.

The zigzag feature was inspired by the Zigzag indicator commonly used in
financial analysis. As shown by an example in Figure 5, the Zigzag indicator
can identify the peak and valley price values in a given price sequence,
commonly used by traders for making trading decisions based on the generalised
trend. Since the distribution of these peak and valley points can provide
important information about future price trend, this paper uses the Zigzag
indicator to generate additional features that identify each point in the
input sequence as a peak, valley, or other value. Following the same procedure
as the similarity feature, the zigzag feature is also one-hot encoded.

As shown in Figure 3, in model training or prediction, the similarity features
and zigzag features are first computed from the input sequences in each data
pair, and then the computed features are stacked with the input sequences to
form an input matrix that serves as the input to the model.

Figure 5: Example of Zigzag indicator

####  IV-A2 Seq2seq Model

The aforementioned input features will be input to the seq2seq model to
forecast the incoming trends. Figure 6 shows the model structure, which
contains an encoder and a decoder. The encoder, which takes in the input
matrix, is consists of bi-directional recurrent neural network layers. We
parameterise this framework into a set of different neural networks, such as
RNN, LSTM or GRU, and with different hidden sizes and layers.

The decoder, which processes the intermediate encoder output, generates the
final output sequence. The decoder consists of multiple layers, including an
attention layer to capture the information from the encoder to enhance the
ability to handle the long sequences. Other layers include multiple
parameterized core bi-directional recurrent neural network layers, Rectified
Linear Unit (ReLU) layers, fully connected layers and $sigmoid$ activation
layers. There is also an optional dropout layer in the encoder and decoder
respectively, which is omitted in Figure 6, to address the problem of over-
fitting.

$y_{n}$…$y_{2}$$y_{1}$SigmoidFully ConnectedFully ConnectedReLUBi-
RNN/LSTM/GRUAttentionBi-RNN/LSTM/GRU$x_{2}$$x_{1}$$x_{n}$…EncoderDecoder
Figure 6: Seq2seq model structure

####  IV-A3 Loss Functions

After the model presented in this paper was designed, some initial experiments
were performed for predicting future price series of length 96 using an input
sequence of length 692. However, we found that the training process was
relatively slow, and the model’s predictions show that it can only predict the
general future price direction, but not the precise price fluctuations. In
order to solve the problem of slow training, we introduce three new customized
loss functions to force the training process to recognise fluctuations in the
target sequence.

The first customized loss function is called Single Peak Valley (SPV) loss. It
adds two penalty terms for the horizontal coordinate difference between peak
and valley values in the prediction and target sequences, based on the regular
RMSE loss. SPV loss is given by

$\begin{array}[]{l}\mathcal{L}_{\textsc{spv}}=\mathcal{L}_{\textsc{rmse}}\times(1+\alpha\times
pd+\beta\times vd),\end{array}$ (1)

where $\mathcal{L}_{\textsc{rmse}}$ is the typical RMSE loss, $pd$ is the
horizontal coordinator difference between the highest point in the predicted
sequence and the target sequence, $vd$ is the horizontal coordinator
difference between the lowest points in the predicted sequence and the target
sequence. Both $\alpha$ and $\beta$ are coefficients, which is initially set
as 0.5.

The second customized loss function is called Multi Peak Valley (MPV) loss .
It is a variant of $\mathcal{L}_{\textsc{spv}}$ that adds penalty terms for
more than one pair of peak and valley points in the target sequence. MPV loss
is given by

$\begin{array}[]{l}\mathcal{L}_{\textsc{mpv}}=\mathcal{L}_{\textsc{rmse}}\times(1+\sum_{i=1}^{k}{(\alpha_{i}\times
pd_{i}+\beta_{i}\times vd_{i})}),\end{array}$ (2)

where $pd_{i}$ is the horizontal coordinator difference between the $i^{th}$
peak points in the predicted sequence and the target sequence, $vd_{i}$ is the
horizontal coordinator difference between the $ith$ valley points in the
predicted sequence and the target sequence. The peak and valley points of the
predicted sequence and the target sequence are selected by using Zigzag
indicator with the same parameter setting. By focusing on more peaks and
troughs, $\mathcal{L}_{\textsc{mpv}}$ is expected to provide more accurate
feedback to the model, thus further accelerating the training process and
helping the model to accurately identify fluctuations in the target sequence.

The last customized loss function is called Weighted RMSE (WRMSE). It differs
from $\mathcal{L}_{\textsc{spv}}$ and $\mathcal{L}_{\textsc{mpv}}$ which add
penalty terms based on horizontal gap between peak and valley point. WRMSE
adds penalty term for each point in the target sequece based on its vertical
distance from the average of the target sequence, given by

$\mathcal{L}_{\textsc{wrmse}}=\sqrt{\frac{\sum_{i=1}^{n}{\left(\hat{y}_{i}-y_{i}\right)^{2}\left(y_{i}-\bar{y}\right)^{2}}}{n}},$
(3)

where $\hat{y}_{i}$ is the $ith$ point in predicted sequence and $y_{i}$ is
the $ith$ point in the target sequence, and $\bar{y}$ is the mean of the
target sequence. The purpose of $\mathcal{L}_{\textsc{spv}}$ and
$\mathcal{L}_{\textsc{mpv}}$ is to more accurately match the horizontal
coordinates of peak and valley points, while the purpose of
$\mathcal{L}_{\textsc{wrmse}}$ is to more rapidly match the predicted sequence
to the vertical value of each point in the target sequence. As such, the
predicted sequence can be more accurately matched to the value of the
corresponding position in the target sequence.

###  IV-B Evaluation Methods

The symmetric mean absolute percentage error (SMAPE) will be used as one of
evaluation metric. Given that this study’s objective is to identify market
price fluctuations, especially the distribution of peak and valley prices in a
short horizon, two extra evaluation metrics are designed based on this
objective.

####  IV-B1 Peak Valley RMSE (PVRMSE)

In order to calculate _PVRMSE_ , the Zigzag indicator is used to identify the
peak and valley points in the predicted and target sequences. Then, the square
errors of corresponding peak and valley points of these two sequences are
calculated according to their occurrence. Let $\mathbf{\hat{p}}$ and
$\mathbf{p}$ denote the predicted and target peak point sequence,
$\mathbf{\hat{v}}$ and $\mathbf{v}$ denote the predicted and target valley
point sequence respectively, where
$\mathbf{\hat{p}}=(\hat{p}_{1},\ldots,\hat{p}_{i},\ldots)$ and vice versa. Let
$M$ and $N$ denote the minimum values among the size of the predicted and
target sequence for $\mathbf{p}$ and $\mathbf{v}$ respectively, i.e.,
$M=min(|\mathbf{\hat{p}}|,|\mathbf{p}|)$ and
$N=min(|\mathbf{\hat{v}}|,|\mathbf{v}|)$ where $|\cdot|$ is the cardinality
operator. Then, PVRMSE is given by

$\textit{PVRMSE}=\sqrt{\frac{\sum_{i=1}^{M}{\left(\hat{p}_{i}-p_{i}\right)^{2}+\sum_{j=1}^{N}{\left(\hat{v}_{j}-v_{j}\right)^{2}}}}{M+N}},$
(4)

which gives us a measure of the errors in predictions among the peak and
valley points.

####  IV-B2 Peak Valley MAE (PVMAE)

Similarly, another commonly used evaluation metric, Mean Absolute Error (MAE),
has been modified for peak valley point matching. Peak Valley MAE (PVMAE) is
given by

$\textit{PVMAE}=\frac{\sum_{i=1}^{M}{\lVert\hat{p}_{i}-p_{i}\rVert+\sum_{j=1}^{N}{\lVert\hat{v}_{j}-v_{j}\rVert}}}{M+N},$
(5)

where the parameters maintain the same mathematical meaning as in PVRMSE.

###  IV-C Experiments Setup

Three major experiments are carried out in this research, including (i)
finding which loss functions are better for accelerating the convergence of
the seq2seq model, (ii) investigating which composition of pattern features
can help the model obtaining the best performance, and (iii) comparing which
RNN variant used in the seq2seq model is the most suitable for the problem of
Forex time-series forecasting and whether the attention mechanism is able to
improve forecasting performance and comparing our model with a baseline ANN
model and the traditional ARIMA model. In the following subsections, the
common settings shared by our seq2seq models with different components, and
the respective parameters for each model variant will be introduced in detail.

####  IV-C1 Common Settings

The settings about the data, shared by all the seq2seq model variants, is
described in Table I.

TABLE I: Common settings

| Parameters
---|---
Data Source[1] | USD to JPY
Time range | 2015-2019[2]
Time Interval | 15 minutes[3]
Input sequence length | 672[4]
Output sequence length | 16[5]
Number of epochs | 150
Learning rate | 1e-4
Number of layers | Encoder: 1
Decoder: 1
Hidden size | 128
Batch size | 128
Dropout rate | 0.0
Teacher forcing ratio[6] | 0.0
Slacked pattern window sizes | 672, 336, 96, 48, 24, 12
Slacked zigzag differences[7] | 0.0063, 0.007, 0.008, 0.0097,
0.012, 0.015, 0.0163, 0.0288

  * •

Notes:

  * •

[1] Forex data [2] 125000 records [3] 4 points/hour [4] 7days [5] 4 hours

  * •

[6] In the decoder, the ratio of feeding the ground truth to the neurons,
otherwise the output of previous point

  * •

[7] Generate average around 3, 6, …, 24 peaks and valleys on the 672-length
sequence

####  IV-C2 Settings for loss function comparison

TABLE II: Loss function parameters Loss | Parameters
---|---
$\mathcal{L}_{\textsc{spv}}$ | $\alpha=\beta=0.5$
$\mathcal{L}_{\textsc{mpv}}$ | $k=3$
$\alpha_{1}=\beta_{1}=0.3$
$\alpha_{2}=\beta_{2}=0.15$
$\alpha_{3}=\beta_{3}=0.05$

To compare the loss functions described in section IV-A3, the model is trained
150 epochs with the same parameter settings in Table I. Certain parameters in
relation to the certain loss functions are listed in Table II.

####  IV-C3 Settings for feature comparison

These model settings refers to our experiments that focus on the generation,
stacking of different granularities of pattern features, and other common
parameters are illustrated in Table I. We study different combinations of
features, including:  1. close price,  2. close price and zigzag,  3. close
price and similarity patterns,  4. close price, zigzag and similarity
patterns.

####  IV-C4 Settings for model component comparison

In this experiment, common parameters in Table I are used. The model is
trained 150 epochs for each RNN variant with or without attention layer, the
compositions are listed below:

  1. 1.

Vallina RNN,

  2. 2.

Vallina RNN and Attention,

  3. 3.

LSTM,

  4. 4.

LSTM and Attention,

  5. 5.

GRU,

  6. 6.

GRU and Attention.

##  V Results

There are three types of experiments conducted, including the comparison of
(i) loss functions, (ii) feature compositions, and (iii) model components and
baseline ANN and ARIMA. By comparing the specifically designed loss functions,
the most appropriate one can be identified and used to accelerate the whole
training process with improved performance.

In addition, the comparison of different compositions of pattern features
allows us to identify the kind of extracted patterns that are useful for time-
series forecasting. Furthermore, by comparing core model components (e.g. RNN
variants and attention mechanism) the optimised internal structure of the
model can be determined for future research.

Additionally, the _PVRMSE_ and _PVMAE_ evaluation metrics (see Section IV-B),
focusing on the most important peak and valley points in the predicting
sequence, are adopted among all the experiments to evaluate the seq2seq model.

###  V-A Loss Function Comparison

This section outlines the loss function that can most effectively helps to
accelerate the model training process. We had compared four different loss
functions—the $\mathcal{L}_{\textsc{rmse}}$, $\mathcal{L}_{\textsc{wrmse}}$,
$\mathcal{L}_{\textsc{spv}}$ and $\mathcal{L}_{\textsc{mpv}}$. Figure 7,
Figure 8 and Table III illustrate that, if trained with the same number of
epochs, the model using $\mathcal{L}_{\textsc{mpv}}$ can provide a prediction
with $\emph{PVRMSE}=4.07\times 10^{-3}$ and $\emph{PVMAE}=2.60\times 10^{-3}$
which outweighs the results of using three other loss functions. It means that
the $\mathcal{L}_{\textsc{mpv}}$ is able to speed up the model training by
providing more specific information in terms of the peak and valley points in
the target sequence. Therefore, this loss function is adopted to be the one
that is used in other experiments and in the process of training the final
model.

Figure 7: Validation loss comparison among different loss functions (PVRMSE)

Figure 8: Validation loss comparison among different loss functions (PVMAE) TABLE III: The comparison of different features Loss functions | PVRMSE | PVMAE | SMAPE
---|---|---|---
|  ($\times 10^{-3}$) |  ($\times 10^{-3}$) |
RMSE ($\mathcal{L}_{\textsc{rmse}}$) | 5.47 | 4.81 | 2.22
WRMSE ($\mathcal{L}_{\textsc{wrmse}}$) | 5.69 | 5.1 | 2.4
Single PeakValley ($\mathcal{L}_{\textsc{spv}}$) | 5.19 | 4.66 | 2.01
Multi PeakValley ($\mathcal{L}_{\textsc{mpv}}$) | 4.07 | 2.60 | 1.95

Additionally, from the result, it can be seen that the performance of the
model with $\mathcal{L}_{\textsc{spv}}$ and $\mathcal{L}_{\textsc{mpv}}$ is
better than that with the other two loss functions. Furthermore, the
performance of $\mathcal{L}_{\textsc{mpv}}$ with the penalty term added for
more horizontal coordinates of peak and valley points are better than that of
$\mathcal{L}_{\textsc{spv}}$ with the penalty term added for only one pair of
peak and valley points. These findings confirm the assumption of adding
suitable penalty terms allow the model to pay more attention to the peak and
valley points of the target sequence during the training process. Thus,
producing a prediction sequence that is more consistent with the ground truth.

Finally, the result shows that $\mathcal{L}_{\textsc{wrmse}}$ achieves the
worse performance out of all the other three loss functions. The direct
comparison of numeric values might not be fair comparison because the computed
loss values might not have the same magnitude as that of the others. Since all
hyperparameters provided to the model are the same, it might indicates that
the model might require further tuning, such as learning rate, for the
$\mathcal{L}_{\textsc{wrmse}}$ loss function. Since it is potential to be the
one has the ability outweighs the others for the problem of forecasting time-
series due to its specific algorithm design (See section IV-A3), future
research can focus on doing more experiments about this loss function by using
a different learning rate.

###  V-B Feature Comparison

The second experiment is to identify which features need to be extracted and
stacked to the input sequence of each window. By comparing the model’s
performance when using combinations of different types of features as model
inputs, the experimental results in this section answer the research question
of _how to identify and extract the patterns that are implicitly embedded in
the Forex trading data_ , which allow us to utilise them via enriching the
model’s input.

As indicated by Table IV, the additional Zigzag features outperform other
features across all metrics, whereas the additional DTW similarity pattern
features had degraded the performance of our baseline model. s In addition,
the result also indicates that combining Zigzag and DTW features might
requires more training time as the dimensionality of the input sequence are
much higher than that of the baseline model.

Figure 9: Validation loss comparison among different features (PVRMSE)

Figure 10: Validation loss comparison among different features (PVMAE) TABLE IV: The comparison of different features Features | PVRMSE | PVMAE | SMAPE
---|---|---|---
|  ($\times 10^{-3}$) |  ($\times 10^{-3}$) |
None | 5.19 | 4.66 | 2.01
Zigzag | 4.96 | 4.45 | 1.63
Similarity Feature | 5.62 | 5.03 | 1.74
Zigzag & Similarity | 6.29 | 5.65 | 2.33

###  V-C Model Component Comparison

The final experiment is about varying the core components of the seq2seq
model, which includes 6 different compositions of Vallina RNN, LSTM, GRU and
attention mechanism. The result, as shown in Figure 11, Figure 12 and Table V,
suggests that attention-based Seq2Seq model with GRU has the highest
performance (PVRMSE: $4.96\times 10^{-3}$, PVMAE: $4.45\times 10^{-3}$) than
the other 5 compositions.This is inline with our initial hypothesis since GRU
is the more advanced version in the RNN family, and attention mechanism can
improve the neural network dealing with a long sequence.

The Diebold Mariano test [5] has been performed to compare our GRU attention
based seq2seq model with a traditional ARIMA, and an ANN model. The ARIMA
model consists of the settings of $p=0$, $d=1$, $q=0$, which refer to the
auto-regressive, differencing, and moving average terms respectively. The ANN
model is composed of 3 dense layers with the dimensionality of (1) input:672,
output:128, (2) input:128, output: 32, (3) input:32, output:16. We use Diebold
Mariano test function with the settings of MSE differential-loss,
$h=|y|^{1/3}+1=16^{1/3}+1=3$, where $h$ is the number of steps ahead of the
prediction, and $y$ is the predicting sequence. As shown in Table VI, the GRU
version achieves lower relative error when compared to ANN. The ARIMA model
has the highest relative error in all of our paired Diebold Mariano test when
compared to the ground truth, with $0.01$ significant level.

Figure 11: Validation loss comparison among different NN models (PVRMSE)

Figure 12: Validation loss comparison among different NN models (PVMAE) TABLE V: The comparison of different model components Model | Attention | PVRMSE | PVMAE | SMAPE
---|---|---|---|---
|  |  ($\times 10^{-3}$) |  ($\times 10^{-3}$) |
Vallina RNN | No | 5.3 | 4.74 | 1.99
Yes | 5.72 | 5.14 | 1.83
LSTM | No | 5.35 | 4.85 | 1.68
Yes | 5.76 | 5.32 | 1.59
GRU | No | 5.35 | 4.85 | 1.63
Yes | 4.96 | 4.45 | 1.63
ARIMA | N/A | 492.28 | 492.23 | 199.97
ANN | N/A | 12.04 | 11.47 | 2.86
TABLE VI: Diebold-Mariano loss-differential test

Model |  GRU+Attention |  ARIMA |  ANN
---|---|---|---
GRU+Attention |  - |  123.93∗ |  0.53
ARIMA |  -123.93∗ |  - |  -120.88∗
ANN |  -0.53 |  120.88∗ |  -

  * •

Notes: Negative sign of the statistics implies the second model (row) has
bigger forecasting errors. ∗Significant at the $0.01$ level.

###  V-D Trend Forecasting Result

Based on the above comparison results, the combination of
$\mathcal{L}_{\textsc{mpv}}$ as the loss function, attention-based Seq2Seq
model with GRU, and the input feature with Zigzag feature is selected as the
optimal model. It is trained for 500 epochs and achieves a performance with
the $\textit{PVRMSE}=4.34\times 10^{-3}$ and $\textit{PVMAE}=3.8\times
10^{-3}$. Figure 13 shows a sample forecasting result of the optimal model. It
indicates that the prediction represented by the green line has a high
similarity compared to the ground truth represented by the blue line in both
scale and shape, which indicates that this model has the ability in
forecasting a relative longer future sequence.

Figure 13: A sample forecasting result of the model

##  VI Discussion and Conclusion

In this research, two types of patterns, including the zigzag peak valley
indicator and the time-series change patterns, are recognized, extracted,
stacked on the 7-day normalized close price data, and then concatenated as the
input for an optimized attention and GRU-based seq2seq model to forecast the
price trends over the next 4 hours. Several comparative experiments are
conducted to select the optimal combination of loss function, feature type,
and model structure. In addition, three customized loss functions are designed
to help speed up the model training process. Finally, in order to meet the
objective of identifying the peak and valley points in the future price curve,
two peak and valley-based evaluation metrics are designed to evaluate model
performance.

Based on the results (refer to section V), the two research objectives
discussed in section I are achieved. In terms of the pattern, the aim is to
verify if adding patterns can improve the forecasting performance or to
identify the approaches to add these patterns. Results in section V-B suggest
that adding patterns of zigzag indicator has a positive effect for improving
the forecast performance. On the other hand, the negative results from the DTW
similarity feature indicates that such a feature set should be utilised in the
model with a different approach, e.g. stacking with other indicators, using in
different granularities, or stacking multiple layers with varied combination
of granularities.

In addition, practically the peak valley prices are the most important points
for financial analysts to make trade decisions. We had persuaded this
direction by designing two methods for our investigation. One method is to use
extra specific evaluation metrics, including PVRMSE, and PVMAE (section IV-B),
which focuses on peaks and valleys during the model evaluation. The other one
is to use customized loss functions, including $\mathcal{L}_{\textsc{spv}}$,
$\mathcal{L}_{\textsc{mpv}}$ and $\mathcal{L}_{\textsc{wrmse}}$ (section
IV-A3), to accelerate the model converging speed by explicitly encouraging the
model learning towards the direction of finding peak and valley points. The
results in section V-A indicate that the loss functions designed in this
research do help for speeding up the model learning and improving the model
performance. These two methods are recommended to be used by research with the
same purpose.

Finally, the result in section V-C illustrates that the combination of GRU
with the attention outperforms other combinations of RNN variants and
attention being used in our seq2seq model in solving the problem of time-
series trend forecasting. Its performance also exceeds that of the traditional
ARIMA, the benchmark ANN model. Thus, it illustrates that GRU with attention
mechanism would be the preferred solution in future studies.

##  VII Future Works

Although there are many achievements in this research as described in section
VI, due to the limitations of time and computational resources, we did not
investigate every aspects on this line of work. Future works can focus on the
following aspects:

  * •

From the financial perspective, find more features that can be extracted to
enrich the price data, and investigate their effects on model performance.

  * •

For the zigzag and similarity pattern features, perform more experiments on
extracting them with different granularities and test their effects on the
model with different set of combinations.

  * •

Convolutional Neural Networks with 1-dimensional convolutional layers can be
adopted to examine if they can automatically extract features, and their
effort on model performance.

  * •

Two loss functions, including the SPV loss $\mathcal{L}_{\textsc{spv}}$ and
MPV loss $\mathcal{L}_{\textsc{mpv}}$, are parameterized with the allowance of
changing the weight coefficients of the error on peak and valley points. More
experiments can be carried out to investigate the set of parameters for the
loss functions that would allow us to achieve the highest performance .

  * •

The prblem with overfitting tends to emerge when we continued to train the
model with more epochs. The model can be modified to add more measures to
mitigate this problem. Another measure that can be done is to train the model
with more data.

  * •

In this research, we only focus on using a fixed interval for the training and
testing data. Future studies can expand this to the data with different
intervals and other currency pairs. Transfer learning can also be carried out
to verify if the model can be improved by training different data.

## Acknowledgment

Many thanks to Dr. Matloob Khushi and Dr. Tin Lai for their valuable guidance
on technical and financial knowledge.

## References

  * Alvarez et al. [2010] Francisco Martinez Alvarez, Alicia Troncoso, Jose C Riquelme, and Jesus S Aguilar Ruiz.  Energy time series forecasting based on pattern sequence similarity.  _IEEE Transactions on Knowledge and Data Engineering_ , 23(8):1230–1243, 2010.
  * Bishop [2006] Christopher M Bishop.  _Pattern recognition and machine learning_.  springer, 2006.
  * Cho et al. [2014] Kyunghyun Cho, Bart Van Merriënboer, Dzmitry Bahdanau, and Yoshua Bengio.  On the properties of neural machine translation: Encoder-decoder approaches.  _arXiv preprint arXiv:1409.1259_ , 2014.
  * Das [1992] S. Das.  Time series analysis.  1992\.
  * Diebold and Mariano [1995] FX Diebold and RS Mariano.  Comparing predictive accuracy’, journal of business and economic statistics, 13 (3), july, 253-263.  _INTERNATIONAL LIBRARY OF CRITICAL WRITINGS IN ECONOMICS_ , 108:253–263, 1995.
  * FXCM [2020] FXCM.  fxcmpy – algorithmic trading with fxcm, 2020.  URL https://www.fxcm.com/fxcmpy/.
  * Ho et al. [2002] Sui-Lau Ho, Min Xie, and Thong Ngee Goh.  A comparative study of neural network and box-jenkins arima modeling in time series prediction.  _Computers & Industrial Engineering_, 42(2-4):371–375, 2002.
  * Hochreiter and Schmidhuber [1997] Sepp Hochreiter and Jürgen Schmidhuber.  Long short-term memory.  _Neural computation_ , 9(8):1735–1780, 1997.
  * Huang and Shih [2003] Shyh-Jier Huang and Kuang-Rong Shih.  Short-term load forecasting via arma model identification including non-gaussian process considerations.  _IEEE Transactions on power systems_ , 18(2):673–679, 2003.
  * Hyndman and Athanasopoulos [2018] Rob J Hyndman and George Athanasopoulos.  _Forecasting: principles and practice_.  OTexts, 2018.
  * Kang et al. [2018] Lei Kang, J Ignacio Toledo, Pau Riba, Mauricio Villegas, Alicia Fornés, and Marçal Rusinol.  Convolve, attend and spell: An attention-based sequence-to-sequence model for handwritten word recognition.  In _German Conference on Pattern Recognition_ , pages 459–472. Springer, 2018.
  * Kim et al. [2018] Sang Hyuk Kim, Hee Soo Lee, Han Jun Ko, Seung Hwan Jeong, Hyun Woo Byun, and Kyong Joo Oh.  Pattern matching trading system based on the dynamic time warping algorithm.  _Sustainability_ , 10(12):4641, 2018.
  * Michael et al. [2019] Johannes Michael, Roger Labahn, Tobias Grüning, and Jochen Zöllner.  Evaluating sequence-to-sequence models for handwritten text recognition.  In _2019 International Conference on Document Analysis and Recognition (ICDAR)_ , pages 1286–1293. IEEE, 2019.
  * Nochai and Nochai [2006] Rangsan Nochai and Titida Nochai.  Arima model for forecasting oil palm price.  In _Proceedings of the 2nd IMT-GT Regional Conference on Mathematics, Statistics and applications_ , pages 13–15, 2006.
  * Petrică et al. [2016] Andreea-Cristina Petrică, Stelian Stancu, and Alexandru Tindeche.  Limitation of arima models in financial and monetary economics.  _Theoretical & Applied Economics_, 23(4), 2016.
  * Qi and Khushi [2020] Ling Qi and Matloob Khushi.  Event-driven lstm for forex price prediction.  IEEE, 2020.
  * Rumelhart et al. [1986] David E Rumelhart, Geoffrey E Hinton, and Ronald J Williams.  Learning representations by back-propagating errors.  _nature_ , 323(6088):533–536, 1986.
  * Shen et al. [2017] Jingyi Shen, Weiping Huang, Dongyang Zhu, and Jun Liang.  A novel similarity measure model for multivariate time series based on lmnn and dtw.  _Neural Processing Letters_ , 45(3):925–937, 2017\.
  * Sherstinsky [2020] Alex Sherstinsky.  Fundamentals of recurrent neural network (rnn) and long short-term memory (lstm) network.  _Physica D: Nonlinear Phenomena_ , 404:132306, 2020.
  * Siddiqui and Shamsi [2018] Tooba Siddiqui and Jawwad Ahmed Shamsi.  Generating abstractive summaries using sequence to sequence attention model.  pages 212–217. IEEE, 2018.
  * Sutskever et al. [2014] Ilya Sutskever, Oriol Vinyals, and Quoc V Le.  Sequence to sequence learning with neural networks.  In _Advances in neural information processing systems_ , pages 3104–3112, 2014.
  * Torres et al. [2005] Jose Luis Torres, Almudena Garcia, Marian De Blas, and Adolfo De Francisco.  Forecast of hourly average wind speed with arma models in navarre (spain).  _Solar energy_ , 79(1):65–77, 2005.
  * Whittle [1951] Peter Whittle.  _Hypothesis testing in time series analysis_ , volume 4.  Almqvist & Wiksells boktr., 1951.
  * Yang et al. [2016] Zichao Yang, Xiaodong He, Jianfeng Gao, Li Deng, and Alex Smola.  Stacked attention networks for image question answering.  pages 21–29, 2016.
  * Zhang et al. [2009] Xiaohang Zhang, Jun Wu, Xuecheng Yang, Haiying Ou, and Tingjie Lv.  A novel pattern extraction method for time series classification.  _Optimization and Engineering_ , 10(2):253–271, 2009.
  * Zhang et al. [2010] Zhe Zhang, Jian Jiang, Xiaoyan Liu, Ricky Lau, Huaiqing Wang, and Rui Zhang.  A real time hybrid pattern matching scheme for stock time series.  In _Proceedings of the Twenty-First Australasian Conference on Database Technologies-Volume 104_ , pages 161–170, 2010.
  * Zhou et al. [2005] Bo Zhou, Dan He, Zhili Sun, and Wee Hock Ng.  Network traffic modeling and prediction with arima/garch.  In _Proc. of HET-NETs Conference_ , pages 1–10, 2005.