after every new SDE is feasible in principle, but not very useful and can also
produce results that are misleading. By their very nature, CEs are relatively
rare within a stream of input SDEs. As a result, if we emit a forecast after
every new SDE, some of these forecasts (possibly even the vast majority) will
have a significant temporal distance from the CE to which they refer. As an
example, consider a pattern from the maritime domain which detects the
entrance of a vessel in the port of Tangiers. We can also try to use this
pattern for forecasting, with the goal of predicting when the vessel will
arrive at the port of Tangiers. However, the majority of the vessel’s messages
may lie in areas so distant from the port (e.g., in the Pacific ocean) that it
would be practically useless to emit forecasts when the vessel is in these
areas. Moreover, if we do emit forecasts from these distant areas, the scores
and metrics that we use to evaluate the quality of the forecasts will be
dominated by these, necessarily low-quality, distant forecasts.

For these reasons, before running a regression experiment, we must first go
through a preprocessing step. We must find the timepoints within a stream
where it is “meaningful” to emit forecasts. We call these timepoints the
_checkpoints_ of the stream. To do so, we must first perform recognition on
the stream to find the timepoints where CEs are detected. We then set a
required distance $d$ that we want to separate a forecast from its CE, in the
sense that we require each forecast to be emitted $d$ events before the CE.
After finding all the CEs in a stream and setting a value for $d$, we set as
checkpoints all the SDEs which occur $d$ events before the CEs. This typically
means that we end up with as many checkpoints as CEs for a given value of $d$
(unless the distance between two consecutive CEs is smaller than $d$, in which
case no checkpoint for the second CE exists). We can then show results for
various values of $d$, starting from the smallest possible value of $1$ (i.e.,
emitting forecasts from the immediately previous SDE before the CE).

At each checkpoint, a forecast interval is produced, as per Section 4.4. Some
of the metrics we can use to assess the quality of the forecasts assume that
forecasts are in the form of points. Such point metrics are the following: the
Root Mean Squared Error (RMSE) and the Mean Absolute Error (MAE) (the latter
is less sensitive than RMSE to outliers). If we denote by $C$ the set of all
checkpoints, by $y_{c}$ the actual distance (in number of events) between a
checkpoint $c$ and its CE (which is always $d$) and by $\hat{y}_{c}$ our
forecast, then the definitions for RMSE and MAE are as follows:

$\mathit{RMSE}=\sqrt{\frac{1}{\lvert C\rvert}\sum_{c\in
C}{{\lvert\hat{y}_{c}-y_{c}\rvert}}^{2}}$ (7) $\mathit{MAE}=\frac{1}{\lvert
C\rvert}\sum_{c\in C}{\lvert\hat{y}_{c}-y_{c}\rvert}$ (8)

When these metrics are used, we need to impose an extra constraint on the
forecasts, requiring that the maximum spread of each forecast is $0$, i.e., we
produce point (instead of interval) forecasts.

Besides these points metrics, we can also use an interval metric, the so-
called _negatively oriented interval score_ [28]. If
$\hat{y}_{c}=(l_{c},u_{c})$ is an interval forecast produced with confidence
$b=1-a$ and $y_{c}$ the actual observation (distance), then the negatively
oriented interval score (NOIS) for this forecast is given by:

$\mathit{NOIS}_{c}=(u_{c}-l_{c})+\frac{2}{a}(l_{c}-y_{c})I_{x<l_{c}}+\frac{2}{a}(y_{c}-u_{c})I_{x>u_{c}}$
(9)

We can then estimate the average negatively oriented interval score (ANOIS) as
follows:

$\mathit{ANOIS}=\frac{1}{\lvert C\rvert}\sum_{c\in C}{\mathit{NOIS}_{c}}$ (10)

The best possible value for ANOIS is $0$ and is achieved only when all
forecasts are point forecasts (so that $u_{c}-l_{c}$ is always $0$) and all of
them are also correct (so that the last two terms in Eq. 9 are always $0$). In
every other case, a forecast is penalized if its interval is long (so that
focused intervals are promoted), regardless of whether it is correct. If it is
indeed correct, no other penalty is added. If it is not correct, then an extra
penalty is added, which is essentially the deviation of the forecast from the
actual observation, weighted by a factor that grows with the confidence of the
forecast. For example, if the confidence is $100\%$, then $b=1.0$, $a=0.0$ and
the extra penalty, according to Eq. 9, grows to infinity. Incorrect forecasts
produced with high confidence are thus severely penalized. Note that if we
emit only point forecasts ($\hat{y}_{c}=l_{c}=u_{c}$), then NOIS and ANOIS
could be written as follows:

$\mathit{NOIS}_{c}=\frac{2}{a}\lvert\hat{y}_{c}-y_{c}\rvert$ (11)
$\mathit{ANOIS}=\frac{1}{\lvert C\rvert}\sum_{c\in
C}{\frac{2}{a}\lvert\hat{y}_{c}-y_{c}\rvert}$ (12)

ANOIS then becomes a weighted version of MAE.

###  6.3 Classification CE Forecasting

The last forecasting task is the most challenging. In contrast to regression
experiments, where we emit forecasts only at a specified distance before each
CE, in classification experiments we emit forecasts regardless of whether a CE
occurs or not. The goal is to predict the occurrence of CEs within a short
future window or provide a “negative” forecast if our model predicts that no
CE is likely to occur within this window. In practice, this task could be the
first one performed at every new event arrival. If a positive forecast is
emitted, then the regression task could also be performed in order to pinpoint
more accurately when the CE will occur.

One issue with classification experiments is that it is not so straightforward
to establish checkpoints in the stream. In regression experiments, CEs provide
a natural point of reference. In classification experiments, we do not have
such reference points, since we are also required to predict the absence of
CEs. As a result, instead of using the stream to find checkpoints, we can use
the structure of the automaton itself. We may not know the actual distance to
a CE, but the automaton can provide us with an “expected” or “possible”
distance, as follows. For an automaton that is in a final state, it can be
said that the distance to a CE is $0$. More conveniently, we can say that the
“process” which it describes has been completed or, equivalently, that there
remains $0\%$ of the process until completion. For an automaton that is in a
non-final state but separated from a final state by $1$ transition, it can be
said that the “expected” distance is $1$. We use the term “expected” because
we are not interested in whether the automaton will actually take the
transition to a final state. We want to establish checkpoints both for the
presence and the absence of CEs. When the automaton fails to take the
transition to a final state (and we thus have an absence of a CE), this
“expected” distance is not an actual distance, but a “possible” one that
failed to materialize. We also note that there might also exist other walks
from this non-final state to a final one whose length could be greater than
$1$ (in fact, there might exist walks with “infinite length”, in case of
loops). In order to estimate the “expected” distance of a non-final state, we
only use the shortest walk to a final state.

After estimating the expected distances of all states, we can then express
them as percentages by dividing them by the greatest among them. A $0\%$
distance will thus refer to final states, whereas a $100\%$ distance to the
state(s) that are the most distant to a final state, i.e., the automaton has
to take the most transitions to reach a final state. These are the start
states. We can then determine our checkpoints by specifying the states in
which the automaton is permitted to emit forecasts, according to their
“expected” distance. For example, we may establish checkpoints by allowing
only states with a distance between $40\%$ and $60\%$ to emit forecasts. The
intuition here is that, by increasing the allowed distance, we make the
forecasting task more difficult. Another option for measuring the distance of
a state to a possible future CE would be to use the waiting-time distribution
of the state and set its expectation as the distance. However, this assumes
that we have gone through the training phase first and learned the
distributions. For this reason, we avoid using this way to estimate distances.

The evaluation task itself consists of the following steps. At the arrival of
every new input event, we first check whether the distance of the new
automaton state falls within the range of allowed distances, as explained
above. If the new state is allowed to emit a forecast, we use its waiting-time
distribution to produce the forecast. Two parameters are taken into account:
the length of the future window $w$ within which we want to know whether a CE
will occur and the confidence threshold $\theta_{fc}$. If the probability of
the first $w$ points of the distribution exceeds the threshold $\theta_{fc}$,
we emit a positive forecast, essentially affirming that a CE will occur within
the next $w$ events; otherwise, we emit a negative forecast, essentially
rejecting the hypothesis that a CE will occur. We thus obtain a binary
classification task.

As a consequence, we can make use of standard classification measures, like
precision and recall. Each forecast is evaluated: a) as a _true positive_ (TP)
if the forecast is positive and the CE does indeed occur within the next $w$
events from the forecast; b) as a _false positive_ (FP) if the forecast is
positive and the CE does not occur; c) as a _true negative_ (TN) if the
forecast is negative and the CE does not occur and d) as a _false negative_
(FN) if the forecast is negative and the CE does occur; Precision is then
defined as $\mathit{Precision}=\frac{TP}{TP+FP}$ and recall (also called
sensitivity or true positive rate) as $\mathit{Recall}=\frac{TP}{TP+FN}$. As
already mentioned, CEs are relatively rare in a stream. It is thus important
for a forecasting engine to be as specific as possible in identifying the true
negatives. For this reason, besides precision and recall, we also use
_specificity_ (also called true negative rate), defined as
$\mathit{Specificity}=\frac{TN}{TN+FP}$.

A classification experiment is performed as follows. For various values of the
“expected” distance and the confidence threshold $\theta_{fc}$, we estimate
precision, recall and specificity on a test dataset. For a given distance,
$\theta_{fc}$ acts as a cut-off parameter. For each value of $\theta_{fc}$, we
estimate the recall (sensitivity) and specificity scores and we plot these
scores as a ROC curve. For each distance, we then estimate the area under
curve (AUC) for the ROC curves. The higher the AUC value, the better the model
is assumed to be.

The setting described above is the most suitable for evaluation purposes, but
might not be the most appropriate when such a system is actually deployed. For
deployment purposes, another option would be to simply set a best, fixed
confidence threshold (e.g., by selecting, after evaluation, the threshold with
the highest F1-score or Matthews correlation coefficient) and emit only
positive forecasts, regardless of their distance. Forecasts with low
probabilities (i.e., negative forecasts) will thus be ignored/suppressed. This
is justified by the fact that a user would typically be more interested in
positive forecasts. For evaluation purposes, this would not be an appropriate
experimental setting, but it would suffice for deployment purposes, where we
would then be focused on fine-tuning the confidence threshold. In this paper,
we focus on evaluating our system and thus do not discuss further any
deployment solution.

##  7 Empirical Evaluation

We now present experimental results on two datasets, a synthetic one (Section
7.3) and a real-world one (Section 7.4). We first briefly discuss in Section
7.1 the models that we evaluated and present our software and hardware
settings in Section 7.2. We show results for SDE forecasting and
classification CE forecasting. Due to space limitations, we omit results for
regression CE forecasting . We intend to present them in future work.

###  7.1 Models Tested

In the experiments that we present, we evaluated the variable-order Markov
model that we have presented in this paper in its two versions: the memory
efficient one that bypasses the construction of a Markov chain and makes
direct use of the $\mathit{PST}$ learned from a stream (Section 4.4.2) and the
computationally efficient one that constructs a $\mathit{PSA}$ (Section
4.4.1). We compared these against four other models inspired by the relevant
literature.

The first, described in [6, 7], is the most similar in its general outline to
our proposed method. It is a previous version of our system presented in this
paper and is also based on automata and Markov chains. The main difference is
that it attempts to construct full-order Markov models of order $m$ and is
thus typically restricted to low values for $m$.

The second model is presented in [40], where automata and Markov chains are
used once again. However, the automata are directly mapped to Markov chains
and no attempt is made to ensure that the Markov chain is of a certain order.
Thus, in the worst case, this model essentially makes the assumption that SDEs
are i.i.d. and $m=0$.

As a third alternative, we evaluated a model that is based on Hidden Markov
Models (HMM), similar to the work presented in [42]. That work uses the Esper
event processing engine [1] and attempts to model a business process as a HMM.
For our purposes, we use a HMM to describe the behavior of an automaton,
constructed from a given symbolic regular expression. The observation variable
of the HMM corresponds to the states of the automaton. Thus, the set of
possible values of the observation variable is the set of automaton states. An
observation sequence of length $l$ for the HMM consists of the sequence of $l$
states visited by the automaton after consuming $l$ SDEs. The $l$ SDEs
(symbols) are used as values for the hidden variable. The last $l$ symbols are
the last $l$ values of the hidden variable. Therefore, this HMM always has $l$
hidden states, whose values are taken from the SDEs, connected to $l$
observations, whose values are taken from the automaton states. We can train
such a HMM with the Baum-Welch algorithm, using the automaton to generate a
training observation sequence from the original training stream. We can then
use this learned HMM to produce forecasts on a test dataset. We produce
forecasts in an online manner as follows: as the stream is consumed, we use a
buffer to store the last $l$ states visited by the pattern automaton. After
every new event, we “unroll” the HMM using the contents of the buffer as the
observation sequence and the transition and emission matrices learned during
the training phase. We can then use the forward algorithm to estimate the
probability of all possible future observation sequences (up to some length),
which, in our case, correspond to future states visited by the automaton.
Knowing the probability of every future sequence of states allows us to
estimate the waiting-time distribution for the current state of the automaton
and thus build a forecast, as already described. Note that, contrary to the
previous approaches, the estimation of the waiting-time distribution via a HMM
must be performed online. We cannot pre-compute the waiting-time distributions
and store the forecasts in a look-up table, due to the possibly large number
of entries. For example, assume that $l=5$ and the size of the “alphabet” (SDE
symbols) of our automaton is $10$. For each state of the automaton, we would
have to pre-compute $10^{5}$ entries. In other words, as with Markov chains,
we still have a problem of combinatorial explosion. We try to “avoid” this
problem by estimating the waiting-time distributions online.

Our last model is inspired by the work presented in [52]. This method comes
from the process mining community and has not been previously applied to CEF.
However, due to its simplicity, we use it here as a baseline method. We again
use a training dataset to learn the model. In the training phase, every time
the pattern automaton reaches a certain state $q$, we simply count how long
(how many transitions) we have to wait until it reaches a final state. After
the training dataset has been consumed, we end up with a set of such “waiting
times” for every state. The forecast to be produced by each state is then
estimated simply by calculating the average “waiting time”.

As far as the Markov models are concerned, we try to increase their order to
the highest possible value, in order to determine if and how high-order values
offer an advantage. We have empirically discovered that our system can
efficiently handle automata and Markov chains that have up to about $1200$
states. Beyond this point, it becomes almost prohibitive (with our hardware)
to create and handle transition matrices with more than $1200^{2}$ elements.
We have thus set this number as an upper bound and increased the order of a
model until this number is reached. This restriction is applied both to full-
order models and variable-order models that use a $\mathit{PSA}$ and an
embedding, since in both of these cases we need to construct a Markov chain.
For the variable-order models that make direct use of a $\mathit{PST}$, no
Markov chain is constructed. We thus increase their order until their
performance scores seem to reach a stable number or a very high number, beyond
which it makes little sense to continue testing.

###  7.2 Hardware and Software Settings

All experiments were run on a 64-bit Debian 10 (buster) machine with Intel
Core i7-8700 CPU @ 3.20GHz × 12 processors and 16 GB of memory. Our framework
was implemented in Scala 2.12.10. We used Java 1.8, with the default values
for the heap size. For the HMM models, we relied on the Smile machine learning
library [2]. All other models were developed by us. No attempt at
parallelization was made.

###  7.3 Credit Card Fraud Management

The first dataset used in our experiments is a synthetic one, inspired by the
domain of credit card fraud management [11]. We start with a synthetically
generated dataset in order to investigate how our method performs under
conditions that are controlled and produce results more readily interpretable.
The data generator was developed in collaboration with
Feedzai222https://feedzai.com, our partner in the SPEEDD
project333http://speedd-project.eu.

In this dataset, each event is supposed to be a credit card transaction,
accompanied by several arguments, such as the time of the transaction, the
card ID, the amount of money spent, the country where the transaction took
place, etc. In the real world, a very small proportion of such transactions
are fraudulent and the goal of a CER system would be to detect, with very low
latency, fraud instances. To do so, a set of fraud patterns must be provided
to the engine. For typical cases of such patterns in a simplified form, see
[11]. In our experiments, we use one such pattern, consisting of a sequence of
consecutive transactions, where the amount spent at each transaction is
greater than that of the previous transaction. Such a trend of steadily
increasing amounts constitutes a typical fraud pattern. The goal in our
forecasting experiments is to predict if and when such a pattern will be
completed, even before it is detected by the engine (if in fact a fraud
instance occurs), so as to possibly provide a wider margin for action to an
analyst.

We generated a dataset consisting of 1,000,000 transactions in total from 100
different cards. About $20\%$ of the transactions are fraudulent. Not all of
these instances of fraud belong to the pattern of increasing amounts. We
actually inject seven different types of known fraudulent patterns in the
dataset, including, for instance, a decreasing trend. Each fraudulent sequence
for the increasing trend consists of eight consecutive transactions with
increasing amounts, where the amount is increased each time by $100$ monetary
units or more. We additionally inject sequences of transactions with
increasing amounts, which are not fraudulent. In those cases, we randomly
interrupt the sequence before it reaches the eighth transaction. In the
legitimate sequences the amount is increased each time by $0$ or more units.
With this setting, we want to test the effect of long-term dependencies on the
quality of the forecasts. For example, a sequence of six transactions with
increasing amounts, where all increases are $100$ or more units is very likely
to lead to a fraud detection. On the other hand, a sequence of just two
transactions with the same characteristics, could still possibly lead to a
detection, but with a significantly reduced probability. We thus expect that
models with deeper memories will perform better. We used $75\%$ of the dataset
for training and the rest for testing. No k-fold cross validation is
performed, since each fold would have exactly the same statistical properties.

Formally, the symbolic regular expression that we use to capture the pattern
of an increasing trend in the amount spent is the following:

$\displaystyle R:=$ $\displaystyle\
(\mathit{amountDiff>0})\cdot(\mathit{amountDiff>0})\cdot(\mathit{amountDiff>0})\cdot(\mathit{amountDiff>0})\cdot$
(13) $\displaystyle\
(\mathit{amountDiff>0})\cdot(\mathit{amountDiff>0})\cdot(\mathit{amountDiff>0})$

$\mathit{amountDiff}$ is an extra attribute (besides the card ID, the amount
spent, the transaction country and the other standard attributes) with which
we enrich each event and is equal to the difference between the amount spent
by the current transaction and that spent by the immediately previous
transaction from the same card. The expression consists of seven terminal sub-
expressions, in order to capture eight consecutive events. The first terminal
sub-expression captures an increasing amount between the first two events in a
fraudulent pattern.

If we attempted to perform forecasting based solely on Pattern (13), then the
minterms that would be created would be based only on the predicate
$\mathit{amountDiff}>0$: namely, the predicate itself, along with its negation
$\neg(\mathit{amountDiff}>0)$. As expected, such an approach does not yield
good results, as the language is not expressive enough to differentiate
between fraudulent and legitimate transaction sequences. In order to address
this lack of informative (for forecasting purposes) predicates, we have
incorporated a mechanism in our system that allows us to incorporate extra
predicates when building a probabilistic model, without affecting the
semantics of the initial expression (exactly the same matches are detected).
We do this by using any such extra predicates during the construction of the
minterms. For example, if $\mathit{country}=\mathit{MA}$ is such an extra
predicate that we would like included, then we would construct the following
minterms for Pattern (13): a)
$m_{1}=(\mathit{amountDiff>0})\wedge(\mathit{country}=\mathit{MA})$; b)
$m_{2}=(\mathit{amountDiff>0})\wedge\neg(\mathit{country}=\mathit{MA})$; c)
$m_{3}=\neg(\mathit{amountDiff>0})\wedge(\mathit{country}=\mathit{MA})$; d)
$m_{4}=\neg(\mathit{amountDiff>0})\wedge\neg(\mathit{country}=\mathit{MA})$).
We can then use these enhanced minterms as guards on the automaton transitions
in a way that does not affect the semantics of the expression. For example, if
an initial transition has the guard $\mathit{amountDiff>0}$, then we can split
it into two derived transitions, one for $m_{1}$ and one for $m_{2}$. The
derived transitions would be triggered exactly when the initial one is
triggered, the only difference being that the derived transitions also have
information about the country. For our experiments and for Pattern (13), if we
include the extra predicate $\mathit{amountDiff}>100$, we expect the model to
be able to differentiate between sequences involving genuine transactions
(where the difference in the amount can by any value above $0$) and fraudulent
sequences (where the difference in the amount is always above $100$ units).

We now present results for SDE forecasting. As already mentioned in Section
7.1, for this type of experiments we do not use the automaton created by
Pattern (13). We instead use only its minterms which will constitute our
“alphabet”. In our case, there are four minterms: a)
$\mathit{amountDiff}>0\wedge\mathit{amountDiff}>100$; b)
$\mathit{amountDiff}>0\wedge\neg(\mathit{amountDiff}>100)$; c)
$\neg(\mathit{amountDiff}>0)\wedge\mathit{amountDiff}>100$; d)
$\neg(\mathit{amountDiff}>0)\wedge\neg(\mathit{amountDiff}>100)$. Thus, the
goal is to predict, as the stream is consumed, which one of these minterms
will be satisfied. Notice that, for every possible event, exactly one minterm
is satisfied (the third one,
$\neg(\mathit{amountDiff}>0)\wedge\mathit{amountDiff}>100$, is actually
unsatisfiable). We use $75\%$ of the original dataset (which amounts to
750.000 transactions) for training and the rest (250.000) for testing. We do
not employ cross-validation, as the dataset is synthetic and the statistical
properties of its folds would not differ.

(a) Average log-loss.

(b) Number of states/nodes.

Figure 10: Results for SDE forecasting from the domain of credit card fraud
management. Fx stands for a Full-order Markov Model of order $x$. Vx stands
for a Variable-order Markov Model (a prediction suffix tree) of maximum order
$x$.

(a) ROC curves for the variable-order model using the $\mathit{PST}$ for
various values of the maximum order $m$. $\mathit{distance}\in[0.2,0.4]$.

(b) ROC curves for the variable-order model using the $\mathit{PST}$ for
various values of the maximum order. $\mathit{distance}\in[0.4,0.6]$.

(c) AUC for ROC curves for all models.

Figure 11: Results for CE forecasting from the domain of credit card fraud
management. Fx stands for a Full-order Markov Model of order $x$. Ex stands
for a Variable-order Markov Model of maximum order $x$ that uses a
$\mathit{PSA}$ and creates an embedding. Tx stands for a Variable-order Markov
Model of maximum order $x$ that is constructed directly from a $\mathit{PST}$.
MEAN stands for the method of estimating the mean of “waiting-times”. HMM
stands for Hidden Markov Model. IID stands for the method assuming (in the
worst case) that SDEs are i.i.d. Ex and Tx models are the ones proposed in
this paper.

Figure 10(a) shows the average log-loss obtained for various models and orders
$m$ and Figure 10(b) shows the number of states for the full-order models or
nodes for the variable-order models, which are prediction suffix trees. The
best result is achieved with a variable-order Markov model of maximum order 7.
The full-order Markov models are slightly better than their equivalent (same
order) variable-order models. This is an expected behavior, since the
variable-order models are essentially approximations of the full-order ones.
We increase the order of the full-order models until $m=3$, in which case the
Markov chain has $\approx\ 750$ states. We avoid increasing the order any
further, because Markov chains with more than 1000 states become significantly
difficult to manage in terms of memory usage (and in terms of the
computational cost of estimating the waiting-time distributions for the
experiments of CE forecasting that follow). Note that a Markov chain with
$1000$ states would require a transition matrix with $1000^{2}$ entries. On
the contrary, we can increase the maximum order of the variable-order model
until we find the best one, i.e., the order after which the average log-loss
starts increasing again. The size of the prediction suffix tree can be allowed
to increase to more than $1000$ nodes, since we are not required to build a
transition matrix.

(a) Throughput.

(b) Training time.

(c) Number of states/nodes.

Figure 12: Results for throughput, training time and number of automaton
states/tree nodes for classification CE forecasting from the domain of credit
card fraud management. modelTime = time to construct the model. wtTime = time
to estimate the waiting-time distributions for all states. inTime = time to
estimate the forecast of all states from their waiting-time distributions.
extraTime = time to determinize an automaton (+ disambiguation time for full-
order models).

We now move on to the classification experiments. Figure 11 shows the ROC
curves of the variable-order model that directly uses a $\mathit{PST}$. We
show results for two different “expected” distance ranges:
$\mathit{distance}\in[0.2,0.4]$ in Figure 11(a) and
$\mathit{distance}\in[0.4,0.6]$ in Figure 11(b). The ideal operating point in
the ROC is the top-left corner and thus, the closer to that point the curve
is, the better. Thus, the first observation is that by increasing the maximum
order we obtain better results. Notice, however, that in Figure 11(b), where
the distance is greater and the forecasting problem is harder, increasing the
order from 6 to 7 yields only a marginal improvement.

Figure 11(c) displays ROC results for different distances and all models, in
terms of the Area Under the ROC Curve (AUC), which is a measure of the models’
classification accuracy. The first observation is that the MEAN and HMM
methods consistently underperform, compared to the Markov models. Focusing on
the Markov models, as expected, the task becomes more challenging and the ROC
scores decrease, as the distance increases. It is also evident that higher
orders lead to better results. The advantage of increasing the order becomes
less pronounced (or even non-existent) as the distance increases. The
variable-order models that use an embedding are only able to go as far as
$m=4$, due to increasing memory requirements, whereas the tree-based versions
can go up to $m=7$ (and possibly even further, but we did not try to extend
the order beyond this point). Although the embedding ($\mathit{PSA}$) can
indeed help achieve better scores than full-order models by reaching higher
orders, this is especially true for the tree-based models which bypass the
embedding. We can thus conclude that full-order models are doing well up to
the order that they we can achieve with them. $\mathit{PSA}$ models can reach
roughly the same levels, as they are also practically restricted. The
performance of $\mathit{PST}$ models is similar to that of the other models
for the same order, but the fact that they can use higher orders allows them
to finally obtain superior performance.

We show performance results in Figure 12, in terms of computation and memory
efficiency. Figure 12(a) displays throughput results. We can observe the
trade-off between the high forecasting accuracy of the tree-based high-order
models and the performance penalty that these models incur. The models based
on $\mathit{PST}$ have a throughput figure that is almost half that of the
full-order models and the embedding-based variable-order ones. In order to
emit a forecast, the tree-based models need to traverse a tree after every new
event arrives at the system, as described in Section 4.4.2. The automata-based
full- and variable-order models, on the contrary, only need to evaluate the
minterms on the outgoing transitions of their current state and simply jump to
the next state. It would be possible to improve the throughput of the tree-
based models, by using caching techniques, so that we can reuse some of the
previously estimated forecasts, but we reserve such optimizations for future
work. By far the worst throughput, however, is observed for the HMM models.
The reason is that the waiting-time distributions and forecasts are always
estimated online, as explained in Section 7.1.

Figure 12(b) shows training times as a stacked, bar plot. For each model, the
total training time is broken down into 4 different components, each
corresponding to a different phase of the forecast building process.
_modelTime_ is the time required to actually construct the model from the
training dataset. _wtTime_ is the time required to estimate the waiting-time
distributions, once the model has been constructed. _inTime_ measures the time
required to estimate the forecast of each waiting-time distribution. Finally,
_extraTime_ measures the time required to determinize the automaton of our
initial pattern. For the full-order Markov models, it also includes the time
required to convert the deterministic automaton into its equivalent,
disambiguated automaton. We observe that the tree-based models exhibit
significantly higher times than the rest, for high orders. The other models
have similar training times, almost always below $5$ seconds. Thus, if we need
high accuracy, we again have to pay a price in terms of training time. Even in
the case of high-order tree-based models though, the training time is almost
half a minute for a training dataset composed of 750,000 transactions, which
allows us to be confident that training could be performed online.

Figure 12(c) shows the memory footprint of the models in terms of the size of
their basic data structures. For automata-based methods, we show the number of
states, whereas for the tree-based methods we show the number of nodes. We see
that variable-order models, especially the tree-based ones, are significantly
more compact than the full-order ones, for the same order. We also observe
that the tree-based methods, for the same order, are much more compact (fewer
nodes) than the ones based on the embedding (more states). This allows us to
increase the order up to $7$ with the tree-based approach, but only up to $4$
with the embedding.

###  7.4 Maritime Situational Awareness

The second dataset that we used in our experiments is a real–world dataset
coming from the field of maritime monitoring. It is composed of a set of
trajectories from ships sailing at sea, emitting AIS (Automatic Identification
System) messages that relay information about their position, heading, speed,
etc., as described in the running example of Section 1.1. These trajectories
can be analyzed, using the techniques of Complex Event Recognition, in order
to detect interesting patterns in the behavior of vessels [44]. The dataset
that we used is publicly available, contains AIS kinematic messages from
vessels sailing in the Atlantic Ocean around the port of Brest, France, and
spans a period from 1 October 2015 to 31 March 2016 [46]. We used a derivative
dataset that contains clean and compressed trajectories, consisting only of
critical points [43]. Critical points are the important points of a trajectory
that indicate a significant change in the behavior of a vessel. Using critical
points, one can reconstruct quite accurately the original trajectory [44]. We
further processed the dataset by interpolating between the critical points in
order to produce trajectories where two consecutive points have a temporal
distance of exactly 60 seconds. The reason for this pre-processing step is
that AIS messages typically arrive at unspecified time intervals. These
intervals can exhibit a very wide variation, depending on many factors (e.g.,
human operators may turn on/off the AIS equipment), without any clear pattern
that could be encoded by our probabilistic model. Consequently, our system
performs this interpolation as a first step.

Figure 13: Trajectories of the vessel with the most matches for Pattern (14)
around the port of Brest. Green points denote ports. Red, shaded circles show
the areas covered by each port. They are centered around each port and have a
radius of 5 km.

(a) Average log-loss.

(b) Number of states/nodes.

Figure 14: Results for SDE forecasting from the domain of maritime monitoring.
Fx stands for a Full-order Markov Model of order $x$. Vx stands for a
Variable-order Markov Model (a prediction suffix tree) of maximum order $x$.

The pattern that we used in the experiments is a movement pattern in which a
vessel approaches the main port of Brest. The goal is to forecast when a
vessel will enter the port. This way, port traffic management may be
optimized, in order to reduce the carbon emissions of vessels waiting to enter
the port. The symbolic regular expression for this pattern is the following:

$\displaystyle R:=$ $\displaystyle\
(\neg\mathit{InsidePort(Brest)})^{*}\cdot(\neg\mathit{InsidePort(Brest)})\cdot$
(14) $\displaystyle\
(\neg\mathit{InsidePort(Brest)})\cdot(\mathit{InsidePort(Brest)})$

The intention is to detect the entrance of a vessel in the port of Brest. The
predicate $\mathit{InsidePort(Brest)}$ evaluates to TRUE whenever a vessel has
a distance of less than 5 km from the port of Brest (see Figure 13). In fact,
the predicate is generic and takes as arguments the longitude and latitude of
any point, but we show here a simplified version, using the port of Brest, for
reasons of readability. The pattern defines the entrance to the port as a
sequence of at least 3 consecutive events, only the last of which satisfies
the $\mathit{InsidePort(Brest)}$ predicate. In order to detect an entrance, we
must first ensure that the previous event(s) indicated that the vessel was
outside the port. For this reason, we require that, before the last event,
there must have occurred at least 2 events where the vessel was outside the
port. We require 2 or more such events to have occurred (instead of just one),
in order to avoid detecting “noisy” entrances.

(a) ROC curves for the variable-order model using the $\mathit{PST}$ for
various values of the maximum order. $\mathit{distance}\in[0.0,0.5]$.

(b) ROC curves for the variable-order model using the $\mathit{PST}$ for
various values of the maximum order. $\mathit{distance}\in[0.5,1.0]$.

(c) AUC for ROC curves for all models.

Figure 15: Results for CE forecasting in the domain of maritime situational
awareness. Fx stands for a Full-order Markov Model of order $x$. Vx stands for
a Variable-order Markov Model of maximum order $x$ using an embedding. Tx
stands for a Variable-order Markov Model of maximum order $x$ using a
prediction suffix tree. MEAN stands for the method of estimating the mean of
“waiting-times”. HMM stands for Hidden Markov Model. Ex and Tx models are the
ones proposed in this paper.

In addition to the $\mathit{InsidePort(Brest)}$ predicate, we included 5 extra
ones providing information about the distance of a vessel from a port when it
is outside the port. Each of these predicates evaluates to TRUE when a vessel
lies within a specified range of distances from the port. The first returns
TRUE when a vessel has a distance between 5 and 6 km from the port, the second
when the distance is between 6 and 7 km and the other three extend similarly 1
km until 10 km. We investigated the sensitivity of our models to the presence
of various extra predicates in the recognition pattern.

For all experimental results that follow, we always present average values
over 4 folds of cross-validation. We start by analyzing the trajectories of a
single vessel and then move to multiple, selected vessels. There are two
issues that we tried to address by separating our experiments into single-
vessel and multiple-vessel ones. First, we wanted to have enough data for
training. For this reason, we only retained vessels for which we can detect a
significant number of matches for Pattern (14). Second, our system can work in
two modes: a) it can build a separate model for each monitored object and use
this collection of models for personalized forecasting; b) it can build a
global model out of all the monitored objects. We thus wanted to examine
whether building a global model from multiple vessels could produce equally
good results, as these obtained for a single vessel with sufficient training
data.

We first used Pattern (14) to perform recognition on the whole dataset in
order to find the number of matches detected for each vessel. The vessel with
the most matches was then isolated and we retained only the events emitted
from this vessel. In total, we detected 368 matches for this vessel and the
number of SDEs corresponding to it is $\approx$ 30.000. Figure 13 shows the
isolated trajectories for this vessel, seemingly following a standard route
between various ports around the Brest area.

(a) Throughput.

(b) Training time.

Figure 16: Throughput and training time results for classification CE
forecasting for maritime situational awareness. modelTime = time to construct
the model. wtTime = time to estimate the waiting-time distributions for all
states. inTime = time to estimate the forecast interval of all states from
their waiting-time distributions. extraTime = time to determinize an automaton
(+ disambiguation time for full-order models).

Figure 14 shows results for SDE forecasting. The best average log-loss is
achieved with a full-order Markov model, with $m=3$, and is $\approx 0.51$.
For the best hyper-parameter values out of those that we tested for the
variable-order model, with $m=10$, we can achieve an average log-loss that is
$\approx 0.57$. Contrary to the case of credit card data, increasing the order
of the variable-order model does not allow us to achieve a better log-loss
score than the best one achieved with a full-oder model. However, as we will
show, this does not imply that the same is true for CE forecasting.

Using the vessel of Figure 13, we obtained the results shown in Figures 15 and
16. Since the original $\mathit{DSFA}$ is smaller in this case (one start and
one final state plus two intermediate states), we have fewer distance ranges
(e.g., there no states in the range $[0.4,0.6]$). Thus, we use only two
distance ranges: $[0,0.5]$ and $[0.5,1]$. We observe the importance of being
able to increase the order of our models for distances smaller than $50\%$.
For distances greater than $50\%$, the area under curve is $\approx$ 0.5 for
all models. This implies that they cannot effectively differentiate between
positives and negatives. Their forecasts are either all positive, where we
have $\mathit{Recall}=100\%$ and $\mathit{Specificity}=0\%$, or all negative,
where we have $\mathit{Recall}=0\%$ and $\mathit{Specificity}=100\%$ (see
Figure 15(b)). Notice that the full-order Markov models can now only go up to
$m=2$, since the existence of multiple extra predicates makes it prohibitive
to increase the order any further. Achieving higher accuracy with higher-order
models comes at a computational cost, as shown in Figure 16. The results are
similar to those in the credit card experiments. The training time for
variable-order models tends to increase as we increase the order, but is
always less than 8 seconds. The effect on throughput is again significant for
the tree-based variable-order models. Throughput figures are also lower here
compared to the credit card fraud experiments, since the predicates that we
need to evaluate for every new input event (like $\mathit{InsidePort(Brest)}$)
involve more complex calculations (the $\mathit{amountDiff}>0$ predicate is a
simple comparison).

As a next step, we wanted to investigate the effect of the optimization
technique mentioned at the end of Section 4.4.2 on the accuracy and
performance of our system. The optimization prunes future paths whose
probability is below a given cutoff threshold. We re-run the experiments
described above for distances between $0\%$ and $50\%$ for various values of
the cutoff threshold, starting from $0.0001$ up to $0.2$. Figure 17 shows the
relevant results. We observe that the accuracy is affected only for high
values of the cutoff threshold, above $0.1$ (Figure 17(a)). We can also see
that throughput remains essentially unaffected (Figure 17(b)). This result is
expected, since the cutoff threshold is only used in the estimation of the
waiting-time distributions. Throughput reflects the online performance of our
system, after the waiting-time distributions have been estimated, and is thus
not affected by the choice of the cutoff threshold. However, the training time
is indeed significantly affected (Figure 17(c)). As expected, the result of
increasing the value of the cutoff threshold is a reduction of the training
time, as fewer paths are retained. Beyond a certain point though, further
increases of the cutoff threshold affect the accuracy of the system.
Therefore, the cutoff threshold should be below $0.01$ so as not to compromise
the accuracy of our forecasts.

(a) AUC-ROC.

(b) Throughput.

(c) Training time.

Figure 17: Effect of cutoff threshold on accuracy, throughput and training
time.

We additionally investigated the sensitivity of our approach to the extra
predicates that are used. Figure 18(a) shows results when the extra 5
predicates referring to the distance of a vessel from the port are modified so
that each “ring” around the port has a width of 3 km, instead of 1 km. With
these extra features, increasing the order does indeed make an important
difference, but only when the order becomes high (5 and beyond), which is
possible only by using the tree-based variable-order models. Moreover, the
best score achieved is still lower than the best score achieved with “rings”
of 1 km (Figure 15(c)). “Rings” of 1 km are thus more appropriate as
predictive features. We also wanted to investigate whether other information
about the vessel’s movement could affect forecasting. In particular, we kept
the 1 km “rings” and we also added a predicate to check whether the heading of
a vessel points towards the port. More precisely, we used the vessel’s speed
and heading to project its location 1 hour ahead in the future and then
checked whether this projected segment and the circle around the port
intersect. The intuition for adding this feature is that the knowledge of
whether a vessel is heading towards the port has predictive value. As shown in
Figure 18(b), this additional information led to higher scores even with low
full-order orders (compare to Figure 15(c)). The heading feature is indeed
important. On the other hand, the high-order models that did not use this
feature seemed to be able to compensate for the missing information about the
vessel’s heading by going into higher orders. A plateau is thus reached, which
cannot be “broken” with the heading information. Notice that here we can only
go up to $m=1$ for full-order models. The inclusion of the heading predicate
leads to an increase of the number of states beyond $1200$.

(a) AUC for ROC curves. Extra features included: concentric rings around the
port every 3 km. Single vessel.

(b) AUC for ROC curves. Extra features included: concentric rings around the
port every 1 km and heading. Single vessel.

(c) AUC for ROC curves. Extra features included: concentric rings around the
port every 1 km. Model constructed for the 9 vessels that have more than 100
matches.

Figure 18: Results for classification CE forecasting from the domain of
maritime monitoring for various sets of extra features and for multiple
vessels.

Finally, we also tested our method when more than one vessel need to be
monitored. Instead of isolating the single vessel with the most matches, we
isolated all vessels which had more than 100 matches. There are in total 9
such vessels in the dataset. The resulting dataset has $\approx$ 222.000
events. Out of the 9 retained vessels, we constructed a global probabilistic
model and produced forecasts. An alternative option would be to build a single
model for each vessel, but in this scenario we wanted to test the robustness
of our aprroach when a global model is built from multiple entities. Figure
18(c) presents the corresponding results. Interestingly, the scores of the
global model remain very close to the scores of the experiments for the single
vessel with the most matches (Figure 15(c)). This is an indication of the
ability of the global model to capture the peculiarities of individual
vessels.

##  8 Summary & Future Work

We have presented a framework for Complex Event Forecasting (CEF), based on a
variable-order Markov model. It allows us to delve deeper into the past and
capture long-term dependencies, not feasible with full-order models. Our
comprehensive evaluation on two application domains has illustrated the
advantages of being able to use such high-order models. Namely, the use of
higher-order modeling allows us to achieve higher accuracy than what is
possible with full-order models or other state-of-the-art solutions. We have
described two alternative ways in which variable-order models may be used,
depending on the imposed requirements. One option is to use a highly efficient
but less accurate model, when online performance is a top priority. We also
provide an option that achieves high accuracy scores, but with a performance
cost. Another important feature of our proposed framework is that it requires
minimal intervention by the user. A given Complex Event pattern is
declaratively defined and subsequently automatically translated to an
automaton and then to a Markov model, without requiring domain knowledge that
should guide the modeling process.

Still, the user needs to set up the model and there seems to be room for
further automation. In particular, the user needs to set the maximum order
allowed by the probabilistic model. Additionally, we have started
investigating ways to handle concept drift by continuously training and
updating the probabilistic model of a pattern.

With respect to the expressive power of our framework, there is one
functionality that we do not currently support and whose incorporation would
require us to move to a more advanced automaton model. This is the
functionality of applying $n$-ary (with $n>1$) predicates to two or more sub-
expressions of an expression, instead of only unary predicates, as is allowed
in symbolic automata. Note that $n$ refers to the number of “terminal symbols”
/ events that a predicate may reference. Each transition predicate of a
symbolic automaton may refer only to one event, the last event consumed. It
cannot apply a predicate to the last event and other earlier events. As an
example, consider the pattern $R:=x\cdot y\ \textsf{\footnotesize WHERE}\
y.\mathit{speed}>x.\mathit{speed}$, detecting an increase in the speed of a
vessel, where we now need to use the variables $x$ and $y$. Such patterns
cannot be captured with $\mathit{SFA}$ since they would require a memory
structure to store some of the past events of a stream, as is possible with
extended symbolic automata [20]. We intend to present in future work an
automaton model which can support patterns with memory, suitable for CER. Some
results towards this direction may be found in [8].

Finally, our framework could also be used for a task that is not directly
related to Complex Event Forecasting. Since predictive modeling and
compression are two sides of the same coin, our framework could be used for
pattern-driven lossless stream compression, in order to minimize the
communication cost, which is a severe bottleneck for geo-distributed CER [27].
The probabilistic model that we construct with our approach could be pushed
down to the event sources, such as the vessels in the maritime domain, in
order to compress each individual stream and then these compressed streams
could be transmitted to a centralized CER engine to perform recognition.

## References

  * [1] Esper.  http://www.espertech.com/esper.
  * [2] Smile - statistical machine intelligence and learning engine.  http://haifengl.github.io/.
  * [3] Naoki Abe and Manfred K. Warmuth.  On the computational complexity of approximating distributions by probabilistic automata.  Machine Learning, 9:205–260, 1992.
  * [4] Rakesh Agrawal, Tomasz Imielinski, and Arun N. Swami.  Mining association rules between sets of items in large databases.  In SIGMOD Conference, pages 207–216. ACM Press, 1993.
  * [5] Adnan Akbar, François Carrez, Klaus Moessner, and Ahmed Zoha.  Predicting complex events for pro-active iot applications.  In WF-IoT, pages 327–332. IEEE Computer Society, 2015.
  * [6] Elias Alevizos, Alexander Artikis, and George Paliouras.  Event forecasting with pattern markov chains.  In DEBS, pages 146–157. ACM, 2017.
  * [7] Elias Alevizos, Alexander Artikis, and George Paliouras.  Wayeb: a tool for complex event forecasting.  In LPAR, volume 57 of EPiC Series in Computing, pages 26–35. EasyChair, 2018.
  * [8] Elias Alevizos, Alexander Artikis, and Georgios Paliouras.  Symbolic automata with memory: a computational model for complex event processing.  CoRR, abs/1804.09999, 2018.
  * [9] Elias Alevizos, Anastasios Skarlatidis, Alexander Artikis, and Georgios Paliouras.  Probabilistic complex event recognition: A survey.  ACM Comput. Surv., 50(5):71:1–71:31, 2017.
  * [10] Alexander Artikis, Chris Baber, Pedro Bizarro, Carlos Canudas-de-Wit, Opher Etzion, Fabiana Fournier, Paul Goulart, Andrew Howes, John Lygeros, Georgios Paliouras, Assaf Schuster, and Izchak Sharfman.  Scalable proactive event-driven decision making.  IEEE Technol. Soc. Mag., 33(3):35–41, 2014.
  * [11] Alexander Artikis, Nikos Katzouris, Ivo Correia, Chris Baber, Natan Morar, Inna Skarbovsky, Fabiana Fournier, and Georgios Paliouras.  A prototype for credit card fraud management: Industry paper.  In DEBS, pages 249–260. ACM, 2017.
  * [12] Ron Begleiter, Ran El-Yaniv, and Golan Yona.  On prediction using variable order markov models.  J. Artif. Intell. Res., 22:385–421, 2004.
  * [13] Peter Bühlmann, Abraham J Wyner, et al.  Variable length markov chains.  The Annals of Statistics, 27(2):480–513, 1999.
  * [14] Buru Chang, Yonggyu Park, Donghyeon Park, Seongsoon Kim, and Jaewoo Kang.  Content-aware hierarchical point-of-interest embedding model for successive POI recommendation.  In IJCAI, pages 3301–3307. ijcai.org, 2018.
  * [15] Chung-Wen Cho, Yi-Hung Wu, Show-Jane Yen, Ying Zheng, and Arbee L. P. Chen.  On-line rule matching for event prediction.  VLDB J., 20(3):303–334, 2011.
  * [16] Maximilian Christ, Julian Krumeich, and Andreas W. Kempa-Liehr.  Integrating predictive analytics into complex event processing by using conditional density estimations.  In EDOC Workshops, pages 1–8. IEEE Computer Society, 2016.
  * [17] John G. Cleary and Ian H. Witten.  Data compression using adaptive coding and partial string matching.  IEEE Trans. Communications, 32(4):396–402, 1984.
  * [18] Don Coppersmith and Shmuel Winograd.  Matrix multiplication via arithmetic progressions.  J. Symb. Comput., 9(3):251–280, 1990.
  * [19] Gianpaolo Cugola and Alessandro Margara.  Processing flows of information: From data stream to complex event processing.  ACM Comput. Surv., 44(3):15:1–15:62, 2012.
  * [20] Loris D’Antoni and Margus Veanes.  Extended symbolic finite automata and transducers.  Formal Methods in System Design, 47(1):93–119, 2015.
  * [21] Loris D’Antoni and Margus Veanes.  The power of symbolic automata and transducers.  In CAV (1), volume 10426 of Lecture Notes in Computer Science, pages 47–67. Springer, 2017.
  * [22] Yagil Engel and Opher Etzion.  Towards proactive event-driven computing.  In DEBS, pages 125–136. ACM, 2011.
  * [23] Opher Etzion and Peter Niblett.  Event Processing in Action.  Manning Publications Company, 2010.
  * [24] Chiara Di Francescomarino, Chiara Ghidini, Fabrizio Maria Maggi, and Fredrik Milani.  Predictive process monitoring methods: Which one suits me best?  In BPM, volume 11080 of Lecture Notes in Computer Science, pages 462–479. Springer, 2018.
  * [25] James C Fu and WY Wendy Lou.  Distribution theory of runs and patterns and its applications: a finite Markov chain imbedding approach.  World Scientific, 2003.
  * [26] Lajos Jeno Fülöp, Árpád Beszédes, Gabriella Toth, Hunor Demeter, László Vidács, and Lóránt Farkas.  Predictive complex event processing: a conceptual framework for combining complex event processing and predictive analytics.  In BCI, pages 26–31. ACM, 2012.
  * [27] Nikos Giatrakos, Elias Alevizos, Alexander Artikis, Antonios Deligiannakis, and Minos N. Garofalakis.  Complex event recognition in the big data era: a survey.  VLDB J., 29(1):313–352, 2020.
  * [28] Tilmann Gneiting and Adrian E Raftery.  Strictly proper scoring rules, prediction, and estimation.  Journal of the American Statistical Association, 102(477):359–378, 2007.
  * [29] Alejandro Grez, Cristian Riveros, and Martín Ugarte.  A formal framework for complex event processing.  In ICDT, volume 127 of LIPIcs, pages 5:1–5:18. Schloss Dagstuhl - Leibniz-Zentrum fuer Informatik, 2019.
  * [30] Ulrich Hedtstück.  Complex event processing: Verarbeitung von Ereignismustern in Datenströmen.  Springer Vieweg, Berlin, 2017.
  * [31] John E. Hopcroft, Rajeev Motwani, and Jeffrey D. Ullman.  Introduction to automata theory, languages, and computation, 3rd Edition.  Pearson international edition. Addison-Wesley, 2007.
  * [32] Srivatsan Laxman, Vikram Tankasali, and Ryen W. White.  Stream prediction using a generative model based on frequent episodes in event sequences.  In KDD, pages 453–461. ACM, 2008.
  * [33] Yan Li, Tingjian Ge, and Cindy Chen.  Data stream event prediction based on timing knowledge and state transitions.  Proceedings of the VLDB Endowment, 13(10), 2020.
  * [34] Zhongyang Li, Xiao Ding, and Ting Liu.  Constructing narrative event evolutionary graph for script event prediction.  In IJCAI, pages 4201–4207. ijcai.org, 2018.
  * [35] David C. Luckham.  The Power of Events: An Introduction to Complex Event Processing in Distributed Enterprise Systems.  Addison-Wesley, 2001.
  * [36] Spyros Makridakis, Evangelos Spiliotis, and Vassilios Assimakopoulos.  Statistical and machine learning forecasting methods: Concerns and ways forward.  PloS one, 13(3):e0194889, 2018.
  * [37] Heikki Mannila, Hannu Toivonen, and A. Inkeri Verkamo.  Discovery of frequent episodes in event sequences.  Data Min. Knowl. Discov., 1(3):259–289, 1997.
  * [38] Alfonso Eduardo Márquez-Chamorro, Manuel Resinas, and Antonio Ruiz-Cortés.  Predictive monitoring of business processes: A survey.  IEEE Trans. Services Computing, 11(6):962–977, 2018.
  * [39] Douglas C Montgomery, Cheryl L Jennings, and Murat Kulahci.  Introduction to time series analysis and forecasting.  John Wiley & Sons, 2015.
  * [40] Vinod Muthusamy, Haifeng Liu, and Hans-Arno Jacobsen.  Predictive publish/subscribe matching.  In DEBS, pages 14–25. ACM, 2010.
  * [41] Jonathan Ozik, Nicholson Collier, Randy Heiland, Gary An, and Paul Macklin.  Learning-accelerated discovery of immune-tumour interactions.  Molecular systems design & engineering, 4(4):747–760, 2019.
  * [42] Suraj Pandey, Surya Nepal, and Shiping Chen.  A test-bed for the evaluation of business process prediction techniques.  In CollaborateCom, pages 382–391. ICST / IEEE, 2011.
  * [43] K Patroumpas, D Spirelis, E Chondrodima, H Georgiou, Petrou P, Tampakis P, Sideridis S, Pelekis N, and Theodoridis Y.  Final dataset of Trajectory Synopses over AIS kinematic messages in Brest area (ver. 0.8) [Data set], 10.5281/zenodo.2563256, 2018.  URL: http://doi.org/10.5281/zenodo.2563256, doi:10.5281/zenodo.2563256.
  * [44] Kostas Patroumpas, Elias Alevizos, Alexander Artikis, Marios Vodas, Nikos Pelekis, and Yannis Theodoridis.  Online event recognition from moving vessel trajectories.  GeoInformatica, 21(2):389–427, 2017.
  * [45] Lawrence R Rabiner.  A tutorial on hidden markov models and selected applications in speech recognition.  Proceedings of the IEEE, 77(2):257–286, 1989.
  * [46] C Ray, R Dreo, E Camossi, and AL Jousselme.  Heterogeneous Integrated Dataset for Maritime Intelligence, Surveillance, and Reconnaissance, 10.5281/zenodo.1167595, 2018.  URL: https://doi.org/10.5281/zenodo.1167595, doi:10.5281/zenodo.1167595.
  * [47] Dana Ron, Yoram Singer, and Naftali Tishby.  The power of amnesia.  In NIPS, pages 176–183. Morgan Kaufmann, 1993.
  * [48] Dana Ron, Yoram Singer, and Naftali Tishby.  The power of amnesia: Learning probabilistic automata with variable memory length.  Machine Learning, 25(2-3):117–149, 1996.
  * [49] Jacques Sakarovitch.  Elements of Automata Theory.  Cambridge University Press, 2009.
  * [50] Andrew James Stothers.  On the complexity of matrix multiplication.  PhD thesis, The University of Edinburgh, 2010.
  * [51] Wil Van Der Aalst.  Process mining: discovery, conformance and enhancement of business processes.  Springer, 2011.
  * [52] Wil M. P. van der Aalst, M. H. Schonenberg, and Minseok Song.  Time prediction based on process mining.  Inf. Syst., 36(2):450–475, 2011.
  * [53] Margus Veanes, Peli de Halleux, and Nikolai Tillmann.  Rex: Symbolic regular expression explorer.  In ICST, pages 498–507. IEEE Computer Society, 2010.
  * [54] Ricardo Vilalta and Sheng Ma.  Predicting rare events in temporal domains.  In ICDM, pages 474–481. IEEE Computer Society, 2002.
  * [55] George A. Vouros, Akrivi Vlachou, Georgios M. Santipantakis, Christos Doulkeridis, Nikos Pelekis, Harris V. Georgiou, Yannis Theodoridis, Kostas Patroumpas, Elias Alevizos, Alexander Artikis, Christophe Claramunt, Cyril Ray, David Scarlatti, Georg Fuchs, Gennady L. Andrienko, Natalia V. Andrienko, Michael Mock, Elena Camossi, Anne-Laure Jousselme, and Jose Manuel Cordero Garcia.  Big data analytics for time critical mobility forecasting: Recent progress and research challenges.  In EDBT, pages 612–623. OpenProceedings.org, 2018.
  * [56] Frans M. J. Willems, Yuri M. Shtarkov, and Tjalling J. Tjalkens.  The context-tree weighting method: basic properties.  IEEE Trans. Information Theory, 41(3):653–664, 1995.
  * [57] Cheng Zhou, Boris Cule, and Bart Goethals.  A pattern based predictor for event streams.  Expert Syst. Appl., 42(23):9294–9306, 2015.

##  Appendix A Appendix

###  A.1 Complete proof of Proposition 3.4

######  Proposition.

For every symbolic regular expression $R$ there exists a symbolic finite
automaton $M$ such that $\mathcal{L}(R)=\mathcal{L}(M)$.

######  Proof A.1.

Except for the first case, for the other three cases the induction hypothesis
is that the theorem holds for the sub-expressions of the initial expression.

(a) Base case of a single predicate. $R=\psi$.

(b) OR. $R=R_{1}+R_{2}$.

(c) Concatenation. $R=R_{1}\cdot R_{2}$.

(d) Iteration. $R^{{}^{\prime}}=R^{*}$.

Figure 19: The four cases for constructing a SFA from a SRE

Case where $R=\psi,\psi\in\Psi$.
We construct a SFA as in Figure 19(a). If $w\in\mathcal{L}(R)$, then $w$ is a
single character and $w\in\llbracket\psi\rrbracket$, i.e., $\psi$ evaluates to
TRUE for $w$. Thus, upon seeing $w$, the SFA of Figure 19(a) moves from
$q^{s}$ to $q^{f}$ and since $q^{f}$ is a final state, then $w$ is accepted by
this SFA. Conversely, if a string $w$ is accepted by this SFA then it must
again be a single character and $\psi$ must evaluate to TRUE since the SFA
moved to its final state through the transition equipped with $\psi$. Thus,
$w\in\llbracket\psi\rrbracket$ and $w\in\mathcal{L}(R)$.

Case where $R=R_{1}+R_{2}$.
We construct a SFA as in Figure 19(b). If $w\in\mathcal{L}(R)$, then either
$w\in\mathcal{L}(R_{1})$ or $w\in\mathcal{L}(R_{2})$ (or both). Without loss
of generality, assume $w\in\mathcal{L}(R_{1})$. From the induction hypothesis,
it also holds that $w\in\mathcal{L}(M_{R_{1}})$. Thus, from Figure 19(b), upon
reading $w$, $M_{R_{1}}$ will have reached $q_{1}^{f}$. Therefore, $M_{R}$
will have reached $q^{f}$ throught the $\epsilon$-transition connecting
$q_{1}^{f}$ to $q^{f}$ and thus $w$ is accepted by $M_{R}$. Conversely, if
$w\in\mathcal{L}(M_{R})$, then the SFA $M_{R}$ of Figure 19(b) must have
reached $q^{f}$ and therefore also $q_{1}^{f}$ or $q_{2}^{f}$ (or both).
Assume it has reached $q_{1}^{f}$. Then $w\in\mathcal{L}(M_{R_{1}})$ and, from
the induction hypothesis $w\in\mathcal{L}(R_{1})$. Similarly, if its has
reached $q_{2}^{f}$, then $w\in\mathcal{L}(R_{2})$. Therefore,
$w\in\mathcal{L}(R_{1})\cup\mathcal{L}(R_{2})=\mathcal{L}(R)$.

Case where $R=R_{1}\cdot R_{2}$.
We construct a SFA as in Figure 19(c). If $w\in\mathcal{L}(R)$, then
$w\in\mathcal{L}(R_{1})\cdot\mathcal{L}(R_{2})$ or $w=w_{1}\cdot w_{2}$ such
that $w_{1}\in\mathcal{L}(R_{1})$ and $w_{2}\in\mathcal{L}(R_{2})$. Therefore,
from the induction hypothesis, upon reading $w_{1}$, $M_{R}$ will have reached
$q_{1}^{f}$ and $q_{2}^{s}$. Upon reading the rest of $w$ ($w_{2}$), again
from the induction hypothesis, $M_{R}$ will have reached $q_{2}^{f}$. As a
result, $w\in\mathcal{L}(M_{R})$. Conversely, if $w\in\mathcal{L}(M_{R})$,
$M_{R}$ will have reached $q_{2}^{f}$ upon reading $w$ and therefore will have
also passed through $q_{1}^{f}$ upon reading a prefix $w_{1}$ of $w$. Thus,
$w=w_{1}\cdot w_{2}$ with $w_{1}\in\mathcal{L}(M_{R_{1}})$ and
$w_{2}\in\mathcal{L}(M_{R_{2}})$. From the induction hypothesis, it also holds
that $w_{1}\in\mathcal{L}(R_{1})$ and $w_{2}\in\mathcal{L}(R_{2})$ and
therefore that $w\in\mathcal{L}(R)$.

Case where $R^{{}^{\prime}}=R^{*}$.
We construct a SFA as in Figure 19(d). If $w\in\mathcal{L}(R^{{}^{\prime}})$,
then $w\in(\mathcal{L}(R))^{*}$ or, equivalently, $w=w_{1}\cdot
w_{2}\cdot\cdots\cdot w_{k}$ such that $w_{i}\in\mathcal{L}(R)$ for all
$w_{i}$. From the induction hypothesis and Figure 19(d), upon reading $w_{1}$,
$M_{R^{{}^{\prime}}}$ will have reached $q_{1}^{f}$ and $q_{1}^{s}$.
Therefore, the same will be true after reading $w_{2}$ and all other $w_{i}$,
including $w_{k}$. Thus, $w\in\mathcal{L}(M_{R^{{}^{\prime}}})$. Note that if
$w=\epsilon$, the $\epsilon$-transition from $q^{s}$ to $q^{f}$ ensures that
$w\in\mathcal{L}(M_{R^{{}^{\prime}}})$. Conversely, assume
$w\in\mathcal{L}(M_{R^{{}^{\prime}}})$. If $w=\epsilon$, then by the
definition of the ∗ operator, $w\in(\mathcal{L(R)})^{*}$. In every other case,
$M_{R^{{}^{\prime}}}$ must have reached $q_{1}^{f}$ and must have passed
through $q_{1}^{s}$. Therefore, $w$ may be written as $w=w_{1}\cdot w_{2}$
where $w_{2}\in\mathcal{M_{R}}$ and, for $w_{1}$, upon reading it,
$M_{R^{{}^{\prime}}}$ must have reached $q_{1}^{s}$. There are two cases then:
either $w_{1}=\epsilon$ and $q_{1}^{s}$ was reached from $q^{s}$ or
$w_{1}\neq\epsilon$ and $q_{1}^{s}$ was reached from $q_{1}^{f}$. In the
former case, $w=\epsilon\cdot w_{2}=w_{2}$ and thus
$w\in(\mathcal{L}(R))^{*}$. In the latter case, we can apply a similar
reasoning recursively to $w_{1}$ in order to split it to sub-strings $w_{i}$
such that $w_{i}\in\mathcal{L}(R)$. Therefore, $w\in(\mathcal{L}(R))^{*}$ and
$w\in\mathcal{L}(R^{{}^{\prime}})$.

###  A.2 Proof of Proposition 3.7

######  Proposition.

If $S=t_{1},t_{2},\cdots$ is a stream of domain elements from an effective
Boolean algebra $\mathcal{A}=(\mathcal{D}$, $\Psi$, $\llbracket\\_\rrbracket$,
$\bot$, $\top$, $\vee$, $\wedge$, $\neg$), where $t_{i}\in\mathcal{D}$, and
$R$ is a symbolic regular expression over the same algebra, then, for every
$S_{m..k}$, $S_{m..k}\in\mathcal{L}(R)$ iff $S_{1..k}\in\mathcal{L}(R_{s})$
(and $S_{1..k}\in\mathcal{L}(M_{R_{s}})$).

######  Proof A.2.

First, assume that $S_{m..k}\in\mathcal{L}(R)$ for some $m,1\leq m\leq k$ (we
set $S_{1..0}=\epsilon$). Then, for $S_{1..k}=S_{1..(m-1)}\cdot S_{m..k}$,
$S_{1..(m-1)}\in\mathcal{L}(\top^{*})$, since $\top^{*}$ accepts every string
(sub-stream), including $\epsilon$. We know that $S_{m..k}\in\mathcal{L}(R)$,
thus
$S_{1..k}\in\mathcal{L}(\top^{*})\cdot\mathcal{L}(R)=\mathcal{L}(\top^{*}\cdot
R)=\mathcal{L}(R_{s})$. Conversely, assume that
$S_{1..k}\in\mathcal{L}(R_{s})$. Therefore,
$S_{1..k}\in\mathcal{L}(\top^{*}\cdot
R)=\mathcal{L}(\top^{*})\cdot\mathcal{L}(R)$. As a result, $S_{1..k}$ may be
split as $S_{1..k}=S_{1..(m-1)}\cdot S_{m..k}$ such that
$S_{1..(m-1)}\in\mathcal{L}(\top^{*})$ and $S_{m..k}\in\mathcal{L}(R)$. Note
that $S_{1..(m-1)}=\epsilon$ is also possible, in which case the result still
holds, since $\epsilon\in\mathcal{L}(\top^{*})$.

###  A.3 Proof of Theorem 4.11

######  Theorem.

Let $\boldsymbol{\Pi}$ be the transition probability matrix of a homogeneous
Markov chain $Y_{t}$ in the form of Equation (1) and $\boldsymbol{\xi}_{init}$
its initial state distribution. The probability for the time index $n$ when
the system first enters the set of states $F$, starting from a state in $F$,
can be obtained from

$P(Y_{n}\in F,Y_{n-1}\notin F,\cdots,Y_{1}\in
F\mid\boldsymbol{\xi_{init}})=\begin{cases}\boldsymbol{\xi_{F}}^{T}\boldsymbol{F}\boldsymbol{1}&\quad\text{if
}n=2\\\
\boldsymbol{\xi_{F}}^{T}\boldsymbol{F_{N}}\boldsymbol{N}^{n-2}(\boldsymbol{I}-\boldsymbol{N})\boldsymbol{1}&\quad\text{otherwise}\\\
\end{cases}$

where $\xi_{F}$ is the vector consisting of the elements of $\xi_{init}$
corresponding to the states of $F$.

######  Proof A.3.

Case where $n=2$.
In this case, we are in a state $i\in F$ and we take a transition that leads
us back to $F$ again. Therefore, $P(Y_{2}\in F,Y_{1}=i\in
F\mid\boldsymbol{\xi_{init}})=\boldsymbol{\xi}(i)\sum_{j\in F}\pi_{ij}$, i.e.,
we first take the probability of starting in $i$ and multiply it by the sum of
all transitions from $i$ that lead us back to $F$. This result folds for a
certain state $i\in F$. If we start in any state of $F$, $P(Y_{2}\in
F,Y_{1}\in F\mid\boldsymbol{\xi_{init}})=\sum_{i\in
F}\boldsymbol{\xi}(i)\sum_{j\in F}\pi_{ij}$. In matrix notation, this is
equivalent to $P(Y_{2}\in F,Y_{1}\in
F\mid\boldsymbol{\xi_{init}})=\boldsymbol{\xi_{F}}^{T}\boldsymbol{F}\boldsymbol{1}$.

Case where $n>2$.
In this case, we must necessarily first take a transition from $i\in F$ to
$j\in N$, then, for multiple transitions we remain in $N$ and we finally take
a last transition from $N$ to $F$. We can write

$\displaystyle P(Y_{n}\in F,Y_{n-1}\notin F,...,Y_{1}\in
F\mid\boldsymbol{\xi_{init}})=$ $\displaystyle P(Y_{n}\in F,Y_{n-1}\notin
F,...,Y_{2}\notin F\mid\boldsymbol{\xi^{{}^{\prime}}_{N}})$ (15)
$\displaystyle=$ $\displaystyle P(Y_{n-1}\in F,Y_{n-2}\notin F,...,Y_{1}\notin
F\mid\boldsymbol{\xi^{{}^{\prime}}_{N}})$

where $\boldsymbol{\xi^{{}^{\prime}}_{N}}$ is the state distribution (on
states of $N$) after having taken the first transition from $F$ to $N$. This
is given by
$\boldsymbol{\xi^{{}^{\prime}}_{N}}=\boldsymbol{\xi_{F}}^{T}\boldsymbol{F_{N}}$.
By using this as an initial state distribution in Eq. 2 and running the index
$n$ from $1$ to $n-1$, as in Eq. 15, we get

$P(Y_{n}\in F,Y_{n-1}\notin F,...,Y_{1}\in
F\mid\boldsymbol{\xi_{init}})=\boldsymbol{\xi_{F}}^{T}\boldsymbol{F_{N}}\boldsymbol{N}^{n-2}(\boldsymbol{I}-\boldsymbol{N})\boldsymbol{1}$

###  A.4 Proof of Correctness for Algorithm 2

######  Proposition.

Let $P$ be a waiting-time distribution with horizon $h$ and let
$\theta_{fc}<1.0$ be a confidence threshold. Algorithm 2 correctly finds the
smallest interval whose probability exceeds $\theta_{fc}$.

######  Proof A.4 (Proof of correctness for Algorithm 2).

Figure 20: One iteration of Algorithm 2. After the $k^{th}$ iteration, $(i,j)$
is either the smallest interval or there already exists a smaller one. In the
$(k+1)^{th}$ iteration, we need to check intervals for which
$j_{k}<j^{\prime}\leq j_{k+1}$, .i.e., whose right limit is greater than
$j_{k}$ and smaller (but including) $j_{k+1}$. Out of those intervals, the
ones that start before $i_{k}$ need not be checked because they are greater
than $(i_{k},j_{k})$. The ones that start after $i_{k+1}$ are also ignored,
since we already know, from the shrinking phase, that their super-interval
$(i_{k+1}+1,j_{k+1})$ does not exceed $\theta_{fc}$. Finally, we can also
ignore the intervals that start between $i_{k}$ and $i_{k+1}$ (not including
$i_{k}$) and end before $j_{k+1}$ (not including $j_{k+1}$), because they
cannot exceed $\theta_{fc}$. If such an interval actually existed, then the
expansion phase would have stopped before $j_{k+1}$.

We only need to prove the loop invariant. Assume that after the $k^{th}$
iteration of the outer while loop $i=i_{k}$ and $j=j_{k}$ and that after the
$(k+1)^{th}$ iteration $i=i_{k+1}$ and $j=j_{k+1}$. If the invariant holds
after the $k_{th}$ iteration, then all intervals with $e\leq j_{k}$ have been
checked and we know that $(s,e)$ is the best interval up to $j_{k}$. It can be
shown that, during the $(k+1)^{th}$ iteration, the intervals up to $j_{k+1}$
that are not explicitly checked are intervals which cannot possibly exceed
$\theta_{fc}$ or cannot be better than the currently held best interval
$(s,e)$. There are three such sets of unchecked intervals (see also Figure
20):

  * •

All intervals $(i^{\prime},j^{\prime})$ such that $i^{\prime}<i_{k}$ and
$j_{k}\leq j^{\prime}\leq j_{k+1}$, i.e., we ignore all intervals that start
before $i_{k}$. Even if these intervals exceed $\theta_{fc}$, they cannot
possibly be smaller than $(s,e)$, since we know that $(s,e)=(i_{k},j_{k})$ or
that $(s,e)$ is even smaller than $(i_{k},j_{k})$.

  * •

All intervals $(i^{\prime},j^{\prime})$ such that $i^{\prime}>i_{k+1}$ and
$j_{k}\leq j^{\prime}\leq j_{k+1}$, i.e., we ignore all intervals that start
after $i_{k+1}$. These intervals cannot possibly exceed $\theta_{fc}$, since
$(i_{k+1}+1,j_{k+1})$ is below $\theta_{fc}$ and all these intervals are sub-
intervals of $(i_{k+1}+1,j_{k+1})$.

  * •

We are thus left with intervals $(i^{\prime},j^{\prime})$ such that $i_{k}\leq
i^{\prime}\leq i_{k+1}$ and $j_{k}\leq j^{\prime}\leq j_{k+1}$. Out of all the
interval that can be constructed from combining these $i^{\prime}$ and
$j^{\prime}$, the algorithm checks the intervals
$(i^{\prime}=i_{k},j^{\prime})$ and $(i^{\prime},j^{\prime}=j_{k+1})$. The
intervals that are thus left unchecked are the intervals
$(i^{\prime},j^{\prime})$ such that $i_{k}<i^{\prime}\leq i_{k+1}$ and
$j_{k}\leq j^{\prime}<j_{k+1}$. The question is: is it possible for such an
interval to exceed $\theta_{fc}$. The answer is negative. Assume that there is
such an interval $(i^{\prime},j^{\prime})$. If this were the case, then the
algorithm, during its expansion phase, would have stopped at $j^{\prime}$,
because $(i_{k},j^{\prime})$ would exceed $\theta_{fc}$. Therefore, these
intervals cannot exceed $\theta_{fc}$.

A similar reasoning allows us to show that the loop invariant holds after the
first iteration. It thus holds after every iteration.

###  A.5 Proofs of Complexity Results

####  A.5.1 Proof of Proposition 5.1

######  Proposition (Step 1 in Figure 9).

Let $S_{1..k}$ be a stream and $m$ the maximum depth of the Counter Suffix
Tree $T$ to be constructed from $S_{1..k}$. The complexity of constructing $T$
is $O(m(k-m))$.

######  Proof A.5.

There are three operations that affect the cost: incrementing the counter of a
node by $1$, with constant cost $i$; inserting a new node, with constant cost
$n$; visiting an existing node with constant cost $v$; We assume that $n>v$.
For every $S_{l-m+1..l}$, $m\leq l\leq k$ of length $m$, there will be $m$
increment operations and $m$ nodes will be “touched”, i.e., either visited if
already existing or created. Therefore, the total number of increment
operations is $(k-m+1)m=km-m^{2}+m=m(k-m)+m$. The same result applies for the
number of node “touches”. It is always true that $m<k$ and typically $m\ll k$.
Therefore, the cost of increments is $O(m(k-m))$ and the cost of
visits/creations is also $O(m(k-m))$. Thus, the total cost is
$O(m(k-m))+O(m(k-m))=O(m(k-m))$. In fact, the worst case is when all
$S_{l-m+1..l}$ are different and have no common suffixes. In this case, there
are no visits to existing nodes, but only insertions, which are more expensive
than visits. Their cost would again be $O(nm(k-m))=O(m(k-m))$, ignoring the
constant $n$.

####  A.5.2 Proof of Proposition 5.3

######  Proposition (Step 3a in Figure 9).

Let $T$ be a $\mathit{PST}$ of maximum depth $m$, learned with the $t$
minterms of a $\mathit{DSFA}$ $M_{R}$. The complexity of constructing a
$\mathit{PSA}$ $M_{S}$ from $T$ is $O(t^{m+1}\cdot m)$.

######  Proof A.6.

We assume that the cost of creating new states and transitions for $M_{S}$ is
constant. In the worst case, all possible suffixes of length $m$ have to be
added to $T$ as leaves. $T$ will thus have $t^{m}$ leaves. The main idea of
the algorithm for converting a $\mathit{PST}$ $T$ to a $\mathit{PSA}$ $M_{S}$
is to use the leaves of $T$ as states of $M_{S}$ and for every symbol
(minterm) $\sigma$ find the next state/leaf and set the transition probability
to be equal to the probability of $\sigma$ from the source leaf. If we assume
that the cost of accessing a leaf is constant (e.g., by keeping separate
pointers to the leaves), the cost for constructing $M_{S}$ is dominated by the
cost of constructing the $k^{m}$ states of $M_{S}$ and the $t$ transitions
from each such state. For each transition, finding the next state requires
traversing a path of length $m$ in $T$. The total cost is thus $O(t^{m}\cdot
t\cdot m)$ = $O(t^{m+1}\cdot m)$.

####  A.5.3 Proof of Proposition 5.7

######  Proposition (Step 4 in Figure 9).

Let $M_{R}$ be a $\mathit{DSFA}$ with $t$ minterms and $M_{S}$ a
$\mathit{PSA}$ learned with the minterms of $M_{R}$. The complexity of
constructing an embedding $M$ of $M_{S}$ in $M_{S}$ with Algorithm 1 is
$O(t\cdot\lvert M_{R}.Q\times M_{S}.Q\rvert)$.

######  Proof A.7.

We assume that the cost of constructing new states and transitions for $M$ is
constant. We also assume that the cost of finding a given state in both
$M_{R}$ and $M_{S}$ is constant, e.g., by using a linked data structure for
representing the automaton with a hash table on its states (or an array), and
the cost of finding the next state from a given state is also constant. In the
worst case, even with the incremental algorithm 1, we would need to create the
full Cartesian product $M_{R}.Q\times M_{S}.Q$ to get the states of $M$. For
each of these states, we would need to find the states of $M_{R}$ and $M_{S}$
from which it will be composed and to create $t$ outgoing transitions.
Therefore, the complexity of creating $M$ would be $O(t\cdot\lvert
M_{R}.Q\times M_{S}.Q\rvert)$.

####  A.5.4 Proof of Proposition 5.9

######  Proposition (Step 5 in Figure 9).

Let $M$ be the embedding of a $\mathit{PSA}$ $M_{S}$ in a $\mathit{DSFA}$
$M_{R}$. The complexity of estimating the waiting-time distribution for a
state of $M$ and a horizon of length $h$ using Theorem 4.10 is
$O((h-1)k^{2.37})$ where $k$ is the dimension of the square matrix
$\boldsymbol{N}$.

######  Proof A.8.

We want to use Equation 2 to estimate the distribution of a state. The
equation is repeated below:

$P(Y_{n}\in F,Y_{n-1}\notin F,...,Y_{1}\notin
F\mid\boldsymbol{\xi_{init}})=\boldsymbol{\xi_{N}}^{T}\boldsymbol{N}^{n-1}(\boldsymbol{I}-\boldsymbol{N})\boldsymbol{1}$

We want to estimate the distribution for the $h$ points of the horizon, i.e.,
for $n=2$, $n=3$ up to $n=h+1$. For $n=2$, we have

$P(Y_{2}\in F,Y_{1}\notin
F\mid\boldsymbol{\xi_{init}})=\boldsymbol{\xi_{N}}^{T}\boldsymbol{N}(\boldsymbol{I}-\boldsymbol{N})\boldsymbol{1}$

For $n=3$, we have

$P(Y_{3}\in F,Y_{2}\notin F,Y_{1}\notin
F\mid\boldsymbol{\xi_{init}})=\boldsymbol{\xi_{N}}^{T}\boldsymbol{N}^{2}(\boldsymbol{I}-\boldsymbol{N})\boldsymbol{1}$

In general, for $n=i$, we can use the power of $\boldsymbol{N}$ that we have
estimated in the previous step for $n=i-1$, $\boldsymbol{N}^{i-2}$, in order
to estimate the next power $\boldsymbol{N}^{i-1}$ via a multiplication by
$\boldsymbol{N}$ so as to avoid estimating this power from scratch. Then
$\boldsymbol{N}^{i-1}$ can be multiplied by
$(\boldsymbol{I}-\boldsymbol{N})\boldsymbol{1}$, which remains fixed for all
$i$ and can thus be calculated only once.

The cost of estimating $(\boldsymbol{I}-\boldsymbol{N})$ is $k^{2}$ due to the
$k^{2}$ subtractions. Multiplying the matrix $(\boldsymbol{I}-\boldsymbol{N})$
by the vector $\boldsymbol{1}$ results in a new vector with $k$ elements. Each
of these elements requires $k$ multiplications and $k-1$ additions or $2k-1$
operations. Thus, the estimation of
$(\boldsymbol{I}-\boldsymbol{N})\boldsymbol{1}$ has a cost of
$k^{2}+k(2k-1)=3k^{2}-k$.

Now, for $n=i$, estimating the power $\boldsymbol{N}^{i-1}$ from
$\boldsymbol{N}^{i-2}$ has a cost of $k^{2.37}$ using an efficient
multiplication algorithm such as the Coppersmith–Winograd algorithm [18] or
the improvement proposed by Stothers [50].

Additionally, $\boldsymbol{N}^{i-1}$ must then be multiplied by the vector
$(\boldsymbol{I}-\boldsymbol{N})\boldsymbol{1}$, with a cost of $2k^{2}-k$,
resulting in a new vector with $k$ elements.

This vector must then be multiplied by $\boldsymbol{\xi_{N}}^{T}$ to produce
the final probability value with a cost of $2k-1$ for the $k$ multiplications
and the $k-1$ additions.

We thus have a fixed initial cost of $3k^{2}-k$ and then for every iteration
$i$ a cost of
$k^{2.37}I_{\\{i>2\\}}+2k^{2}-k+2k-1=k^{2.37}I_{\\{i>2\\}}+2k^{2}+k-1$, where
$I_{\\{i>2\\}}$ is an indicator function ($1$ for $i>2$ and $0$ otherwise).
Note that the cost $k^{2.37}$ is not included for $i=2$ because in this case
we do not need to raise $\boldsymbol{N}$ to a power. The total cost would thus
be:

$\displaystyle 3k^{2}-k+$ $\displaystyle k^{2.37}\cdot 0+2k^{2}+k-1$
$\displaystyle\text{for }n=2$ $\displaystyle+$ $\displaystyle k^{2.37}\cdot
1+2k^{2}+k-1$ $\displaystyle\text{for }n=3$ $\displaystyle\cdots$
$\displaystyle+$ $\displaystyle k^{2.37}\cdot 1+2k^{2}+k-1$
$\displaystyle\text{for }n=h+1$ $\displaystyle=$
$\displaystyle(h-1)k^{2.37}+(2h+3)k^{2}+(h-1)k-h$
$\displaystyle=O((h-1)k^{2.37})$

####  A.5.5 Proof of Proposition 5.11

######  Proposition (Step 6 in Figure 9).

For a waiting-time distribution with a horizon of length $h$, the complexity
of finding the smallest interval that exceeds a confidence threshold
$\theta_{fc}$ with Algorithm 2 is $O(h)$.

######  Proof A.9.

Indexes $i$ and $j$ of Algorithm 2 scan the distribution only once. The cost
for $j$ is the cost of $h$ points of the distribution that need to be accessed
plus $h-1$ additions. Similarly, the cost for $i$ is the cost of (at most) $h$
accessed points plus the cost of (at most) $h-1$ subtractions. Thus the total
cost is $O(h)$.

####  A.5.6 Proof of Proposition 5.5

######  Proposition (Step 3b in Figure 9).

Let $T$ be a $\mathit{PST}$ of maximum depth $m$, learned with the $t$
minterms of a $\mathit{DSFA}$ $M_{R}$. The complexity of estimating the
waiting-time distribution for a state of $M_{R}$ and a horizon of length $h$
directly from $T$ is $O((m+3)\frac{t-t^{h+1}}{1-t})$.

######  Proof A.10.

After every new event arrival, we first have to construct the tree of future
states, as shown in Figure 7(b). In the worst case, no paths can be pruned and
the tree has to be expanded until level $h$. The total number of nodes that
have to be created is thus a geometric progress:
$t+t^{2}+\cdots+t^{h}=\sum_{i=1}^{h}t^{i}=\frac{t-t^{h+1}}{1-t}$. Assuming
that it takes constant time to create a new node, this formula gives the cost
of creating the nodes of the trees. Another cost that is involved concerns the
time required to find the proper leaf of the $\mathit{PST}$ $T$ before the
creation of each new node. In the worst case, all leaves will be at level $m$.
The cost of each search will thus be $m$. The total search cost for all nodes
will be $mt+mt^{2}+\cdots+mt^{h}=\sum_{i=1}^{h}mt^{i}=m\frac{t-t^{h+1}}{1-t}$.
The total cost (node creation and search) for constructing the tree is
$\frac{t-t^{h+1}}{1-t}+m\frac{t-t^{h+1}}{1-t}=(m+1)\frac{t-t^{h+1}}{1-t}$.
With the tree of future states at hand, we can now estimate the waiting-time
distribution. In the worst case, the tree will be fully expanded and we will
have to access all its paths until level $h$. We will first have to visit the
$t$ nodes of level 1, then the $t^{2}$ nodes of level 2, etc. The access cost
will thus be $t+t^{2}+\cdots+t^{h}=\sum_{i=1}^{h}t^{i}=\frac{t-t^{h+1}}{1-t}$.
We also need to take into account the cost of estimating the probability of
each node. For each node, one multiplication is required, assuming that we
store partial products and do not have to traverse the whole path to a node to
estimate its probability. As a result, the number of multiplications will also
be $\frac{t-t^{h+1}}{1-t}$. The total cost (path traversal and
multiplications) will thus be $2\frac{t-t^{h+1}}{1-t}$, where we ignore the
cost of summing the probabilities of final states, assuming it is constant. By
adding the cost of constructing the tree ($(m+1)\frac{t-t^{h+1}}{1-t}$) and
the cost of estimating the distribution ($2\frac{t-t^{h+1}}{1-t}$), we get a
complexity of $O((m+3)\frac{t-t^{h+1}}{1-t})$.