∎

22institutetext: HSE University
International Laboratory of Stochastic Analysis and its Applications
Pokrovsky boulevard 11, 109028 Moscow, Russia

22email<EMAIL_ADDRESS>22institutetext<EMAIL_ADDRESS>
# Modelling the Bitcoin prices and the media attention to Bitcoin via the
jump-type processes ††thanks:

Ekaterina Morozova and Vladimir Panov

(Received: )

###### Abstract

In this paper, we present a new bivariate model for the joint description of
the Bitcoin prices and the media attention to Bitcoin. Our model is based on
the class of the Lévy processes and is able to realistically reproduce the
jump-type dynamics of the considered time series. We focus on the low-
frequency setup, which is for the Lévy - based models essentially more
difficult than the high-frequency case. We design a semiparametric estimation
procedure for the statistical inference on the parameters and the Lévy
measures of the considered processes. We show that the dynamics of the market
attention can be effectively modelled by the Lévy processes with finite Lévy
measures, and propose a data-driven procedure for the description of the
Bitcoin prices.

###### Keywords:

Bitcoin, media attention, Lévy process, low-frequency data, deconvolution

Mathematics Subject Classiffication (2020): 60G51; 62M99; 91G30

##  1 Introduction

Stochastic modelling of the Bitcoin prices has been attracted the attention of
researchers for the last 10 years. The models, which are able to realistically
describe the dynamics of the cryptocurrencies, are quite useful for many
interesting research areas such as the prediction of bubbles, the development
of methods to raise individual capital (in particular, portfolio
diversification) and many others, see the survey by Härdle et al. (2020).

As usual in the price modelling, the hugest research stream is concentrated on
the GARCH-type models and related stochastic volatility models, see, e.g.,
Dyhrberg (2016), Katsiampa (2017), Naimy and Hayek (2018), Tiwari et al.
(2019). It is important to note that the dynamics of the Bitcoin price
essentially differs from the dynamics of many other financial instruments, due
to rapid changes (jumps) in the dynamics of returns and volatility. This
observation serves as the main motivation of the application of more
complicated approaches for modelling the Bitcoin prices.

One of the most promising ideas is to use various characteristics of the media
attention to the Bitcoin. In fact, as it was mentioned by many researches,
cryptocurrencies are prominently addressed by media, see, e.g., Glaser et al.
(2014) and many further references until the end of the article. In this area,
the media attention is quite often measured by the volume of Google searches
(in terms of Search Volume Index, SVI) and Wikipedia requests, see Da et al.
(2011), Kristoufek (2013), Urquahart (2018), Figa-Talamanca and Patacca
(2019).

The stochastic factor $M_{t},t\geq 0$, representing the attention index, can
be included in the model in different ways. For instance, in Cretarola et al.
(2020), the process $M_{t}$ is described by a separate stochastic differential
equation, which influences the price process $S_{t},t\geq 0,$ with some shift
in time. Namely, the SDE for $S_{t}$ includes $M_{t-\tau}$ with the delay
parameter $\tau,$ which can be estimated by the profile likelihood approach.
Nevertheless, this parameter can change over time - in particular, $\tau$ can
be both positive and negative (changes in media attention can influence the
prices of Bitcoin and vice versa).

In our paper, we introduce a new model for the media attention and the Bitcoin
prices. Our approach is based on the class of Lévy process, which is a well-
used tool for the construction of the asset-pricing models, see the brilliant
books by Schoutens (2003), Cont and Tankov (2004), and the recent papers by
Panov and Samarin (2019), Gardini et al. (2021).

In the statistical part of this research, we assume that the processes $M_{t}$
and $S_{t}$ are observed on the equidistant grid $t=0,\Delta,...,n\Delta$ with
some fixed $\Delta>0.$ For the Lévy - based models this case is essentially
more difficult than the case of high - frequency data (i.e., the case
$\Delta\to 0$ as $n\to\infty$), see, e.g., the papers by Belomestny and his
coauthors (2010, 2013, 2015, 2019), Neumann and Reiss (2009), Comte and Genon-
Catalot (2010), Gugushvili (2012). Nevertheless, this kind of data allows to
avoid the difficulties raised by the unclear behaviour of the delay parameter
(see above), since we consider the increments of prices over some fixed
intervals, and it is not important whether the changes in the Bitcoin prices
cause the changes in media attention, or vice versa.

In what follows, we consider the daily data on logarithmic returns of the
Bitcoin price and media attention collected over a five-year period, from 2017
to 2021. Our approach to modelling is based on two key observations. First,
let us note that the Bitcoin prices and the attention index to the Bitcoin
market (represented by the Google Trends to the word “Bitcoin”) jump on the
same intervals of the fixed length $\Delta=1$ corresponding to the daily data.
For example, Figure 1 represents the dynamics of logarithmic returns of the
media attention and the Bitcoin price in six consecutive months of 2021. It
can be seen that the spikes in Bitcoin price are well-aligned with those of
the media attention: at each time moment, there is typically a rapid change in
both values, with the spikes sometimes being headed in opposite directions.
This observation is further supported by nonparametric correlation
coefficients between the logarithmic returns of media attention and the
absolute values for the Bitcoin prices, which are presented in Table 1. As
will be shown later, the values of coefficients are slightly higher for the
observations that can be classified as jumps.

Table 1: Correlation coefficients for the logarithmic returns of media attention and absolute values of logarithmic returns of the Bitcoin prices. | Pearson | Kendall’s $\tau$ | Spearman’s $\rho$
---|---|---|---
Coefficient | 0.35357 | 0.173711 | 0.25628
p-value | 1.0027$\cdot 10^{-54}$ | 1.595$\cdot 10^{-28}$ | 1.0954$\cdot 10^{-28}$

Another key observation is that the class of Lévy processes with finite Lévy
measures can be efficiently used for the description of the media attention,
but is no longer appropriate for modelling the Bitcoin prices. Hence, in order
to describe the latter, one has to construct another model, taking into
account the relationship between the Bitcoin prices and media attention.

These two observations yield the following procedure:

  1. 1.

to build the model for the media attention;

  2. 2.

to determine the intervals with jumps in the media attention (which coincide
with the intervals containing jumps in the Bitcoin prices);

  3. 3.

to use this information for the construction of the model for the Bitcoin
prices.

The paper is organised as follows. In the next section, we focus on the Step 1
and present the jump model for the media attention index. Subsection 2.2 is
devoted to the presentation of the estimation procedure, which is summarised
on pages 2.2 and 1. The efficiency of the described procedure is illustrated
by numerical examples both on the simulated and real data (subsections 2.3 and
2.4 resp.). Later, in subsection 2.5, we show that this kind of models is
inappropriate for describing the Bitcoin prices. Section 3 deals with the
second step of the aforementioned plan and shows the possibility to determine
the jump times. Finally, the model for the Bitcoin prices is constructed in
Section 4.

Figure 1: The dynamics of logarithmic returns of media attention and Bitcoin
prices in June-December 2021.

##  2 Stochastic model for media attention

###  2.1 Set-up

For modelling the attention index to the Bitcoin market (denoted below by
$M_{t},t\geq 0$), we use the class of Lévy processes with the jump component
represented by a compound Poisson process. More precisely,

$\displaystyle M_{t}=\exp\Bigl{\\{}\mu t+\sigma
W_{t}+\sum_{i=1}^{N_{t}}\xi_{i}\Bigr{\\}},\qquad t\geq 0,$ (1)

where $\mu\in\mathbb{R},\sigma\in\mathbb{R}_{+},$ $N_{t}$ is a (homogeneous)
Poisson process with intensity $\lambda$ and $\xi_{1},\xi_{2},...$ are i.i.d.
random variables with absolutely continuous distribution having the density
$p:\mathbb{R}\to\mathbb{R}_{+}.$ This class of processes includes the Merton
model (corresponding to the case when $\xi_{1}$, $\xi_{2}$,.. have normal
distribution) and the Kou model (double-exponential distribution).

It is known that any Lévy process $L_{t},t\geq 0,$ can be represented in the
form $L_{t}=\mu t+\sigma W_{t}+J_{t}$, where $J_{t}$ is a pure-jump process,
see Sato, K. (1999). As it is shown in many papers (see the references in the
introduction), the statistical inference for this general case is rather
complicated, provided that the Lévy measure of $J_{t}$ may be infinite and
only the low-frequency data are available. Indeed, in this case it is
difficult to distinguish the small jumps of $J_{t}$ from the continuous part
$\mu t+\sigma W_{t}$, since the changes over time interval of the fixed length
$\Delta$ can be caused both by a jump or by the continuous part. Taking into
account this fact, we concentrate on the Lévy processes with finite Lévy
measure, that is, on the processes in the form (1).

###  2.2 Statistical inference

Assuming that the the process $M_{t}$ is observed on the discrete grid
$\\{0,\Delta,2\Delta,...\\}$ with fixed $\Delta>0$ (low-frequency data), we
aim to estimate the parameters $(\mu,\sigma,\lambda)$ and the density function
$p(\cdot)$. For the statistical inference we use the method described in
Belomestny, D., and Reiss, M. (2015). Note that

$D_{k}:=\log M_{k\Delta}-\log M_{(k-1)\Delta},\qquad k=1,2,...$

with $\log M_{0}=0,$ form a sequence of i.i.d. random variables, and due to
the Lévy-Khinchine formula, the characteristic function of their distribution
is equal to

$\displaystyle\phi_{\Delta}(u)={\mathbb{E}}[e^{\mathrm{i}u\log M_{\Delta}}]$
$\displaystyle=$ $\displaystyle\exp\Bigl{\\{}\Delta\bigl{(}\mathrm{i}\mu
u-\frac{1}{2}\sigma^{2}u^{2}+\int_{\mathbb{R}}\bigl{(}e^{\mathrm{i}ux}-1\bigr{)}\nu(dx)\bigr{)}\Bigr{\\}}$
(2) $\displaystyle=$ $\displaystyle\exp\Bigl{\\{}\Delta\bigl{(}\mathrm{i}\mu
u-\frac{1}{2}\sigma^{2}u^{2}-\lambda+\lambda\mathcal{F}[p](u)\bigr{)}\Bigr{\\}},\qquad
u\in\mathbb{R},\qquad$

where $\nu(B)=\int_{B}p(u)du,B\in\mathcal{B}(\mathbb{R}),$ is the Lévy measure
of $M_{t}$ and $\mathcal{F}[p](u)=\int_{\mathbb{R}}e^{\mathrm{i}ux}p(x)dx$ is
the Fourier transform of $p(\cdot).$ In what follows, we use the function

$\varphi_{\Delta}(u):=\frac{1}{\Delta}\log(\phi_{\Delta}(u)),\qquad
u\in\mathbb{R},$

which can be naturally estimated by

$\displaystyle\widehat{\varphi}_{\Delta}(u)$ $\displaystyle:=$
$\displaystyle\frac{1}{\Delta}\log\Bigl{(}\frac{1}{n}\sum_{k=1}^{n}e^{\mathrm{i}uD_{k}}\Bigr{)}$
$\displaystyle=$ $\displaystyle-\frac{1}{\Delta}\log
n+\frac{1}{2\Delta}\log\Bigl{[}\Bigl{(}\sum_{k=1}^{n}\cos(uD_{k})\Bigr{)}^{2}+\Bigl{(}\sum_{k=1}^{n}\sin(uD_{k})\Bigr{)}^{2}\Bigr{]}$
$\displaystyle\hskip
113.81102pt+\mathrm{i}\frac{1}{\Delta}\arctan\left(\frac{\sum_{k=1}^{n}\sin(uD_{k})}{\sum_{k=1}^{n}\cos(uD_{k})}\right),\qquad
u\in\mathbb{R}.$

Due to the Riemann-Lebesque lemma, $\mathcal{F}[p](u)\to 0$ as $|u|\to\infty$.
Therefore, from (2) we get

$\displaystyle\operatorname{Re}\bigl{(}\varphi_{\Delta}(u)\bigr{)}$
$\displaystyle=$ $\displaystyle-\frac{1}{2}\sigma^{2}u^{2}-\lambda+o(1),$
$\displaystyle\operatorname{Im}\bigl{(}\varphi_{\Delta}(u)\bigr{)}$
$\displaystyle=$ $\displaystyle\mu u+o(1),$

as $|u|\to\infty.$ The first relation yield that for large $u$ the real part
of $\varphi_{\Delta}(u)$ is linear in $u^{2},$ and the corresponding
coefficients are proportional to $\sigma^{2}$ and $\lambda.$ This observation
leads to the following definition of the estimates for $\sigma^{2}$ and
$\lambda$:

$\displaystyle\bigl{(}\hat{\sigma}^{2},\hat{\lambda}\bigr{)}$ $\displaystyle=$
$\displaystyle\operatornamewithlimits{arg\,min}_{\sigma^{2},\lambda}\int_{\mathbb{R}_{+}}w^{U_{n}}(u)\Bigl{[}\operatorname{Re}\bigl{(}\widehat{\varphi}_{\Delta}(u)\bigr{)}+\frac{1}{2}\sigma^{2}u^{2}+\lambda\Bigr{]}^{2}du$
(3) $\displaystyle=$
$\displaystyle\operatornamewithlimits{arg\,min}_{\sigma^{2},\lambda}\int_{\varepsilon}^{1}w(u)\Bigl{[}\operatorname{Re}\bigl{(}\widehat{\varphi}_{\Delta}(uU_{n})\bigr{)}+\frac{1}{2}\sigma^{2}u^{2}U_{n}^{2}+\lambda\Bigr{]}^{2}du,$

where $w^{U_{n}}(u)=U_{n}^{-1}w(u/U_{n})$ with some unbounded increasing
sequence of positive numbers $U_{n}$ and a function
$w:[\varepsilon,1]\to\mathbb{R}$, $\varepsilon>0.$

Similarly, for the imaginary part of $\varphi_{\Delta}(u)$, we also deal with
a kind of asymptotic linear regression model. Define the estimate of $\mu$ as

$\displaystyle\hat{\mu}$ $\displaystyle=$
$\displaystyle\operatornamewithlimits{arg\,min}_{\mu}\int_{\mathbb{R}_{+}}w^{V_{n}}(u)\Bigl{[}\operatorname{Im}\bigl{(}\widehat{\varphi}_{\Delta}(u)\bigr{)}-\mu
u\Bigr{]}^{2}du$ (4) $\displaystyle=$
$\displaystyle\operatornamewithlimits{arg\,min}_{\mu}\int_{\varepsilon}^{1}w(u)\Bigl{[}\operatorname{Im}\bigl{(}\widehat{\varphi}_{\Delta}(uV_{n})\bigr{)}-\mu
uV_{n}\Bigr{]}^{2}du$

with some unbounded increasing sequence of positive numbers $V_{n}$, which may
coincide with $U_{n}$.

Let us summarise the ideas presented above into the following algorithm.

[couleur=blue!15, logo=] Algorithm 1: Estimation of $(\mu,\sigma^{2},\lambda)$

Data: $n$ observations of the process $M_{t}$ at the equidistant grid
$t=\Delta,2\Delta,...,n\Delta.$

Initiate: Fix some positive numbers $U_{n},V_{n}$ and $\varepsilon\in(0,1)$.

Fix a function $w(\cdot)\geq 0$ supported on $[\varepsilon,1]$,

e.g., $w(x)={\mathbb{I}}\\{x\in[\varepsilon,1]\\}$.

Fix some $N\in\mathbb{N}$ and divide $[\varepsilon,1]$ into $N$ equidistant
intervals $I_{j}=[u_{j-1},u_{j}],j=1..N$ with $u_{0}=\varepsilon,u_{N}=1.$

Let $\tilde{u}_{j}\in I_{j},j=1..N$.

Algorithm:

  1. 1.

Estimate the function
$\varphi_{\Delta}(u):=\frac{1}{\Delta}\log({\mathbb{E}}[e^{\mathrm{i}uM_{\Delta}}])$

at the points $u=\tilde{u}_{j}U_{n},j=1..N,$ and
$u=\tilde{u}_{j}V_{n},j=1..N,$ by
$\displaystyle\hat{\varphi}_{\Delta}(u)=\frac{1}{\Delta}\log\Bigl{(}\frac{1}{n}\sum_{k=1}^{n}e^{\mathrm{i}uD_{k}}\Bigr{)},$
where $D_{k}:=M_{k\Delta}-M_{(k-1)\Delta},\;k=1..n,$ and $M_{0}=0.$

* 2.

Compute

$\displaystyle\Lambda_{d}$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{N}w(\tilde{u}_{j})\tilde{u}_{j}^{2d},\qquad
d=0,1,2,$ $\displaystyle\Psi_{d}$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{N}w(\tilde{u}_{j})\operatorname{Re}\bigl{(}\widehat{\varphi}_{\Delta}(\widetilde{u}_{j}U_{n})\bigr{)}(\tilde{u}_{j}U_{n})^{2d},\qquad
d=0,1.$ $\displaystyle\Upsilon$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{N}w(\tilde{u}_{j})\operatorname{Im}\bigl{(}\widehat{\varphi}_{\Delta}(\widetilde{u}_{j}V_{n})\bigr{)}\tilde{u}_{j}V_{n}.$

* 3.

Estimate $\sigma^{2}$ and $\lambda$ by

$\displaystyle\bigl{(}\sigma_{n}^{2},\lambda_{n}\bigr{)}$ $\displaystyle:=$
$\displaystyle\operatornamewithlimits{arg\,min}_{\sigma^{2},\lambda}\sum_{j=1}^{N}w(\widetilde{u}_{j})\Bigl{[}\operatorname{Re}\bigl{(}\widehat{\varphi}_{\Delta}(\widetilde{u}_{j}U_{n})\bigr{)}+\frac{1}{2}\sigma^{2}\widetilde{u}_{j}^{2}U_{n}^{2}+\lambda\Bigr{]}^{2}$
$\displaystyle=$
$\displaystyle\Bigl{(}2\frac{\Psi_{0}\Lambda_{1}U_{n}^{2}-\Psi_{1}\Lambda_{0}}{(\Lambda_{2}\Lambda_{0}-\Lambda_{1}^{2})U_{n}^{4}},\frac{\Psi_{1}\Lambda_{1}-\Psi_{0}\Lambda_{2}U_{n}^{2}}{(\Lambda_{2}\Lambda_{0}-\Lambda_{1}^{2})U_{n}^{2}}\Bigr{)}$

* 4.

Estimate $\mu$ by

$\displaystyle\hat{\mu}_{n}$ $\displaystyle:=$
$\displaystyle\operatornamewithlimits{arg\,min}_{\mu}\sum_{j=1}^{N}w(\widetilde{u}_{j})\Bigl{[}\operatorname{Im}\bigl{(}\widehat{\varphi}_{\Delta}(\widetilde{u}_{j}V_{n})\bigr{)}-\mu\widetilde{u}_{j}V_{n}\Bigr{]}^{2}$
$\displaystyle=$ $\displaystyle\Upsilon/(\Lambda_{1}V_{n}^{2}).$

###### Remark 1

In practice, the choice of $\varepsilon$ and $V_{n},U_{n}$ is based on the
idea that for any $u\in[\varepsilon U_{n},U_{n}]$ and any $u\in[\varepsilon
V_{n},V_{n}]$, the function $\varphi_{\Delta}(u)$ should be well approximated
by the following sum

$\displaystyle\varphi_{\Delta}(u)\approx\mathrm{i}\mu
u-\frac{1}{2}\sigma^{2}u^{2}-\lambda,$

and at the same time by its natural empirical estimate

$\displaystyle\varphi_{\Delta}(u)\approx\widehat{\varphi}_{\Delta}(u):=\frac{1}{\Delta}\log\Bigl{(}\frac{1}{n}\sum_{k=1}^{n}e^{\mathrm{i}uD_{k}}\Bigr{)}.$

Finally, we turn towards nonparametric estimation of the function $p.$ From
(2), we get that the natural estimator of $p$ is

$\displaystyle\hat{p}(x)$ $\displaystyle=$
$\displaystyle\widehat{F}^{-1}\Bigl{(}\bigl{(}\frac{1}{\hat{\lambda}}\bigl{(}\widehat{\varphi}_{\Delta}(u)-\mathrm{i}\hat{\mu}u+\frac{1}{2}\hat{\sigma}^{2}u^{2}\bigr{)}+1\bigr{)}K\bigl{(}u/T_{n}\bigr{)}\Bigr{)}$
(5) $\displaystyle=$
$\displaystyle\frac{1}{2\pi\hat{\lambda}}\int_{\mathbb{R}}e^{-\mathrm{i}ux}\bigl{(}\widehat{\varphi}_{\Delta}(u)-\mathrm{i}\hat{\mu}u+\frac{1}{2}\hat{\sigma}^{2}u^{2}+\hat{\lambda}\bigr{)}K\left(u/T_{n}\right)du,$

where $K:\mathbb{R}\to\mathbb{R}_{+}$ is a smoothing kernel defined as

$\displaystyle K(x)=\begin{cases}1,&|x|\leq 0.05,\\\
\exp\left\\{-\frac{e^{-1/(|x|-0.05)}}{1-|x|}\right\\},&0.05<|x|<1,\\\
0,&|x|\geq 1,\end{cases}$ (6)

and $T_{n}$ is an unbounded increasing sequence of positive numbers.

The particular choice of the kernel in the form (6) is not important and any
function from a broad class can be chosen. Nevertheless, as it is shown in
Belomestny, D. (2011), this kind of kernels (known as flat-top kernels) is
quite useful in semiparametric inference, since it often guarantees the
optimality of the estimators in the minimax sense. Our empirical studies
confirm the efficiency of this choice. [couleur=blue!15, logo=] Algorithm 2:
Estimation of the Lévy density $p:\mathbb{R}\to\mathbb{R}_{+}$

Data: $m$ points $x_{1},...,x_{m}$ (can be chosen by the researcher)

Initiate: Fix a positive numbers $T_{n}$.

Divide $[-1,1]$ into $N$ equidistant intervals

$I_{j}=[u_{j-1},u_{j}],j=1..N$ of length $\delta=2/N$

with $u_{0}=-1,u_{N}=1.$

Let $\tilde{u}_{j}\in I_{j},j=1..N$.

Algorithm:

  1. 1.

Estimate the function
$\varphi_{\Delta}(u):=\frac{1}{\Delta}\log({\mathbb{E}}[e^{\mathrm{i}uM_{\Delta}}])$

at the points $u=\breve{u}_{j}:=u_{j}T_{n},j=1..N,$ by
$\displaystyle\hat{\varphi}_{\Delta}(u)=\frac{1}{\Delta}\log\Bigl{(}\frac{1}{n}\sum_{k=1}^{n}e^{\mathrm{i}uD_{k}}\Bigr{)},$
where $D_{k}:=M_{k\Delta}-M_{(k-1)\Delta},\;k=1..n,$ and $M_{0}=0.$

* 2.

Estimate the density function $p(x)$ at the points $x_{1},..,x_{m}$ by

$\displaystyle\hat{p}_{n}(x_{s})$ $\displaystyle:=$
$\displaystyle\frac{T_{n}\delta}{2\pi\hat{\lambda}_{n}}\sum_{j=1}^{N}e^{-\mathrm{i}\breve{u}_{j}x_{s}}\bigl{(}\widehat{\varphi}_{\Delta}(\breve{u}_{j})-\mathrm{i}\hat{\mu}_{n}\breve{u}_{j}\bigr{.}$
$\displaystyle\bigl{.}\hskip
113.81102pt+\frac{1}{2}\hat{\sigma}_{n}^{2}\breve{u}_{j}^{2}+\hat{\lambda}_{n}\bigr{)}K\left(\breve{u}_{j}\right).$

###  2.3 Simulation study

The aim of the current section is to analyse the behaviour of the proposed
estimators on the simulated data. In what follows, we consider the Merton
jump-diffusion model, which is defined as (1) with normally distributed jump
sizes $\xi$. For illustrative purposes, we consider the same specification of
the model as in Belomestny, D., and Reiss, M. (2015); namely, the parameters
are chosen to be $\mu=0$, $\sigma=1$, $\lambda=10$ and $\Delta=0.1$.

Figure 2 represents the estimates
$\operatorname{Re}(\widehat{\varphi}_{\Delta}(u))$ based on 25 independent
samples of size $n=1000$ superimposed with the true values of
$\operatorname{Re}(\varphi_{\Delta}(u))$ and its approximation due to the
Riemann-Lebesgue lemma. It can be seen that the three lines are fairly close
in the interval $u\in[3,6]$. Hence, as suggested by Remark 1, the parameters
$\varepsilon$ and $U_{n},V_{n}$ can be chosen as $\varepsilon=0.5$ and
$U_{n}=V_{n}=6$.

Figure 2: The real part of $\varphi_{\Delta}(u)$ for the Merton jump-diffusion
model, its approximation and the estimates based on 25 repetitions.

The boxplots for the estimates of the parameters of the model for 25 samples
of size $n=1000,5000$ and 10000 are presented in Figure 3. It can be seen that
the estimates get closer to the true value with the growth of $n$.

Figure 3: Boxplots for the estimates of $\sigma^{2}$ (left), $\lambda$
(middle) and $\mu$ (right) for different sample sizes based on 25 repetitions.

Finally, Figure 4 depicts the true density $p(u)=(2\pi)^{-1/2}e^{-u^{2}/2}$ of
$\xi$ superimposed with its estimates (2) based on 25 independent samples. The
last plot represents the boxplots for the values of the mean-squared errors
computed as

$MSE(\hat{p}_{n}(x))=\frac{1}{m}\sum\limits_{s=1}^{m}\left(\hat{p}_{n}(x_{s})-p(x_{s})\right)^{2},$

where the points $x_{m}$ are chosen on an equidistant grid from $-5$ to 5,
with $m=1000$. Again, it can be observed that the estimate becomes more
accurate with the growth of $n$.

Figure 4: The true density $p(u)=\frac{1}{\sqrt{2\pi}}e^{-u^{2}/2}$ of $\xi$
superimposed with 25 density estimates (2) and the corresponding mean-squared
errors.

###  2.4 Real data

The current study considers the daily data over 5 most recent years, from
January 1st, 2017 to December 31st, 2021, the total of 1820 observations. The
media attention is represented by the Google
Trends111URL:https://trends.google.com/trends, which measure the popularity of
a topic relative to the maximum search interest in a given time frame. Namely,
the value of 100 corresponds to the most popular topic, while the value of 10
means that it is 10 times less popular. Since the chi-squared test for the
logarithmic returns of this characteristic is resulted into the p-value of
0.2056, the observations within the sample can be assumed to be independent,
justifying the use of the models with independent increments and the proposed
estimation procedure.

To begin with, we fit two simple partial cases of the model (1) — the Brownian
motion with a drift ($\xi\equiv 0$), which is purely continuous, and the pure-
jump Cauchy model ($\mu=\sigma=0$ and $\xi$ having the Cauchy distribution).
In both cases the parameters are estimated using the maximum likelihood
approach. The Kolmogorov-Smirnov test results into the p-values of $3.997\cdot
10^{-15}$ and $8.602\cdot 10^{-6}$ for two models, respectively. As can be
seen from Figure 9, presented in the appendix, the theoretical distributions
indeed do not provide an efficient fit to the true density of logarithmic
returns of media attention, since the scale parameter is evidently
underestimated in the normal model and overestimated in the Cauchy case. This
observation suggests the use of the model (1) in the general setting, which
includes both the continuous and the jump components.

Now we apply Algorithm 1 for fitting the parameters of the model (1). The
values of parameters $U_{n}$, $V_{n}$ and $\varepsilon$ were chosen as those
minimising the average quadratic error of the true and simulated densities
based on 25 samples of size 1000. The resulting parameter estimates are
$\hat{\sigma}_{n}=0.099$, $\hat{\lambda}_{n}=0.201$ and
$\hat{\mu}_{n}=-0.022$. Figure 5 represents the true density superimposed with
densities of simulated data for different sample sizes (500, 1000 and 5000)
and the optimal values of the parameters equal to $U_{n}=17$, $V_{n}=15$ and
$\varepsilon=0.1$. It can be seen that the density curves are very close, and
the quadratic errors become smaller as the sample size grows. In addition, the
Wilcoxon test does not reject the null hypothesis of zero median between the
true and simulated samples with an average p-value of 0.26141 obtained based
on 25 simulation runs. Hence, we can conclude that the proposed model indeed
provides a reasonable fit to the considered data.

Figure 5: The true density of logarithmic returns of media attention (red)
superimposed with densities of 25 simulated samples of different sizes, and
the corresponding average quadratic errors.

###  2.5 Modelling the bitcoin prices: first attempts

At the end of this section, let us mention that formally the model (1) can
also be applied for the description of the Bitcoin prices. Nevertheless, below
we provide some empirical evidence that the models of this kind are not
appropriate for the considered data.

As before, we consider the daily data on the same time period, which was
obtained from Finam222URL: https://finam.ru. Since the chi-squared test
resulted into the p-value of 0.1461, the logarithmic returns of the Bitcoin
prices were assumed to be independent.

On the first step, the normal and the Cauchy models were fitted to the Bitcoin
data. As before, Figure 10 (see Appendix) suggests that the first two models
do not provide a reasonable fit, with scale parameter being under- and
overestimated in the normal and Cauchy models, respectively. The Kolmogorov-
Smirnov test also rejects the null of the normal and Cauchy distributions,
with p-values of $8.815\cdot 10^{-14}$ and 0.00341. However, unlike the case
of media attention, it turns out that the Lévy process defined as in (1) is
also inappropriate for the Bitcoin price description. Figure 11 depicts the
true and simulated densities for values of $U_{n}=V_{n}=44$ and
$\varepsilon=0.5$ chosen as those minimising the discrepancy between the true
and simulated densities. It can be observed that even for these parameter
values the curves do not seem coherent, and the average quadratic errors both
take unreasonably large values and does not decrease with the growth of sample
size. Thus, there is an evidence that another model is needed for the
description of the Bitcoin prices.

##  3 Determination of the jump times

As we have already mentioned in the Introduction (see Figure 1), the jumps in
media attention and the jumps in the Bitcoin prices occur on the time
intervals. Let us use the constructed model for the media attention for the
determination of the jump times.

Assume that on each time interval of length $\Delta$, only one jump in the
media attention may occur. This assumption is partially confirmed to our
numerical studies: in fact, since $\hat{\lambda}\approx 0.2$ and $\Delta=1,$
one would expect approximately 1 jump in 5 intervals. Therefore,

$\displaystyle\log M_{k\Delta}-\log M_{(k-1)\Delta}\stackrel{{\scriptstyle
Law}}{{=}}\mu\Delta+\sigma\sqrt{\Delta}\zeta+\xi_{k}J_{k},\qquad k=1,...,K$

where $\zeta$ is a standard normal random variable and $J_{k}$ is an indicator
of the jump at time $k=1,2,...$ From now on, we will think of $J_{k}$ as of a
non-random value, which is determined by external factors (e.g., external
shocks in the market).

At each time moment $k=1,2,...$ one should decide between the following 2
possibilities:

$\displaystyle\log M_{k\Delta}-\log M_{(k-1)\Delta}$
$\displaystyle\stackrel{{\scriptstyle Law}}{{=}}$
$\displaystyle\widehat{\mu}_{n}\Delta+\widehat{\sigma}_{n}\sqrt{\Delta}\zeta$
(7) $\displaystyle\log M_{k\Delta}-\log M_{(k-1)\Delta}$
$\displaystyle\stackrel{{\scriptstyle Law}}{{=}}$
$\displaystyle\widehat{\mu}_{n}\Delta+\widehat{\sigma}_{n}\sqrt{\Delta}\zeta+\xi,$
(8)

where the estimates $\widehat{\mu}_{n}$ and $\widehat{\sigma}_{n}$ are defined
in Section 2, $\xi$ is a random variable with a distribution determined by the
density function $\widehat{p}_{n},$ $\xi$ and $\zeta$ are independent. The
densities of the r.v. in the right-hand sides of (7) and (8) are depicted in
Figure 6.

This figure suggests the following data-driven decision rule. All values of
media attention between $x_{1}\approx-0.2175$ and $x_{2}\approx 0.1715$ are
assumed to belong to the first distribution, and for these values we set
$J_{k}=0$. To the contrary, observations that are smaller than $x_{1}$ or
greater than $x_{2}$, which correspond to the atypically large or small values
for the model (7), are assumed to come from the distribution with the jump
component, and the value of $J_{k}$ is set equal to 1. Hence, the set

$\mathcal{K}=\Bigl{\\{}k=1..K:\qquad J_{k}=1\Bigr{\\}},$

corresponds to the jump times, while its complement
$\mathcal{K}^{c}=\\{1..K\\}\backslash\mathcal{K}$ represents the time moments
without jumps.

Figure 7 represents the jump times from June till December 2021 determined by
the procedure described above. It can be observed that the estimated jump
times indeed correspond to the significant changes both in media attention and
the Bitcoin prices. The total number of jump times constituted 337, which is
coherent with the results of Section 2, demonstrating that the jumps should
occur approximately once in four days. The values of correlation coefficients
computed for the observations classified as jumps are presented in Table 2. It
can be noticed that the values of these coefficients, all being significant,
are greater than those presented in Table 1 and computed for all observations
in total. Hence, there is indeed a significant dependence between the jumps of
two processes.

Figure 6: The densities of $\widehat{\mu}\Delta+\widehat{\sigma}\sqrt{\Delta}\zeta$ (black) and of $\widehat{\mu}\Delta+\widehat{\sigma}\sqrt{\Delta}\zeta+\xi$ (red)  Figure 7: The jumps in logarithmic returns of media attention and Bitcoin prices in June-December 2021  Table 2: Correlation coefficients for the logarithmic returns of media attention and absolute values of logarithmic returns of the Bitcoin price, computed only for the time moments classified as jumps. | Pearson | Kendall’s $\tau$ | Spearman’s $\rho$
---|---|---|---
Coefficient | 0.44357 | 0.3026 | 0.43869
p-value | 1.1204$\cdot 10^{-17}$ | 1.2176$\cdot 10^{-16}$ | 2.77314$\cdot 10^{-17}$

##  4 Modelling the Bitcoin prices

In view of the results of the previous two sections, we consider the following
model for the Bitcoin prices $S_{k\Delta},\;k=1..K,$

$\displaystyle\log\bigl{(}S_{k\Delta}\bigr{)}-\log\bigl{(}S_{(k-1)\Delta}\bigr{)}=\widetilde{\mu}\Delta+\widetilde{\sigma}\bigl{(}\widetilde{W}_{k\Delta}-\widetilde{W}_{(k-1)\Delta}\bigr{)}+\eta_{k}{\mathbb{I}}\bigl{\\{}k\in\mathcal{K}\bigr{\\}},\,$

where $\widetilde{\mu}\in\mathbb{R},\widetilde{\sigma}\in\mathbb{R}_{+}$ are
two constant parameters, $\widetilde{W}$ is a Brownian motion, and
$\eta_{1},\eta_{2},...$ are i.i.d. random variables independent on
$\widetilde{W}$. The set $\mathcal{K}$ is determined above in Section 3. So,
the jumps in the Bitcoin price occur at the same times as those of the media
attention.

###  4.1 Estimation of the continuous part

Due to our procedure, the sample is divided into two parts, corresponding to
the indices from the set $\mathcal{K}$ and to its compliment
$\mathcal{K}^{c}$. In what follows, we will use the second part of the sample
for the estimation of the parameters $\widetilde{\mu}$ and
$\widetilde{\sigma}$, and the first part — for the estimation of the density
of $\eta_{1}.$

Since the sample

$\log\bigl{(}S_{k\Delta}\bigr{)}-\log\bigl{(}S_{(k-1)\Delta}\bigr{)},\qquad
k\in\mathcal{K}^{c},$

is normally distributed with the mean $\widetilde{\mu}\Delta$ and the variance
$\widetilde{\sigma}^{2}\Delta,$ we estimate the parameters $\widetilde{\mu}$
and $\widetilde{\sigma}$ by the maximum likelihood approach, which leads to
$\widehat{\widetilde{\sigma}}=0.037$ and $\widehat{\widetilde{\mu}}=0.002$.
Table 3 presents the real values of the first three moments, standard
deviation and median for the continuous part of the Bitcoin returns along with
the confidence intervals obtained for the normal distribution. It can be
observed that the confidence intervals have a small length and cover the true
values. The null hypothesis of the Wilcoxon test is not rejected with a
p-value of 0.32404, obtained as an average over 25 simulated samples of size
1000 as before.

Table 3: 99% confidence intervals for the moments, standard deviation and median of the continuous part of the logarithmic returns of Bitcoin prices. | 1st moment$\cdot 10^{3}$ | 2d moment$\cdot 10^{3}$ | 3d moment$\cdot 10^{3}$ | SD$\cdot 10^{3}$ | Median$\cdot 10^{3}$
---|---|---|---|---|---
Lower CI | $-0.5017$ | 1.3515 | $-0.0027$ | 35.2067 | $-1.0035$
True value | 1.9634 | 1.3612 | 0.0104 | 36.854 | 2.0861
Upper CI | 4.4285 | 1.3708 | 0.0234 | 38.7042 | 5.1756

###  4.2 Estimation of the jump part

Similar to the case of media attention, we assume that the jump part of the
Bitcoin price can be described by a Lévy process. That is, for each
$k\in\mathcal{K}$ the characteristic function $\phi_{\eta_{k}-\eta_{k-1}}(u)$
admits a representation

$\phi_{\eta_{k}-\eta_{k-1}}(u)=\exp\left\\{\Delta\left(\mathrm{i}\mu^{\circ}u-\frac{1}{2}(\sigma^{\circ})^{2}u^{2}-\lambda^{\circ}+\lambda^{\circ}\mathcal{F}[p^{\circ}](u)\right)\right\\},$

where as before, $\sigma^{\circ}>0,\mu^{\circ}\in\mathbb{R}$ are some
parameters, and $p^{\circ}$ is the density of pure jumps. Thus, we can apply
the estimation procedure described in Section 2.2. As before, the values of
the technical parameters $U_{n}$, $V_{n}$ and $\varepsilon$, defined as in
Section 2.2, are determined as those minimising the mean-squared error between
the real and estimated densities, and turn out to be $U_{n}=15$, $V_{n}=23$
and $\varepsilon=0.6$. The parameter estimates are $\sigma^{\circ}=0.0576$,
$\lambda^{\circ}=0.01362$ and $\mu^{\circ}=0.0071$. Figure 8 represents the
true density of jumps in logarithmic returns of the Bitcoin price and
simulated densities obtained from 25 independent samples of size 1000. It can
be seen that the densities are very close. The Wilcoxon test does not reject
the null hypothesis with an average p-value of 0.611. In addition, Table 4
indicates that the confidence intervals for moments, standard deviation and
median include the true parameter values. Hence, we conclude that the proposed
model provides a good fit to the jump part of the Bitcoin price.

Figure 8: The true density of jumps in logarithmic returns of the Bitcoin price superimposed with densities of 25 simulated samples of size 1000.  Table 4: 99% confidence intervals for the moments, standard deviation and median of the jump part of the logarithmic returns of Bitcoin price | 1st moment | 2d moment | 3d moment | SD | Median
---|---|---|---|---|---
Lower CI | $-0.0074$ | 0.0044 | $-0.0991$ | 0.0658 | $-0.0001$
True value | 0.0028 | 0.0044 | $-0.0003$ | 0.0662 | 0.0062
Upper CI | 0.0119 | 0.0381 | 0.0442 | 0.1951 | 0.0111

##  5 Discussion

The results presented in this paper follow the promising and rather popular in
the literature idea to use the indicators of the media attention for the
prediction — or at least description — of the Bitcoin prices. We consider the
stochastic models with independent increments, namely, the Lévy processes.
Using the recent data, we provide some empirical evidence that the Lévy
processes with finite Lévy measure (e.g., as in the Kou and the Merton models)
describe well the daily returns of media attention measured by Google Trends,
but are insufficient for modelling the returns of the Bitcoin price. Hence, we
develop a model for the latter indicator using the observation that the jumps
in media attention and the Bitcoin price are closely linked. Determining the
jump times from the model constructed for media attention, we describe the
Bitcoin price by modelling the jump and the continuous parts separately. As
shown by various statistical tests, the proposed approach provides a good fit
to the considered data, allowing to successfully capture its various
statistical characteristics.

## Appendix

Figure 9: The true density of logarithmic returns of media attention (black)
superimposed with the normal density (left) and Cauchy density (right) with
estimated parameters  Figure 10: The true density of logarithmic returns of
Bitcoin price (black) superimposed with the normal density (left) and Cauchy
density (right) with estimated parameters  Figure 11: The true density of
logarithmic returns of Bitcoin price (red) superimposed with densities of 25
simulated samples of different size, and the corresponding average quadratic
errors

## References

  * Belomestny, D. (2010) Belomestny, D (2010) Spectral estimation of the fractional order of a Lévy process. The Annals of Statistics 38(1):317–351
  * Belomestny, D. (2011) Belomestny, D (2011) Statistical inference for time-changed Lévy processes via composite characteristic function estimation. The Annals of Statistics 39(4):2205–2242
  * Belomestny, D. and Panov, V. (2013) Belomestny, D and Panov, V (2013) Estimation of the activity of jumps in time-changed lévy models. Electron J Statist 7:2970–3003
  * Belomestny, D., and Reiss, M. (2015) Belomestny, D, and Reiss, M (2015) Estimation and Calibration of Lévy Models via Fourier Methods. In: Belomestny, D, Comte, F, Genon-Catalot, V, Masuda, H, Reiss, M (ed) Lévy Matters IV. Estimation for Discretly Observed Lévy processes, Springer, pp 1–76
  * Belomestny, D., Panov, V. and Woerner, J. (2019) Belomestny, D, Panov, V and Woerner, J (2019) Low-frequency estimation of continuous-time moving average Lévy processes. Bernoulli 25(2):902–931
  * Comte, F., and Genon-Catalot, V. (2010) Comte, F, and Genon-Catalot, V (2010) Nonparametric adaptive estimation for pure jump Lévy processes. Annales de l’Institut Henri Poincaré - Probabilités et Statistiques 46(3):595–617
  * Cont, R. and Tankov, P. (2004) Cont, R and Tankov, P (2004) Financial modelling with jump process. Chapman & Hall, CRC Press UK
  * Cretarola, A., Figà-Talamanca and Patacca (2020) Cretarola, A, Figà-Talamanca and Patacca (2020) Market attention and Bitcoin price modeling: Theory, estimation and option pricing. Decisions in Economics and Finance 43(1):187–228
  * Da, Z., Engelberg, J. and Gao, P. (2011) Da, Z, Engelberg, J and Gao, P (2011) In search of attention. The Journal of Finance 66(5):1461–1499
  * Dyhrberg, A. (2016) Dyhrberg, A (2016) Hedging capabilities of bitcoin. Is it the virtual gold? Finance Research Letters 16:139–144
  * Figa-Talamanca, G. and Patacca, M. (2019) Figa-Talamanca, G and Patacca, M (2019) Does market attention affect Bitcoin returns and volatility? Decisions in Economics and Finance 42(1):135–155
  * Gardini, M., Sabino, P. and Sasso, E. (2021) Gardini, M, Sabino, P and Sasso, E (2021) Correlating Lévy processes with self-decomposability: applications to energy markets. Decisions in Economics and Finance 44(2):1253–1280
  * Glaser, F., Zimmermann, K., Haferkorn, M., Weber, M. and Siering, M. (2014) Glaser, F, Zimmermann, K, Haferkorn, M, Weber, M and Siering, M (2014) Bitcoin-asset or currency? revealing users’ hidden intentions. Revealing Users’ Hidden Intentions (April 15, 2014) ECIS
  * Gugushvili, S. (2012) Gugushvili, S (2012) Nonparametric inference for discretely sampled Lévy processes. Annales de l’Institut Henri Poincaré - Probabilités et Statistiques 48(1):282–307
  * Härdle, W., Harvey, C., and Reule, R. (2020) Härdle, W, Harvey, C, and Reule, R (2020) Understanding cryptocurrencies. Journal of Financial Econometrics 18(2):181–208
  * Katsiampa, P. (2017) Katsiampa, P (2017) Volatility estimation for Bitcoin: A comparison of GARCH models. Economics letters 158:3–6
  * Kristoufek, L. (2013) Kristoufek, L (2013) BitCoin meets Google Trends and Wikipedia: Quantifying the relationship between phenomena of the Internet era. Scientific reports 3(1):1–7
  * Naimy, V. and Hayek, M. (2018) Naimy, V and Hayek, M (2018) Modelling and predicting the Bitcoin volatility using GARCH models. Int J Math Model Numer Optimisation 8(3):197–215
  * Neumann, M., and Reiss, M. (2009) Neumann, M, and Reiss, M (2009) Nonparametric estimation for Lévy processes from low-frequency observations. Bernoulli 15(1):223–248
  * Panov, V. and Samarin, E. (2019) Panov, V and Samarin, E (2019) Multivariate asset-pricing model based on subordinated stable processes. Applied Stochastic Models in Business and Industry 35(4):1060–1076
  * Sato, K. (1999) Sato, K (1999) Lévy processes and infinitely divisible distributions. Cambridge University Press
  * Schoutens, W. (2003) Schoutens, W (2003) Lévy processes in finance: pricing financial derivatives. Wiley Online Library
  * Tiwari, A., Kumar, S. and Pathak, R. (2019) Tiwari, A, Kumar, S and Pathak, R (2019) Modelling the dynamics of Bitcoin and Litecoin: GARCH versus stochastic volatility models. Applied Economics 51(37):4073–4082
  * Urquhart, A. (2018) Urquhart, A (2018) What causes the attention of Bitcoin? Economics Letters 166:40–44


11institutetext: Biomedical Image Analysis Group, Imperial College London,
London, UK
22institutetext: Data Science Institute, Imperial College London, London, UK
33institutetext: Department of Medicine, Imperial College London, London, UK
33email: {chen.chen15<EMAIL_ADDRESS>
# Unsupervised Multi-modal Style Transfer for Cardiac MR Segmentation

Chen Chen (✉) Equal contribution11 Cheng Ouyang* (✉) 11 Giacomo Tarroni 11 Jo
Schlemper 11
Huaqi Qiu 11 Wenjia Bai 2233 Daniel Rueckert 11

###### Abstract

In this work, we present a fully automatic method to segment cardiac
structures from late-gadolinium enhanced (LGE) images without using labelled
LGE data for training, but instead by transferring the anatomical knowledge
and features learned on annotated balanced steady-state free precession
(bSSFP) images, which are easier to acquire. Our framework mainly consists of
two neural networks: a multi-modal image translation network for style
transfer and a cascaded segmentation network for image segmentation. The
multi-modal image translation network generates realistic and diverse
synthetic LGE images conditioned on a single annotated bSSFP image, forming a
synthetic LGE training set. This set is then utilized to fine-tune the
segmentation network pre-trained on labelled bSSFP images, achieving the goal
of unsupervised LGE image segmentation. In particular, the proposed cascaded
segmentation network is able to produce accurate segmentation by taking both
shape prior and image appearance into account, achieving an average Dice score
of 0.92 for the left ventricle, 0.83 for the myocardium, and 0.88 for the
right ventricle on the test set.

##  1 Introduction

Cardiac segmentation from late-gadolinium enhanced (LGE) cardiac magnetic
resonance (CMR) which highlights myocardial infarcted tissue is of great
clinical importance, enabling quantitative measurements useful for treatment
planning and patient management. To this end, the segmentation of the
myocardium is an important first step for myocardial infarction analysis.

Since manual segmentation is tedious and likely to suffer from inter-observer
variability, it is of great interest to develop an accurate automated
segmentation method. However, this is a challenging task due to the fact that
1) the infarcted myocardium presents an enhanced and heterogeneous intensity
distribution different from the normal myocardium region and 2) the border
between infarcted myocardium and blood pool appears blurry and ambiguous [1].
While the borders of the myocardium can be difficult to delineate on LGE
images, they are clear and easy to identify on the balanced steady-state free
precession (bSSFP) CMR images, which have high signal-to-noise ratio and whose
contrast is less sensitive to pathology (see red arrows in Fig. 1 (a)).
Conventional methods [2, 3] use the segmentation result from the bSSFP CMR of
the same patient as prior knowledge to assist the segmentation on LGE CMR
images. These methods generally require accurate registration between the
bSSFP and LGE images, which can be challenging as the imaging field-of-view
(FOV), image contrast and resolution between the two acquisitions can vary
significantly [1, 4]. Fig. 1 (b) visualizes the discrepancy between the
intensity distributions of the two imaging modalities in the cardiac
structures (specifically, left ventricle (LV), myocardium (MYO), and right
ventricle (RV)).

Figure 1: The differences of image appearance (a) and intensity distributions
(b) in the cardiac region (the union of LV, MYO, RV) between LGE images and
bSSFP images

Most recently, a deep neural network-based method has been proposed to segment
the three cardiac structures directly from LGE images [5], reporting superior
performance. However, this supervised segmentation method requires a large
amount of labelled LGE data. Because of the heterogeneous intensity
distribution of the myocardium in LGE images and the scarcity of experienced
image analysts, it is difficult to perform accurate manual segmentations on
LGE images and collect a large training set, compared to that on bSSFP images.

In this paper, we present a fully automatic framework that addresses the above
mentioned issues by training a segmentation model without using manual
annotations on LGE images. This is achieved by transferring the anatomical
knowledge and features learned on annotated bSSFP images, which are easier to
acquire. Our framework mainly consists of two neural networks:

  * •

A multi-modal image translation network: this network is used for translating
annotated bSSFP images into LGE images through style transfer. Of note, the
network is trained in an unsupervised fashion where the training bSSFP images
and LGE images are unpaired. In addition, unlike common one-to-one translation
networks, this network allows the generation of multiple synthetic LGE images
conditioned on a single bSSFP image;

  * •

A cascaded segmentation network for LGE images consisting of two U-net [6]
models (Cascaded U-net): this network is first trained using the labelled
bSSFP images and then fine-tuned using the synthetic LGE data generated by the
image translation network.

The main contributions of our work are the following: 1) we employ a
translation network that can generate realistic and diverse synthetic LGE
images given a single bSSFP image. This network enables generative model-based
data augmentation for unsupervised domain adaptation, which not only closes
the domain gap between the two modalities, but also improves the
generalization properties of the following segmentation network by increasing
data variety; 2) we demonstrate that the proposed two-stage cascaded network,
which takes both anatomical shape information and image appearance information
into account, produces accurate segmentation on LGE images, greatly
outperforming baseline methods; 3) the proposed framework can be easily
extended to other unsupervised cross-modality domain adaptation applications
where labels of one modality are not available.

##  2 Methodology

The proposed method aims at learning an LGE image segmentation model using
labelled bSSFP $\\{({\bm{{x}}}_{b},{\bm{{y}}}_{b})\\}$ and unlabelled LGE
$\\{{\bm{{x}}}_{l}\\}$ only. Specifically, the proposed method is a two-stage
framework. In the first stage, an unsupervised image translation network is
trained to translate each bSSFP image ${\bm{{x}}}_{b}$ into multiple instances
of LGE-like images, noted as $\\{{\bm{{x}}}_{bl}\\}$. In the second stage,
these LGE-stylized bSSFP images are used together with their original labels
$\\{({\bm{{x}}}_{bl},{\bm{{y}}}_{b})\\}$ to adapt an image segmentation
network pre-trained on labelled bSSFP images to segment LGE images.

###  2.1 Image Translation

We employ the state-of-the-art multi-modal unsupervised image-to-image
translation network (MUNIT) [7] as our multi-modal image translator. Let
$\\{{\bm{{x}}}_{l}\\}$ and $\\{{\bm{{x}}}_{b}\\}$ denote unpaired images from
the two different imaging modalities (domains): LGE and bSSFP, given an image
drawn from one domain as input, the network is able to change the appearance
(i.e. image style) of the image to that of the other domain while preserving
the underlying anatomical structure [8]. This is achieved by learning
disentangled image representations.

Figure 2: Overview of the multi-modal image translation network. The network
employs the structure of MUNIT [7], which consists of two encoder-decoder
pairs for the two domains: bSSFP and LGE, respectively.

As shown in Fig. 2, each image ${\bm{{x}}}$ is disentangled into (a) a domain-
invariant content code $\bm{c}$: $\bm{c}=E^{c}{({\bm{{x}}})}$ and (b) a
domain-specific style code $\bm{s}$: $\bm{s}=E^{s}{({\bm{{x}}})}$ using the
content encoder $E^{c}$ and the style encoder $E^{s}$ relative to its domain
where the content code captures the anatomical structure and the style code
carries the information for rendering the structure which is determined by the
imaging modality. The image-to-image translation from one domain to the other
is achieved by swapping latent codes in two domains. For example, translating
a bSSFP image $\bm{x_{b}}$ to be stylized as LGE, is achieved by feeding the
content code $\bm{c_{b}}$ for the bSSFP image and the style code $\bm{s_{l}}$
into the LGE decoder $D_{l}$:
${\bm{{x}}}_{bl}=D_{l}({\bm{{c}}}_{b},{\bm{{s}}}_{l})$.

Of note, during training, each style encoder is trained to embed images into a
latent space that matches the standard Gaussian distribution
${\mathcal{{N}}}(0,I)$, minimizing the Kullback-Leibler (KL) divergence
between the two. This allows to generate an arbitrary number of synthetic LGE
images given a single bSSFP image during inference, by repeatedly sampling the
style code from the prior distribution ${\mathcal{{N}}}(0,I)$. Of note,
although this prior distribution is unimodal, the distribution of translated
images in the output space is multi-modal thanks to the nonlinearity of the
decoder[7]. We apply this translation network to translate annotated bSSFP
images, resulting in a synthetic labelled LGE dataset, which will then be used
to finetune a segmentation network. For more details about training the
translation network, readers are referred to the original work by Huang et al.
[7].

###  2.2 Image Segmentation

Figure 3: Overview of the two-stage cascaded segmentation network. The
architecture of each U-net is the same as the one of the vanilla U-net, except
for two main differences: (1) batch normalization is applied after each
convolutional layer; (2) a dropout layer (dropout rate=0.1) is applied after
each concatenation operation in the network’s expanding path to encourage
model generalizability. Of note, in this diagram, we simplify the training
procedure by omitting the pre-training procedure using labelled bSSFP images.

Let $\bm{x}_{l}$ be an observed LGE image, the aim of the segmentation task is
to estimate label maps $\bm{y}_{l}$ having observed $\bm{x}_{l}$ by modeling
the posterior $p(\bm{y}_{l}|\bm{x}_{l})$. Inspired by curriculum learning [9]
and transfer learning, we first train a segmentation network using annotated
bSSFP images (source domain; easy examples) and then fine-tune it to segment
LGE images (target domain; hard examples). Since labelled LGE images
$\\{({\bm{{x}}}_{l},{\bm{{y}}}_{l})\\}$ are not available for finetuning, we
use a synthetic dataset
$\mathcal{X}_{bl}:\\{({\bm{{x}}}_{bl},{\bm{{y}}}_{b})\\}_{1..N}$ generated by
the aforementioned multi-modal image translator. Ideally, the posterior
modelled by the network $p(\bm{y}_{b}|\bm{x}_{bl})$ matches
$p(\bm{y}_{l}|\bm{x}_{l})$ when image space and label space are shared. For
simplicity, we use $\bm{x}$ and $\bm{y}$ to denote an image and its
corresponding label map from the synthetic dataset in the following
paragraphs.

The segmentation network is a two-stage cascaded network which consists of two
U-nets [6], see Fig. 3. Specifically, given an image $\bm{x}$ as input, the
first U-net (U-net 1) aims at predicting four-class pixel-wise probabilistic
maps $\bm{p_{1}}=f_{\text{U-net}}^{1}(\bm{x};\theta)$ for the three cardiac
structures (i.e. LV, MYO, RV) and the background class (BG). Inspired by the
auto-context architecture [10], we combine these learned probabilistic maps
$\bm{p_{1}}$ from the first network with the raw image $\bm{x}$ to form a
5-channel input to train the second U-net (U-net 2) for fine-grained
segmentation: $\bm{{p_{2}}}={{f_{\text{U-net}}^{2}}}(\bm{x},\bm{p_{1}};\phi)$.
By combining the appearance information from the image ${\bm{{x}}}$ with the
shape prior information from the initial segmentation $\bm{p_{1}}$ as input,
the cascaded network has the potential to produce more precise and robust
segmentations even in the presence of unclear boundaries for the different
cardiac structures.

To train the network, we use a composite segmentation loss function
$\mathcal{L}_{seg}$ which consists of two loss terms:
$\mathcal{L}_{seg}=\mathcal{L}_{wce}+\lambda\mathcal{L}_{edge}.$ The first
term $\mathcal{L}_{wce}$ is a weighted cross entropy loss:
$\mathcal{L}_{wce}=-\sum_{m}\omega^{m}\bm{y}^{m}\log\left(\bm{p}^{m}\right)$
where $w^{m}$ denotes the weight for class $m$ and $\bm{p}^{m}$ is the
corresponding predicted probability map. We set the weight for myocardium
$\omega^{MYO}$ to be higher than the weights for the other three classes to
address class imbalance problem since there is a lower percentage of pixels
that corresponds to the myocardium class in CMR images. The second term
$\mathcal{L}_{edge}$ is an edge-based loss which penalizes the disagreement on
the contours of the cardiac structures. Specifically, we apply two 2D $3\times
3$ Sobel filters [11] $S_{k}$ (k=1,2) to the soft prediction maps $\bm{p}$ as
well as the one-hot heatmaps $\bm{y}$ of the ground truth to extract edge
information along horizontal and vertical directions. The edge loss is then
computed by calculating the $l_{2}$ distance between the predicted edge maps
and the ground truth edge maps: $\mathcal{L}_{edge}=\sum_{m,m\neq
BG}\sum_{k=1,2}\left\|f_{S_{k}}({\bm{p}^{m}})-f_{S_{k}}(\bm{y}^{m})\right\|_{2}$,
where $f_{S_{k}}(\bm{p}^{m})$ is the edge map extracted by applying the sobel
filter $S_{k}$ to the predicted probabilistic map $\bm{p}^{m}$ for foreground
class $m$.

By using the edge loss together with the weighted cross entropy for
optimization, the network is encouraged to focus more on the contours of the
three structures and the myocardium, which are usually more difficult to
delineate. In our experiments, we set $\lambda$ = 0.5 to balance the
contribution of the two losses.

###  2.3 Post-processing

At inference time, each slice from a previously unseen LGE stack is fed to the
cascaded network to get the probabilistic maps for the four classes. Dense
conditional random field (CRF) [12] is then applied to refine the 2D predicted
segmentation mask slice by slice. After that, 3D morphological dilation and
erosion operations are applied to the whole segmentation stack to further
improve the global smoothness. In particular, we perform the operations in a
hierarchical order: first we apply them to the binary map covering all the
three structures, then to the MYO and the LV labels, separately.

##  3 Experiments and Results

###  3.1 Data

The framework was trained and evaluated on the Multi-sequence Cardiac MR
Segmentation Challenge (MS-CMRSeg 2019)
dataset111https://zmiclab.github.io/mscmrseg19/. We used a subset of 40 bSSFP
and 40 LGE images to train the image translation network. Then, we created a
synthetic dataset by applying the learned translation network to 30 labelled
bSSFP images. Specifically, for each bSSFP image, we randomly sampled the
style code from ${\mathcal{{N}}}(0,I)$ five times, resulting in a set of 150
synthetic LGE images in total. This synthetic dataset and the original 30
bSSFP images with corresponding labels formed the training set for the
segmentation network. Exemplar results of these synthetic LGE images are
provided in the supplemental material. For validation, we used a subset of 5
annotated LGE images provided by the challenge organizers.

###  3.2 Implementation Details

Image preprocessing. To deal with the different image size and heterogeneous
pixel spacing between different imaging modalities, all images were resampled
to a pixel spacing of $1.25~{}mm\times 1.25~{}mm$ and then cropped to
$192\times 192$ pixels, with the heart roughly at the center of each image.
This spatial normalization would reduce the computational cost and task
complexity in the following training procedure of image translation and
segmentation, making the networks focus on the relevant regions. To identify
the heart, we trained a localization network based on U-net using the 30
annotated bSSFP images in the training set to produce rough segmentations for
the three structures. The localization network employs instance normalization
layers which perform style normalization [13], encouraging the network
invariance to image style changes (e.g. image contrast). As a result, the
network is able to produce coarse masks localizing the heart on all bSSFP
images and most LGE images even though it was trained on bSSFP images only. In
case that this network might fail to locate the heart on certain LGE slices,
we summed the segmentation masks across slices in each volume and then cropped
them according to the center of the aggregated mask. After cropping, each
image was intensity normalized.

Network training. (1) For the image translation network, we used the official
implementation222https://github.com/NVlabs/MUNIT of [7]. Network configuration
and hyper-parameters were kept the same as in [7] except the input and output
images are 2D, single-channel. It was trained for 20k iterations with a batch
size of 1. (2) For the segmentation network, we first trained the first U-net
with the labelled bSSFP images and then fine-tuned it with synthetic LGE
images. This procedure was replicated to train the second U-net with the
parameters of the first U-net being fixed. Both networks were optimized using
the composite loss $\mathcal{L}_{seg}$ where adam was used for stochastic
gradient descent. The learning rate was initially set to 0.001 and was then
decreased to $1\times 10^{-5}$ for fine-tuning. The weights for BG, LV, MYO,
and RV in $\mathcal{L}_{wce}$ were empirically set to $0.2:0.25:0.3:0.25$.
During training, we applied data augmentation on the fly. Specifically,
elastic deformations, random scaling and random rotations as well as gamma
augmentation [14] were used. The algorithm was implemented using python and
PyTorch and was trained for 1000 epochs in total on an
NVIDIA${}^{\tiny{\text{\textregistered}}}$ Tesla P40 GPU.

###  3.3 Results

To evaluate the accuracy of segmentation results, the Dice metric and the
average surface distance (ASD) between the automatic segmentation and the
corresponding manual segmentation for each volume were calculated.

We compare the proposed method with two baseline methods: (1) a registration-
based method and (2) a single U-net. Specifically, for the registration-based
method, each LGE segmentation result was obtained by directly registering the
corresponding bSSFP labels to the LGE image using MIRTK toolkit
333https://mirtk.github.io/ for ease of comparison. The transformation matrix
was learned by applying mutual information-based registration
(Rigid+Affine+FFD) between the two images. For U-net, we trained it with two
settings: a) U-net: trained on labelled bSSFP images only; b) U-net with fine-
tuning (FT): trained on labelled bSSFP images and then fine-tuned using the
synthetic LGE data, which is the same training procedure of the proposed
method. Quantitative and qualitative results are shown in Table 1 and Fig. 4.

While the registration-based method (MIRTK) outperforms the U-net (see row 1
and row 2 in Table 1), it still fails to produce accurate segmentation on the
myocardium (see the red number in row 1), indicating the limitation of this
registration-based method. However, by contrast, neural network-based methods
(row 3-5) fine-tuned using the _synthetic LGE dataset_ significantly improves
the segmentation accuracy, increasing the Dice score for MYO by $\sim 15\%$.
This improvement demonstrates the learned translation network is capable of
generating realistic LGE images while preserving the domain-invariant
structural information that is informative to optimize the segmentation
network. In particular, compared to U-net (FT), the proposed Cascaded U-net
(FT) achieves more accurate segmentation performance with improvement in terms
of both Dice and ASD (see blue numbers). The model even produces robust
segmentation results on the challenging apical and basal slices (please see
the last column in Fig. 4). This demonstrates the benefit of integrating the
high-level shape knowledge and low-level image appearance to guide the
segmentation procedure. In addition, the proposed post-processing further
refines the segmentation results through smoothing, reducing the average ASD
from 1.37 to 1.26 (see the last row in Table 1).

Table 1: Dice scores and ASD (mm) of the proposed segmentation method
(Cascaded U-net) and baseline methods on the validation set. Blue numbers
indicate the best scores among the results obtained by those methods before
post-processing (PP) whereas red numbers are those mean Dice scores under
0.700. FT: fine-tuning using the synthetic LGE dataset. N/A means that the ASD
value cannot be calculated due to missing predictions for that cardiac
structure.

| Dice | ASD
---|---|---
Method | LV | MYO | RV | AVG* | LV | MYO | RV | AVG*
MIRTK | 0.819 | 0.665 | 0.831 | 0.772 | 2.56 | 1.65 | 2.11 | 2.11
U-net | 0.624 | 0.441 | 0.577 | 0.547 | 10.03 | 6.07 | N/A | N/A
U-net (FT) | 0.874 | 0.781 | 0.896 | 0.850 | 1.78 | 1.50 | 1.28 | 1.52
Cascaded U-net (FT) | 0.895 | 0.812 | 0.898 | 0.868 | 1.41 | 1.46 | 1.23 | 1.37
Cascaded U-net (FT) + PP | 0.897 | 0.816 | 0.895 | 0.869 | 1.17 | 1.42 | 1.18 | 1.26

  * •

* For ease of comparison, we calculate the average (AVG) Dice score and the average ASD score over the three structures for each method.

Figure 4: Segmentation results for the proposed Cascaded U-net and the
baseline approaches. Our proposed method (the right-most column) produces more
anatomically plausible segmentation results on the images, greatly
outperforming the baseline methods, especially in the challenging cases: the
apical (the top row) and the basal slices (the bottom row).

Finally, we applied ensemble learning to improve our model’s performance in
the test phase. Specifically, we trained the proposed segmentation network for
multiple times, each time regenerating a new synthetic LGE dataset for fine-
tuning. We trained four models in total. Our final submission result for each
test image was obtained by averaging the probabilistic maps from these models
and then assigning to each pixel the class with the highest score. In the
testing stage of the competition, the method achieves very promising
segmentation performance on a relative large test set (40 subjects), with an
average Dice score of 0.92 for LV, 0.83 for MYO, and 0.88 for RV; an ASD of
1.66 for LV, 1.76 for MYO, and 2.16 for RV.

##  4 Conclusion

In this paper, we showed that synthesizing multi-modal LGE images from
labelled bSSFP images to finetune a pre-trained segmentation network shows
impressive segmentation performance on LGE images even though the network has
not seen _real_ labelled LGE images before. We also demonstrated that the
proposed segmentation network (Cascaded U-net) outperformed the baseline
methods by a significant margin, suggesting the benefit of integrating the
high-level shape knowledge and low-level image appearance to guide the
segmentation procedure. More importantly, our cascaded segmentation network is
independent of the particular architecture of underlying convolutional neural
networks. In other words, the basic neural network (U-net) in our work can be
replaced with any of the state-of-the-art segmentation network to potentially
improve the prediction accuracy and robustness. Moreover, the proposed
solution based on unsupervised multi-modal style transfer is not only limited
to the cardiac image segmentation but can be extended to other multi-modal
image analysis tasks where the manual annotations of one modality are not
available. Future work will focus on the application of the method to the
problems such as domain adaptation for multi-modality brain segmentation.

## References

  * [1] Xiahai Zhuang.  Multivariate mixture model for cardiac segmentation from multi-sequence MRI.  In MICCAI, pages 581–588, 2016.
  * [2] YingLi Lu et al.  Automatic myocardium segmentation of LGE MRI by deformable models with prior shape data.  JCMR, 15(1):P14, 2013.
  * [3] Qian Tao et al.  Automated left ventricle segmentation in late gadolinium-enhanced MRI for objective myocardial scar assessment.  JMRI, 42(2):390–399, 2015.
  * [4] Xiahai Zhuang.  Multivariate mixture model for myocardium segmentation combining multi-source images.  PAMI, 2018.
  * [5] Qian Yue et al.  Cardiac Segmentation from LGE MRI Using Deep Neural Network Incorporating Shape and Spatial Priors.  MICCAI, 2019.
  * [6] Ronneberger et al.  U-net: Convolutional networks for biomedical image segmentation.  In MICCAI, 2015.
  * [7] Xun Huang et al.  Multimodal unsupervised image-to-image translation.  In ECCV, 2018.
  * [8] Chen Qin, Bibo Shi, Rui Liao, Tommaso Mansi, Daniel Rueckert, and Ali Kamen.  Unsupervised deformable registration for multi-modal images via disentangled representations.  In Information Processing in Medical Imaging, pages 249–261. Springer, Cham, June 2019.
  * [9] Yoshua Bengio, Jérôme Louradour, Ronan Collobert, and Jason Weston.  Curriculum learning.  In Proceedings of the 26th Annual International Conference on Machine Learning, ICML ’09, pages 41–48, New York, NY, USA, 2009. ACM.
  * [10] Zhuowen Tu and Xiang Bai.  Auto-context and its application to high-level vision tasks and 3D brain image segmentation.  PAMI, 2010.
  * [11] Irwin Sobel and G Feldman.  A 3x3 isotropic gradient operator for image processing.  Pattern Classification and Scene Analysis, pages 271–272, 01 1973\.
  * [12] Philipp Krähenbühl et al.  Efficient inference in fully connected CRFs with gaussian edge potentials.  In NeuralIPS, 2011.
  * [13] Xun Huang et al.  Arbitrary style transfer in real-time with adaptive instance normalization.  In ICCV, 2017.
  * [14] Chen Chen, Wenjia Bai, and Daniel Rueckert.  Multi-task learning for left atrial segmentation on GE-MRI.  In Statistical Atlases and Computational Models of the Heart. Atrial Segmentation and LV Quantification Challenges, 2019.

## Supplemental Material

Figure 5: Exemplar synthetic LGE images generated from bSSFP images using the
multi-modal image translation network. Given one bSSFP image (column 1), the
translation network translates the image into multi-modal LGE-like images
(column 2 to 4). These translated images differ in image brightness and
contrast as well as the intensity distribution in the cardiac region, while
preserving the same cardiac anatomy. These synthetic images, in together with
the annotations on the original bSSFP images (the last column) contribute to
the synthetic dataset which is used to fine-tune the proposed segmentation
network.