(Version: 2023-02-24)

Uniqueness Problem for the Backward Differential

Equation of a Continuous-State Branching Process111This work was supported by
the National Key R&D Program of China (No. 2020YFA0712900) and the National
Natural Science Foundation of China (No.12271029).

Pei-Sen Li 222School of Mathematics and Statistics, Beijing Institute of
Technology, Beijing 100872, People’s Republic of China<EMAIL_ADDRESS>and Zenghu Li 333School of Mathematical Sciences, Beijing Normal University,
Beijing 100875, People’s Republic of China<EMAIL_ADDRESS>

Abstract: The distributional properties of a multi-dimensional continuous-
state branching process are determined by its cumulant semigroup, which is
defined by the backward differential equation. We provide a proof of the
assertion of Rhyzhov and Skorokhod (Theory Probab. Appl., 1970) on the
uniqueness of the solutions to the equation, which is based on a
characterization of the process as the pathwise unique solution to a system of
stochastic equations.

Key words and phrases. Continuous-state branching process, multi-dimensional,
backward differential equation, stochastic equation, generator.

MSC 2010. Primary 60J80; secondary 60H20

##  1 Introduction

Let $\mathbb{R}_{+}=[0,\infty)$ and let $D=\mathbb{R}_{+}^{m}$ be its $m$th
power for an integer $m\geq 1$. Suppose that $\xi=\\{\xi(t):t\geq 0\\}$ is a
stochastically continuous Markov process in $D$ with homogeneous transition
semigroup $\\{P(t):t\geq 0\\}$. We call $\xi$ a continuous-state branching
process (CB-process) if its transition semigroup satisfies the branching
property:

(1.1) $\displaystyle P(t,x+y,\cdot)=P(t,x,\cdot)*P(t,y,\cdot),\quad x,y\in D,$

where “$*$” denotes the convolution operation. Here $\xi$ may not be
conservative, so it is in general a Markov process in the one-point
compactification $\bar{D}:=D\cup\\{\infty\\}$ with $\infty$ being the
cemetery.

The CB-process is a special form of the measure-valued branching processes
studied in Dawson [7], Dynkin [12], Le Gall [20] and Li [22]. A more general
class of finite-dimensional Markov processes, the affine processes, were
investigated in Duffie et al. [10]. The process has played important roles in
the study of different areas including biological populations, statistical
physics, mathematical finance and so on. The problem of characterizing the
most general branching Markov processes in finite- or infinite-dimensional
state spaces has been studied by Duffie et al. [10], Dynkin et al. [13],
Rhyzhov and Skorokhod [25], Silverstein [28] and Watanabe [29].

The basic structures of the multi-dimensional CB-process were investigated by
Rhyzhov and Skorokhod [25] and Watanabe [29]. Let
$\mathbb{C}_{-}=\\{x+\mathrm{i}y:x\leq 0,y\in\mathbb{R}\\}$, where
$\mathrm{i}=\sqrt{-1}$, be the closed left half complex plane and let
$\mathbb{C}_{--}=\\{x+\mathrm{i}y:x<0,y\in\mathbb{R}\\}$ be its interior. The
$m$th powers of them are denoted by $\mathbb{C}_{-}^{m}$ and
$\mathbb{C}_{--}^{m}$, respectively. From the branching property (1.1), it
follows that $P(t,x,\cdot)$ is an infinitely divisible distribution on $D$. A
more careful analysis shows that the Laplace transform of $P(t,x,\cdot)$ is
characterized by, for $\lambda\in\mathbb{C}_{-}^{m}$,

(1.2)
$\displaystyle\int_{D}\mathrm{e}^{\langle\lambda,z\rangle}P(t,x,\mathrm{d}z)=\exp\big{\\{}\langle
x,K(t,\lambda)\rangle\big{\\}},$

where $\langle\cdot,\cdot\rangle$ denotes the inner product and
$K(t,\lambda)=(K_{1}(t,\lambda),\cdots,K_{m}(t,\lambda))\in\mathbb{C}_{-}^{m}$.
From (1.2) we have, for $i=1,\cdots,m$,

(1.3) $\displaystyle
K_{i}(t,\lambda)=\log\int_{D}\mathrm{e}^{\langle\lambda,z\rangle}P(t,\delta_{i},\mathrm{d}z),$

where $\delta_{i}=(1_{\\{i=1\\}},\cdots,1_{\\{i=m\\}})\in D$. The
Chapman–Kolmogorov equation of $\\{P(t):t\geq 0\\}$ implies

(1.4) $\displaystyle K(s+t,\lambda)=K(s,K(t,\lambda)),\quad s,t\geq
0,\lambda\in\mathbb{C}_{-}^{m}.$

We call $\\{K(t):t\geq 0\\}$ the cumulant semigroup of $\xi$ or $\\{P(t):t\geq
0\\}$. In view of (1.2), we have

(1.5) $\displaystyle P(t,x,D)=\exp\big{\\{}\langle
x,K(t,0)\rangle\big{\\}},\quad t\geq 0,x\in D.$

From the stochastic continuity of $\\{\xi(t):t\geq 0\\}$ it follows that
$t\mapsto K(t,\lambda)$ is continuous in $t\geq 0$ for each
$\lambda\in\mathbb{C}_{-}^{m}$ and the restriction of $\lambda\mapsto
K(t,\lambda)$ to $\mathbb{C}_{--}^{m}$ is a holomorphic self-map. In fact, the
infinite divisibility of the transition probabilities yields a Lévy–Khintchine
type representation of the function $\lambda\mapsto K(t,\lambda)$. For $m=1$
this is known as a Bernstein function; see, e.g., Schilling et al. [27]. It
was shown by Rhyzhov and Skorokhod [25] that for any
$\lambda\in\mathbb{C}_{--}^{m}$ the following derivative exists:

(1.6) $\displaystyle
H(\lambda)=\frac{\mathrm{d}}{\mathrm{d}t}K(t,\lambda)\Big{|}_{t=0}.$

Rhyzhov and Skorokhod [25] also derived from (1.4) and (1.6) the differential
equation, for $\lambda\in\mathbb{C}_{--}^{m}$,

(1.7)
$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}K(t,\lambda)=H(K(t,\lambda)),\quad
K(0,\lambda)=\lambda.$

This is referred to as the backward differential equation of the cumulant
semigroup since it corresponds to Kolmogorov’s backward differential equation
of the transition semigroup $\\{P(t):t\geq 0\\}$. A similar differential
equation for the generating functions of integer-valued branching processes
was given by Athreya and Ney [1, p.106].

We call $H=(H_{1},\cdots,H_{m})$ the branching mechanism of $\xi$. In complex
analysis, one usually refers to (1.7) as the Loewner–Kufarev equation; see,
e.g., Bracci et al. [6], Duren [11] and Gumenyuk et al. [16, 17]. It is also
called the generalized Riccati equation; see Duffie et al. [10] and the
references therein. For $x=(x_{1},\cdots,x_{m})\in D$ let
$|x|=|x_{1}|+\cdots+|x_{m}|$. In Rhyzhov and Skorokhod [25], the following
Lévy–Khintchine type representation of $H=(H_{1},\cdots,H_{m})$ was
identified: for $i=1,\cdots,m$,

(1.8) $\displaystyle H_{i}(\lambda)=\langle
a_{i},\lambda\rangle+\frac{1}{2}\beta_{i}\lambda_{i}^{2}+\int_{D^{\circ}}\Big{(}\mathrm{e}^{\langle\lambda,z\rangle}-1-\frac{\lambda_{i}z_{i}}{1+|z|^{2}}\Big{)}\pi_{i}(\mathrm{d}z),$

where $\beta_{i}\geq 0$, $a_{i}=(a_{i1},\cdots,a_{im})\in\mathbb{R}^{m}$
satisfies $a_{ij}\geq 0$ for $j\neq i$ and $\pi_{i}(\mathrm{d}z)$ is a
$\sigma$-finite measure on $D^{\circ}:=D\setminus\\{0\\}$ satisfying

(1.9) $\displaystyle\int_{D^{\circ}}\Big{[}\Big{(}\sum_{j=1,j\neq
i}^{m}z_{j}+z_{i}^{2}\Big{)}\land 1\Big{]}\pi_{i}(\mathrm{d}z)<\infty.$

The property $a_{ij}\geq 0$ for $j\neq i$ follows essentially from (1.6) and
the fact that $\lambda\mapsto K(t,\lambda)$ is the cumulant of an infinitely
divisible distribution on $\mathbb{R}_{+}^{m}$. It was derived rigorously in
Duffie et al. [10, Theorem 2.7]. Note that the property was stated incorrectly
as $a_{ij}>0$ for $j\neq i$ in [25, p.706]. The integrability condition (1.9)
also corrects the condition given in [25, p.706], where it was only claimed
that

(1.10) $\displaystyle\int_{\\{|z|<1\\}}\Big{(}\sum_{j=1,j\neq
i}^{m}z_{j}+z_{i}^{2}\Big{)}\pi_{i}(\mathrm{d}z)<\infty.$

It is easy to see that (1.9) is satisfied if and only if (1.10) holds and

(1.11) $\displaystyle\pi_{i}(\\{z\in D^{\circ}:|z|\geq 1\\})<\infty,$

which is a standard assumption for the Lévy measure; see, e.g., Sato [26,
Theorem 8.1]. For instance, if $\lambda=-(1,\cdots,1)$, then the integrand in
(1.8) becomes

$\displaystyle f(z)=\mathrm{e}^{-|z|}-1+\frac{z_{i}}{1+|z|^{2}},$

which satisfies

$\displaystyle-1\leq
f(z)\leq\mathrm{e}^{-|z|}-1+\frac{|z|}{1+|z|^{2}}\leq-\Big{(}\frac{1}{2}-\frac{1}{\mathrm{e}}\Big{)},\quad|z|\geq
1.$

Then the integral on the right-hand side of (1.8) diverges when (1.11) is not
satisfied. This shows that (1.9) is the correct assumption.

For the convenience of our exploration, we need another representation of the
branching mechanism $H=(H_{1},\cdots,H_{m})$. Let
$\alpha_{i}=(\alpha_{i1},\cdots,\alpha_{im})$, where $\alpha_{ij}=a_{ij}$ for
$j\neq i$ and

$\displaystyle\alpha_{ii}=a_{ii}+\int_{D^{\circ}}z_{i}\Big{(}1_{\\{|z|\leq
1\\}}-\frac{1}{1+|z|^{2}}\Big{)}\pi_{i}(\mathrm{d}z).$

Then, for $i=1,\cdots,m$,

(1.12) $\displaystyle
H_{i}(\lambda)=\langle\alpha_{i},\lambda\rangle+\frac{1}{2}\beta_{i}\lambda_{i}^{2}+\int_{D^{\circ}}\big{(}\mathrm{e}^{\langle\lambda,z\rangle}-1-\lambda_{i}z_{i}1_{\\{|z|\leq
1\\}}\big{)}\pi_{i}(\mathrm{d}z).$

It is known that the one-dimensional CB-process is conservative if and only if
the following condition is satisfied:

(1.13)
$\displaystyle\int_{0+}\frac{\mathrm{d}\lambda}{0\vee[-H(-\lambda)]}=\infty;$

see, e.g., [15, 18, 19]. There seems no necessary and sufficient condition in
such a simple form for the conservativeness of multi-dimensional CB-processes.

###### Example 1.1

Let us consider the one-dimensional case, i.e., $m=1$. For $\sigma>0$ and
$0<\alpha\leq 1$, an $\alpha$-stable branching mechanism is defined by
$H(\lambda)=-\sigma(-\lambda)^{\alpha}$. In this case, for
$\lambda\in\mathbb{C}_{--}$ the unique solution to (1.7) is given by

$\displaystyle
K(t,\lambda)=\bigg{\\{}\begin{array}[]{ll}-[\sigma(1-\alpha)t+(-\lambda)^{1-\alpha}]^{1/(1-\alpha)}\\!\\!\\!&~{}~{}\mbox{for
$0<\alpha<1$,}\cr\lambda\mathrm{e}^{\sigma t}\\!\\!\\!&~{}~{}\mbox{for
$\alpha=1$.}\end{array}$

This is a holomorphic function in $\mathbb{C}_{--}$ and has a unique
continuous extension on $\mathbb{C}_{-}$ with

$\displaystyle
K(t,0)=\bigg{\\{}\begin{array}[]{ll}-[\sigma(1-\alpha)t]^{1/(1-\alpha)}\\!\\!\\!&~{}~{}\mbox{for
$0<\alpha<1$,}\cr 0\\!\\!\\!&~{}~{}\mbox{for $\alpha=1$,}\end{array}$

which solves (1.7) for $\lambda=0$. The solution to (1.7) is unique for
$\lambda=0$ if and only if $\alpha=1$, as $K^{\infty}(t,0):=0$ also solves the
equation for $\lambda=0$.

###### Example 1.2

Let $m=1$ and consider the $\frac{1}{2}$-stable branching mechanism
$H(\lambda)=-2\sqrt{-\lambda}$. In this case, for $\lambda\in\mathbb{C}_{--}$
we have

$\displaystyle K(t,\lambda)=-(t+\sqrt{-\lambda})^{2},$

which maps $\mathbb{C}_{--}$ to the domain
$\\{x+\mathrm{i}y:x<-t\sqrt{t^{2}+2|y|}\\}$. The uniqueness of solution to
(1.7) breaks down for $\lambda=0$. In fact, for any $r\in[0,\infty]$ the map
$t\mapsto K^{r}(t,0):=-(t-r)^{2}1_{\\{t>r\\}}$ is a solution to (1.7) for
$\lambda=0$.

###### Remark 1.1

Let $H$ be a branching mechanism defined by (1.8) or (1.12). From the result
of Duffie et al. [10, Theorem 2.7 and Proposition 6.1] it follows that for any
$\lambda\in\mathbb{C}_{--}^{m}$ there is a unique solution to (1.7) and one
can define the transition semigroup of a CB-process by (1.2); see also
Watanabe [29, Theorems 2 and 2′]. By (1.3) we can extend $\lambda\mapsto
K(t,\lambda)$ to a continuous self-map of $\mathbb{C}_{-}^{m}$. For
$\lambda\in\mathbb{C}_{--}^{m}$ we can rewrite (1.7) into its integral form

(1.16) $\displaystyle
K(t,\lambda)=\lambda+\int_{0}^{t}H(K(s,\lambda))\mathrm{d}s,\quad t\geq 0.$

In view of (1.8) or (1.12), it is clear that the function $\lambda\mapsto
H(\lambda)$ admits a unique continuous extension on $\mathbb{C}_{-}^{m}$ with
$H(0)=0$. By approximation from the interior and application of dominated
convergence we see that (1.16) remains true for
$\lambda\in\partial\mathbb{C}_{-}^{m}:=\mathbb{C}_{-}^{m}\setminus\mathbb{C}_{--}^{m}$.
By differentiating both sides of the integral equation we see that (1.7) also
holds for all $\lambda\in\mathbb{C}_{-}^{m}$. That proves the existence of the
solution to (1.7) for all $\lambda\in\mathbb{C}_{-}^{m}$; see also Duffie et
al. [10, Proposition 6.4 and its proof].

Rhyzhov and Skorokhod [25, p.706] wrote: It can be shown that if (1.7) has a
unique solution belonging to $\mathbb{C}_{-}^{m}$ for all
$\lambda\in\mathbb{C}_{-}^{m}$, then $\exp\\{\langle x,K(t,\lambda)\rangle\\}$
will be a transform of type (1.1) of the transition probability of a branching
process. The following assertion is very plausible: in order that the equation
(1.7) have a unique solution for all $\lambda\in\mathbb{C}_{-}^{m}$, it is
necessary and sufficient that it have a unique solution for $\lambda=0$. In
view of Remark 1.1, we can restate their assertion equivalently as: in order
that the equation (1.7) have a unique solution for all
$\lambda\in\mathbb{C}_{-}^{m}$, it is necessary and sufficient that $t\mapsto
K(t,0)\equiv 0$ is the unique solution to (1.7) for $\lambda=0$.

In this paper, we provide a proof of the assertion of Rhyzhov and Skorokhod
[25, p.706]. The proof is based on a characterization of the $m$-dimensional
CB-process as the pathwise unique solution to a system of stochastic
equations.

##  2 The system of stochastic equations

Recall that the vector-valued function $H=(H_{1},\cdots,H_{m})$ is defined by
(1.12) with parameters $(a_{i},\beta_{i},\pi_{i})$, $i=1,\cdots,m$. Suppose
that $({\it\Omega},\mathscr{F},\mathscr{F}_{t},\mathbf{P})$ is a filtered
probability space satisfying the usual hypotheses. For each $k=1,\cdots,m$ let
$\\{B_{k}(t)\\}$ be an $(\mathscr{F}_{t})$-Brownian motion and let
$\\{N_{k}(\mathrm{d}s,\mathrm{d}z,\mathrm{d}u)\\}$ be a time-space
$(\mathscr{F}_{t})$-Poisson random measure on $(0,\infty)\times
D^{\circ}\times(0,\infty)$ with intensity
$\mathrm{d}s\pi_{k}(\mathrm{d}z)\mathrm{d}u$. Let
$\tilde{N}_{k}(\mathrm{d}s,\mathrm{d}z,\mathrm{d}u):=N_{k}(\mathrm{d}s,\mathrm{d}z,\mathrm{d}u)-\mathrm{d}s\pi_{k}(\mathrm{d}z)\mathrm{d}u$
denote the compensated measure. Let $D_{1}=\\{x\in D^{\circ}:|x|\leq 1\\}$ and
$D_{1}^{c}=\\{x\in D^{\circ}:|x|>1\\}$. We consider the system of stochastic
integral equations, for $k=1,\cdots,m$,

(2.3) $\displaystyle\xi_{k}(t)\\!\\!\\!$ $\displaystyle=\\!\\!\\!$
$\displaystyle\xi_{k}(0)+\int_{0}^{t}\sqrt{\beta_{k}\xi_{k}(s)}\mathrm{d}B_{k}(s)+\sum_{i=1}^{m}\int_{0}^{t}\alpha_{ik}\xi_{i}(s)\mathrm{d}s$
$\displaystyle+\int_{0}^{t}\int_{D_{1}}\int_{0}^{\xi_{k}(s-)}z_{k}\tilde{N}_{k}(\mathrm{d}s,\mathrm{d}z,\mathrm{d}u)+\int_{0}^{t}\int_{D_{1}^{c}}\int_{0}^{\xi_{k}(s-)}z_{k}N_{k}(\mathrm{d}s,\mathrm{d}z,\mathrm{d}u)$
$\displaystyle+\sum_{i=1,i\neq
k}^{m}\int_{0}^{t}\int_{D^{\circ}}\int_{0}^{\xi_{i}(s-)}z_{k}N_{i}(\mathrm{d}s,\mathrm{d}z,\mathrm{d}u).$

We say an $(\mathscr{F}_{t})$-adapted càdlàg process $\\{\xi(t):t\geq
0\\}=\\{(\xi_{1}(t),\cdots,\xi_{m}(t)):t\geq 0\\}$ in $D$ with possibly finite
lifetime is a solution to (2.3) if the equation holds almost surely when $t$
is replaced by $t\land\zeta_{n}$ for every $t\geq 0$ and $n\geq 1$, where
$\zeta_{n}=\inf\\{t\geq 0:|\xi(t)|\geq n\\}$ with $\inf\emptyset=\infty$ by
convention. We make the convention that $\xi(t)=\infty$ for
$t\geq\zeta:=\lim_{n\to\infty}\zeta_{n}$. By saying $\\{\xi(t)\\}$ is
conservative we mean $\zeta=\infty$ almost surely.

To prove the existence and uniqueness of the solution to (2.3), we introduce
an approximating sequence of the stochastic equation. For each $n\geq 1$ let
us consider the equations, for $k=1,\cdots,m$,

(2.6) $\displaystyle x_{nk}(t)\\!\\!\\!$ $\displaystyle=\\!\\!\\!$
$\displaystyle\xi_{k}(0)+\int_{0}^{t}\sqrt{\beta_{k}x_{nk}(s)}\mathrm{d}B_{k}(s)+\sum_{i=1}^{m}\int_{0}^{t}\alpha_{ik}x_{ni}(s)\mathrm{d}s$
$\displaystyle+\int_{0}^{t}\int_{D_{1}}\int_{0}^{x_{nk}(s-)}z_{k}\tilde{N}_{k}(\mathrm{d}s,\mathrm{d}z,\mathrm{d}u)+\int_{0}^{t}\int_{D_{1}^{c}}\int_{0}^{x_{nk}(s-)}(z_{k}\land
n)N_{k}(\mathrm{d}s,\mathrm{d}z,\mathrm{d}u)\quad$
$\displaystyle+\sum_{i=1,i\neq
k}^{m}\int_{0}^{t}\int_{D^{\circ}}\int_{0}^{x_{ni}(s-)}(z_{k}\land
n)N_{i}(\mathrm{d}s,\mathrm{d}z,\mathrm{d}u),$

The truncation $z_{k}\land n:=\min\\{z_{k},n\\}$ in (2.6) changes the jump
sizes larger than $n$ into jump size $n$. With this modification of the
equation, we can apply Barczy et al. [2, Theorem 4.6] to see that there is a
pathwise unique solution $\\{x_{n}(t):t\geq
0\\}=\\{(x_{n1}(t),\cdots,x_{nm}(t)):t\geq 0\\}$ to (2.6).

###### Theorem 2.1

For any $\mathscr{F}_{0}$-measurable random variable $\xi(0)\in D$, there is a
pathwise unique solution $\\{\xi(t):t\geq 0\\}$ to (2.3).

Proof. For $n\geq 1$ define the stopping time $\zeta_{n}=\inf\\{t\geq
0:|x_{n}(t)|\geq n\\}$. Then $|x_{n}(t)|<n$ for $0\leq t<\zeta_{n}$. Since the
processes $\\{x_{nk}(t):t\geq 0\\}$, $k=1,\cdots,m$ are nonnegative and can
only have positive jumps, they do not have jumps with jump sizes larger than
$n$ in the time interval $(0,\zeta_{n})$, so the truncation $(z_{k}\land n)$
in (2.6) makes no difference in this time interval. It follows that
$x_{n+1}(t)=x_{n}(t)$ for $0\leq t<\zeta_{n}$ and $\zeta_{n}$ is nondecreasing
in $n\geq 1$. Then we can define a process $\\{\xi(t):t\geq
0\\}=\\{(\xi_{1}(t),\cdots,\xi_{m}(t)):t\geq 0\\}$ in $\bar{D}$ such that
$\xi(t)=x_{n}(t)$ for $0\leq t<\zeta_{n}$ and $\xi(t)=\infty$ for
$t\geq\zeta:=\lim_{n\to\infty}\zeta_{n}$. It is easy to see that
$\\{\xi(t):t\geq 0\\}$ is a solution to (2.3) and $\zeta_{n}=\inf\\{t\geq
0:|\xi(t)|\geq n\\}$. This proves the existence of solution to (2.3). Now
suppose that $\\{y(t):t\geq 0\\}=\\{(y_{1}(t),\cdots,y_{m}(t)):t\geq 0\\}$ is
also a solution to (2.3). Let $\tau_{n}=\inf\\{t\geq 0:|y(t)|\geq n\\}$ and
$\eta_{n}=\tau_{n}\land\zeta_{n}$ for $n\geq 1$. Then $|y(t)|<n$ and
$|\xi(t)|<n$ for $0\leq t<\eta_{n}$, and so the processes $\\{y(t):t\geq 0\\}$
and $\\{\xi(t):t\geq 0\\}$ coincide in the time interval $(0,\eta_{n})$. In
view of (2.3) and (2.6) we have $\xi(\eta_{n})=y(\eta_{n})$, which implies
$\eta_{n}=\zeta_{n}=\tau_{n}$. It follows that
$\zeta=\lim_{n\to\infty}\zeta_{n}=\lim_{n\to\infty}\eta_{n}$, and so the
processes $\\{y(t):t\geq 0\\}$ and $\\{\xi(t):t\geq 0\\}$ are
indistinguishable. Then the pathwise uniqueness of solutions holds for (2.3).
$\square$

The system of stochastic equations (2.3) gives a construction of the sample
path of the $m$-dimensional CB-process. The existence and pathwise uniqueness
of solutions to (2.3) were established in Dawson and Li [8] for $m=1$ under a
stronger moment condition. A flow of discontinuous CB-processes was
constructed in Bertoin and Le Gall [5] by weak solutions to a special form of
the one-dimensional stochastic equation; see also [4]. Their construction was
extended to general flows in Dawson and Li [9] using strong solutions. See,
e.g., Bernis et al. [3], Fang and Li [14], Li [21, 22] and Pardoux [24] for a
number of applications of those stochastic equations. Theorem 2.1 weakens the
moment conditions in Barczy et al. [2] and Ma [23], where multi-dimensional
continuous state branching processes with immigration were constructed in
terms of stochastic equations.

##  3 A time-space martingale problem

Let $C_{0}^{2}(D)$ be the set of twice continuously differentiable functions
$f=f(x)$ on $D$ that together with their derivatives up to the second order
are rapidly decreasing as $|x|\to\infty$. We extend the functions $f=f(x)$ in
$C_{0}^{2}(D)$ and those derivatives trivially to $\infty$ so that they become
bounded continuous functions on $\bar{D}$. For $x\in D$ and $f\in
C_{0}^{2}(D)$ set

(3.2) $\displaystyle Af(x)\\!\\!\\!$ $\displaystyle=\\!\\!\\!$
$\displaystyle\sum_{k=1}^{m}x_{k}\langle\alpha_{k},\nabla
f(x)\rangle+\frac{1}{2}\sum_{k=1}^{m}\beta_{k}x_{k}f_{kk}^{\prime\prime}(x)$
$\displaystyle+\sum_{k=1}^{m}x_{k}\int_{D^{\circ}}\big{[}{\it\Delta}_{z}f(x)-z_{k}f^{\prime}_{k}(x)1_{D_{1}}(z)\big{]}\pi_{k}(\mathrm{d}z),$

where ${\it\Delta}_{z}f(x)=f(x+z)-f(x)$. Then we can understand $Af$ as a
bounded continuous function on $\bar{D}$ with $Af(\infty)=0$.

For a fixed constant $u>0$ let $C^{1,2}([0,u]\times D)$ be the set of
functions $f=f(t,x)$ on $[0,u]\times D$ that are continuously differentiable
in $t\in[0,u]$ and twice continuously differentiable in
$x=(x_{1},\cdots,x_{m})\in D$. For $f\in C^{1,2}([0,u]\times D)$ and
$i,j=1,\cdots,m$ write

$\displaystyle f_{0}^{\prime}(t,x)=\frac{\partial}{\partial
t}f(t,x),~{}f_{i}^{\prime}(t,x)=\frac{\partial}{\partial
x_{i}}f(t,x),~{}f_{ij}^{\prime\prime}(t,x)=\frac{\partial}{\partial
x_{i}x_{j}}f(t,x).$

Let $C_{0}^{1,2}([0,u]\times D)$ be the set of functions $f=f(t,x)$ in
$C^{1,2}([0,u]\times D)$ that together with the above derivatives are rapidly
decreasing as $|x|\to\infty$. We extend all the functions in
$C_{0}^{1,2}([0,u]\times D)$ and those derivatives trivially to
$[0,u]\times\\{\infty\\}$ so that they become bounded continuous functions on
$[0,u]\times\bar{D}$. The theorem below gives the characterization of the
solution $\\{\xi(t):t\geq 0\\}$ to (2.3) by a time-space martingale problem,
which identifies the operator $A$ defined by (3.2) as (a restriction of) its
generator.

###### Theorem 3.1

Let $\\{\xi(t):t\geq 0\\}$ be the pathwise unique solution to (2.3). Then for
any $0\leq t\leq u$ and $f\in C_{0}^{1,2}([0,u]\times D)$ we have

(3.3) $\displaystyle
f(t,\xi(t))=f(0,\xi(0))+\int_{0}^{t}\big{[}f_{0}^{\prime}(s,\xi(s))+Af(s,\xi(s))\big{]}\mathrm{d}s+M(t),$

where $A$ acts on the function $x\mapsto f(s,x)$ and $\\{M(t):0\leq t\leq
u\\}$ is a bounded martingale defined by

(3.5) $\displaystyle M(t)\\!\\!\\!$ $\displaystyle=\\!\\!\\!$
$\displaystyle\sum_{k=1}^{m}\int_{0}^{t}\sqrt{\beta_{k}\xi_{k}(s)}f_{k}^{\prime}(s,\xi(s))\mathrm{d}B_{k}(s)$
$\displaystyle+\sum_{k=1}^{m}\int_{0}^{t}\int_{D^{\circ}}\int_{0}^{\xi_{k}(s-)}{\it\Delta}_{z}f(s,\xi(s-))\tilde{N}_{k}(\mathrm{d}s,\mathrm{d}z,\mathrm{d}u).$

Proof. For $n\geq 1$ let $\\{x_{n}(t):t\geq
0\\}=\\{(x_{n1}(t),\cdots,x_{nm}(t)):t\geq 0\\}$ be the pathwise unique
solution to (2.6). By Itô’s formula, for $0\leq t\leq u$ we have

$\displaystyle
f(t,x_{n}(t))=f(0,\xi(0))+\int_{0}^{t}f_{0}^{\prime}(s,x_{n}(s-))\mathrm{d}s+\sum_{k=1}^{m}\sum_{i=1}^{m}\int_{0}^{t}\alpha_{ik}x_{ni}(s-)f_{k}^{\prime}(s,x_{n}(s-))\mathrm{d}s$
$\displaystyle\qquad\qquad+\sum_{k=1}^{m}\int_{0}^{t}\sqrt{\beta_{k}x_{nk}(s-)}f_{k}^{\prime}(s,x_{n}(s-))\mathrm{d}B_{k}(s)+\frac{1}{2}\sum_{k=1}^{m}\int_{0}^{t}\beta_{k}x_{nk}(s-)f_{kk}^{\prime\prime}(s,x_{n}(s-))\mathrm{d}s$
$\displaystyle\qquad\qquad+\sum_{k=1}^{m}\int_{0}^{t}\int_{D_{1}}\int_{0}^{x_{nk}(s-)}z_{k}f_{k}^{\prime}(s,x_{n}(s-))\tilde{N}_{k}(\mathrm{d}s,\mathrm{d}z,\mathrm{d}u)$
$\displaystyle\qquad\qquad+\sum_{k=1}^{m}\int_{0}^{t}\int_{D_{1}^{c}}\int_{0}^{x_{nk}(s-)}(z_{k}\land
n)f_{k}^{\prime}(s,x_{n}(s-))N_{k}(\mathrm{d}s,\mathrm{d}z,\mathrm{d}u)$
$\displaystyle\qquad\qquad+\sum_{k=1}^{m}\sum_{i=1,i\neq
k}^{m}\int_{0}^{t}\int_{D^{\circ}}\int_{0}^{x_{ni}(s-)}(z_{k}\land
n)f_{k}^{\prime}(s,x_{n}(s-))N_{i}(\mathrm{d}s,\mathrm{d}z,\mathrm{d}u)$
$\displaystyle\qquad\qquad+\sum_{k=1}^{m}\int_{0}^{t}\int_{D^{\circ}}\int_{0}^{x_{nk}(s-)}\big{[}f(s,x_{n}(s-)+(z\land
n))-f(s,x_{n}(s-))$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad-\,\langle
z\land n,\nabla
f(s,x_{n}(s-))\rangle\big{]}N_{k}(\mathrm{d}s,\mathrm{d}z,\mathrm{d}u)$
$\displaystyle\qquad\quad=f(0,\xi(0))+\int_{0}^{t}f_{0}^{\prime}(s,x_{n}(s-))\mathrm{d}s+\sum_{i=1}^{m}\int_{0}^{t}x_{ni}(s-)\langle\alpha_{i},\nabla
f(s,x_{n}(s-))\rangle\mathrm{d}s$
$\displaystyle\qquad\qquad+\sum_{k=1}^{m}\int_{0}^{t}\sqrt{\beta_{k}x_{nk}(s-)}f_{k}^{\prime}(s,x_{n}(s-))\mathrm{d}B_{k}(s)+\frac{1}{2}\sum_{k=1}^{m}\int_{0}^{t}\beta_{k}x_{nk}(s-)f_{kk}^{\prime\prime}(s,x_{n}(s-))\mathrm{d}s$
$\displaystyle\qquad\qquad+\sum_{k=1}^{m}\int_{0}^{t}\int_{D_{1}}\int_{0}^{x_{nk}(s-)}z_{k}f_{k}^{\prime}(s,x_{n}(s-))\tilde{N}_{k}(\mathrm{d}s,\mathrm{d}z,\mathrm{d}u)$
$\displaystyle\qquad\qquad+\sum_{k=1}^{m}\int_{0}^{t}\int_{D^{\circ}}\int_{0}^{x_{nk}(s-)}\big{[}f(s,x_{n}(s-)+(z\land
n))-f(s,x_{n}(s-))$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad-\,z_{k}f_{k}^{\prime}(s,x_{n}(s-))1_{D_{1}}(z)\big{]}N_{k}(\mathrm{d}s,\mathrm{d}z,\mathrm{d}u)$
$\displaystyle\qquad\quad=f(0,\xi(0))+\int_{0}^{t}f_{0}^{\prime}(s,x_{n}(s-))\mathrm{d}s+\sum_{i=1}^{m}\int_{0}^{t}x_{ni}(s-)\langle\alpha_{i},\nabla
f(s,x_{n}(s-))\rangle\mathrm{d}s$
$\displaystyle\qquad\qquad+\sum_{k=1}^{m}\int_{0}^{t}\sqrt{\beta_{k}x_{nk}(s-)}f_{k}^{\prime}(s,x_{n}(s-))\mathrm{d}B_{k}(s)+\frac{1}{2}\sum_{k=1}^{m}\int_{0}^{t}\beta_{k}x_{nk}(s-)f_{kk}^{\prime\prime}(s,x_{n}(s-))\mathrm{d}s$
$\displaystyle\qquad\qquad+\sum_{k=1}^{m}\int_{0}^{t}\int_{D^{\circ}}\int_{0}^{x_{nk}(s-)}{\it\Delta}_{z\land
n}f(s,x_{n}(s-))\tilde{N}_{k}(\mathrm{d}s,\mathrm{d}z,\mathrm{d}u)$
$\displaystyle\qquad\qquad+\sum_{k=1}^{m}\int_{0}^{t}x_{nk}(s-)\mathrm{d}s\int_{D^{\circ}}\big{[}{\it\Delta}_{z\land
n}f(s,x_{n}(s-))-z_{k}f_{k}^{\prime}(s,x_{n}(s-))1_{D_{1}}(z)\big{]}\pi_{k}(\mathrm{d}z),$

where $z\land n=(z_{1}\land n,\cdots,z_{m}\land n)$. Then

$\displaystyle
f(t,x_{n}(t))=f(0,\xi(0))+\int_{0}^{t}\big{[}f_{0}^{\prime}(s,x_{n}(s-))+A_{n}f(s,x_{n}(s-))\big{]}\mathrm{d}s+M_{n}(t),$

where

$\displaystyle M_{n}(t)\\!\\!\\!$ $\displaystyle=\\!\\!\\!$
$\displaystyle\sum_{k=1}^{m}\int_{0}^{t}\sqrt{\beta_{k}x_{nk}(s-)}f_{k}^{\prime}(s,x_{n}(s-))\mathrm{d}B_{k}(s)$
$\displaystyle+\sum_{k=1}^{m}\int_{0}^{t}\int_{D^{\circ}}\int_{0}^{x_{nk}(s-)}{\it\Delta}_{z\land
n}f(s,x_{n}(s-))\tilde{N}_{k}(\mathrm{d}s,\mathrm{d}z,\mathrm{d}u)$

and

$\displaystyle A_{n}f(s,x)\\!\\!\\!$ $\displaystyle=\\!\\!\\!$
$\displaystyle\sum_{k=1}^{m}x_{k}\langle\alpha_{k},\nabla
f(s,x)\rangle+\frac{1}{2}\sum_{k=1}^{m}\beta_{k}x_{k}f_{kk}^{\prime\prime}(s,x)$
$\displaystyle+\sum_{k=1}^{m}x_{k}\int_{D^{\circ}}\big{[}{\it\Delta}_{z\land
n}f(s,x)-z_{k}f_{k}^{\prime}(s,x)1_{D_{1}}(z)\big{]}\pi_{k}(\mathrm{d}z).$

Let $\zeta_{n}=\inf\\{t\geq 0:|x_{n}(t)|\geq n\\}$ for $n\geq 1$. By the proof
of Theorem 2.1, the sequence $\\{\zeta_{n}\\}$ is nondecreasing. Moreover, we
have $\xi(t)=x_{n}(t)$ for $0\leq t<\zeta_{n}$ and $\xi(t)=\infty$ for
$t\geq\zeta:=\lim_{n\to\infty}\zeta_{n}$. It follows that

(3.6) $\displaystyle
f(t,x_{n}(t\land\zeta_{n}))=f(0,\xi(0))+\int_{0}^{t\land\zeta_{n}}\big{[}f_{0}^{\prime}(s,\xi(s-))+A_{n}f(s,\xi(s-))\big{]}\mathrm{d}s+M_{n}(t\land\zeta_{n}),$

where

(3.8) $\displaystyle M_{n}(t\land\zeta_{n})\\!\\!\\!$
$\displaystyle=\\!\\!\\!$
$\displaystyle\sum_{k=1}^{m}\int_{0}^{t\land\zeta_{n}}\sqrt{\beta_{k}\xi_{k}(s-)}f_{k}^{\prime}(s,\xi(s-))\mathrm{d}B_{k}(s)$
$\displaystyle+\sum_{k=1}^{m}\int_{0}^{t\land\zeta_{n}}\int_{D^{\circ}}\int_{0}^{\xi_{k}(s-)}{\it\Delta}_{z\land
n}f(s,\xi(s-))\tilde{N}_{k}(\mathrm{d}s,\mathrm{d}z,\mathrm{d}u).$

The definition of $\zeta_{n}$ implies that
$\lim_{n\to\infty}x_{n}(\zeta_{n})=\infty$. By considering the cases
$\zeta\leq t$ and $\zeta>t$ separately, we see that
$\lim_{n\to\infty}x_{n}(t\land\zeta_{n})=\xi(t\land\zeta)$. Then letting
$n\to\infty$ in (3.6) and (3.8) we obtain

$\displaystyle
f(t,\xi(t\land\zeta))=f(0,\xi(0))+\int_{0}^{t\land\zeta}\big{[}f_{0}^{\prime}(s,\xi(s-))+Af(s,\xi(s-))\big{]}\mathrm{d}s+M_{\zeta}(t),$

where

$\displaystyle M_{\zeta}(t)\\!\\!\\!$ $\displaystyle=\\!\\!\\!$
$\displaystyle\sum_{k=1}^{m}\int_{0}^{t\land\zeta}\sqrt{\beta_{k}\xi_{k}(s-)}f_{k}^{\prime}(s,\xi(s-))\mathrm{d}B_{k}(s)$
$\displaystyle+\sum_{k=1}^{m}\int_{0}^{t\land\zeta}\int_{D^{\circ}}\int_{0}^{\xi_{k}(s-)}{\it\Delta}_{z}f(s,\xi(s-))\tilde{N}_{k}(\mathrm{d}s,\mathrm{d}z,\mathrm{d}u).$

Recall that we understand $f$, $f_{k}^{\prime}$ and $Af$ as bounded continuous
functions on $[0,u]\times\bar{D}$ with
$f(s,\infty)=f_{k}^{\prime}(s,\infty)=Af(s,\infty)=0$ for all $s\in[0,u]$.
Then we have (3.3) and (3.5). From (3.2) it is clear that

$\displaystyle|M(t)|\leq 2\|f\|+(\|f_{0}^{\prime}\|+\|Af\|)u,\quad 0\leq t\leq
u.$

where $\|\cdot\|$ denotes the supremum norm of functions on
$[0,u]\times\bar{D}$. Then $\\{M(t):0\leq t\leq u\\}$ is a bounded martingale
by (3.5). $\square$

###### Corollary 3.2

The pathwise unique solution $\\{\xi(t):t\geq 0\\}$ to (2.3) is a càdlàg
process in $\bar{D}$.

Proof. By applying Theorem 3.1 to the function $f(t,x)\equiv\mathrm{e}^{-x}$
and $u=1,2,\cdots$ we see that $\\{\mathrm{e}^{-\xi(t)}:t\geq 0\\}$ is a
càdlàg process in $[0,1]$. Then $\\{\xi(t):t\geq 0\\}$ is a càdlàg process in
$\bar{D}$. $\square$

###### Corollary 3.3

The pathwise unique solution $\\{\xi(t):t\geq 0\\}$ to (2.3) is an
$m$-dimensional CB-process with transition semigroup $\\{P(t):t\geq 0\\}$
defined by (1.2) and (1.7).

Proof. For $\lambda\in\mathbb{C}_{--}^{m}$ let $t\mapsto
K(t,\lambda)\in\mathbb{C}_{--}^{m}$ be the unique solution to (1.7). Fix $u>0$
and define

$\displaystyle f(t,x)=\mathrm{e}^{\langle K(u-t,\lambda),x\rangle},\quad 0\leq
t\leq u,x\in D.$

Then $f\in C_{0}^{1,2}([0,u]\times D)$. It is easy to show that for
$k=1,\cdots,m$,

$\displaystyle f_{k}^{\prime}(t,x)=K_{k}(u-t,\lambda)\mathrm{e}^{\langle
K(u-t,\lambda),x\rangle},~{}f_{kk}^{\prime\prime}(t,x)=K_{k}(u-t,\lambda)^{2}\mathrm{e}^{\langle
K(u-t,\lambda),x\rangle}$

and

$\displaystyle{\it\Delta}_{z}f(x)=(\mathrm{e}^{\langle
K(u-t,\lambda),z\rangle}-1)\mathrm{e}^{\langle K(u-t,\lambda),x\rangle}.$

By (1.7) and (3.2) we have

$\displaystyle f_{0}^{\prime}(t,x)=-\mathrm{e}^{\langle
K(u-t,\lambda),x\rangle}\langle H(K(u-t,\lambda)),x\rangle=-Af(t,x).$

Note that the equalities above extend trivially to $x\in\bar{D}$. An
application of Theorem 3.1 shows that

$\displaystyle\mathrm{e}^{\langle
K(u-t,\lambda),\xi(t)\rangle}=\mathrm{e}^{\langle
K(u,\lambda),\xi(0)\rangle}+M(t),\quad 0\leq t\leq u,$

where $\\{M(t):0\leq t\leq u\\}$ is a bounded martingale. It follows that

$\displaystyle\mathbf{E}\big{[}\mathrm{e}^{\langle\lambda,\xi(u)\rangle}|\mathscr{F}_{t}\big{]}\\!\\!\\!$
$\displaystyle=\\!\\!\\!$ $\displaystyle\mathbf{E}\big{[}\mathrm{e}^{\langle
K(u,\lambda),\xi(0)\rangle}+M(u)|\mathscr{F}_{t}\big{]}$
$\displaystyle=\\!\\!\\!$ $\displaystyle\mathrm{e}^{\langle
K(u,\lambda),\xi(0)\rangle}+M(t)=\mathrm{e}^{\langle
K(u-t,\lambda),\xi(t)\rangle}.$

Then $\\{\xi(t):t\geq 0\\}$ is an $m$-dimensional CB-process with transition
semigroup $\\{P(t):t\geq 0\\}$ defined by (1.2) and (1.7). $\square$

###### Remark 3.4

In general, the process $\\{\xi(t):t\geq 0\\}$ lives in the extended state
space $\bar{D}$. For this reason, we need the functions $f$ and $Af$ to have
continuous extensions on $[0,u]\times\bar{D}$ in Theorem 3.1 and in the proof
of Corollary 3.3. This is the reason that we assume
$\lambda\in\mathbb{C}_{--}^{m}$ in the proof.

##  4 The backward differential equation

By Remark 1.1, there exists the solution $t\mapsto K(t,\lambda)$ to (1.7) for
every $\lambda\in\mathbb{C}_{-}^{m}$. The next theorem confirms the assertion
of Rhyzhov and Skorokhod [25, p.706]:

###### Theorem 4.1

The following statements are equivalent:

  * (i)

For every $\lambda\in\mathbb{C}_{-}^{m}$ there is a unique solution $t\mapsto
K(t,\lambda)\in\mathbb{C}_{-}^{m}$ to (1.7).

  * (ii)

For $\lambda=0$ the unique solution to (1.7) is given by $t\mapsto
K(t,0)\equiv 0$.

  * (iii)

The CB-process with transition semigroup given by (1.2) and (1.7) is
conservative.

Proof. (i)$\Rightarrow$(ii) This is trivial.

(ii)$\Rightarrow$(iii) Suppose that $t\mapsto K(t,0)\equiv 0$ is the unique
solution to (1.7) for $\lambda=0$. Then (1.5) implies $P(t,x,D)=1$ for all
$t\geq 0$ and $x\in D$, which means the CB-process is conservative.

(iii)$\Rightarrow$(i) Suppose that the CB-process with transition semigroup
given by (1.2) and (1.7) is conservative. Let $\\{\xi(t):t\geq 0\\}$ be the
realization of the process defined by (2.3) with $\xi(0)=x\in D$. Then
$\\{\xi(t):t\geq 0\\}$ is a càdlàg process living in $D$. Suppose that
$t\mapsto J(t,\lambda)$ is also a solution to (1.7) for
$\lambda\in\mathbb{C}_{-}^{m}$. Fix $u>0$ and define

(4.1) $\displaystyle g(t,x)=\mathrm{e}^{\langle J(u-t,\lambda),x\rangle},\quad
0\leq t\leq u,x\in D.$

It is easy to see that $g\in C^{1,2}([0,u]\times D)$. By (2.3) and Itô’s
formula, as in the proof of Theorem 3.1 we get

(4.4) $\displaystyle
g(t,\xi(t))=g(0,\xi(0))+\int_{0}^{t}g_{0}^{\prime}(s,\xi(s-))\mathrm{d}s+\sum_{i=1}^{m}\int_{0}^{t}\xi_{i}(s-)\langle\alpha_{i},\nabla
g(s,\xi(s-))\rangle\mathrm{d}s$
$\displaystyle\qquad\qquad+\frac{1}{2}\sum_{k=1}^{m}\int_{0}^{t}\beta_{k}\xi_{k}(s-)g_{kk}^{\prime\prime}(s,\xi(s-))\mathrm{d}s+Z(t)$
$\displaystyle\qquad\qquad+\sum_{k=1}^{m}\int_{0}^{t}\xi_{k}(s-)\mathrm{d}s\int_{D^{\circ}}\big{[}{\it\Delta}_{z}g(s,\xi(s-))-z_{k}g_{k}^{\prime}(s,\xi(s-))1_{D_{1}}(z)\big{]}\pi_{k}(\mathrm{d}z),$

where

$\displaystyle Z(t)\\!\\!\\!$ $\displaystyle=\\!\\!\\!$
$\displaystyle\sum_{k=1}^{m}\int_{0}^{t}\sqrt{\beta_{k}\xi_{k}(s-)}g_{k}^{\prime}(s,\xi(s-))\mathrm{d}B_{k}(s)$
$\displaystyle+\sum_{k=1}^{m}\int_{0}^{t}\int_{D^{\circ}}\int_{0}^{\xi_{k}(s-)}{\it\Delta}_{z}g(s,\xi(s-))\tilde{N}_{k}(\mathrm{d}s,\mathrm{d}z,\mathrm{d}u).$

Since $t\mapsto J(t,\lambda)$ solves (1.7), one can use (4.1) to see

$\displaystyle g_{0}^{\prime}(s,x)\\!\\!\\!$ $\displaystyle=\\!\\!\\!$
$\displaystyle-\sum_{i=1}^{m}x_{i}\langle\alpha_{i},\nabla
g(s,x)\rangle-\frac{1}{2}\sum_{k=1}^{m}\beta_{k}x_{k}g_{kk}^{\prime\prime}(s,x)$
$\displaystyle-\sum_{k=1}^{m}x_{k}\int_{D^{\circ}}\big{[}{\it\Delta}_{z}g(s,x)-z_{k}g_{k}^{\prime}(s,x)1_{D_{1}}(z)\big{]}\pi_{k}(\mathrm{d}z),$

From (4.4) it follows that

(4.5) $\displaystyle\mathrm{e}^{\langle
J(u-t,\lambda),\xi(t)\rangle}=\mathrm{e}^{\langle
J(u,\lambda),x\rangle}+Z(t).$

Then $\\{Z(t):0\leq t\leq u\\}$ is a bounded martingale. (Here we cannot apply
Theorem 3.1 directly since $g\notin C_{0}^{1,2}([0,u]\times D)$.) By taking
the expectations in both sides of (4.5) for $t=u$ we obtain
$\mathbf{E}[\mathrm{e}^{\langle\lambda,\xi(u)\rangle}]=\mathrm{e}^{\langle
J(u,\lambda),x\rangle}$. Then $J(u,\lambda)=K(u,\lambda)$ by (1.2) and
Corollary 3.3. This shows the implication “(iii)$\Rightarrow$(i)”. $\square$

It would be interesting to have an analytic proof of the assertion of Rhyzhov
and Skorokhod [25, p.706].

Acknowledgements. We would like to thank the two anonymous referees for their
valuable comments and suggestions, which have led to a number of improvements
in the presentation of the work.

## References

  * [1] Athreya, K.B. and Ney, P.E. (1972): Branching Processes. Springer, Berlin.
  * [2] Barczy, M.; Li, Z. and Pap, G. (2015): Stochastic differential equation with jumps for multi-type continuous state and continuous time branching processes with immigration. ALEA Lat. Am. J. Probab. Math. Stat. 12, 1: 129–169.
  * [3] Bernis, G.; Brignone, R.; Scotti, S. and Sgarra, C. (2021): A Gamma Ornstein-Uhlenbeck model driven by a Hawkes process. Math. Financ. Econ. 15, 747–773.
  * [4] Bertoin, J. and Le Gall, J.-F. (2000): The Bolthausen–Sznitman coalescent and the genealogy of continuous-state branching processes. Probab. Theory Related Fields 117, 249–266.
  * [5] Bertoin, J. and Le Gall, J.-F. (2006): Stochastic flows associated to coalescent processes III: Limit theorems. Illinois J. Math. 50, 147–181.
  * [6] Bracci, F.; Contreras, M.D. and Díaz-Madrigal, S. (2020): Continuous semigroups of holomorphic functions in the unit disc. Springer, Switzerland.
  * [7] Dawson, D.A. (1993): Measure-valued Markov processes. In: Ecole d’Eté de Probabilités de Saint-Flour XXI-1991, 1–260. Lecture Notes Math. 1541. Springer, Berlin.
  * [8] Dawson, D.A. and Li, Z. (2006): Skew convolution semigroups and affine Markov processes. Ann. Probab. 34, 1103–1142.
  * [9] Dawson, D.A. and Li, Z. (2012): Stochastic equations, flows and measure-valued processes. Ann. Probab. 40, 813–857.
  * [10] Duffie, D.; Filipović, D. and Schachermayer, W. (2003): Affine processes and applications in finance. Ann. Appl. Probab. 13, 984–1053.
  * [11] Duren, P.L. (1983): Univalent functions. Springer, New York.
  * [12] Dynkin, E.B. (1994): An Introduction to Branching Measure-Valued Processes. Amer. Math. Soc., Providence, RI.
  * [13] Dynkin, E.B., Kuznetsov, S.E. and Skorokhod, A.V. (1994): Branching measure-valued processes. Probab. Theory Related Fields 99, 55–96.
  * [14] Fang, R. and Li, Z. (2022): Construction of continuous-state branching processes in varying environments. Ann. Appl. Probab. 32, 3645–3673.
  * [15] Grey, D.R. (1974): Asymptotic behavior of continuous time, continuous state-space branching processes. J. Appl. Probab. 11, 669–677.
  * [16] Gumenyuk, P.; Hasebe, T. and Pérez, J.-L. (2024): Loewner Theory for Bernstein functions I: evolution families and differential equations. Constr. Approx. Available online: doi.org/10.1007/s00365-023-09675-9.
  * [17] Gumenyuk, P.; Hasebe, T. and Pérez, J.-L. (2022+): Loewner Theory for Bernstein functions II: applications to inhomogeneous continuous-state branching processes. arXiv:2211.12442.
  * [18] Ikeda, N. and Watanabe, S. (1970): On uniqueness and non-uniqueness of solutions for a class of non-linear equations and explosion problem for branching processes. J. Fac. Sci. Univ. Tokyo Sect. I 167, 187–214.
  * [19] Kawazu, K. and Watanabe, S. (1971): Branching processes with immigration and related limit theorems. Theory Probab. Appl. 16, 36–54.
  * [20] Le Gall, J.-F. (1999): Spatial Branching Processes, Random Snakes and Partial Differential Equations. Lectures in Mathematics ETH Zürich, Birkhäuser, Basel.
  * [21] Li, Z. (2020): Continuous-state branching processes with immigration. A Chapter in: From Probability to Finance, pp. 1–69, edited by Y. Jiao. Mathematical Lectures from Peking University. Springer, Singapore.
  * [22] Li, Z. (2022): Measure-Valued Branching Markov Processes. 2nd Ed. Springer, Berlin.
  * [23] Ma, R. (2013): Stochastic equations for two-type continuous-state branching processes with immigration. Acta Math. Sci. English Ed. 29, 287–294.
  * [24] Pardoux, E. (2016). Probabilistic Models of Population Evolution: Scaling Limits, Genealogies and Interactions. Springer, Switzerland.
  * [25] Rhyzhov, Y.M. and Skorokhod, A.V. (1970): Homogeneous branching processes with a finite number of types and continuous varying mass. Theory Probab. Appl. 15, 704–707.
  * [26] Sato, K. (1999): Lévy Processes and Infinitely Divisible Distributions. Cambridge Univ. Press, Cambridge.
  * [27] Schilling, R.L.; Song, R. and Vondraček, Z. (2012): Bernstein Functions: Theory and Applications. 2nd Ed. De Gruyter, Berlin/Boston.
  * [28] Silverstein, M.L. (1968): A new approach to local time. J. Math. Mech. 17, 1023–1054.
  * [29] Watanabe, S. (1969): On two dimensional Markov processes with branching property. Trans. Amer. Math. Soc. 136, 447–466.


11institutetext:  Institute
University of East Anglia 11email<EMAIL_ADDRESS>
# A tale of two toolkits, report the third: on the usage and performance of
HIVE-COTE v1.0

Anthony Bagnall  Michael Flynn  James Large  Jason Lines  Matthew Middlehurst

###### Abstract

The Hierarchical Vote Collective of Transformation-based Ensembles (HIVE-COTE)
is a heterogeneous meta ensemble for time series classification. Since it was
first proposed in 2016, the algorithm has undergone some minor changes and
there is now a configurable, scalable and easy to use version available in two
open source repositories. We present an overview of the latest stable HIVE-
COTE, version 1.0, and describe how it differs to the original. We provide a
walkthrough guide of how to use the classifier, and conduct extensive
experimental evaluation of its predictive performance and resource usage. We
compare the performance of HIVE-COTE to three recently proposed algorithms.

###### Keywords:

Time series Classification Heterogeneous ensembles HIVE-COTE

##  1 Introduction

The Hierarchical Vote Collective of Transformation-based Ensembles (HIVE-COTE)
is a heterogeneous meta ensemble for time series classification. The key
principle behind HIVE-COTE is that time series classification (TSC) problems
are best approached by careful consideration of the data representation, and
that with no expert knowledge to the contrary, the most accurate algorithm
design is to ensemble classifiers built on different representations.

HIVE-COTE was first described in an ICDM paper [14] which was later expanded
to become an ACM Transactions paper [15]. At the time, HIVE-COTE was
significantly more accurate on average than other known approaches [3] on the
85 datasets that were then the complete UCR archive [6]. HIVE-COTE was an
improvement over the 2015 version, called just COTE on publication [4] but
later renamed Flat-COTE to differentiate it from its successor. Flat-COTE is a
standard ensemble of a range of classifiers built on different
representations. It was itself a natural extension of the Elastic Ensemble
[13] which only contains nearest neighbour classifiers using different
distance measures. HIVE-COTE takes a more structured approach than Flat-COTE.
The original HIVE-COTE, which we will henceforth refer to as HIVE-COTE alpha,
contained the following classification modules: Elastic Ensemble (EE) [13];
Shapelet Transform Classifier (STC) [11]; Time Series Forest (TSF) [8]; and
Bag of Symbolic-Fourier Approximation Symbols (BOSS) [19]. Each module is
encapsulated and built on the train data independently of the others. For new
data, each module passes an estimate of class probabilities to the control
unit, which combines them to form a single prediction. It does this by
weighting the probabilities of each module by an estimate of its testing
accuracy formed from the training data.

Our goal with HIVE-COTE alpha was to achieve the highest level of accuracy
without concern for the computational resources. This has lead to the
perception that HIVE-COTE is very slow and does not scale well. Whilst this is
true if used in its basic form, it is in fact very simple to restructure HIVE-
COTE so it achieves the same level of accuracy in orders of magnitude less
time. We have made many small changes to HIVE-COTE with the goal of making it
scalable and more useful. We describe these improvements and encapsulate them
as HIVE-COTE version 1.0. The changes in HIVE-COTE are both algorithmic and
engineering in nature.

The two slowest components of HIVE-COTE alpha are STC (in its old format) and
EE. STC used to conduct a full enumeration of all possible shapelets. We have
found that this enormous computational effort is not only unnecessary, but
often results in over fitting. EE requires cross validation of numerous
nearest neighbour classifiers and is very slow on training and testing. EE
resulted from of a comparative study of nearest neighbour distance measures.
Our hypothesis was there was no significant difference between the numerous
distance measures and dynamic time warping when used in nearest neighbour
classifiers, which is true. We only ensembled as an afterthought. We were
surprised to see significant improvement. Its design was necessarily ad hoc to
avoid over fitting. We have found that dropping EE all together does not make
HIVE-COTE much worse. The main changes are:

  1. 1.

STC no longer fully enumerates the shapelet space.

  2. 2.

EE is dropped altogether from HIVE-COTE.

  3. 3.

The STC, BOSS and RISE components include revisions to improve efficiency.

  4. 4.

All components and HIVE-COTE 1.0 itself are now contractable (you can set a
run time limit), checkpointable (you can save a version to continue building
later) and tuneable (select parameters based on train set performance).

  5. 5.

HIVE-COTE 1.0 can be threaded, built from existing results and easily
configured.

The aim of this report is to describe the changes, showcase the usage of HIVE-
COTE and to present some new benchmark results that should be used in all
future experiments, and to demonstrate the scalability of HIVE-COTE.

##  2 HIVE-COTE 1.0 Design

Figure 1 provides an overview of the HIVE-COTE structure. The top level
ensemble structure and the implementation of each component are described
below.

Figure 1: An overview of the ensemble structure of HIVE-COTE 1.0. Each module
produces an estimate of the probability of membership of each class. The
control unit (CAWPE) combines these probabilities, weighted by an estimate of
the quality of the module found on the train data.

###  2.1 Ensemble Structure

HIVE-COTE adopts the Cross-validation Accuracy Weighted Probabilistic Ensemble
(CAWPE) ensemble structure [12], summmarised in Algorithm 1.

Algorithm 1 HIVE-COTE classify(A test case $\bf{x}$)

0: prediction for case $\bf{x}$

0: A set of classifiers $<M_{1},\ldots,M_{k}>$, an exponent $\alpha$, a set of
weights $w_{i}$, and the number of classes $c$

1: $\\{\hat{p}_{1},\ldots,\hat{p}_{c}\\}=\\{0,\ldots,0\\}$

2: for $i\leftarrow 1$ to $k$ do

3: for $j\leftarrow 1$ to $c$ do

4: $\hat{{q_{j}}}\leftarrow\hat{p}((y=j|M_{i},\bf{x})$

5: $\hat{p}_{j}\leftarrow\hat{p}_{j}+w_{i}^{\alpha}\cdot\hat{q}_{j}$

6: return $argmax_{j=1\ldots c}\hat{p}_{j}$

CAWPE uses an estimate of the accuracy of each classifier to weight the
probability estimates of each component. It constructs a tilted distribution
through exponentiation (using $\alpha$) to extenuate differences in
classifiers. The weight for each component is found either through ten fold
cross validation, or, if the classifier has the capability to estimate its own
performance, internally.

###  2.2 Time Series Forest (TSF)

TSF [8] aims to capture basic summary features from intervals of a time
series. For any given time series of length $m$ there are $m(m-1)/2$ possible
intervals that can be extracted. TSF takes a random forest-like approach to
sampling these intervals. A formal description of TSF is provided in Algorithm
2.

For each tree, $r$ intervals are randomly selected (lines 5-7), each with a
random start position and length. Each interval is summarised by the mean,
standard deviation and slope (lines 8-11), and the summaries of each interval
are concatenated into a single feature vector of length $3r$ for each time
series. A decision tree is built on this concatenated feature vector (line
12). New cases are classified using a majority vote of all trees in the
forest. The version of TSF used in the bake off [3] employed the random tree
used by random forest. However, the decision tree described in [8] has some
minor differences to the random tree. It makes no difference in terms of
accuracy, but the tree from [8], the time series tree, has advantages in terms
of interpretability. Hence, the current version of TSF uses the original
version, which we call TimeSeriesTree. Contracting is enforced by the method
timeRemaining. TSF simply builds until it has the required number of trees or
the time runs out.

Algorithm 2 buildTSF(A list of $n$ cases length $m$, ${\bf T}=({\bf X,y})$)

0: the number of trees, $k$; the minimum interval length, $p$; the number of
intervals per tree, $r$. (default $k\leftarrow 500$, $p\leftarrow 3$, and
$r\leftarrow\sqrt{m}$)

1: Let ${\bf F}=({\bf F_{1}}\ldots{\bf F_{k}})$ be the trees in the forest

2: $i\leftarrow 1$

3: while $i<k$ $\and$ timeRemaining() do

4: Let ${\bf S}$ be a list of $n$ cases $({\bf s_{1}}\ldots{\bf s_{n}})$ with
$3r$ attributes

5: for $j\leftarrow 1$ to $r$ do

6: $b\leftarrow$ randBetween$(1,m-p)$

7: $e\leftarrow$ randBetween$(b+p,m)$

8: for $t\leftarrow 1$ to $n$ do

9: $s_{t,3(j-1)+1}\leftarrow$ mean$({\bf x}_{t},b,e)$

10: $s_{t,3(j-1)+2}\leftarrow$ standardDeviation$(\bf{x}_{t},b,e)$

11: $s_{t,3(j-1)+3}\leftarrow$ slope$({\bf x}_{t},b,e)$

12: $\bf{F_{i}}.$buildTimeSeriesTree(${\bf S},y$)

###  2.3 Random Interval Spectral Ensemble (RISE)

Like TSF, RISE [15] is a tree based interval ensemble. Unlike TSF, it uses a
single interval for each tree, and it uses spectral features rather than
summary statistics. RISE was recently updated to be faster and contractable
[10]. During the build process, summarised in Algorithm 3, a single random
interval is selected for each tree. The first tree is a special case in which
the whole series is used (lines 5 and 6). Otherwise, an interval with length
that is a power of 2 (line 9) is chosen. The same interval for each series is
then transformed using the Fast Fourier Transform (FFT) and Auto Correlation
Function (ACF). This is a change on the original RISE which also used the
partial autocorrelation function (PACF) and autoregressive (AR) model
features. Restriction to the Fourier and ACF coefficients does not decrease
accuracy, but makes the algorithm much faster. The power spectrum coefficients
are concatenated with the first 100 ACF coefficients to form a new training
set. In the tsml implementation of RISE the base classifier used is the
RandomTree classifier used by random forest (line 13). In the test process
class probabilities are assigned as a proportion of base classifier votes.

RISE controls the contract run time by creating an adaptive model of the time
to build a single tree (lines 4 and 14). This is important for long series
(such as audio), where very large intervals can mean very few trees. Details
are in [10].

Algorithm 3 buildRISE(A list of $n$ cases of length $m$, ${\bf T}=({\bf
X,y})$)

0: the number of trees, $k$; the minimum interval length, p. (default
$k\leftarrow 500$, $p\leftarrow$min($16,m/2$))

1: Let ${\bf F}\leftarrow<{\bf F_{1}}\ldots{\bf F_{k}}>$ be the trees in the
forest.

2: $i\leftarrow 1$

3: while $i<k$ $\and$ timeRemaining() do

4: buildAdaptiveTimingModel()

5: if $i=1$ then

6: $r\leftarrow m$

7: else

8: $max\leftarrow$ findMaxIntervalLength()

9: $r\leftarrow$ findPowerOf2Interval($p$, $max$)

10: $b\leftarrow$ randBetween($1,m-r$)

11: $\bf{T^{\prime}}\leftarrow$ removeAttributesOutsideOfRange($\bf{T},b,r$)

12: $\bf{S}\leftarrow$ getSpectralFeatures($\bf{T^{\prime}}$)

13: ${\bf F_{i}}$.buildRandomTreeClassifier($\bf{S}$,$y$)

14: updateAdaptiveModel($r$)

15: $i\leftarrow i+1$

###  2.4 Bag of SFA Symbolic Features (BOSS)

Dictionary based classifiers convert real valued time series into a sequence
of discrete symbol words, then base classification on these words. Commonly, a
sliding window of length $w$ is run across a series. For each window, the real
valued series of length $w$ is converted through approximation and
discretisation processes into a symbolic string of length $l$ (referred to as
a word), which consists of $\alpha$ possible letters. The occurrence in a
series of each word from the dictionary defined by $l$ and $\alpha$ is counted
and, once the sliding window has completed, the series is transformed into a
histogram. Classification is based on the histograms of the words extracted
from the series, rather than the raw data. The Bag of Symbolic-Fourier-
Approximation Symbols (BOSS) [19] ensemble was found to be the most accurate
dictionary based classifier in the bake off [3]. Hence, it forms our benchmark
for new dictionary based approaches.

Algorithm 4 gives a formal description of the bag forming process of an
individual BOSS classifier. Windows may or may not be normalised (lines 6 and
7). Words are created using Symbolic Fourier Approximation (SFA) [20] (lines
8-13). SFA first finds the Fourier transform of the window (line 8), ignoring
the first term if normalisation occurs (lines 9-12). It then discretises the
first $l$ Fourier terms into $\alpha$ symbols to form a word in the method
SFAlookup, using a bespoke supervised discretisation algorithm called Multiple
Coefficient Binning (MCB) (line 13). Lines 14-16 implement the process of not
counting self similar words: if two consecutive windows produce the same word,
the second occurrence is ignored. This is to avoid a slow-changing pattern
relative to the window size being over-represented in the resulting histogram.

BOSS uses a non-symmetric distance function in conjunction with a nearest
neighbour classifier. Only the words contained in the test instance’s
histogram (i.e. the word count is above zero) are used in the distance
calculation, but it is otherwise the Euclidean distance.

Algorithm 4 baseBOSS(A list of $n$ time series of length $m$, ${\bf T}=({\bf
X,y})$)

0: the word length $l$, the alphabet size $\alpha$, the window length $w$,
normalisation parameter $z$

1: Let ${\bf H}$ be a list of $n$ histograms $({\bf h}_{1},\ldots,{\bf
h}_{n})$

2: Let ${\bf B}$ be a matrix of $l$ by $\alpha$ breakpoints found by MCB

3: for $i\leftarrow 1$ to $n$ do

4: for $j\leftarrow 1$ to $m-w+1$ do

5: ${\bf s}\leftarrow x_{i,j}\ldots x_{i,j+w-1}$

6: if $z$ then

7: $s\leftarrow$normalise($s$)

8: ${\bf q}\leftarrow$ DFT(${\bf s},l,\alpha$,$p$) { q is a vector of the
complex DFT coefficients}

9: if $z$ then

10: ${\bf q^{\prime}}\leftarrow(q_{2}\ldots q_{l/2+1})$

11: else

12: ${\bf q^{\prime}}\leftarrow(q_{1}\ldots q_{l/2})$

13: ${\bf r}\leftarrow$ SFAlookup(${\bf q^{\prime},B}$)

14: if ${\bf r}\neq{\bf p}$ then

15: $pos\leftarrow$index(${\bf r}$)

16: ${h}_{i,pos}\leftarrow{h}_{i,pos}+1$

17: ${\bf p}\leftarrow{\bf r}$

The final classifier is an ensemble of individual BOSS classifiers. The
original BOSS ensemble built all models over a pre-defined parameter space for
$w$, $l$, $z$ and $\alpha$ and retained all base classifiers with accuracy of
92% or higher of the best. This introduces instability in memory usage and
carries a time overhead. HIVE-COTE 1.0 uses contractable BOSS (cBOSS) [17] as
its dictionary based classifier. cBOSS changes the method used by BOSS to form
its ensemble to improve efficiency and allow for a number of usability
improvements. cBOSS was shown to be an order of magnitude faster than BOSS on
both small and large datasets from the UCR archive while showing no
significant difference in accuracy [17]. It randomly samples the parameter
space without replacement (line 8), subsamples the data for each base
classifier (line 10), and retains a fixed number of base classifiers. An
exponential weighting scheme based on train accuracy, such as the one used in
HIVE-COTE, is introduced for ensemble members.

Algorithm 5 cBOSS(A list of $n$ cases length $m$, ${\bf T}=({\bf X,y})$)

0: the maximum number of base classifiers, $k$; the number of parameter
samples, $s$; the proportion of cases to sample, $p$. (default $k\leftarrow
50$; $s\leftarrow 250$; $p\leftarrow 0.7$)

1: Let $w$ be window length, $l$ be word length, $z$ be normalise/not
normalise and $\alpha$ be alphabet size.

2: Let ${\bf B}\leftarrow<{\bf B_{1}}\ldots{\bf B_{k}}>$ be a list of $k$ BOSS
classifiers

3: Let ${\bf W}\leftarrow<{w}_{1},\ldots,{w}_{k}>$ be a list of classifier
weights

4: Let ${\bf R}$ be a set of possible BOSS parameter combinations

5: $i\leftarrow 0$

6: $minAcc\leftarrow\infty$

7: while $i<s\;\and\;|{\bf R}|>0\;\and\;$ timeRemaining() do

8: $[l,a,w,p]\leftarrow$ randomSampleParameters(${\bf R}$)

9: ${\bf R}={\bf R}\setminus\\{[l,a,w,p]\\}$

10: ${\bf T^{\prime}}\leftarrow$ subsampleData(${\bf T}$)

11: ${\bf cls}\leftarrow$ baseBOSS(${\bf T^{\prime}},l,a,w,z$)

12: $acc\leftarrow$ estimateAccuracy(${\bf T^{\prime}}$,$cls$) { estimate
accuracy on train data}

13: if $i<k$ then

14: ${\bf B_{i}}\leftarrow cls$, $w_{i}\leftarrow acc^{4}$

15: if $acc<minAcc$ then

16: $minAcc\leftarrow acc$, $idx\leftarrow i$

17: else if $acc>min\\_acc$ then

18: ${\bf B_{idx}}\leftarrow{\bf cls}$, $w_{idx}\leftarrow acc^{4}$

19: $[minAcc,idx]\leftarrow$ findLowestAcc(${\bf B}$)

20: $i\leftarrow i+1$

###  2.5 Shapelet Transform Classifier (STC)

There are two significant changes to the STC used in HC 1.0. Firstly, it only
fully enumerates the shapelet space when there is sufficient time to do so.
Secondly, it uses a Rotation Forest classifier [18] rather than a
heterogeneous ensemble. The shapelet transform is highly configurable: it can
use a range of sampling/search techniques in addition to alternative quality
measures. We present the default settings and direct the interested reader to
the code. The amount of time for the shapelet search is now a parameter. The
algorithm calculates how many possible shapelets there are in a data set, then
estimates how many shapelets it can sample from each series. After searching,
it updates its timing model using simple reinforcement learning. These
operations are encapsulated in operations estimateNumberOfShapelets (line 7)
and updateTimingModel (line 11). Shapelets are randomly sampled in the method
sampleShapelets (line 8). If the algorithm is allowed more shapelets than the
series contains, it evaluates them all. We have experimented with a range of
alternative neighbourhood search algorithms, but nothing is much better than
random search. Once the shapelets are generated, they are evaluated using
information gain (line 9). We use a one vs all evaluation for multi-class
problems [5]. Overlapping shapelets are removed in line 10, before the
candidates are merged into the overall pool, with the weakest members of the
population being deleted. Once the search is complete, the transform is
performed (line 13) and the classifier constructed (line 14).

Algorithm 6 STC(A list of $n$ cases length $m$, ${\bf T}=({\bf X,y})$)

0: the maximum number of shapelets to keep, $k$; the shapelet search time,
$t$. (default $k\leftarrow 1000$, $t\leftarrow 1$ hour.

1: Let ${\bf S}$ be a list of up to $k$ shapelets

2: Let ${\bf R}$ be a rotation forest classifier.

3: $i\leftarrow 0$

4: $minIG\leftarrow 0$

5:

6: while shapeletTimeRemaining($t$) do

7: $p\leftarrow$ estimateNumberOfShapelets($t,m,n$)

8: ${\bf S^{\prime}}\leftarrow$ sampleShapelets(${\bf x_{i}}$,$p$)

9: ${\bf s}\leftarrow$ evaluateShapelets($S^{\prime}$,${\bf T}$)

10: ${\bf S^{\prime}}\leftarrow$ removeSelfSimilar($S^{\prime}$, $s$)

11: updateTimingModel()

12: ${\bf S}\leftarrow$ merge($S$,$S^{\prime}$)

13: ${\bf X^{\prime}}\leftarrow$ shapeletTransform($X$,$S$)

14: ${\bf R}$.buildRotationForest(${\bf X^{\prime}},$y)

##  3 HIVE-COTE 1.0 Usability

We help maintain two toolkits that include time series classification
functionality. sktime111https://github.com/alan-turing-institute/sktime is an
open source, Python based, sklearn compatible toolkit for time series
analysis. sktime is designed to provide a unifying API for a range of time
series tasks such as annotation, prediction and forecasting (see [16] for a
description of the overarching design of sktime and [1] for an experimental
comparison of some of the classification algorithms available). The Java
toolkit for time series machine learning, tsml222https://github.com/uea-
machine-learning/tsml, is Weka compatible and is the descendent of the
codebase used to perform the bake off. The two toolkits will eventually
converge to include all the features described here. Experiments reported in
this paper are conducted with tsml, as it has more functionality.

###  3.1 Java Implementation of HIVE-COTE 1.0 in tsml

The HIVE_COTE class is in the package tsml.classifiers.hybrids and can be used
as any other Weka classifier. The default configuration is that described in
this paper. The code described here is all available in the class
EX07_HIVE_COTE_Examples with more detail and comments. A basic build is
described in Listing 1. It cannot handle missing values, unequal length series
or multivariate problems yet.

⬇

1 HIVE_COTE hc = new HIVE_COTE();

2 //this setup called in default constructor in April 2020

3 hc.setupHIVE_COTE_1_0();

4 Instances[] trainTest =

5 DatasetLoading.sampleItalyPowerDemand(0);

6 hc.buildClassifier(trainTest[0]);

7 int correct=0;

8 for(Instance ins: trainTest[1]){

9 double c=hc.classifyInstance(ins);

10 if(c==ins.classValue())

11 correct++;

Listing 1: A most basic way using HIVE-COTE 1.0 in tsml

We rarely build the classifier in this way. Instead, we build the components
and post process the meta ensemble. This is most easily done using our
Experiments class, which formats the output in a standard way. An example code
snippet is in Listing 2. Details on optional input flags not given below can
be found in the code.

⬇

1 String[] settings=new String[6];

2 //Where to get data

3 settings[0]="-dp=src/main/java/experiments/data/tsc/";

4 //Where to write results

5 settings[1]="-rp=Temp/";

6 //Whether to generate train files or not

7 settings[2]="-gtf=true";

8 //Classifier name: See ClassifierLists for valid options

9 settings[3]="-cn=TSF";

10 //Problem file

11 settings[4]="-dn=Chinatown";

12 //Resample number: 1 gives the default train/test split

13 settings[5]="-f=1";

14 Experiments.ExperimentalArguments expSettings =

15 new Experiments.ExperimentalArguments(settings);

16 Experiments.setupAndRunExperiment(expSettings);

Listing 2: Using Experiments.java to build a single component.

HIVE_COTE can read component results directly from file using syntax of the
form given in Listing 3. It will look in the directory structure created by
Experiments. In this example, when building TSF, it will look for
“C:/Temp/TSF/Predictions/Chinatown/trainFold0.csv”. Currently, this file
loading method requires all the classifier results to be present in order to
build.

⬇

1 HIVE_COTE hc=new HIVE_COTE();

2 hc.setBuildIndividualsFromResultsFiles(true);

3 hc.setResultsFileLocationParameters("C:/Temp/", "Chinatown", 0);

4 String[] components={"TSF","RISE","cBOSS","STC"};

5 hc.setClassifiersNamesForFileRead(components);

Listing 3: Buildind HIVE-COTE from existing results files

HIVE_COTE is configurable for different components, threadable (see Listing 4)
and contractable (see Listing 5). In sequential mode, it simply divides the
time equally between components. When threaded, it gives the full contract
time to each component. It does not yet thread individual components; it is on
our development list. You can set the maximum build time for HIVE_COTE if the
components all implement the TrainTimeContractable interface.

⬇

1 HIVE_COTE hc = new HIVE_COTE();

2 EnhancedAbstractClassifier[] c=new EnhancedAbstractClassifier[2];

3 c[0]=new RISE();

4 c[1]=new TSF();

5 String[] names={"RISE","TSF"};

6 hc.setClassifiers(c,names,null);

7 hc.enableMultiThreading(2);

Listing 4: Threaded build of HIVE-COTE with bespoke classifiers

⬇

1//Ways of setting the contract time

2 HIVE_COTE hc = new HIVE_COTE();

3//Minute, hour or day limit

4 hc.setMinuteLimit(1);

5 hc.setHourLimit(2);

6 hc.setDayLimit(1);

7//Specify units

8 hc.setTrainTimeLimit(30, TimeUnit.SECONDS);

9 hc.setTrainTimeLimit(1, TimeUnit.MINUTES);

10//Or just give it in nanoseconds

11 hc.setTrainTimeLimit(10000000000L);

12 hc.buildClassifier(train);

Listing 5: Contracting HIVE-COTE for a rom existing results files

Finally, HIVE_COTE is tuneable. Our method of implementing tuned classifiers
is to wrap the base classifier in a TunedClassifier object which interacts
through the method setOptions. An example of tuning the $\alpha$ parameter is
given in Figure 6. We advise tuning using results loaded from file. We have
found tuning $\alpha$ makes no significant difference. We have not finished
evaluating tuning which components to use.

⬇

1 HIVE_COTE hc=new HIVE_COTE();

2 hc.setBuildIndividualsFromResultsFiles(true);

3 hc.setResultsFileLocationParameters(resultsPath, dataset, fold);

4 hc.setClassifiersNamesForFileRead(cls);

5 TunedClassifier tuner=new TunedClassifier();

6 tuner.setClassifier(hc);

7 ParameterSpace pc=new ParameterSpace();

8 double[] alphaVals={1,2,3,4,5,6,7,8,9,10};

9 pc.addParameter("A",alphaVals);

10 tuner.setParameterSpace(pc);

Listing 6: Tuning HIVE-COTE $\alpha$ parameter from existing results files

###  3.2 Python Implementation of HIVE-COTE 1.0 in sktime

A Python implementation of HIVE-COTE is under development and available in
sktime. As previously discussed, this version of the algorithm is less mature
in terms of its development and the results reported in this paper are from
the Java version of the algorithm. The sktime implementation is in an alpha
state, and will eventually converge on the same functionally as the more
developed Java implementation, but currently has a number of limitations (such
as building from file and running constituents in parallel). Further, the
Python implementations of the constituent classifiers are less efficient than
the Java implementations, and as such, HIVE-COTE 1.0 in sktime is slower than
the tsml implementation on equivalent inputs.

The interface and basic usage of HIVE-COTE in sktime is very similar to that
of tsml. The terminology is slightly different however as the sktime version
uses fit and predict derived from scikit-learn while the Java version uses
build and classify from Weka/tsml. Notionally the process of constructing and
making predictions with HIVE-COTE are equivalent and a simple example of
fitting and predicting with HIVE-COTE in sktime is given in Listing 7. The
most up-to-date implementation of HIVE-COTE can be found in the sktime toolkit
on GitHub under the hive-cote branch 333https://github.com/alan-turing-
institute/sktime/blob/hive_cote/sktime/contrib/meta/ensembles.py.

⬇

1def basic_hive_cote(data_dir, dataset_name):

2 # using the default constructor for the HIVE-COTE class

3 hc = HIVE-COTE()

4

5 # loading training data

6 train_x, train_y = load_data(

7 data_dir + dataset_name + "_TRAIN.ts")

8

9 # building HIVE-COTE 1.0 sequentially

10 hc.fit(train_x, train_y)

11

12 # loading testing data

13 test_x, test_y = load_data(

14 data_dir + dataset_name + "_TEST.ts")

15

16 # predict class values of the test data

17 preds = hc.predict(train_x)

18

19 # calculate the test accuracy

20 acc = accuracy_score(self.train_y, preds)

Listing 7: A simple example of using HIVE-COTE 1.0 in sktime

##  4 Performance

To measure performance of the new HIVE-COTE, we evaluate each component and
the algorithm itself on 112 of the 128 UCR archive datasets. These 112
datasets are all equal length and have no missing values. Figure 2 shows the
critical difference diagram for HIVE-COTE 1.0 (HC 1.0) and its four
components. This broadly mirrors the performance presented in [15].

Figure 2: Critical difference diagram for HIVE-COTE 1.0 and its four
components on 112 UCR TSC problems. Full results are available from
www.timeseriesclassification/Results/HC1.0.xls

We have compared the results of the four components to the original results
and found there is no significant different. Comparison of HIVE-COTE alpha
version and 1.0 identify a small, but significant, difference. Removing EE
makes HIVE-COTE worse on 48 and better on 33 of the 85 datasets used in the
original experiments. The mean reduction in accuracy is 0.6%. The differences
in accuracy on specific problems identify those where a distance based
approach may be the best. MedicalImages, SonyAIBORobotSurface1, WordSynonyms
and Lightning7 were all more than 5% less accurate with EE removed. We are
trading this small loss in average accuracy for significant gains in run time
and reduction in memory usage. We have explored ways of making EE more
efficient, but as yet none of these approaches have met the criteria for
maintaining accuracy and providing sufficient speed up. Since HIVE-COTE alpha
was proposed, three new algorithms have achieved equivalent accuracy. TS-CHIEF
[21] is a tree ensemble that embeds dictionary, spectral and distance based
representations. InceptionTime [9] is a deep learning ensemble, combining five
homogeneous networks each with random weight initialisations for stability.
ROCKET [7] uses a large number (10,000 by default) of randomly parameterised
convolution kernels in conjunction with a linear ridge regression classifier.
We use the configurations of each classifier described in their respective
publications. Figure 3 shows the critical difference diagram for these three
classifiers and HC 1.0. There is no significant difference between any of
them.

Figure 3: Critical difference diagram for current state of the art on 109 UCR
TSC problems.

The differences between HIVE-COTE and the other three are summarised in Table
1. TS-CHIEF is the most similar to HIVE-COTE, with an average accuracy just
0.25% lower and a high correlation between accuracies. ROCKET has a high
variation in performance, as reflected in the high standard deviation of
differences and the difference between the mean and median of the differences.
It also has the lowest correlation to HC 1.0.

Table 1: Summary of differences in accuracy between HIVE-COTE and the other three algorithms. A negative difference indicates HIVE-COTE is more accurate. | Mean | Median | Std Dev of differences | Correlation
---|---|---|---|---
TS-CHIEF | -0.25% | 0.00% | 3.801 | 95.54%
InceptionTime | -0.65% | 0.00% | 5.82 | 89.46%
ROCKET | -1.41% | 0.05% | 8.64 | 80.28%

Table 2 shows the results for problems with the 10 biggest differences between
HIVE-COTE and ROCKET. ROCKET does very poorly on some problems (hence the
large average difference in Table 1). InceptionTime also shares this
characteristic of occasionally simply failing on a problem. For example,
InceptionTime is over 20% less accurate on the problem Rock than the other
three classifiers. It is not obvious why this happens to both ROCKET and
InceptionTime. It may be the result of over fitting.

Table 2: Accuracy for ten problems with the biggest difference between HIVE-COTE and ROCKET (five negative and five positive differences). | TS-CHIEF | ROCKET | HC 1.0 | InceptionTime
---|---|---|---|---
PigAirwayPressure | 96.01% | 19.55% | 95.77% | 92.21%
SemgHandMovementCh2 | 88.50% | 65.26% | 88.90% | 55.10%
EthanolLevel | 60.56% | 62.53% | 84.90% | 87.55%
CinCECGTorso | 95.34% | 86.41% | 99.37% | 83.28%
ScreenType | 59.42% | 60.90% | 72.42% | 70.56%
FiftyWords | 84.27% | 82.51% | 77.16% | 82.68%
ChlorineConcentration | 66.08% | 79.61% | 73.39% | 86.36%
MedicalImages | 79.91% | 80.51% | 74.04% | 79.63%
WordSynonyms | 79.37% | 76.44% | 69.32% | 75.18%
SonyAIBORobotSurface1 | 88.97% | 95.81% | 82.63% | 95.42%

It is worthwhile considering the run time of these algorithms. However,
comparisons are made more difficult because of the different software. HIVE-
COTE and TS-CHIEF are both built using tsml, so are directly comparable. Table
3 summarises the time taken to train the classifiers on 109 UCR problems.
Three problems (HandOutlines, NonInvasiveFetalECGThorax1 and
NonInvasiveFetalECGThorax2) are omitted because TS-CHIEF did not complete
within 7 days (the job limit on our cluster). Of the HC-1.0 components, STC is
by far the slowest. This is caused by the classifier, Rotation Forest, not the
transform, which is contracted to take at most one hour. The STC design
principle is to choose a large number of shapelets (up to 1000) and let the
classifier sort out their importance. Rotation forest is the best approach for
problems with all real valued attribute [2], and we have developed a
contracted version that can, if necessary, be used to control the build time
(see Listing 6). Most of the extra computation required by HIVE-COTE is in
forming the estimates of the accuracy on the train data. TS-CHIEF is the
slowest algorithm, but it is approximately the same as HIVE-COTE.

Table 3: Time in hours to train a classifier for 112 of the UCR problems | TSF | cBOSS | RISE | STC | HC 1.0 | TS-CHIEF
---|---|---|---|---|---|---
Mean | 0.13 | 0.11 | 0.15 | 2.30 | 4.26 | 5.75
Total | 14.15 | 12.31 | 16.05 | 251.12 | 464.49 | 626.33

Test time can be a factor for deploying classifiers in time critical
situations. Table 4 summarises the time (in minutes) taken to predict the test
cases. STC is the slowest component when testing, again caused by the
classifier not the transform. TS-CHIEF is the slowest in testing.

Table 4: Time in minutes to make predictions on the test data for 109 of the UCR problems | TSF | cBOSS | RISE | STC | HC 1.0 | TS-CHIEF
---|---|---|---|---|---|---
Mean | 0.01 | 0.08 | 0.07 | 0.62 | 0.78 | 4.09
Total | 1.42 | 8.63 | 7.09 | 67.50 | 84.64 | 445.55

We do not have reliable times for ROCKET, which is run in sktime. It is
undoubtedly faster than HIVE-COTE and TS-CHIEF though. We can also measure
maximum memory usage of the classifiers, as this is also often a serious
bottleneck for scalability. Table 5 summarises the memory usage of the six
tsml classifiers. The pattern is similar to that observed with run time. TSF,
RISE and cBOSS have a light memory footprint. STC, and hence HIVE-COTE, have a
larger memory usage. TS-CHIEF is the most memory hungry, and it seemingly does
not scale well in terms of memory. For, example, on the 11 problems where
HIVE-COTE required more than 1 GB, TS-CHIEF required approximately four times
the memory. The highest memory usage on the 109 problem it could complete was
18 GB on FordA and FordB. It seems highly likely that the memory usage of TS-
CHIEF could be improved without loss of accuracy.

Table 5: Memory usage in MB for 109 UCR problems | TSF | cBOSS | RISE | STC | HC 1.0 | TS-CHIEF
---|---|---|---|---|---|---
Mean | 162 | 255 | 263 | 1,464 | 1,618 | 2,004
Max | 1061 | 3432 | 740 | 3,533 | 4,426 | 18,532

##  5 Conclusions

The purpose of this report is to present a more practical version of HIVE-COTE
and compare it to recent advances in the field of TSC. On average, there is no
real difference between InceptionTime, TS-CHIEF, ROCKET and HIVE-COTE in terms
of accuracy. ROCKET is undoubtedly the fastest, but it is prone to fail badly
on the occasional data set. InceptionTime is slow and requires a GPU. It also
fails on some data sets. HIVE-COTE and TS-CHIEF are broadly comparable. HIVE-
COTE is presently more configurable and controllable.

We hope the results presented here and on the accompanying website serve to
act as a baseline comparison for any new research in the field in the future.
We have presented HIVE-COTE as is, without improvements that we have been
working on. Our next goal is to find a more accurate version that is
comparable in run time and memory usage, and to demonstrate scalability by
introducing some much larger problems into the UCR archive. We are also
developing versions that can handle unequal length series and multivariate TSC
problems.

## References

  * [1] A.Bagnall, F. Király, M. Löning, M. Middlehurst, and G. Oastler.  A tale of two toolkits, report the first: benchmarking time series classification algorithms for correctness and efficiency.  ArXiv e-prints, arXiv:1909.05738, 2019.
  * [2] A. Bagnall, A. Bostrom, G. Cawley, M. Flynn, J. Large, and J. Lines.  Is rotation forest the best classifier for problems with continuous features?  ArXiv e-prints, arXiv:1809.06705, 2018.
  * [3] A. Bagnall, J. Lines, A. Bostrom, J. Large, and E. Keogh.  The great time series classification bake off: a review and experimental evaluation of recent algorithmic advances.  Data Mining and Knowledge Discovery, 31(3):606–660, 2017.
  * [4] A. Bagnall, J. Lines, J. Hills, and A. Bostrom.  Time-series classification with COTE: The collective of transformation-based ensembles.  IEEE Transactions on Knowledge and Data Engineering, 27:2522–2535, 2015.
  * [5] A. Bostrom and A. Bagnall.  Binary shapelet transform for multiclass time series classification.  Transactions on Large-Scale Data and Knowledge Centered Systems, 32:24–46, 2017.
  * [6] H. Dau, A. Bagnall, K. Kamgar, M. Yeh, Y. Zhu, S. Gharghabi, C. Ratanamahatana, A. Chotirat, and E. Keogh.  The ucr time series archive.  IEEE/CAA Journal of Automatica Sinica, 6(6):1293–1305, 2019.
  * [7] Angus Dempster, François Petitjean, and Geoffrey I Webb.  Rocket: Exceptionally fast and accurate time series classification using random convolutional kernels.  arXiv preprint arXiv:1910.13051, 2019.
  * [8] H. Deng, G. Runger, E. Tuv, and M. Vladimir.  A time series forest for classification and feature extraction.  Information Sciences, 239:142–153, 2013.
  * [9] HI. Fawaz, B. Lucas, G. Forestier, C. Pelletier, D. Schmidt, J. Weber, G. Webb, L. Idoumghar, P. Muller, and F. Petitjean.  Inceptiontime: Finding alexnet for time series classification.  ArXiv, 2019.
  * [10] M. Flynn, J. Large, and A. Bagnall.  The contract random interval spectral ensemble (c-RISE): The effect of contracting a classifier on accuracy.  In International Conference on Hybrid Artificial Intelligence Systems, volume 11734 of Lecture Notes in Computer Science, pages 381–392. 2019.
  * [11] J. Hills, J. Lines, E. Baranauskas, J. Mapp, and A. Bagnall.  Classification of time series by shapelet transformation.  Data Mining and Knowledge Discovery, 28(4):851–881, 2014.
  * [12] J. Large, J. Lines, and A. Bagnall.  A probabilistic classifier ensemble weighting scheme based on cross validated accuracy estimates.  Data Mining and Knowledge Discovery, 33(6):1674––1709, 2019.
  * [13] J. Lines and A. Bagnall.  Time series classification with ensembles of elastic distance measures.  Data Mining and Knowledge Discovery, 29:565–592, 2015.
  * [14] J. Lines, S. Taylor, and A. Bagnall.  HIVE-COTE: The hierarchical vote collective of transformation-based ensembles for time series classification.  In Proc. 16th IEEE International Conference on Data Mining, 2016\.
  * [15] J. Lines, S. Taylor, and A. Bagnall.  Time series classification with HIVE-COTE: The hierarchical vote collective of transformation-based ensembles.  ACM Trans. Knowledge Discovery from Data, 12(5), 2018.
  * [16] M. Löning, A. Bagnall, S. Ganesh, V. Kazakov, J. Lines, and F. J. Király.  A unified interface for machine learning with time series.  ArXiv e-prints, arXiv:1909.07872, 2019.
  * [17] M. Middlehurst, W. Vickers, and A. Bagnall.  Scalable dictionary classifiers for time series classification.  In International Conference on Intelligent Data Engineering and Automated Learning, Lecture Notes in Computer Science, pages 11–19. 2019.
  * [18] J. Rodriguez, L. Kuncheva, and C. Alonso.  Rotation forest: A new classifier ensemble method.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 28(10):1619–1630, 2006.
  * [19] P. Schäfer.  The BOSS is concerned with time series classification in the presence of noise.  Data Mining and Knowledge Discovery, 29(6):1505–1530, 2015.
  * [20] P. Schäfer and M. Högqvist.  SFA: a symbolic Fourier approximation and index for similarity search in high dimensional datasets.  In Proceedings of the 15th International Conference on Extending Database Technology, pages 516–527, 2012.
  * [21] A. Shifaz, C. Pelletier, F. Petitjean, and G. Webb.  TS-CHIEF: A scalable and accurate forest algorithm for time series classification.  ArXiv e-prints, arXiv:1906.10329, 2019.