11institutetext: Univ. Orléans, INSA Centre Val de Loire, LIFO EA 4022,
F-45067, Orléans, France 22institutetext: Worldline
22email<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>

# Random Dilated Shapelet Transform: A New Approach for Time Series Shapelets

Antoine Guillaume 1122 0000-0001-8975-982X Christel Vrain 11
0000-0003-3307-0753 Wael Elloumi 22 0000-0002-5886-2450

###### Abstract

Shapelet-based algorithms are widely used for time series classification
because of their ease of interpretation, but they are currently outperformed
by recent state-of-the-art approaches. We present a new formulation of time
series shapelets including the notion of dilation, and we introduce a new
shapelet feature to enhance their discriminative power for classification.
Experiments performed on 112 datasets show that our method improves on the
state-of-the-art shapelet algorithm, and achieves comparable accuracy to
recent state-of-the-art approaches, without sacrificing neither scalability,
nor interpretability.

###### Keywords:

Time series Shapelets Classification

##  1 Introduction

Time series occur in a multitude of domains, covering a wide range of
applications, which have impacts in many parts of society. The ever-increasing
quantity of data and the publications of laws regarding models
interpretability are setting new constraints for applications across
industries.

Recent research in time series classification produced highly accurate
classifiers, using either deep learning approaches [10], or meta-ensemble
methods [13, 18]. Despite being the most accurate approaches, they are among
the slowest, which make them hard to apply on use-cases with huge amount of
data. Other methods, based on random approaches, notably the RandOm
Convolutional KErnel Transform (ROCKET) [4], achieve comparable accuracy with
extreme scalability. Even though recent works on post-hoc methods and specific
frameworks [5] improved the interpretability of those approaches, they lack a
"by design" interpretability.

On the other hand, time series shapelets [25] have been widely used in time
series classification for their ease of interpretation, which is a critical
aspect to some application domains such as health and security. The downside
is that shapelet algorithms are often outperformed by recent approaches, both
in terms of accuracy and scalability. Most Shapelet approaches tried to solve
the scalability issues at the expense of some classification accuracy, notably
through the use of symbolic approximation techniques [11], while others used
random shapelets [24]. Recently, a symbolic sequence ensemble learning [19]
method was proposed, which improved the predictive power of approximation-
based methods, while other work focused on finding a new discriminative
feature [7] to consider during the extraction process.

In this work, we present the Random Dilated Shapelet Transform, an adaptation
of time series shapelets that includes the notion of dilation, one of the core
mechanism of the success of convolutional kernel approaches. We also extend on
the work of [7] and introduce a new feature to enhance the discriminative
power of shapelets. Our contributions can be summarized as follows:

  * •

an adaptation of time series shapelets allowing the use of dilation, and a
feature to capture a new discriminative property of shapelets,

  * •

an interpretable, scalable and accurate shapelet algorithm, which allows
shapelet based algorithm to catch-up with the state-of-the-art,

  * •

an experimental study about the sensitivity of our method parameters and a
comparative study against the state-of-the-art algorithms for time series
classification.

##  2 Background

In this section, we present a brief review of time series classification, and
make a focus on shapelet methods. From now on, we use calligraphic font
(${\cal X}$) to denote a collection of elements (e.g. a set of time series),
capital ($X$) for one element (e.g. a time series), and lowercase ($x$) for a
value of this element. In this work we consider supervised classification: the
ensemble of input time series will be denoted by ${\cal
X}=\\{X_{1},...,X_{n}\\}$ with $X_{i}=\\{x_{1},...,x_{m}\\}$ a time series and
$Y=\\{y_{1},...,y_{n}\\}$ their respective classes.

###  2.1 Time series classification

We present a brief overview of the algorithms identified as state-of-the-art
and used in our experimental section, and we report the reader to a recent
review [1] for a more detailed view of the field.

  * •

Shapelet Transform Classifier (STC) [2], is regarded as a state of the art for
shapelet algorithms in terms of accuracy. This algorithm iteratively
initializes new shapelets, assesses their discriminative power, and removes
those that are too similar. The goal being to maximize the discriminative
power of an ensemble of shapelets. A Rotation Forest is then applied as a
classifier.

  * •

Temporal Dictionary Ensemble (TDE) [17] is an ensemble of dictionary-based
classifiers. It uses some variants of the BOSS classifier [20] and WEASEL
[21], as base estimators and optimizes their parameters through a Gaussian
process.

  * •

Diverse Representation Canonical Interval Forest Classifier (DrCIF) [18], is
an extension of the CIF algorithm [16]. After selecting random intervals from
different representations of the time series, it uses the Catch22 [14] method
to extract a feature matrix.

  * •

RandOm Convolutional KErnel Transform (ROCKET) [4], randomly generates a huge
set of convolutional kernels, and extracts as features the maximum and the
proportion of positive values of the convolution for each time series and each
kernel. It is followed by a Ridge classifier. An ensemble version of this
method, called ARSENAL, was introduced by [18].

  * •

Inception-Time [10] is an ensemble of Inception networks, which introduce
Inception modules as replacement for traditional fully convolutional layers,
notably to mitigate the vanishing gradient problem.

  * •

Hierarchical Vote Collective of Transformation-based Ensembles (HC1) [13]is a
meta-ensemble classifier using a variety of time series classifiers, such as
STC, with a novel ensemble learning scheme which estimate the weight of each
base classifier in the final decision. This estimation is based on performance
in a 10-fold validation scheme.

Variants of this method were developed, such as TS-CHIEF [22] and HC2 [18],
that both modify the set of base classifiers to improve accuracy. HC2 also
modified the meta-ensemble procedure, using a Out-Of-Bag estimate instead of a
10-fold validation to estimate the performance of each base classifier, which
improved scalability compared to HC1.

###  2.2 Shapelets

Shapelets [25] were originally defined as time series subsequences
representative of class membership. In the following, we define a shapelet $S$
as a vector $S=\\{s_{1},...,s_{l}\\}$ with $l$ its length. All shapelet-based
algorithms have the same protocol to extract features from a shapelet $S$ and
a time series $X=\\{x_{1},...,x_{m}\\}$, by using a distance vector
$f(S,X)=\\{f_{1},...,f_{m-(l-1)}\\}$ defined as :

$f_{i}=\sqrt{\sum_{j=1}^{l}(X_{i+(j-1)}-s_{j})^{2}}$ (1)

In this definition, a point $f_{i}$ is simply the Euclidean distance between
$S$ and the subsequence of length $l$ starting at index $i$ in $X$. The
minimum value of $f(S,X)$ is then extracted as a feature, which can be
interpreted as an indicator of the presence of the pattern represented by $S$
in $X$. A popular variant of this distance function consists in using a
z-normalized Euclidean distance, where $S$ and all subsequences of $X$ are
z-normalized independently, allowing to add scale invariance to the
translation invariance of the initial formulation. Then, as presented in the
Shapelet Transform [12], by using a set of shapelets $\cal S$, one can then
transform an ensemble of time series $\cal X$ into a feature matrix of shape
$(|{\cal X}|,|{\cal S}|)$, and use it as input in a non-temporal classifier
such as a decision tree.

The step of generating and selecting shapelet candidates is the main
difference between most approaches. In order to speed up the exhaustive
search, Fast Shapelet [11] use input discretization, while Ultra Fast Shapelet
[24] use random sampling. FLAG [8] build shapelets location indicators from
the data to reduce the set of admissible candidates, and GENDIS [23] use an
evolutionary algorithm initialized by a clustering on the set of possible
candidates. Learning Time Series Shapelet [6] use a gradient-descent
optimization that iteratively change the values of a set of shapelets. MrSEQL
[19], while not strictly speaking a shapelet algorithm, searches for
discriminative symbolic sequences in a variety of symbolic representations of
the inputs.

Since the publication of Localized Random Shapelet (LRS) [7], which showed the
benefit of extracting $argmin\ d(S,X)$ to discriminate time series based on
the location of the minimum between $S$ and $X$, it has been included in most
recent approaches. Based on their results, we will also use this feature in
our method.

##  3 Proposed method

In this section, we introduce the main components of our method: the use of
dilation in the shapelet formulation and the features extracted from the
distance vector between a shapelet and a time series. We put emphasis on the
dilation and on the Shapelet Occurrence feature that are new contributions to
shapelet algorithms. We give some simple visual examples to illustrate these
notions, and report the visualization on real data to the experimental
section.

###  3.1 Dilated Shapelets

To introduce the notion of dilation in shapelets, we define now a shapelet $S$
as $S=\\{\\{v_{1},...,v_{l}\\},d\\}$ with $l$ the length parameter and $d$ the
dilation parameter. In practice, the dilation is used in the distance function
$f$, where each value of the shapelet will be compared to a dilated
subsequence of the input time series. More formally, consider a time series
$X=\\{x_{1},...,x_{m}\\}$ and a dilated shapelet $S$, we now define
$f(S,X)=\\{f_{1},...,f_{m-(l-1)\times d}\\}$ as :

$f_{i}=\sqrt{\sum_{j=1}^{l}(X_{i+(j-1)\times d}-s_{j})^{2}}$ (2)

The interest of using dilation in shapelets is to make them non-contiguous
subsequences. It allows a shapelet to either match a non-contiguous pattern,
or a contiguous one, by focusing on key points of the pattern without covering
it entirely, as illustrated in Figure 1. Note that this formulation is
equivalent to the original shapelet formulation when $d=1$.

Figure 1: An example of two possible shapelets (in orange), positioned on a
synthetic pattern (in blue): (a) one without dilation, and (b) a much smaller
one but with dilation

###  3.2 Shapelet Occurrence feature

If we consider a shapelet $S$ and two time series $X_{1}$ and $X_{2}$, we can
imagine multiple ways of discriminating $X_{1}$ and $X_{2}$ using $S$.

  * •

$S$ can be present in $X_{1}$ but not in $X_{2}$. This is captured by $min\
f(S,X_{i})$, with smaller distances indicating better matches between the
shapelet and the series.

  * •

$S$ can be present in both series, but not at the same place. This is captured
by the $argmin$ feature introduced by LRS [7].

  * •

$S$ can be present in both series, but not at the same scale. In this case, a
normalized distance would not be able to discriminate the series.

  * •

$S$ can be present in both series, but occurs a different number of times in
$X_{1}$ compared to $X_{2}$. This is captured by a new feature, called
Shapelet Occurrence (SO).

Those points are illustrated in Figure 2. Deciding whether scaling is
important or not is highly dependent on the application, but without prior
knowledge, one cannot know which to choose. For this reason, we introduce a
parameter in Section 3.3 allowing to tune the amount of normalized shapelets.

Figure 2: Synthetic examples of possible discriminative properties between two
classes

To the best of our knowledge, the number of occurrences of a shapelet has
never been considered as a feature. This requires another modification to the
definition of $S$ as $S=\\{\\{v_{1},...,v_{l}\\},d,\lambda\\}$, with $\lambda$
a threshold allowing us to compute the Shapelet Occurrence ($SO$) feature as
$SO=|\\{i|f(S,X)_{i}<\lambda\\}|$. Although the parameter $\lambda$ could be
set randomly, we discuss in Section 3.3 a method to set the value of this
threshold.

###  3.3 Random Dilated Shapelet Transform (RDST)

Our objective for this algorithm is to produce an accurate but scalable
approach. As our shapelet formulation adds attributes compared to the initial
formulation [25], optimizing a set of dilated shapelets with a threshold
$\lambda$ will be costly, and this explains why we choose a random approach.

For simplicity, we present our approach in the context of univariate and even
length time series, with ${\cal X}=\\{X_{1},...,X_{n}\\}$ a set of time series
($X_{i}=\\{x_{1},...,x_{m}\\}$) and $Y=\\{y_{1},...,y_{n}\\}$ their respective
classes. Our method takes as input four parameters that are: $n\\_shapelets$
the number of shapelets to generate, $L$ a set of possible lengths for the
shapelets, $p\\_norm$ the proportion of shapelets that will use
z-normalization, and $(P_{1},P_{2})\in[0,100]$ a pair used as percentile
bounds for the sampling of the threshold $\lambda$.

Given the definition of a shapelet $S=\\{\\{v_{1},...,v_{l}\\},d,\lambda\\}$,
we initialize each parameter as follows:

  * •

the length $l$ is uniformly drawn from $L$,

  * •

the dilation $d$, in the same way as ROCKET [4], is set to $d=\left\lfloor
2^{x}\right\rfloor$ with $x$ uniformly drawn in $[0,log_{2}\frac{m}{l}]$,

  * •

we randomly choose whether the shapelet will use a z-normalized distance with
probability $p\\_norm$,

  * •

for setting the values, a sample $X$ is uniformly drawn from ${\cal X}$, and
an admissible start point $i$ (given $l,d$) is randomly selected. Then, values
are set to $[X_{i},...,X_{i+(l-1)\times d}]$.

  * •

finally, given a shapelet $S$, to fix the value of $\lambda$, we take a sample
$X$ from the same class as the one used for extracting the shapelet value, and
uniformly draw a value between the two percentiles $(P_{1},P_{2})$ of
$f(S,X)$.

The strategy employed to find $\lambda$ is a classic trade-off between time
and accuracy. If scalability was not a focus, we could compute the distance
vector for more samples in ${\cal X}$, and optimize the value of $\lambda$
based on an information measure. After computing the distance vector between
all pairs of time series and shapelets, the output of our method is a feature
matrix of size $(|{\cal X}|,3\times n\\_shapelets)$, with the three features
extracted from the distance vector $f(S,X)$ being the $min,argmin,SO(S,X)$.

Following the arguments of the authors of ROCKET [4], we use a Ridge
Classifier after the transformation of ${\cal X}$, as the L2 regularization
used in Ridge is of critical importance due to the high number of features
that are generated, while being scalable and interpretable.

##  4 Experiments

Our focus in this section is to study the influence of the four parameters of
our method on classification accuracy, as well as comparing its performance to
recent state-of-the-art approaches. All the experiments were run on a DELL
PowerEdge R730 on Debian 9 with 2 XEON E5-2630 Corei7 (92 cores) and 64GB of
RAM. We provide a python package111https://github.com/baraline/convst using
community standards to run the method and the interpretability tool on any
dataset, along with all result tables, and reproducibility instructions for
our experiments.

In the following, we use the 112 univariate datasets from the UCR archive [3]
and when comparing to state-of-the-art results, we use the same resamples
scheme as the one used in their experiments. We use critical difference
diagrams to display the mean ranks of objects, with cliques (formed by
horizontal bars) computed using the Wilcoxon-Holm post-hoc analysis [9], with
a p-value of $0.05$. A clique indicates that the accuracy difference between
objects is not statistically significant.

###  4.1 Sensitivity Analysis

We conduct a sensitivity analysis on the four input parameters of our
algorithm and their effect on classification accuracy on 40 datasets selected
randomly, with raw results and selected datasets in a specific file in the
online repository. For each parameter analysis, all other parameters remain
fixed at the following default values : $n\\_shapelets=10000$, $p\\_norm=0.9$,
$L=[7,9,11]$, $P1=5$ and $P2=15$. Figure 3 and Figure 4 give the mean accuracy
ranks of each method over the 40 datasets, with the accuracy of each method
and each dataset computed as the mean of the same 10 resamples. Given the
tested set of values, the most impactful parameter is the number of shapelets,
with a noticeable increase in performance above 10000 shapelets. All other
parameters only display minor gains and thus seem to be stable. Based on those
results, for all further experiments we set as default parameters
$n\\_shapelets=10000$, $p\\_norm=0.8$, $L=[11]$ and $P1=5$, $P2=10$, and
report results for datasets used in sensitivity analysis and the others.

Figure 3: Accuracy ranks for (a) different number of shapelets, and (b)
different shapelet lengths Figure 4: Accuracy ranks for (a) different
percentiles bounds, and (b) proportion of z-normalized shapelets

###  4.2 Scalability

We perform a comparison of the scalability of our approach against Hive-Cote
1.0 (HC1), Hive-Cote 2.0 (HC2), DrCIF, ROCKET, and the Shapelet Transform
Classifier (STC). Note that when used as a component in HC1 and HC2, STC is by
default subject to a time contract of two hours. Except from this default
configuration in HC1 and HC2, we are not setting any time contract in other
algorithms. Both STC and RDST are by default sampling 10000 shapelets, and
ROCKET use 10000 kernels.

We are aware that the runtime of HC1, HC2 and STC could be reduced with time
contracts. But, as our goal in this section is to contextualize the gain in
classification accuracy against the time complexity of each method, we present
the results with the time contracts used to generate the accuracy results of
the next section.

We use the Crop Dataset and the Rock Dataset of the UCR archive for evaluating
the scalability respectively on the number of time series and their length. As
all competing algorithms implemented in the sktime package of [15] can use
parallel processing, we set each algorithm to use 90 cores. Figure 5 reports
the mean training time over 10 resamples, showing the very competitive
scalability of RDST. Note that due to job time limitation on our machine and
the computational cost of HC2, we could not consider all samples for the Crop
dataset. We report the reader interested in the implementation details of our
algorithm to the web page of the project.

Figure 5: Result of the scalability study of the competing algorithms for
current state-of-the-art, for (a) number of time series and (b) time series
length. Y-axis use log-scale.

###  4.3 Comparative study

We present the results of our comparative study using the mean accuracy over
the same 30 resamples for each of the 112 datasets as HC2 [18] used in their
study, and compare our approach against their experimental result. Figure 6
gives the mean accuracy rank of each method over the 40 datasets used for
setting the defaults parameters in sensitivity analysis, and for the 72
others. The full result tables including standard deviation per dataset and
more visualizations of the results are available online as supplementary
materials.

Figure 6: Mean accuracy ranks of each method for the 40 dataset used in
sensitivity analysis and the 72 others.

Given the scalability and simplicity of our method, having an accuracy
comparable to the prior developments of HC2 and to deep learning approaches is
a very promising result. Notably for future developments where focus would
shift to accuracy rather than scalability. For reference, using RDST without
any distance normalization is equivalent to STC in terms of mean accuracy
rank, with the same protocol as above.

###  4.4 Interpretability

Given a set of $M$ shapelets, RDST generates $3M$ features. Each feature is
linked to a weight for each class in the Ridge classifier, as it is trained in
a one-vs-all fashion. Given a class, we can then visualize either global or
local information. Locally, we can inspect a shapelet to show how it
discriminates the current class, and where the shapelet is positioned with
either training or testing data, as shown in Figure 7. Globally, we can
display the distribution of weights for each feature type ($min$,
$\operatorname*{arg\,min}$ and $SO$) or by shapelet characteristics such as
length, dilation, or use of normalization as shown in Figure 8. While this
only provides a basic interpretation of the results, we believe a more formal
framework could be developed to extract explanations from this data.

Figure 7: The most important shapelet for class 0 of the Coffee dataset,
according to weights of the Ridge classifier, with distribution displayed on
the testing data, and two testing samples for visualization. Figure 8: A
global interpretation of RDST, with (a) distribution of weights for each type
of feature, and (b) distribution of weights per dilation.

##  5 Conclusions and future work

The Random Dilated Shapelet Transform introduces new ways of increasing the
global performance of shapelet algorithms, notably through the use of
dilation, allowing the use of small non-contiguous subsequences as shapelets,
efficiently covering areas of interest in the data. We have shown in our
experiments that this new method improves on the state-of-the-art for shapelet
algorithms with a good scalability compared to most of the approaches. This
work offers many perspectives for future work, notably a generalized version
to process uneven length or multivariate time series, as well as modifications
of the shapelet generation process to better leverage class information. A
more formal explainability framework is also one of our main priorities with
this work, since being able to extract clear and visual explanations for
domain experts is an extremely desirable property.

#### Acknowledgements

This work is supported by the ANRT CIFRE grant n°2019/0281 in partnership with
Worldline and the University of Orléans.

## References

  * [1] Bagnall, A., Lines, J., Bostrom, A., Large, J., Keogh, E.: The great multivariate time series classification bake off: a review and experimental evaluation of recent algorithmic advances. Data Mining and Knowledge Discovery 35, 401–449 (08 2020)
  * [2] Bostrom, A., Bagnall, A.: Binary shapelet transform for multiclass time series classification. In: Madria, S., Hara, T. (eds.) Big Data Analytics and Knowledge Discovery. pp. 257–269. Springer International Publishing, Cham (2015)
  * [3] Dau, H.A., Bagnall, A., Kamgar, K., Yeh, C.C.M., Zhu, Y., Gharghabi, S., Ratanamahatana, C.A., Keogh, E.: The ucr time series archive (2019)
  * [4] Dempster, A., Petitjean, F., Webb, G.: Rocket: exceptionally fast and accurate time series classification using random convolutional kernels. Data Mining and Knowledge Discovery 34 (09 2020)
  * [5] Fauvel, K., Masson, V., Fromont, E.: A Performance-Explainability Framework to Benchmark Machine Learning Methods: Application to Multivariate Time Series Classifiers. In: IJCAI-PRICAI 2020 - Workshop on Explainable Artificial Intelligence (XAI). pp. 1–8. Yokohama, Japan (Jan 2021)
  * [6] Grabocka, J., Schilling, N., Wistuba, M., Schmidt-Thieme, L.: Learning time-series shapelets. In: Proceedings of the 20th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. p. 392–401. KDD ’14, Association for Computing Machinery, New York, NY, USA (2014)
  * [7] Guilleme, M., Malinowski, S., Tavenard, R., Renard, X.: Localized Random Shapelets. In: International Workshop on Advanced Analysis and Learning on Temporal Data. pp. 85–97. Wurzburg, Germany (2019)
  * [8] Hou, L., Kwok, J.T., Zurada, J.M.: Efficient learning of timeseries shapelets. In: Proceedings of the Thirtieth AAAI Conference on Artificial Intelligence. p. 1209–1215. AAAI’16, AAAI Press (2016)
  * [9] Ismail Fawaz, H., Forestier, G., Weber, J., Idoumghar, L., Muller, P.A.: Deep learning for time series classification: a review. Data Mining and Knowledge Discovery 33(4), 917–963 (2019)
  * [10] Ismail Fawaz, H., Lucas, B., Forestier, G., Pelletier, C., Schmidt, D., Weber, J., Webb, G., Idoumghar, L., Muller, P.A., Petitjean, F.: Inceptiontime: Finding alexnet for time series classification. Data Mining and Knowledge Discovery 34, 1–27 (11 2020)
  * [11] Keogh, E.J., Rakthanmanon, T.: Fast shapelets: A scalable algorithm for discovering time series shapelets. In: SDM (2013)
  * [12] Lines, J., Davis, L.M., Hills, J., Bagnall, A.: A shapelet transform for time series classification. In: Proceedings of the 18th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. p. 289–297. KDD ’12, Association for Computing Machinery, New York, NY, USA (2012)
  * [13] Lines, J., Taylor, S., Bagnall, A.: Hive-cote: The hierarchical vote collective of transformation-based ensembles for time series classification. In: 2016 IEEE 16th International Conference on Data Mining (ICDM). pp. 1041–1046 (2016)
  * [14] Lubba, C.H., Sethi, S., Knaute, P., Schultz, S., Fulcher, B., Jones, N.: catch22: Canonical time-series characteristics: Selected through highly comparative time-series analysis. Data Mining and Knowledge Discovery 33 (08 2019)
  * [15] Löning, M., Bagnall, T., Middlehurst, M., Ganesh, S., Oastler, G., Király, F., Lines, J., ViktorKaz, Walter, M., RNKuhns, Rockenschaub, P., Owoseni, T., jesellier, Bulatova, G., Lovkush, Meyer, S.M., AidenRushbrooke, Schäfer, P., oleskiewicz, Xu, Y.X., Ansari, A., Hongyi, Sakshi, A., Arelo, Large, J., Orduz, J., Kumar, B., Gorelli, M.E., Wang, A., Pong, J.: alan-turing-institute/sktime: v0.8.2 (Oct 2021)
  * [16] Middlehurst, M., Large, J., Bagnall, A.: The canonical interval forest (cif) classifier for time series classification. In: 2020 IEEE International Conference on Big Data (Big Data). pp. 188–195 (2020)
  * [17] Middlehurst, M., Large, J., Cawley, G., Bagnall, A.: The temporal dictionary ensemble (tde) classifier for time series classification. In: Hutter, F., Kersting, K., Lijffijt, J., Valera, I. (eds.) Machine Learning and Knowledge Discovery in Databases. pp. 660–676. Springer International Publishing, Cham (2021)
  * [18] Middlehurst, M., Large, J., Flynn, M., Lines, J., Bostrom, A., Bagnall, A.: Hive-cote 2.0: a new meta ensemble for time series classification. Machine Learning 110 (12 2021)
  * [19] Nguyen, T., Gsponer, S., Ilie, I., O’reilly, M., Ifrim, G.: Interpretable time series classification using linear models and multi-resolution multi-domain symbolic representations. Data Min. Knowl. Discov. 33(4), 1183–1222 (jul 2019)
  * [20] Schäfer, P.: The boss is concerned with time series classification in the presence of noise. Data Mining and Knowledge Discovery 29, 1505–1530 (2014)
  * [21] Schäfer, P., Leser, U.: Fast and accurate time series classification with weasel. In: Proceedings of the 2017 ACM on Conference on Information and Knowledge Management. p. 637–646. CIKM ’17, Association for Computing Machinery, New York, NY, USA (2017)
  * [22] Shifaz, A., Pelletier, C., Petitjean, F., Webb, G.: Ts-chief: a scalable and accurate forest algorithm for time series classification. Data Mining and Knowledge Discovery 34 (05 2020)
  * [23] Vandewiele, G., Ongenae, F., De Turck, F.: Gendis: Genetic discovery of shapelets. Sensors 21(4) (2021)
  * [24] Wistuba, M., Grabocka, J., Schmidt-Thieme, L.: Ultra-fast shapelets for time series classification. ArXiv abs/1503.05018 (2015)
  * [25] Ye, L., Keogh, E.: Time series shapelets: A new primitive for data mining. In: Proceedings of the 15th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. p. 947–956. KDD ’09, Association for Computing Machinery, New York, NY, USA (2009)


On Approximating Ruin Probability of Double Stochastic Compound Poisson
Processes111
Amir T. Payandeh Najafabadia,222Corresponding author<EMAIL_ADDRESS>Phone No. +98-21-29903011; Fax No. +98-21-22431649 & Dan Kucerovskyb

a Department of Mathematical Sciences, Shahid Beheshti University, G.C. Evin,
1983963113, Tehran, Iran.

b Department of Mathematics and Statistics, University of New Brunswick,
Fredericton, N.B. CANADA E3B 5A3.

Abstract

Consider a surplus process which both of collected premium and payed claim
size are two independent compound Poisson processes. This article derives two
approximated formulas for the ruin probability of such surplus process, say
double stochastic compound poisson process. More precisely, it provides two
mixture exponential approximations for ruin probability of such double
stochastic compound poisson process. Applications to long–term Bonus–Malus
systems and a heavy-tiled claim size distribution have been given. Improvement
of our findings compare to the Cramér-Lundberg upper bound has been given.
keywords: Ruin probability; Double Stochastic Compound Poisson Processes;
Bonus–Malus system; Heavy-tailed distributions; Laplace transforms; Cramér-
Lundberg upper bound.

##  1 Introduction

Consider double stochastic compound Poisson process

$\displaystyle U_{t}$ $\displaystyle=$ $\displaystyle
u+\sum_{i=1}^{N_{1}(t)}C_{i}-\sum_{j=1}^{N_{2}(t)}X_{j},$ (1)

where $C_{1},C_{2},\cdots$ and $X_{1},X_{2},\cdots,$ respectively, are two
independent i.i.d. random samples from independent random premium $C$ and
random claim size $X$, two independent Poisson processes $N_{1}(t)$ and
$N_{2}(t)$ (with intensity rates $\lambda_{1}$ and $\lambda_{2}$) are,
respectively, stand for claims and purchase request processes, and $u$
represents initial wealth/reserve $u$ of the process. Moreover, suppose that
non-negative and continuous random premium $C$ and claim $X,$ respectively,
have density functions $f_{C}$ and $f_{X}$ with some additional properties to
be discussed later (see Assumption 1 and Assumption 2).

The ruin probability for such process can be defined by

$\displaystyle\psi(u)$ $\displaystyle=$ $\displaystyle
p(T_{u}<\infty|U_{0}=u),$ (2)

where $T_{u}$ is the hitting time, i.e., $T_{u}:=\inf\\{t:~{}U_{t}\leq 0\\}.$

Several authors studied the ruin probability under a stochastic income
assumption such as the surplus process 1. For instance, Lappo (2004) provided
an equation for the ruin probability of a surplus process which contains two
compound Poisson processes. Rongming et al. (2007) considered a surplus
process with random premium and geometric Lévy investments return process.
They obtained an integro-differential equation for the ruin probability of
such process. Lim & Qi (2009) considered a discrete-time surplus process and
established an equation for the ultimate ruin probability. Moreover, they
obtained an upper bound for the ruin probability and studied its properties
via a simulation study. Under the compound binomial processes, Bao & Wang
(2013) provided a difference equation and a defective renewal equation
satisfied respectively by the expected discounted penalty function. Another
version of the surplus compound binomial process has been studied by Yu
(2013a). Yu (2013b) considered surplus process that both stochastic premium
income and stochastic claims occurrence are driven by the Markovian regime-
switching process. Then, he derived the Laplace transform of the time of ruin.
Temnov (2014) provided a recursive formula similar to the Beekman convolution
formula to evaluate the ruin probability for a stochastic premium surplus
process. Landriault & Shi (2014) studied the finite-time ruin probability for
double stochastic compound poisson processes whenever busy period follows a
fluid flow model. Gatto & Baumgartner (2014a) using the saddlepoint method to
derive an approximation for ruin probability of compound poisson risk process
perturbed by diffusion. Gatto & Baumgartner (2014b) extended Gatto &
Baumgartner (2014)’s findings to compound poisson risk process perturbed by a
Wiener process with infinite time horizon. Cai & Yang (2014) provided an
integro-differential equation for ruin probability of a compound Poisson
surplus process perturbed by diffusion with debit interest.

It is well-known that, in the situation that, random premium $C$ and random
claim size $X$ are two independent exponential random variables survival
probability $\widetilde{\psi}(\cdot)$ can be found explicitly as an
exponential function (Melnikov, 2011, Proposition 8.1). By an induction
argument this fact can be readily generalized to situation that $C$ and $X$
are two independent mixture exponential random variables. In this situation,
survival probability $\widetilde{\psi}(\cdot)$ appears as a mixture of some
exponential functions. This article utilized this fact and approximates
survival probability $\widetilde{\psi}(\cdot)$ appears as a mixture of some
exponential functions. The rest of this article is organized as follows. Some
mathematical background for the problem has been collected in Section 2.
Section 3 provides the main contribution of this article. Applications of the
results along with a and comparison with the Cramér-Lundberg upper bound have
been given in Section 4.

##  2 Preliminaries

The following recalls the exponential type $T$ functions which plays a vital
role in the rest of this article.

######  Definition 1.

An $L_{1}({\mathbb{R}})\cap L_{2}({\mathbb{R}})$ function $f$ is said to be of
exponential type $T$ on ${\mathbb{C}}$ if there are positive constants $M$ and
$T$ such that $|f(\omega)|\leq M\exp\\{T|\omega|\\},$ for
$\omega\in{\mathbb{C}}.$

The Fourier transforms of exponential type functions are continuous functions
which are infinitely differentiable everywhere and are given by a Taylor
series expansion over every compact interval, see Champeney (1987, page 77)
and Walnut (2002, page 81). The Fourier transform to exponential type
functions are also called band-limited functions, see Bracewell (2000, page
119) for more details on band-limited functions. The well-known Paley-Wiener
theorem states that the Fourier transform of an $L_{2}({\mathbb{R}})$ function
vanishes outside of an interval $[-T,T],$ if and only if the function is an
entire function of exponential type $T$, see Dym & McKean (1972, page 158) for
more details. We also need one form of the half-line version of the Paley-
Wiener theorem, namely that if a function in $f\in L_{2}({\mathbb{R}})$
vanishes on the left half-line, then its corresponding Fourier transform
$\hat{f}$ is holomorphic and uniformly bounded in the upper half plane.
Ablowitz & Fokas (1990 §4).

We also have need of the theory of residues for meromorphic functions.
Briefly, $\mbox{Res}(f,z_{i})$ is the coefficient of $\frac{1}{z-z_{i}}$ in
the Laurent series expansion of $f$ around $z_{i}.$ It is possible to define
Laurent series about the point $z_{i}=\infty$, and hence it is possible to
define a residue at infinity, $\mbox{Res}(f,\infty).$ The importance of the
theory of residues probably rests on two facts: many integrals can be
evaluated in terms of a sum of residues, and at poles of finite order residues
can be evaluated efficiently by formulas:

$\mbox{Res}(f,z_{0})=\frac{1}{(m-1)!}\lim_{s\rightarrow
z_{0}}\left(\frac{d}{ds}\right)^{m-1}\\!\\!\\!\\!\\!\\!\\!\\!\\!(s-z_{0})^{m}f(s),$

where $m$ is the order of the pole, see for example Ablowitz & Fokas (1990
§4).

The following explores the Laplace transform of a function that behaves like
an exponential function about zero.

######  Lemma 1.

Suppose $f(\cdot)$ is an exponential type function which behaves like an
exponential function about zero, i.e.,
$f(\epsilon)=\alpha_{0}e^{-\beta_{0}\epsilon},$ where $\alpha_{0}$ and
$\beta_{0}$ are two given positive numbers. Then, the Laplace transform
$f(t+a)$ is

$\displaystyle\mathcal{L}(f(t+a);t;s)$ $\displaystyle=$ $\displaystyle
e^{sa}\mathcal{L}(f(t);t;s)-\frac{\alpha_{0}}{s+\beta_{0}}\left(e^{as}-e^{-a\beta_{0}}\right).$
(3)

_Proof._ Suppose $f$ is of exponential type $T,$ the Taylor series expansion
for $f(t+a)$ is

$f(t+a)=\sum_{0}^{\infty}\frac{a^{k}}{k!}f^{(k)}(t),$

where $f^{(k)}(\cdot)$ stands for the $k^{\hbox{th}}$ derivative of
$f(\cdot).$

Using the Weierstrass M-test along with properties of exponential type
function (i.e., $|f^{(k)}(t)|\leq|T|^{k}\|f\|_{\infty}.$). One may conclude
that the above Taylor series converges uniformly with respect to both
$a\in{\mathbb{C}}$ and $t\in{\mathbb{R}}.$ The uniform convergence justifies
the following taking term-by-term Laplace transform.

$\displaystyle\mathcal{L}(f(t+a);t;s)$ $\displaystyle=$
$\displaystyle\mathcal{L}(f(t);t;s)+\sum_{k=1}^{\infty}\frac{a^{k}}{k!}\mathcal{L}(f^{(k)}(t);t;s)$
$\displaystyle=$
$\displaystyle\mathcal{L}(f(t);t;s)\left[1+\sum_{k=1}^{\infty}\frac{(as)^{k}}{k!}\right]-\sum_{k=1}^{\infty}\frac{a^{k}}{k!}\sum_{j=0}^{k-1}s^{j}f^{(k-j-1)}(0)$
$\displaystyle=$ $\displaystyle
e^{as}\mathcal{L}(f(t);t;s)-\sum_{k=1}^{\infty}\frac{a^{k}}{k!}\sum_{j=0}^{k-1}s^{j}\alpha_{0}(-\beta_{0})^{k-j-1}$
$\displaystyle=$ $\displaystyle
e^{as}\mathcal{L}(f(t);t;s)-\frac{\alpha_{0}}{s+\beta_{0}}\left(e^{as}-e^{-a\beta_{0}}\right),$

where the second equality arrived from an application of the Laplace transform
of a derivative, see Schiff (1999) for more details. $\square$

We consider a class of random variables with their corresponding
density/probability functions and moment generating functions satisfy the
following two assumptions:

Assumption 1.

Their density/probability and survival functions are exponential type
functions, twice (continuously) differentiable, and bounded on
${\mathbb{R}}^{+};$

Assumption 2.

Their moment generating functions are or can be extended to a meromorphic
function on ${\mathbb{C}}.$

Exponential, mixture exponential, Erlang, Pareto (with finite mean), normal,
mixture normal, lognormal, etc are some distributions that satisfy Assumption1
(Cai, 2004) while Assumption 2 holds for most of the common distributions,
such as binomial, negative binomial, poisson, uniform, normal, chi-squared
(with an even degree of freedom), gamma (with integer shape parameter),
laplace, etc. On the other hand, several nonmeromorphic moment generating
functions can be analytically continued to a meromorphic functions. Geometric
and Multinomial distributions are two examples which have nonmeromorphic
moment generating functions but their moment generating functions can be
analytically continued to a meromorphic function, see Kucerovsky & Payandeh
(2014) and Sasv’ri (2013, §3.3) for more details.

By conditioning of the survival probability $\widetilde{\psi}(\cdot)$ of
surplus process 1 on the first arriving premium and claim along with
properties of compound Poisson process. One may derive the following integro-
differential equation for survival probability $\widetilde{\psi}(\cdot)$ of
surplus process 1, proof may be found in Rongming et al. (2007), Melnikov
(2011, Theorem 8.1), among others.

$-(\lambda_{1}+\lambda_{2})\widetilde{\psi}(u)+\lambda_{1}E(\widetilde{\psi}(u+C))+\lambda_{2}\int_{0}^{u}\widetilde{\psi}(u-x)f_{X}(x)dx=0,$
(4)

where $\lim_{u\rightarrow\infty}\widetilde{\psi}(u)=1$ and two random premium
$C$ and claim size $X$ satisfy assumptions Assumption 1 and Assumption 2.

In the situation that, $C$ and $X$ are two independent exponential random
variables with rates $\mu_{1}$ and $\mu_{2},$ respectively. Then, survival
probability can be found explicitly as $\widetilde{\psi}(u)=Ae^{-Bu},$ where
$A:=(\mu_{1}+\mu_{2})\lambda_{2}/(\lambda_{1}+\lambda_{2})$ and
$B:=(\mu_{1}\lambda_{1}-\mu_{2}\lambda_{2})/(\lambda_{1}+\lambda_{2}),$ see
Melnikov (2011, Proposition 8.1) for more details. This fact can be
generalized for situation that $C$ and $X$ are two independent mixture
exponential random variables.

Hereafter now, we assume for very small surplus/reserve, say $u=\epsilon,$
survival probability behaves like an exponential function, i.e.,

Assumption 3.

For very small surplus/reserve $u=\epsilon$ survival probability behaves like
$\widetilde{\psi}(\epsilon)=\alpha_{0}e^{-\beta_{0}u},$

where two positive numbers $\alpha_{0}$ and $\beta_{0}$ have to be determined.
The above assumption can be explained by the above observation along with the
fact that a given density/probability function can be approximated, with some
degree of accuracy, by a mixture exponential distribution.

Coefficient $\beta_{0}$ in Assumption 3 can be found _either_ by integro-
differential Equation 4 and letting $u\rightarrow\infty$ _or_ the fact that
$\lim_{u\rightarrow\infty}\psi(u)e^{Ru}=c\in(0,1),$ where adjustment
coefficient $R$ is positive solution of
$\lambda_{1}M_{C}(R)+\lambda_{2}M_{X}(-R)=\lambda_{1}+\lambda_{2}$ (Kaas, et
al., Page 112), see Proposition 2 for more details. While coefficient
$\alpha_{0}$ cannot be found using integro-differential Equation 4 or other
theoretical assumptions on ruin probability. Therefore, analogue to situation
that $C$ and $X$ are two independent exponential random variables, we set
$\alpha_{0}:=(1/E(C)+1/E(X))\lambda_{2}/(\lambda_{1}+\lambda_{2}).$

##  3 Main results

This section utilizes integro-differential Equation 4 to derive two
approximate formulas for the ruin probability of a double stochastic compound
Poisson process. 4. We seek an analytical solution $\widetilde{\psi}(\cdot)$
which is an exponential type function. In the other word, we assume:

Assumption 4.

$|\widetilde{\psi}(\omega)|\leq Me^{T|\omega|},~{}\omega\in{\mathbb{C}},$ for
some real numbers $M$ and $T$ in ${\mathbb{R}}.$

If this assumption is not met, as might be the case if, for example, there are
point masses in $\psi(\cdot),$ our methods may still result in a formal
solution that can be verified by substitution into integro-differential
Equation 4.

The following from provides an explicit expression for the ruin probability of
a double stochastic compound Poisson process in an analytical model.

######  Theorem 1.

Suppose $U_{t}$ represents double stochastic compound Poisson process 1 which
their random claim size $X$ and random premium $C$ satisfy Assumption 1 and
Assumption 2. Moreover, suppose that survival probability
$\widetilde{\psi}(\cdot)$ satisfies Assumption 3 and Assumption 4. Then,

i)

the Laplace transform of the ruin probability $\psi(\cdot)$ satisfies

$\displaystyle\mathcal{L}(\psi(u);u,s)$ $\displaystyle=$
$\displaystyle\frac{N(s)}{s(s+\beta_{0})D(s)};$ (5)

ii)

poles of $\frac{N(s)}{s(s+\beta_{0})D(s)},$ say $z_{j},$ can be represented as
$z_{j}:=-a_{j}\pm ib_{j},$ where $a_{j}>0;$

iii)

There is at most one simple pole on the negative real axis,

where
$D(s):=-\lambda_{1}-\lambda_{2}+\lambda_{1}M_{C}(s)+\lambda_{2}M_{X}(-s),$
$N(s):=(s+\beta_{0})D(s)-s\lambda_{1}\alpha_{0}(M_{C}(s)-M_{C}(-\beta_{0})),$
and $M_{C}(\cdot)$ and $M_{X}(\cdot)$ represent moment generating functions of
$C$ and $X,$ respectively.

_Proof._ Taking a Laplace transform from integro-differential Equation 4 with
an application of Lemma 1 complete part (i). For part (ii) observe that: by
Assumption 3 and the Paley-Wiener theorem, we may conclude that an expression
$g(s):=\frac{N(s)}{s(s+\beta_{0})D(s)}$ is holomorphic and uniformly bounded
in the right-half plane. Therefore, $g(s)$ has no poles at any
$z\in{\mathbb{C}}$ with non-negative real part. On the other hand, since
$\psi(\cdot)$ is real-valued function, its corresponding Laplace transform is
real everywhere on the real line (except possibly at poles) and hence by the
Schwartz reflection principle,

$\overline{\left(\frac{N(s)}{s(s+\beta_{0})D(s)}\right)}=\frac{N(\overline{s})}{\overline{s}(\overline{s}+\beta_{0})D(\overline{s})}.$

Therefore, poles of $g(s)$ are located at $z_{j}:=-a_{j}\pm ib_{j}.$ To
establish $a_{j}>0,$ observe that if there is a pole at some
$z_{i}\in{\mathbb{C}},$ then there is also a pole at the complex conjugate,
$\overline{z_{i}}.$ On the other hand, since $D(s)$ is a convex (concave up)
function on the real line, it can have at most two zeros on the real line.
Moreover, $D(s)$ has a zero at the origin, so there can be at most one other
zero of $D(s)$ on the real line, necessarily located on the negative real
line. For part (iii) observe that a simple pole for $g(s)$ on the negative
real axis can be $-\beta_{0}.$ $\square$

Using the above theorem, one need only consider the roots of $D(s)=0,$ which
are completely determined by the moment generating functions, random claim
size $X,$ and random random premium $C,$ say respectively $M_{X}(\cdot)$ and
$M_{C}(\cdot).$ Residues at simple poles are particularly easy to evaluate,
and the case of simple poles is the generic case. If there are finitely many
poles and all are of finite order, they may be made to be simple poles by an
infinitesimal perturbation of the problem. Thus the following result is of
practical importance. The condition that appears below, of the derivatives
being non-zero, is more or less equivalent to all the zeros of $D(s)$ being
simple zeros:

######  Corollary 1.

Under the same conditions as in Theorem 1, we then have that the ruin
probability $\psi(\cdot)$ is

$\displaystyle\psi(u)$ $\displaystyle\approx$
$\displaystyle\sum_{D(z_{i})=0,\,\Re(z_{i})<0}\frac{N(z_{i})}{z_{i}(z_{i}+\beta_{0})D^{\prime}(z_{i})}e^{z_{i}u}.$

whenever the derivatives that appear are all non-zero.

_Proof._ The residue of $N(s)/(s(s+\beta_{0})D(s))$ at a simple zero $z_{i}$
of $D(s)$ is

$\lim_{s\rightarrow
z_{i}}\frac{(s-z_{i})N(s)}{s(s+\beta_{0})D(s)}=\frac{N(z_{i})}{z_{i}(z_{i}+\beta_{0})D^{\prime}(z_{i})}.$

Conversely, if $D^{\prime}(s)$ is nonzero at a zero of $D(s)$ then this zero
is a simple zero, and $N(s)/(s(s+\beta_{0})D(s))$ has at most a simple pole
(as the zero is assumed to have negative real part). $\square$

Since the ruin probability is real and non-negative, and since the poles come
in pairs, we can draw some further conclusions about the structure of the
expression for $\psi(\cdot).$

######  Proposition 1.

Under the same conditions as in Theorem 1. The ruin probability $\psi(u)$ can
be simplified by a finite sum of the form

$\displaystyle\psi(u)$ $\displaystyle\approx$
$\displaystyle\alpha_{0}e^{-c_{0}u}+\sum_{i}A_{i}e^{-c_{i}u}\cos(d_{i}u)+\sum_{i}B_{i}e^{c_{i}u}\sin(d_{i}u),$
(6)

where $c_{i}$ and $d_{i}$ are, respectively, real and imaginary parts of the
roots of $D(s)=0$ and $A_{i},$ $B_{i}$ are real numbers evaluated by
substituting 6 in Equation 4.

_Proof._ We have the above general form of a sum involving trigonometric
functions and exponentials comes directly from Theorem 1, plus the fact that
the poles come in pairs (Theorem 1.) The $b_{i}$ are the imaginary parts of
the roots, and the $a_{i}$ are the real parts. Theorem 1 moreover tells us
that there is at most one root on the negative real axis. Since the ruin
probability must remain non-negative for large $u,$ but must vanish at
infinity, it follows that the exponential corresponding to this purely real
root must be the dominant term for large $u,$ and thus the purely real root is
negative, nonzero, and is to the right of all the complex roots. $\square$

######  Remark 1.

It is a corollary of the proof above that under our conditions there does
always exist a unique negative and nonzero real root of $D(s).$

######  Corollary 2.

Under the same conditions as in Theorem 1, and also assuming $D(s)/s=0$ has
only a simple non-zero root at $z_{1},$ then, the ruin probability is given by

$\displaystyle\psi(u)$ $\displaystyle\approx$
$\displaystyle\frac{\lambda_{2}}{\lambda_{1}+\lambda_{2}-\lambda_{1}M_{C}(z_{1})}e^{z_{1}u},$

where $M_{C}(z_{1})$ stands for the moment generating function of the random
variable $C$ at $z_{1}.$

_Proof._ Substituting $u=0$ in Equation 4, one may conclude that

$\displaystyle\alpha_{0}$ $\displaystyle:=$ $\displaystyle\psi(0)$
$\displaystyle=$
$\displaystyle\frac{\lambda_{1}}{\lambda_{1}+\lambda_{2}}E(\psi(C))+\frac{\lambda_{2}}{\lambda_{1}+\lambda_{2}}$
$\displaystyle=$
$\displaystyle\frac{\lambda_{1}}{\lambda_{1}+\lambda_{2}}\frac{N(z_{1})}{-z_{1}D^{\prime}(z_{1})}M_{C}(z_{1})+\frac{\lambda_{2}}{\lambda_{1}+\lambda_{2}}$
$\displaystyle=$
$\displaystyle\frac{\lambda_{1}}{\lambda_{1}+\lambda_{2}}\alpha_{0}M_{C}(z_{1})+\frac{\lambda_{2}}{\lambda_{1}+\lambda_{2}},$

where two last equalities arrive from double applications of Corollary 1, with
simple non-zero root $z_{1}$ at $u=0.$ $\square$

In general, the moment generating function of a distribution, when it exists,
will grow rapidly as we move to the right on the real axis, but decreases
along the negative real axis, and decreases along the imaginary axis (because
there it is the characteristic function). Thus, it is reasonable to expect
that $M_{C}(z_{k})$ will be small for the off-axis roots of $D(s)$, since
these roots are located in the second and third quadrants of the complex
plane. Furthermore, since the moment generating function oscillates in the
off-axis direction, sums involving the off-axis roots are likely to be quite
small due to cancellations. If this is the case, then we can hope to
approximate by neglecting the sum involving off-axis roots roots in the above
calculation. We thus obtain the following approximation result.

######  Proposition 2.

. Under the same conditions as in Theorem 1, then the ruin probability
$\psi(\cdot)$ about the origin can be approximated by one of the following
exponential functions.

i)

If equation $D(s)=0$ has unique strictly negative real root, say $z_{1},$

$\displaystyle\psi(\epsilon)$ $\displaystyle\approx$
$\displaystyle\frac{\lambda_{2}}{\lambda_{1}+\lambda_{2}-\lambda_{1}M_{C}(z_{1})}e^{z_{1}\epsilon};$

ii)

If $\lambda_{1}E(C)\geq\lambda_{2}E(X)$

$\displaystyle\psi(\epsilon)$ $\displaystyle\approx$
$\displaystyle\frac{E(C)+E(X)}{E(C)E(X)(\lambda_{1}+\lambda_{2})}e^{-\beta_{0}\epsilon};$

iii)

$\displaystyle\psi(\epsilon)$ $\displaystyle\approx$
$\displaystyle\frac{E(C)+E(X)}{E(C)E(X)(\lambda_{1}+\lambda_{2})}e^{-R\epsilon},$

where $\beta_{0}$ is positive solution of
$\lambda_{1}M_{C}(-\beta_{0})+\lambda_{2}M_{X}(\beta_{0})=\lambda_{1}+\lambda_{2}$
and adjustment coefficient $R$ is positive solution of
$\lambda_{1}M_{C}(R)+\lambda_{2}M_{X}(-R)=\lambda_{1}+\lambda_{2}.$

_Proof._ Part (i) arrives by an application of Corollary 2. For part (ii)
observe that substituting an exponential function $\alpha_{0}e^{-\beta_{0}u}$
into integro-differential Equation 4 and letting $u\rightarrow\infty.$ One may
obtain equation
$\lambda_{1}M_{C}(-\beta_{0})+\lambda_{2}M_{X}(\beta_{0})=\lambda_{1}+\lambda_{2}.$
Now, set
$h(\beta_{0}):=\lambda_{1}M_{C}(-\beta_{0})+\lambda_{2}M_{X}(\beta_{0}).$
Using assumption on random premium $C$ and random claim size $X,$ one may
conclude that $h(0)=\lambda_{1}+\lambda_{2},$
$\lim_{\beta_{0}\rightarrow\infty}h(\beta_{0})=\infty,$
$\lim_{\beta_{0}\rightarrow\infty}h^{\prime}(\beta_{0})=-\lambda_{1}E(C)+\lambda_{2}E(X)\leq
0,$ and $h^{\prime\prime}(\beta_{0})\geq 0.$ Therefore, the above equation has
a unique positive solution in $\beta_{0}.$ Part (iii) arrives by an
application of the fact that
$\lim_{u\rightarrow\infty}\psi(u)e^{Ru}=c\in(0,1),$ (Kaas, et al., Page 112).
Coefficient $\alpha_{0}$ for both parts (ii & iii) cannot be found using
integro-differential Equation 4 or other theoretical assumptions on ruin
probability. Therefore, analogue to situation that $C$ and $X$ are two
independent exponential random variables, we set
$\alpha_{0}:=(1/E(C)+1/E(X))\lambda_{2}/(\lambda_{1}+\lambda_{2}).$ $\square$

Residues at simple poles are particularly easy to evaluate, and the case of
simple poles is the generic case. If there are finitely many (finite order)
poles. They may be made to be simple poles by an infinitesimal perturbation of
the problem. Therefore, an approximation for the ruin probability can be
obtained by assuming all roots of $D(s)=0$ are simple. The following theorem
improves accuracy of approximation result given by Theorem 1 whenever $D(s)=0$
has finitely many simple zeros.

######  Theorem 2.

Supposing that under the conditions of the above result, $D(s)$ has finitely
many roots in the left half-plane. Then, $\psi(u)\approx\sum A_{i}e^{z_{i}u},$
where $A_{i}$ can be found out by the following linear system of equations

$A_{j}=\frac{1}{\lambda_{2}z_{j}M^{\prime}_{X}(-z_{j})}\left(\lambda_{2}M_{X}(-z_{j})-\lambda_{2}+z_{j}\sum_{i\neq
j}\frac{A_{i}}{z_{j}-z_{i}}\left(\lambda_{1}+\lambda_{2}-\lambda_{1}M_{C}(z_{i})-\lambda_{2}M_{X}(-z_{j})\right)\right).$
(7)

_Proof._ Roots of $D(s)=0$ may be made to be simple poles by an infinitesimal
perturbation of the problem. Therefore, $D(\cdot)$ has finitely many simple
roots, $z_{i}$, we have by Lemma 1 and Proposition 1 that $\psi(u)\approx\sum
A_{i}e^{z_{i}u}.$ Substituting $\psi(u)\approx\sum A_{i}e^{z_{i}u}$ into
Equation 4 and taking the Laplace transform, we obtain

$\displaystyle-(\lambda_{1}+\lambda_{2})\left(\frac{1}{s}-\sum_{i}\frac{A_{i}}{z_{j}-z_{i}}\right)+\lambda_{1}\left(\frac{1}{s}-\sum_{i}\frac{A_{i}}{z_{j}-z_{i}}M_{C}(z_{i})\right)+\lambda_{2}\left(\frac{1}{s}-\sum_{i}\frac{A_{i}}{z_{j}-z_{i}}\right)M_{X}(-s)$
$\displaystyle=$ $\displaystyle 0.$

The desired proof arrives by setting $s=z_{j}$ along with an application of
the $l^{\prime}H\hat{o}pital^{\prime}$s rule. $\square$

It is clear from the above that the solution obtained is in fact unique. The
only further point to mention with respect to the above theorem that since we
now are dealing with approximations, it could happen that the sum $\sum
A_{i}e^{z_{i}u}$ obtained as in Theorem 2 is not quite real or not quite non-
negative, so we adjust it in the above by taking the real part and then the
floor with zero. The amount by which the sum $\sum A_{i}e^{z_{i}u}$ fails to
be real and non-negative can be taken as a practical indication of the amount
of error in the approximation. Normally, one would pick the complex zeros that
are closest to the real root of $D(s),$ aiming to ensure that the complex
moment function $M_{C}(z)$ has quite small absolute value at the zeros that
have been neglected. Moreover, the set of roots chosen should be closed under
complex conjugation, meaning that if $a+ib$ is in the set of roots, then
$a-ib$ should also be in the set of roots.

It would be worthwhile mentioning that, the ruin probability may be impacted
by different ways of adjustment (i.e., adjusting either $\psi(u)$ or complex
roots). Sensitivity of such adjustment should be investigated and managed by
the amount of error which obtained by substituting adjusted $\psi(\cdot)$ into
Equation 4.

The following section is devoted to applications of Corollary 2 and Theorem 2.

##  4 Applications

In the most of practical application the ruin probability of surplus process 1
has been approximated by the Cramér-Lundberg upper bound $e^{-Ru},$ where
adjustment coefficient $R$ is positive solution of
$\lambda_{1}M_{C}(R)+\lambda_{2}M_{X}(-R)=\lambda_{1}+\lambda_{2}.$ Certainly,
$-R$ is one of solution of our equation $D(s)=0.$ Practical implementation of
the above findings along with a comparison with the well-known Cramér-Lundberg
upper bound have been given in the following examples. Example 1 considers a
situation that equation $D(s)=0$ has just one simple non-positive solution. In
this situation our method improves the Cramér-Lundberg upper bound by a
constant coefficient. While in the Example 2 that equation $D(s)=0$ has more
than one simple non-positive root, our method provides a significant
improvement on the Cramér-Lundberg upper bound.

######  Example 1.

Suppose an insurance company plans to charge its policyholders based upon
their risks (such insurance systems well known as Bonus–Malus systems).
Moreover suppose that the insurance company considers three different
scenarios with 5, 10, and 20 premiums’ values such that under all of these
scenarios expected revenue of the company stay the same by selling an
insurance contract. Based upon a statistical investigation an actuary
suggested the following premium values and probability that a given
policyholder falls in a given level, say ${\bf\pi},$ of each scenario333These
three different scenarios also can be viewed as three different Bonus–Malus
systems which stabilized, in the long run, around their corresponding
equilibrium distributions ${\bf\pi},$ well known as long–term Bonus–Malus
systems..

Table 1: Premiums values ${\bf C}$ and probability that a given policyholder
falls in a given level for such three scenarios.
Scenario (number of its levels) ${\bf\pi}$ Premium value ${\bf C}$ S1 (with 5
levels) $\pi_{i}\equiv 0.2$ (0.6, 1, 1.4, 1.8, 2.2) S2 (with 10 levels)
$\pi_{i}\equiv 0.1$ (0.4, 0.8, 1, 1.2, 1.4, 1.5, 1.7, 1.8, 2, 2.2) S3 (with 20
levels) $\pi_{i}\equiv 0.05$ (0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 1.1, 1.2, 1.3,
– – 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.1, 2.2, 2.6, 2.7)

Moreover, suppose that the number of sold contracts, $N_{1}(t),$ and number of
arrived claims, $N_{2}(t),$ are two independent Poisson processes with
intensity $\lambda_{1}=18$ and $\lambda_{2}=11$ and claim size distributions
given by Table 2.

It is easy to show that, for all cases, equation $D(s)=0$ has just one simple
non-positive root. Using results of Corollary 2, one may evaluate the ruin
probability in an exact mode. The ruin probability (as well as the Cramér-
Lundberg upper bound) for such three different scenarios have been given by
Table 2.

Table 2: The ruin probability and Cramér-Lundberg upper bound for the above
three scenarios under different claim size distributions.
Scenario (number of its levels) Claim size distribution Scenario 1 (5)
Scenario 2 (10) Scenario 3 (20) Gamma(1,3) $\psi(u)=0.313e^{-1.760u}$
$\psi(u)=0.312e^{-1.759u}$ $\psi(u)=0.280e^{-1.746u}$ $(\psi(u)\leq
e^{-1.760u})$ $(\psi(u)\leq e^{-1.759u})$ $(\psi(u)\leq e^{-1.746u})$
Gamma(1,1) $\psi(u)=0.072e^{-0.408u}$ $\psi(u)=0.072e^{-0.407u}$
$\psi(u)=0.071e^{-0.402u}$ $(\psi(u)\leq e^{-0.408u})$ $(\psi(u)\leq
e^{-0.407u})$ $(\psi(u)\leq e^{-0.402u})$ Gamma(3,2)
$\psi(u)=0.057e^{-0.229u}$ $\psi(u)=0.057e^{-0.227u}$
$\psi(u)=0.056e^{-0.223u}$ $(\psi(u)\leq e^{-0.229u})$ $(\psi(u)\leq
e^{-0.227u})$ $(\psi(u)\leq e^{-0.223u})$ Gamma(5,3)
$\psi(u)=0.053e^{-0.174u}$ $\psi(u)=0.053e^{-0.173u}$
$\psi(u)=0.053e^{-0.170u}$ $(\psi(u)\leq e^{-0.174u})$ $(\psi(u)\leq
e^{-0.173u})$ $(\psi(u)\leq e^{-0.170u})$

As one may observe, our method just improve the Cramér-Lundberg upper bound by
a constant coefficient.

Figure 1 illustrates several comparison regarding to the ruin probability of
such three scenarios under different Claim size distributions.

Figure 1: Figures (a) Density function of claim size distributions; (b) to (d)
comparison between the ruin probability of such claim size distributions with
respect to 5, 10, and 20 premium values scenarios, respectively; (e) to (h)
comparison between the ruin probability of 5, 10, and 20 premium values
scenarios with respect to different claim size distributions (i.e, Gamma(1,1),
Gamma(1,3), Gamma(3,2), Gamma(5,3)).

From the above figures, one may conclude that the ruin probability $\psi(u)$
as function of initial wealth, $u,$ goes to zero and its decay rate depends on
tail of its corresponding claim size distribution. Moreover, the ruin
probability decreases as the number of premium levels increases. The second
observation in different context has been pointed by several authors, see
Denuit, et al. (2007) for more detail.

######  Example 2.

Suppose the claim size is a considerable heavy-tailed distribution with
density function

$\displaystyle f_{X}(x)$ $\displaystyle=$ $\displaystyle
0.150e^{-(x-4)^{2}/2}+\frac{0.075x}{(1+x^{2})^{2}}+0.074e^{-3x}+\frac{0.224e^{-1/(2x)}}{x^{1.5}},x\geq
0$ (8)

and random premium $C$ has the Gamma distribution with parameters
$\kappa=\theta=3.$ Moreover, suppose that the number of sold contracts,
$N_{1}(t),$ and number of arrived claims, $N_{2}(t),$ are two independent
Poisson processes with intensity $\lambda_{1}=18$ and $\lambda_{2}=11.$ Using
a mathematical software such as Maple, one may show that an equation $D(s)=0$
has infinity many simple non-positive roots, see Figure 2(b) for an
illustration.

For simplicity, we just considered roots which their both real and imaginary
parts are fall within the interval $[-1,2].$ In this subset of complex plane
${\mathbb{C}},$three roots $z_{1}:=-0.5639909305-1.521665917I,$
$z_{2}:=-0.4233614172,$ and $z_{3}:=-0.5639909305+1.521665918I$ have been
found. Now using Theorem 2, one may verify that the ruin probability can be
evaluated using $\psi(u)=A_{1}e^{z_{1}u}+A_{2}e^{z_{2}u}+A_{3}e^{z_{3}u},$
where

$\displaystyle A_{1}$ $\displaystyle=$
$\displaystyle-\frac{-M_{1}z_{3}\beta_{23}z_{2}\beta_{32}+M_{1}T_{2}T_{3}+z_{1}\beta_{21}z_{2}\beta_{32}M_{3}+z_{1}\beta_{21}M_{2}T_{3}+z_{1}\beta_{31}T_{2}M_{3}+z_{1}\beta_{31}z_{3}\beta_{23}M_{2}}{T_{1}z_{3}\beta_{23}z_{2}\beta_{32}-T_{1}T_{2}T_{3}+T_{2}z_{3}\beta_{13}z_{1}\beta_{31}+z_{2}\beta_{21}^{2}z_{1}T_{3}+z_{3}\beta_{23}z_{2}\beta_{21}z_{1}\beta_{31}+z_{3}\beta_{13}z_{1}\beta_{21}z_{2}\beta_{32}}$
$\displaystyle A_{2}$ $\displaystyle=$
$\displaystyle\frac{-\beta_{21}\beta_{31}M_{3}z_{2}z_{1}-\beta_{21}z_{2}M_{1}T_{3}+\beta_{31}z_{3}\beta_{13}M_{2}z_{1}-z_{2}T_{1}\beta_{32}M_{3}-z_{2}M_{1}z_{3}\beta_{13}\beta_{32}-T_{1}M_{2}T_{3}}{T_{1}z_{3}\beta_{23}z_{2}\beta_{32}-T_{1}T_{2}T_{3}+T_{2}z_{3}\beta_{13}z_{1}\beta_{31}+z_{2}\beta_{21}^{2}z_{1}T_{3}+z_{3}\beta_{23}z_{2}\beta_{21}z_{1}\beta_{31}+z_{3}\beta_{13}z_{1}\beta_{21}z_{2}\beta_{32}}$
$\displaystyle A_{3}$ $\displaystyle=$
$\displaystyle\frac{-z_{3}\beta_{23}z_{2}\beta_{21}M_{1}+T_{1}T_{2}M_{3}+T_{1}z_{3}\beta_{23}M_{2}-z_{2}\beta_{21}^{2}z_{1}M_{3}+T_{2}z_{3}\beta_{13}M_{1}+z_{3}\beta_{13}z_{1}\beta_{21}M_{2})}{(T_{1}z_{3}\beta_{23}z_{2}\beta_{32}-T_{1}T_{2}T_{3}+T_{2}z_{3}\beta_{13}z_{1}\beta_{31}+z_{2}\beta_{21}^{2}z_{1}T_{3}+z_{3}\beta_{23}z_{2}\beta_{21}z_{1}\beta_{31}+z_{3}\beta_{13}z_{1}\beta_{21}z_{2}\beta_{32}},$

where $T_{j}:=\lambda_{2}z_{j}M^{\prime}_{X}(-z_{j})),$
$M_{j}:=\lambda_{2}M_{X}(-z_{j})-\lambda_{2},$ and
$\beta_{ij}=(\lambda_{1}+\lambda_{2}-\lambda_{1}M_{C}(z_{i})-\lambda_{2}M_{X}(-z_{j}))/(z_{j}-z_{i})$
for $i,j=1,2,3.$ Therefore, the ruin probability can be approximated by

$\displaystyle\psi(u_{0})$ $\displaystyle=$
$\displaystyle\frac{1}{100}||(-0.26+0.52I)e^{(-0.5640-1.5217I)u_{0}}-2.87e^{-0.4234u_{0}}+(-0.26-0.52I)e^{(-0.5640+1.5217I)u_{0}}||.$

The Cramér-Lundberg upper bound for this situation is $e^{-0.4234u_{0}}$ which
is significantly improved by our method, see Figure 2(c).

Figure 2: Figures (a) Density function of mixture claim size density 8
accompanied with Gamma (1,1) density function; (b) Roots of equation $D(z)=0$
where $|z|\leq 20\sqrt{2};$ and (c) $\frac{1}{20}\times$ the ruin probability
and the Cramér-Lundberg upper bound of the surplus process 1 under mixture
claim size density 8.

##  5 Conclusion and suggestions

The Ruin theory provides some potential tools to study today’s solvency
assessments for non-life insurance companies, see Wüthrich (2014) for a
practical application. In most practical applications the ruin probability of
surplus process 1 has been approximated by the Cramér-Lundberg upper bound
$e^{-Ru_{0}},$ where adjustment coefficient $R.$

This article assumed the ruin probability $\psi(\cdot)$ is an exponential type
$T$ and $L_{2}({\mathbb{R}})$ function. Then, it derived several approximate
formulas for the ruin probability of a double stochastic compound Poisson
process. The approximated ruin probability constructed based upon roots of
equation $D(s)=0.$ Certainly, $-R$ is one of such roots. Therefore, the
Cramér-Lundberg upper bound, multiplied by a constant, will be a part of our
approximated ruin probability. Since $-R$ is the largest roots of the equation
$D(s)=0,$ our approximation method improved the Cramér-Lundberg upper bound.
Such improvement is significant whenever equation $D(s)=0$ has more than one
simple non-positive root, see Example 2.

The exponential type assumption can be dropped and a formal solution can be
verified by substitution into Equation 4. But, the generally realistic
$L_{2}({\mathbb{R}})$ assumption cannot be dropped, because our methods always
produce solutions with finite $L_{2}({\mathbb{R}})$ norm. Furthermore, it is
essential to our methods that the Laplace transform of $f_{X}$ does not have
any branch points in the complex plane.

## Acknowledgements

The support of Shahid Beheshti University and Natural Sciences and Engineering
Research Council (NSERC) of Canada are gratefully acknowledged by authors.

## References

  * [1] Ablowitz, M. J. & Fokas, A. S. (1990). Complex variable, introduction and application. Springer-Verlag, Berlin.
  * [2] Bao, Z., & Wang, J. (2013). On The Compound Binomial Risk Mondel With Stochastic Income. International Journal of Pure and Applied Mathematics, 82(3), 377-389.
  * [3] Bracewell, R. N. (2000). The Fourier transform and its applications. 3rd. ed. McGraw-Hill publisher, New York.
  * [4] Cai, J. (2004). Ruin probabilities and penalty functions with stochastic rates of interest. Stochastic process Appl., 112, 53–78.
  * [5] Cai, J., & Yang, H. (2014). On the decomposition of the absolute ruin probability in a perturbed compound Poisson surplus process with debit interest. Annals of Operations Research, 212(1), 61–77.
  * [6] Champeney, D. C. (1987). A handbook of Fourier theorems. Cambridge University press, New York.
  * [7] Denuit, M., Marechal, X., Pitrebois, S., & Walhim, J. F. (2007). Actuarial Modelling of Claim Counts: Risk Classification, Credibility and Bonus–Malus systems. John Wiley & Sons, New York.
  * [8] Dym, H. & Mckean, H. P. (1972). Fourier series and integrals. Probability and Mathematical Statistics. Academic Press, New York-London.
  * [9] Duffin, R. J. & Schaeffer, A. C. (1937). Some Inequalities concerning functions of Exponential Type. Bull. AMS, 43, 554–556.
  * [10] Kaas, R., Goovaerts, M., Dhaene, J., & Denuit, M. (2008). Modern actuarial risk theory: using R.
  * [11] Kucerovsky, D. Z. & Payandeh, A. T. (2014). On Ruin Probability Of Long–Term Bonus–Malus Systems. Submitted.
  * [12] Gatto, R., & Baumgartner, B. (2014a). Saddlepoint Approximations to the Probability of Ruin in Finite Time for the Compound Poisson Risk Process Perturbed by Diffusion. Methodology and Computing in Applied Probability, (ahead-of-print), 1–19.
  * [13] Gatto, R., & Baumgartner, B. (2014b). Value at ruin and tail value at ruin of the compound Poisson process with diffusion and efficient computational methods. Methodology and Computing in Applied Probability, (ahead-of-print), 1–22.
  * [14] Landriault, D., & Shi, T. (2014). First passage time for compound Poisson processes with diffusion: ruin theoretical and financial applications. Scandinavian Actuarial Journal, 4, 368–382.
  * [15] Lappo, P. M. (2004). Probability of ruin in collective risk model with random premiums. Vestn, Beloruss. Gos. Univ. Ser. 1 Fiz. Mat. Inform, 121, 103–105.
  * [16] Lim, S. P. & Qi, L. (2009). A discrete-time risk model with investment and interference under random premium. Applied mathematic: Journal of Chinese university series A., 24, 9–14.
  * [17] Melnikov, A. (2011). Risk analysis in finance and insurance. CRC Press.
  * [18] Needham, T. (1999). Visual complex analysis. Clarendon press. Oxford.
  * [19] Rongming, W., Lin, X., & Dingjun, Y. (2007). Ruin problems with stochastic premium stochastic return on investments. Front. Math. China, 2, 467–490.
  * [20] Sasv’ri, Z. (2013). Multivariate Characteristic and Correlation Functions.Walter de Gruyter GmbH, Berlin
  * [21] Schiff, J. L. (1999). _The Laplace transform: Theory and applications._ Springer, New York.
  * [22] Stein, E. M. (1957). Functions of Exponential Type. Annals of Math, 65, 582–592.
  * [23] Temnov, G. (2014). Risk Models with Stochastic Premium and Ruin Probability Estimation. Journal of Mathematical Sciences, 196(1), 84–96.
  * [24] Walnut, D. F. (2002). An introduction to wavelet analysis. 2nd ed. Birkhäuser publisher, New York.
  * [25] Wüthrich, M. V. (2014). From Ruin Theory To Solvency In Non-life Insurance. Scandinavian Actuarial Journal, (ahead-of-print), 1–11.
  * [26] Yu, W. (2013a). Randomized Dividends in a Discrete Insurance Risk Model with Stochastic Premium Income. Mathematical Problems in Engineering, 2013.
  * [27] Yu, W. (2013b). On the Expected Discounted Penalty Function for a Markov Regime-Switching Insurance Risk Model with Stochastic Premium Income. Discrete Dynamics in Nature and Society, 2013, 1–9.